REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_advance_3d.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2
3using namespace amrex;
4
5/** Corresponds to step3d_uv.F and step3d_t.F in ROMS
6 *
7 * @param[in ] lev level of refinement (coarsest level is 0)
8 * @param[inout] mf_cons scalar variables (temp, salt, scalar, etc)
9 * @param[inout] mf_u u-velocity
10 * @param[inout] mf_v v-velocity
11 * @param[inout] mf_sstore intermediate-step scalar variables
12 * @param[inout] mf_ru RHS of total u-velocity
13 * @param[inout] mf_rv RHS of total v-velocity
14 * @param[in ] mf_DU_avg1 time-averaged u-flux for 2D equations
15 * @param[in ] mf_DU_avg2 time-averaged u-flux for 3D equation coupling
16 * @param[in ] mf_DV_avg1 time-averaged v-flux for 2D equations
17 * @param[in ] mf_DV_avg2 time-averaged v-flux for 3D equation coupling
18 * @param[inout] mf_ubar vertically integrated u-momentum
19 * @param[inout] mf_vbar vertically integrated v-momentum
20 * @param[inout] mf_Akv vertical viscosity coefficient
21 * @param[inout] mf_Akt vertical diffusivity coefficient
22 * @param[inout] mf_Hz vertical height of cells
23 * @param[inout] mf_Huon u-volume flux
24 * @param[inout] mf_Huon v-volume flux
25 * @param[inout] mf_z_w vertical coordinates on w points
26 * @param[in ] mf_h bathymetry
27 * @param[in ] mf_pm 1 / dx
28 * @param[in ] mf_pn 1 / dy
29 * @param[in ] mf_msku land-sea mask at u-points
30 * @param[in ] mf_mskv land-sea mask at v-points
31 * @param[in ] N number of vertical levels
32 * @param[in ] dt_lev time step at this refinement level
33 */
34void
35REMORA::advance_3d (int lev, MultiFab& mf_cons,
36 MultiFab& mf_u , MultiFab& mf_v ,
37 MultiFab* mf_sstore,
38 MultiFab* mf_ru , MultiFab* mf_rv,
39 std::unique_ptr<MultiFab>& mf_DU_avg1,
40 std::unique_ptr<MultiFab>& mf_DU_avg2,
41 std::unique_ptr<MultiFab>& mf_DV_avg1,
42 std::unique_ptr<MultiFab>& mf_DV_avg2,
43 std::unique_ptr<MultiFab>& mf_ubar,
44 std::unique_ptr<MultiFab>& mf_vbar,
45 std::unique_ptr<MultiFab>& mf_Akv,
46 std::unique_ptr<MultiFab>& mf_Akt,
47 std::unique_ptr<MultiFab>& mf_Hz,
48 std::unique_ptr<MultiFab>& mf_Huon,
49 std::unique_ptr<MultiFab>& mf_Hvom,
50 std::unique_ptr<MultiFab>& mf_z_w,
51 MultiFab const* mf_h,
52 MultiFab const* mf_pm,
53 MultiFab const* mf_pn,
54 MultiFab const* mf_mskr,
55 MultiFab const* mf_msku,
56 MultiFab const* mf_mskv,
57 const int N, Real dt_lev)
58{
59 BL_PROFILE("REMORA::advance_3d()");
60 const int nrhs = 0;
61 int nnew = 0;
62
63 int iic = istep[lev];
64 int ntfirst = 0;
65
66 // ROMS accumulates the finer level's perimeter tracer flux over its substeps and
67 // subtracts the mismatch against the coarse flux from the coarse ring (correct_tracer,
68 // nesting.F). YAFluxRegister is the same construction.
69 const bool refluxing = do_reflux && do_substep && finest_level > 0 &&
70 (solverChoice.coupling_type == CouplingType::two_way);
71
72 // Because zeta may have changed
73 stretch_transform(lev);
74
75 // These temporaries used to be made in advance_3d_ml and passed in;
76 // now we make them here
77
78 const BoxArray& ba = mf_cons.boxArray();
79 const DistributionMapping& dm = mf_cons.DistributionMap();
80
81 //Only used locally, probably should be rearranged into FArrayBox declaration
82 MultiFab mf_AK (convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0)); //2d missing j coordinate
83 MultiFab mf_DC (ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
84 MultiFab mf_Hzk(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
85
86 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
87 {
88 Array4<Real > const& u = mf_u.array(mfi);
89 Array4<Real > const& v = mf_v.array(mfi);
90
91 Array4<Real > const& ru = mf_ru->array(mfi);
92 Array4<Real > const& rv = mf_rv->array(mfi);
93
94 Array4<Real > const& AK = mf_AK.array(mfi);
95 Array4<Real > const& DC = mf_DC.array(mfi);
96
97 Array4<Real > const& Hzk = mf_Hzk.array(mfi);
98 Array4<Real > const& Akv = mf_Akv->array(mfi);
99
100 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
101
102 Array4<Real const> const& DU_avg1 = mf_DU_avg1->const_array(mfi);
103 Array4<Real const> const& DV_avg1 = mf_DV_avg1->const_array(mfi);
104
105 Array4<Real const> const& pm = mf_pm->const_array(mfi);
106 Array4<Real const> const& pn = mf_pn->const_array(mfi);
107
108 Array4<Real const> const& msku = mf_msku->const_array(mfi);
109 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
110
111 Box bx = mfi.tilebox();
112 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
113 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
114
115 Box xbx = mfi.nodaltilebox(0);
116 Box ybx = mfi.nodaltilebox(1);
117
118 Box gbx2D = gbx2;
119 gbx2D.makeSlab(2,0);
120
121 Box tbxp1 = bx;
122 Box tbxp11 = bx;
123 Box tbxp2 = bx;
124 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
125 tbxp2.grow(IntVect(NGROW,NGROW,0));
126 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
127
128 FArrayBox fab_FC(convert(gbx2,IntVect(0,0,1)),1,amrex::The_Async_Arena());
129 FArrayBox fab_BC(convert(gbx2,IntVect(0,0,1)),1,amrex::The_Async_Arena());
130 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
131 FArrayBox fab_W(tbxp2,1,amrex::The_Async_Arena());
132
133 auto FC = fab_FC.array();
134 auto BC = fab_BC.array();
135 auto CF = fab_CF.array();
136
137 Real cff;
138 if (iic==ntfirst) {
139 cff=Real(0.25)*dt_lev;
140 } else if (iic==ntfirst+1) {
141 cff=Real(0.25)*dt_lev*Real(3.0)/Real(2.0);
142 } else {
143 cff=Real(0.25)*dt_lev*Real(23.0)/Real(12.0);
144 }
145
146 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
147 {
148 u(i,j,k) += cff * (pm(i,j,0)+pm(i-1,j,0)) * (pn(i,j,0)+pn(i-1,j,0)) * ru(i,j,k,nrhs);
149 u(i,j,k) *= two / (Hz(i-1,j,k) + Hz(i,j,k));
150 });
151
152 ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
153 {
154 v(i,j,k) += cff * (pm(i,j,0)+pm(i,j-1,0)) * (pn(i,j,0)+pn(i,j-1,0)) * rv(i,j,k,nrhs);
155 v(i,j,k) *= two / (Hz(i,j-1,k) + Hz(i,j,k));
156 });
157
158 // NOTE: DC is only used as scratch in vert_visc_3d -- no need to pass or return a value
159 // NOTE: may not actually need to set these to zero
160
161 // Reset to zero on the box on which they'll be used
162 mf_DC[mfi].template setVal<RunOn::Device>(zero,xbx);
163 fab_CF.template setVal<RunOn::Device>(zero,xbx);
164
165 vert_visc_3d(xbx,1,0,u,Hz,Hzk,AK,Akv,BC,DC,FC,CF,nnew,N,dt_lev);
166
167 // Reset to zero on the box on which they'll be used
168 mf_DC[mfi].template setVal<RunOn::Device>(zero,ybx);
169 fab_CF.template setVal<RunOn::Device>(zero,ybx);
170
171 vert_visc_3d(ybx,0,1,v,Hz,Hzk,AK,Akv,BC,DC,FC,CF,nnew,N,dt_lev);
172
173 // Reset to zero on the box on which they'll be used
174 mf_DC[mfi].template setVal<RunOn::Device>(zero,xbx);
175 fab_CF.template setVal<RunOn::Device>(zero,xbx);
176
177 vert_mean_3d(xbx,1,0,u,Hz,DU_avg1,DC,CF,pn,msku,nnew,N);
178
179 // Reset to zero on the box on which they'll be used
180 mf_DC[mfi].template setVal<RunOn::Device>(zero,ybx);
181 fab_CF.template setVal<RunOn::Device>(zero,ybx);
182
183 vert_mean_3d(ybx,0,1,v,Hz,DV_avg1,DC,CF,pm,mskv,nnew,N);
184 }
185
186 // Apply physical boundary conditions to u and v
187 (*physbcs[lev])(mf_u,*mf_msku,0,1,mf_u.nGrowVect(),t_old[lev],xvel_bc(),0,*xvel_old[lev]);
188 (*physbcs[lev])(mf_v,*mf_mskv,0,1,mf_v.nGrowVect(),t_old[lev],yvel_bc(),0,*yvel_old[lev]);
189
190#ifdef REMORA_USE_NETCDF
191 // Fill the data which is stored in the boundary data read from netcdf files
192 if (solverChoice.boundary_from_netcdf)
193 {
194 fill_from_bdyfiles(lev, mf_u,*mf_msku,t_old[lev],xvel_bc(),BdyVars::u,0,0,*xvel_old[lev],dt_lev);
195 fill_from_bdyfiles(lev, mf_v,*mf_mskv,t_old[lev],yvel_bc(),BdyVars::v,0,0,*yvel_old[lev],dt_lev);
196 }
197
198 if (solverChoice.do_rivers) {
199 river_source_transport->update_interpolated_to_time(model_time(t_old[lev]));
200 }
201#endif
202 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
203 {
204 Array4<Real > const& u = mf_u.array(mfi);
205 Array4<Real > const& v = mf_v.array(mfi);
206
207 Array4<Real > const& DC = mf_DC.array(mfi);
208
209 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
210
211 Array4<Real const> const& DU_avg1 = mf_DU_avg1->const_array(mfi);
212 Array4<Real const> const& DV_avg1 = mf_DV_avg1->const_array(mfi);
213
214 Array4<Real const> const& DU_avg2 = mf_DU_avg2->const_array(mfi);
215 Array4<Real const> const& DV_avg2 = mf_DV_avg2->const_array(mfi);
216
217 Array4<Real> const& ubar = mf_ubar->array(mfi);
218 Array4<Real> const& vbar = mf_vbar->array(mfi);
219
220 Array4<Real> const& Huon = mf_Huon->array(mfi);
221 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
222
223 Array4<Real const> const& pm = mf_pm->const_array(mfi);
224 Array4<Real const> const& pn = mf_pn->const_array(mfi);
225 Array4<Real const> const& msku = mf_msku->const_array(mfi);
226 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
227
228 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
229
230 FArrayBox fab_FC(gbx2,1,amrex::The_Async_Arena());
231 auto FC = fab_FC.array();
232
233#ifdef REMORA_USE_NETCDF
234 if (solverChoice.do_rivers) {
235 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
236 Array4<Real const> const& z_w = mf_z_w->const_array(mfi);
237 Array4<int const> const& river_pos = vec_river_position[lev]->const_array(mfi);
238 Array4<Real const> const& river_transport = river_source_transport->fab_interp->array();
239 int* river_direction_d = river_direction.data();
240 ParallelFor(gbx1, [=] AMREX_GPU_DEVICE(int i, int j, int k)
241 {
242 int iriver = river_pos(i,j,0);
243 if (iriver >= 0) {
244 if (river_direction_d[iriver] == 0) {
245 Real on_u = two / (pn(i,j,0)+pn(i-1,j,0));
246 Real cff = one / (on_u * Real(0.5) * (z_w(i-1,j,k+1) - z_w(i-1,j,k) + z_w(i,j,k+1) - z_w(i,j,k)));
247 u(i,j,k) = cff * river_transport(iriver,0,k);
248 } else {
249 Real om_v = two / (pm(i,j,0)+pm(i,j-1,0));
250 Real cff = one / (om_v * Real(0.5) * (z_w(i,j-1,k+1) - z_w(i,j-1,k) + z_w(i,j,k+1) - z_w(i,j,k)));
251 v(i,j,k) = cff * river_transport(iriver,0,k);
252 }
253 }
254 });
255 }
256#endif
257
258#if 0
259 // Reset to zero on the box on which they'll be used
260 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(xbx,IntVect(0,0,1)));
261 fab_CF.template setVal<RunOn::Device>(zero,grow(xbx,IntVect(0,0,1)));
262
263 update_massflux_3d(xbx,1,0,u,ubar,Huon,Hz,pn,DU_avg1,DU_avg2,DC,FC,nnew);
264
265 // Reset to zero on the box on which they'll be used
266 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(ybx,IntVect(0,0,1)));
267 fab_CF.template setVal<RunOn::Device>(zero,grow(ybx,IntVect(0,0,1)));
268
269 update_massflux_3d(ybx,0,1,v,vbar,Hvom,Hz,pm,DV_avg1,DV_avg2,DC,FC,nnew);
270
271#else
272
273 // Reset to zero on the box on which they'll be used
274 fab_FC.template setVal<RunOn::Device>(zero,gbx2);
275 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(gbx2,IntVect(0,0,1)));
276 update_massflux_3d(lev,gbx2,1,0,u,ubar,Huon,Hz,pn,DU_avg1,DU_avg2,DC,FC,msku,nnew);
277
278 // Reset to zero on the box on which they'll be used
279 fab_FC.template setVal<RunOn::Device>(zero,gbx2);
280 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(gbx2,IntVect(0,0,1)));
281 update_massflux_3d(lev,gbx2,0,1,v,vbar,Hvom,Hz,pm,DV_avg1,DV_avg2,DC,FC,mskv,nnew);
282
283#endif
284 }
285
286 // WE BELIEVE THESE VALUES SHOULD ALREADY BE FILLED
287 // mf_Huon->FillBoundary(geom[lev].periodicity());
288 // mf_Hvom->FillBoundary(geom[lev].periodicity());
289
290 // ************************************************************************
291 // This should fill both temp and salt with temp/salt currently in cons_old
292 // ************************************************************************
293
294 MultiFab mf_W(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW+1,NGROW+1,0));
295 mf_W.setVal(zero);
296 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
297 {
298 Array4<Real> const& Huon = mf_Huon->array(mfi);
299 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
300
301 Array4<Real const> const& z_w = mf_z_w->const_array(mfi);
302 Array4<Real const> const& h = mf_h->const_array(mfi);
303
304 Box bx = mfi.tilebox();
305 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
306 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
307 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
308
309 Box tbxp1 = bx;
310 Box tbxp11 = bx;
311 Box tbxp2 = bx;
312 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
313 tbxp2.grow(IntVect(NGROW,NGROW,0));
314 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
315
316 FArrayBox fab_FC(surroundingNodes(gbx2,2),1,amrex::The_Async_Arena());
317 FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
318 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
319
320 auto W = mf_W.array(mfi);
321
322 //
323 //------------------------------------------------------------------------
324 // Vertically integrate horizontal mass flux divergence.
325 //------------------------------------------------------------------------
326 //
327 //Should really use gbx3uneven
328 //TODO: go over these boxes and compare to other spots where we do the same thing
329 Box gbx1D = gbx1;
330 gbx1D.makeSlab(2,0);
331
332 ParallelFor(gbx1D, N+1,
333 [=] AMREX_GPU_DEVICE (int i, int j, int , int kk)
334 {
335 // Starting with zero vertical velocity at the bottom, integrate
336 // from the bottom (k=0) to the free-surface (k=N). The w(:,:,N(ng))
337 // contains the vertical velocity at the free-surface, d(zeta)/d(t).
338 // Notice that barotropic mass flux divergence is not used directly.
339 //
340 int k = kk + 1;
341 W(i,j,k) = W(i,j,k-1) - (Huon(i+1,j,k-1)-Huon(i,j,k-1)) - (Hvom(i,j+1,k-1)-Hvom(i,j,k-1));
342 });
343 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
344 {
345 W(i,j,N+1)=W(i,j,N+1)/(z_w(i,j,N+1)+h(i,j,0,0)); // wrk_i
346 });
347 ParallelFor(gbx1D, N,
348 [=] AMREX_GPU_DEVICE (int i, int j, int , int kk)
349 {
350 int k = kk + 1;
351 W(i,j,k) = W(i,j,k)- W(i,j,N+1)*(z_w(i,j,k)+h(i,j,0,0));
352 });
353 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
354 {
355 W(i,j,N+1) = zero;
356 });
357 }
358
359 const int nstp = (iic) % 2;
360 nnew = 1-nstp;
361 if (solverChoice.vert_mixing_type == VertMixingType::GLS) {
362 gls_corrector(lev, vec_gls[lev].get(), vec_tke[lev].get(), mf_W, vec_Akv[lev].get(),
363 vec_Akt[lev].get(),vec_Akk[lev].get(), vec_Akp[lev].get(), vec_mskr[lev].get(),
364 vec_msku[lev].get(), vec_mskv[lev].get(),
365 nstp, nnew, N, dt_lev);
366 }
367 advance_biology(lev, *cons_old[lev], mf_cons, N, dt_lev);
368 nnew = 0;
369
370 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
371 {
372 Array4<Real> const& Hz = mf_Hz->array(mfi);
373
374 Array4<Real> const& Huon = mf_Huon->array(mfi);
375 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
376
377 Array4<Real const> const& pm = mf_pm->const_array(mfi);
378 Array4<Real const> const& pn = mf_pn->const_array(mfi);
379 Array4<Real const> const& mskr = mf_mskr->const_array(mfi);
380 Array4<Real const> const& msku = mf_msku->const_array(mfi);
381 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
382
383 Box bx = mfi.tilebox();
384 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
385 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
386
387 Box tbxp1 = bx;
388 Box tbxp11 = bx;
389 Box tbxp2 = bx;
390 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
391 tbxp2.grow(IntVect(NGROW,NGROW,0));
392 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
393
394 FArrayBox fab_FC(surroundingNodes(gbx2,2),1,amrex::The_Async_Arena());
395 FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
396 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
397
398 auto FC = fab_FC.array();
399 auto W = mf_W.array(mfi);
400
401 // Advective tracer fluxes for the coarse-fine correction, all components at once.
402 FArrayBox fab_fx_reg, fab_fy_reg, fab_fz_reg;
403 if (refluxing) {
404 fab_fx_reg.resize(surroundingNodes(bx,0), ncons, The_Async_Arena());
405 fab_fy_reg.resize(surroundingNodes(bx,1), ncons, The_Async_Arena());
406 fab_fz_reg.resize(surroundingNodes(bx,2), ncons, The_Async_Arena());
407 fab_fx_reg.template setVal<RunOn::Device>(zero);
408 fab_fy_reg.template setVal<RunOn::Device>(zero);
409 // No refinement in the vertical, so nothing crosses a z interface.
410 fab_fz_reg.template setVal<RunOn::Device>(zero);
411 }
412 //
413 //-----------------------------------------------------------------------
414 // rhs_t_3d
415 //-----------------------------------------------------------------------
416 //
417 for (int i_comp=0; i_comp < ncons; i_comp++)
418 {
419#ifdef REMORA_USE_NETCDF
420 FArrayBox* fab_river_source;
421 if (solverChoice.do_rivers_cons[i_comp]) {
422 fab_river_source = river_source_cons[i_comp]->fab_interp;
423 }
424 const Array4<const int>& river_pos = (solverChoice.do_rivers) ? vec_river_position[lev]->const_array(mfi) : Array4<const int>();
425 const Array4<const Real>& river_source = (solverChoice.do_rivers_cons[i_comp]) ? fab_river_source->const_array() : Array4<const Real>();
426#else
427 const Array4<const int >& river_pos = Array4<const int>();
428 const Array4<const Real>& river_source = Array4<const Real>();
429#endif
430 Array4<Real> const& sstore = mf_sstore->array(mfi, i_comp);
431 rhs_t_3d(lev,bx, mf_cons.array(mfi,i_comp), sstore, Huon, Hvom,
432 Hz, pn, pm, W, FC, mskr, msku, mskv, river_pos, river_source, nrhs, nnew, N,dt_lev,
433 refluxing ? fab_fx_reg.array(i_comp) : Array4<Real>(),
434 refluxing ? fab_fy_reg.array(i_comp) : Array4<Real>());
435 }
436
437 if (refluxing) {
438 std::array<FArrayBox const*, AMREX_SPACEDIM> flux
439 { &fab_fx_reg, &fab_fy_reg, &fab_fz_reg };
440 const Real* dx = geom[lev].CellSize();
441 if (lev < finest_level) {
442 getAdvFluxReg(lev+1)->CrseAdd(mfi, flux, dx, dt_lev, RunOn::Device);
443 }
444 if (lev > 0) {
445 getAdvFluxReg(lev)->FineAdd(mfi, flux, dx, dt_lev, RunOn::Device);
446 }
447 // The flux FABs are temporaries scoped to this tile, so don't move on until the
448 // register has finished reading them. ERF does the same at its own CrseAdd.
449 Gpu::streamSynchronize();
450 }
451 } // mfi
452
453 FillPatch(lev, t_old[lev], mf_cons, cons_new, BCVars::cons_bc, BdyVars::t,0,true,false,0,0,dt_lev,*cons_old[lev]);
454
455 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
456 {
457 Array4<Real> const& AK = mf_AK.array(mfi);
458 Array4<Real> const& DC = mf_DC.array(mfi);
459
460 Array4<Real> const& Hzk = mf_Hzk.array(mfi);
461 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
462
463 Box bx = mfi.tilebox();
464
465 // Copy the tilebox
466 Box tbxp1 = bx;
467 Box tbxp11 = bx;
468 Box tbxp2 = bx;
469 Box tbxp21 = bx;
470 //make only gbx be grown to match multifabs
471 tbxp21.grow(IntVect(NGROW,NGROW,NGROW-1));
472 tbxp2.grow(IntVect(NGROW,NGROW,0));
473 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
474 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
475
476 FArrayBox fab_FC(tbxp2,1,amrex::The_Async_Arena());
477 FArrayBox fab_BC(tbxp2,1,amrex::The_Async_Arena());
478 FArrayBox fab_CF(tbxp21,1,amrex::The_Async_Arena());
479 FArrayBox fab_W(tbxp2,1,amrex::The_Async_Arena());
480
481 auto FC = fab_FC.array();
482 auto BC = fab_BC.array();
483 auto CF = fab_CF.array();
484
485 for (int i_comp=0; i_comp < ncons; i_comp++) {
486 vert_visc_3d(bx,0,0,mf_cons.array(mfi,i_comp),Hz,Hzk,
487 AK,mf_Akt->array(mfi,akt_comp(i_comp)),BC,DC,FC,CF,nnew,N,dt_lev);
488 }
489 } // MFiter
490 FillPatch(lev, t_old[lev], *cons_new[lev], cons_new, BCVars::cons_bc, BdyVars::t,0,true,false,0,0,dt_lev,*cons_old[lev]);
491
492#ifdef REMORA_USE_NETCDF
493 // Nudge every tracer that was given a climatology toward it. temp and salt are just
494 // components 0 and 1 of this loop, so their behavior is unchanged.
495 for (int icomp = 0; icomp < ncons; ++icomp) {
496 if (!solverChoice.do_cons_clim_nudg[icomp]) { continue; }
497
498 cons_clim_data_from_file[icomp]->update_interpolated_to_time(model_time(t_old[lev]), lev, cons_new[lev], geom, ref_ratio);
499
500 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
501 {
502 Box bx = mfi.growntilebox(IntVect(1,1,0));
503 Array4<const Real> const& cons_nudg_coeff = vec_nudg_coeff[BdyVars::cons(icomp)][lev]->const_array(mfi);
504 Array4<const Real> const& cons_clim = cons_clim_data_from_file[icomp]->get_interpolated_mf(lev)->const_array(mfi);
505 Array4< Real> const& cons = cons_new[lev]->array(mfi, icomp);
506 Array4<const Real> const& Hz = vec_Hz[lev]->const_array(mfi);
507 Array4<const Real> const& pm = vec_pm[lev]->const_array(mfi);
508 Array4<const Real> const& pn = vec_pn[lev]->const_array(mfi);
509
510 apply_clim_nudg(bx, 0, 0, cons, cons, cons_clim, cons_nudg_coeff, Hz, pm, pn, dt_lev);
511 }
512 }
513#endif
514}
two
constexpr amrex::Real two
Definition REMORA_Constants.H:8
one
constexpr amrex::Real one
Definition REMORA_Constants.H:7
zero
constexpr amrex::Real zero
Definition REMORA_Constants.H:6
CouplingType::two_way
@ two_way
VertMixingType::GLS
@ GLS
NGROW
#define NGROW
Definition REMORA_IndexDefines.H:30
akt_comp
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int akt_comp(int icomp) noexcept
Definition REMORA_IndexDefines.H:28
setVal
mf_h setVal(geomdata.ProbHi(2))
dx
const Real dx
Definition REMORA_InitAnalyticMasks_DogboneAnalytic.H:6
REMORA::ncons
int ncons
Number of conserved scalars in the state (temperature + salt + passive scalars + biology tracers)
Definition REMORA.H:1797
REMORA::update_massflux_3d
void update_massflux_3d(int lev, const amrex::Box &bx, const int ioff, const int joff, const amrex::Array4< amrex::Real > &phi, const amrex::Array4< amrex::Real > &phibar, const amrex::Array4< amrex::Real > &Hphi, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &pm_or_pn, const amrex::Array4< amrex::Real const > &Dphi1, const amrex::Array4< amrex::Real const > &Dphi2, const amrex::Array4< amrex::Real > &DC, const amrex::Array4< amrex::Real > &FC, const amrex::Array4< amrex::Real const > &msk, const int nnew)
Correct mass flux.
Definition REMORA_update_massflux_3d.cpp:23
REMORA::xvel_bc
int xvel_bc() const noexcept
Definition REMORA.H:1448
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:603
REMORA::cons_new
amrex::Vector< amrex::MultiFab * > cons_new
multilevel data container for current step's scalar data: temperature, salinity, passive tracer
Definition REMORA.H:393
REMORA::stretch_transform
void stretch_transform(int lev)
Calculate vertical stretched coordinates.
Definition REMORA_DepthStretchTransform.H:92
REMORA::vec_mskr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr
land/sea mask at cell centers (2D)
Definition REMORA.H:584
REMORA::advance_3d
void advance_3d(int lev, amrex::MultiFab &mf_cons, amrex::MultiFab &mf_u, amrex::MultiFab &mf_v, amrex::MultiFab *mf_sstore, amrex::MultiFab *mf_ru, amrex::MultiFab *mf_rv, std::unique_ptr< amrex::MultiFab > &mf_DU_avg1, std::unique_ptr< amrex::MultiFab > &mf_DU_avg2, std::unique_ptr< amrex::MultiFab > &mf_DV_avg1, std::unique_ptr< amrex::MultiFab > &mf_DV_avg2, std::unique_ptr< amrex::MultiFab > &mf_ubar, std::unique_ptr< amrex::MultiFab > &mf_vbar, std::unique_ptr< amrex::MultiFab > &mf_Akv, std::unique_ptr< amrex::MultiFab > &mf_Akt, std::unique_ptr< amrex::MultiFab > &mf_Hz, std::unique_ptr< amrex::MultiFab > &mf_Huon, std::unique_ptr< amrex::MultiFab > &mf_Hvom, std::unique_ptr< amrex::MultiFab > &mf_z_w, amrex::MultiFab const *mf_h, amrex::MultiFab const *mf_pm, amrex::MultiFab const *mf_pn, amrex::MultiFab const *mf_mskr, amrex::MultiFab const *mf_msku, amrex::MultiFab const *mf_mskv, const int N, const amrex::Real dt_lev)
Advance the 3D variables.
Definition REMORA_advance_3d.cpp:35
REMORA::vec_tke
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_tke
Turbulent kinetic energy.
Definition REMORA.H:667
REMORA::do_substep
int do_substep
Whether to substep fine levels in time.
Definition REMORA.H:1831
REMORA::vec_gls
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_gls
Turbulent generic length scale.
Definition REMORA.H:669
REMORA::rhs_t_3d
void rhs_t_3d(int lev, const amrex::Box &bx, const amrex::Array4< amrex::Real > &t, const amrex::Array4< amrex::Real const > &tempstore, const amrex::Array4< amrex::Real const > &Huon, const amrex::Array4< amrex::Real const > &Hvom, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &pn, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &W, const amrex::Array4< amrex::Real > &FC, const amrex::Array4< amrex::Real const > &mskr, const amrex::Array4< amrex::Real const > &msku, const amrex::Array4< amrex::Real const > &mskv, const amrex::Array4< int const > &river_pos, const amrex::Array4< amrex::Real const > &river_source, int nrhs, int nnew, int N, const amrex::Real dt_lev, const amrex::Array4< amrex::Real > &fx_reg=amrex::Array4< amrex::Real >(), const amrex::Array4< amrex::Real > &fy_reg=amrex::Array4< amrex::Real >())
RHS terms for tracer.
Definition REMORA_rhs_t_3d.cpp:28
REMORA::yvel_bc
int yvel_bc() const noexcept
Definition REMORA.H:1449
REMORA::fill_from_bdyfiles
void fill_from_bdyfiles(int lev, amrex::MultiFab &mf_to_fill, const amrex::MultiFab &mf_mask, const amrex::Real time, const int bccomp, const int bdy_var_type, const int icomp_to_fill, const int icomp_calc=0, const amrex::MultiFab &mf_calc=amrex::MultiFab(), const amrex::Real=zero)
Fill boundary data from netcdf file.
Definition REMORA_BoundaryConditions_netcdf.cpp:20
REMORA::vert_mean_3d
void vert_mean_3d(const amrex::Box &bx, const int ioff, const int joff, const amrex::Array4< amrex::Real > &phi, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &Dphi_avg1, const amrex::Array4< amrex::Real > &DC, const amrex::Array4< amrex::Real > &CF, const amrex::Array4< amrex::Real const > &pm_or_pn, const amrex::Array4< amrex::Real const > &msk, const int nnew, const int N)
Adjust 3D momentum variables based on vertical mean momentum.
Definition REMORA_vert_mean_3d.cpp:19
REMORA::vec_Hz
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Hz
Width of cells in the vertical (z-) direction (3D, Hz in ROMS)
Definition REMORA.H:424
REMORA::vec_Akt
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akt
Vertical diffusion coefficient (3D)
Definition REMORA.H:444
REMORA::vec_msku
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:586
REMORA::xvel_old
amrex::Vector< amrex::MultiFab * > xvel_old
multilevel data container for last step's x velocities (u in ROMS)
Definition REMORA.H:386
REMORA::vec_mskv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:588
REMORA::physbcs
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1710
REMORA::getAdvFluxReg
amrex::YAFluxRegister * getAdvFluxReg(int lev)
flux register between lev and lev-1
Definition REMORA.H:1716
REMORA::istep
amrex::Vector< int > istep
which step?
Definition REMORA.H:1689
REMORA::yvel_old
amrex::Vector< amrex::MultiFab * > yvel_old
multilevel data container for last step's y velocities (v in ROMS)
Definition REMORA.H:388
REMORA::apply_clim_nudg
void apply_clim_nudg(const amrex::Box &bx, int ioff, int joff, const amrex::Array4< amrex::Real > &var, const amrex::Array4< amrex::Real const > &var_old, const amrex::Array4< amrex::Real const > &var_clim, const amrex::Array4< amrex::Real const > &clim_coeff, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &pn, const amrex::Real dt_lev=zero)
Apply climatology nudging.
Definition REMORA_apply_clim_nudg.cpp:19
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1949
REMORA::vec_river_position
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > vec_river_position
iMultiFab for river positions; contents are indices of rivers
Definition REMORA.H:1627
REMORA::vec_Akk
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akk
Turbulent kinetic energy vertical diffusion coefficient.
Definition REMORA.H:673
REMORA::do_reflux
int do_reflux
correct the coarse tracer with the finer level's accumulated advective flux at their interface....
Definition REMORA.H:1847
REMORA::vert_visc_3d
void vert_visc_3d(const amrex::Box &bx, const int ioff, const int joff, const amrex::Array4< amrex::Real > &phi, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real > &Hzk, const amrex::Array4< amrex::Real > &AK, const amrex::Array4< amrex::Real const > &Akv, const amrex::Array4< amrex::Real > &BC, const amrex::Array4< amrex::Real > &DC, const amrex::Array4< amrex::Real > &FC, const amrex::Array4< amrex::Real > &CF, const int nnew, const int N, const amrex::Real dt_lev)
Calculate effects of vertical viscosity or diffusivity.
Definition REMORA_vert_visc_3d.cpp:23
REMORA::advance_biology
void advance_biology(int lev, amrex::MultiFab const &mf_cons_old, amrex::MultiFab &mf_cons_new, int N, amrex::Real dt_lev)
Apply biological tracer source/sink terms.
Definition REMORA_Biology.cpp:498
REMORA::river_direction
amrex::Gpu::DeviceVector< int > river_direction
Vector over rivers of river direction: 0: u-face; 1: v-face; 2: w-face.
Definition REMORA.H:1629
REMORA::FillPatch
void FillPatch(int lev, amrex::Real time, amrex::MultiFab &mf_to_be_filled, amrex::Vector< amrex::MultiFab * > const &mfs, const int bccomp, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=true, const bool fill_set=false, const int n_not_fill=0, const int icomp_calc=0, const amrex::Real dt=zero, const amrex::MultiFab &mf_calc=amrex::MultiFab(), amrex::Vector< amrex::MultiFab * > const &mfs_crse_old={}, amrex::Vector< amrex::MultiFab * > const &mfs_crse_new={})
Fill a new MultiFab by copying in phi from valid region and filling ghost cells.
Definition REMORA_FillPatch.cpp:30
REMORA::cons_old
amrex::Vector< amrex::MultiFab * > cons_old
multilevel data container for last step's scalar data: temperature, salinity, passive tracer
Definition REMORA.H:384
REMORA::cons_clim_data_from_file
amrex::Vector< std::unique_ptr< NCTimeSeries > > cons_clim_data_from_file
Vector over cons components of climatology data read from file.
Definition REMORA.H:1613
REMORA::river_source_cons
amrex::Vector< std::unique_ptr< NCTimeSeriesRiver > > river_source_cons
Vector of data containers for scalar data in rivers.
Definition REMORA.H:1616
REMORA::river_source_transport
std::unique_ptr< NCTimeSeriesRiver > river_source_transport
Data container for momentum transport in rivers.
Definition REMORA.H:1618
REMORA::vec_nudg_coeff
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:678
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:605
REMORA::vec_Akv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akv
Vertical viscosity coefficient (3D)
Definition REMORA.H:442
REMORA::t_old
amrex::Vector< amrex::Real > t_old
old time at each level, in seconds since start_time
Definition REMORA.H:1702
REMORA::model_time
double model_time(amrex::Real elapsed) const noexcept
Time on the model clock, in seconds, of an elapsed time such as t_new.
Definition REMORA.H:2170
REMORA::vec_Akp
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akp
Turbulent length scale vertical diffusion coefficient.
Definition REMORA.H:675
REMORA::gls_corrector
void gls_corrector(int lev, amrex::MultiFab *mf_gls, amrex::MultiFab *mf_tke, amrex::MultiFab &mf_W, amrex::MultiFab *mf_Akv, amrex::MultiFab *mf_Akt, amrex::MultiFab *mf_Akk, amrex::MultiFab *mf_Akp, amrex::MultiFab *mf_mskr, amrex::MultiFab *mf_msku, amrex::MultiFab *mf_mskv, const int nstp, const int nnew, const int N, const amrex::Real dt_lev)
Corrector step for GLS calculation.
Definition REMORA_gls.cpp:247
BCVars::cons_bc
static constexpr int cons_bc
Definition REMORA_IndexDefines.H:36
BdyVars::t
static constexpr int t
cons component Temp_comp
Definition REMORA_IndexDefines.H:73
BdyVars::u
static constexpr int u
Definition REMORA_IndexDefines.H:71
BdyVars::v
static constexpr int v
Definition REMORA_IndexDefines.H:72
BdyVars::cons
int cons(int icomp) noexcept
Definition REMORA_IndexDefines.H:76
SolverChoice::do_rivers_cons
amrex::Vector< int > do_rivers_cons
Definition REMORA_DataStruct.H:916
SolverChoice::vert_mixing_type
VertMixingType vert_mixing_type
Definition REMORA_DataStruct.H:956
SolverChoice::do_cons_clim_nudg
amrex::Vector< int > do_cons_clim_nudg
Definition REMORA_DataStruct.H:1105
SolverChoice::coupling_type
CouplingType coupling_type
Definition REMORA_DataStruct.H:928