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 const BoxArray& ba = mf_cons.boxArray();
76 const DistributionMapping& dm = mf_cons.DistributionMap();
77
78 //Only used locally, probably should be rearranged into FArrayBox declaration
79 MultiFab mf_AK (convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0)); //2d missing j coordinate
80 MultiFab mf_DC (ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
81 MultiFab mf_Hzk(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
82
83 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
84 {
85 Array4<Real > const& u = mf_u.array(mfi);
86 Array4<Real > const& v = mf_v.array(mfi);
87
88 Array4<Real > const& ru = mf_ru->array(mfi);
89 Array4<Real > const& rv = mf_rv->array(mfi);
90
91 Array4<Real > const& AK = mf_AK.array(mfi);
92 Array4<Real > const& DC = mf_DC.array(mfi);
93
94 Array4<Real > const& Hzk = mf_Hzk.array(mfi);
95 Array4<Real > const& Akv = mf_Akv->array(mfi);
96
97 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
98
99 Array4<Real const> const& DU_avg1 = mf_DU_avg1->const_array(mfi);
100 Array4<Real const> const& DV_avg1 = mf_DV_avg1->const_array(mfi);
101
102 Array4<Real const> const& pm = mf_pm->const_array(mfi);
103 Array4<Real const> const& pn = mf_pn->const_array(mfi);
104
105 Array4<Real const> const& msku = mf_msku->const_array(mfi);
106 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
107
108 Box bx = mfi.tilebox();
109 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
110 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
111
112 Box xbx = mfi.nodaltilebox(0);
113 Box ybx = mfi.nodaltilebox(1);
114
115 Box gbx2D = gbx2;
116 gbx2D.makeSlab(2,0);
117
118 Box tbxp1 = bx;
119 Box tbxp11 = bx;
120 Box tbxp2 = bx;
121 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
122 tbxp2.grow(IntVect(NGROW,NGROW,0));
123 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
124
125 FArrayBox fab_FC(convert(gbx2,IntVect(0,0,1)),1,amrex::The_Async_Arena());
126 FArrayBox fab_BC(convert(gbx2,IntVect(0,0,1)),1,amrex::The_Async_Arena());
127 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
128 FArrayBox fab_W(tbxp2,1,amrex::The_Async_Arena());
129
130 auto FC = fab_FC.array();
131 auto BC = fab_BC.array();
132 auto CF = fab_CF.array();
133
134 Real cff;
135 if (iic==ntfirst) {
136 cff=Real(0.25)*dt_lev;
137 } else if (iic==ntfirst+1) {
138 cff=Real(0.25)*dt_lev*Real(3.0)/Real(2.0);
139 } else {
140 cff=Real(0.25)*dt_lev*Real(23.0)/Real(12.0);
141 }
142
143 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
144 {
145 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);
146 u(i,j,k) *= two / (Hz(i-1,j,k) + Hz(i,j,k));
147 });
148
149 ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
150 {
151 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);
152 v(i,j,k) *= two / (Hz(i,j-1,k) + Hz(i,j,k));
153 });
154
155 // NOTE: DC is only used as scratch in vert_visc_3d -- no need to pass or return a value
156 // NOTE: may not actually need to set these to zero
157
158 // Reset to zero on the box on which they'll be used
159 mf_DC[mfi].template setVal<RunOn::Device>(zero,xbx);
160 fab_CF.template setVal<RunOn::Device>(zero,xbx);
161
162 vert_visc_3d(xbx,1,0,u,Hz,Hzk,AK,Akv,BC,DC,FC,CF,nnew,N,dt_lev);
163
164 // Reset to zero on the box on which they'll be used
165 mf_DC[mfi].template setVal<RunOn::Device>(zero,ybx);
166 fab_CF.template setVal<RunOn::Device>(zero,ybx);
167
168 vert_visc_3d(ybx,0,1,v,Hz,Hzk,AK,Akv,BC,DC,FC,CF,nnew,N,dt_lev);
169
170 // Reset to zero on the box on which they'll be used
171 mf_DC[mfi].template setVal<RunOn::Device>(zero,xbx);
172 fab_CF.template setVal<RunOn::Device>(zero,xbx);
173
174 vert_mean_3d(xbx,1,0,u,Hz,DU_avg1,DC,CF,pn,msku,nnew,N);
175
176 // Reset to zero on the box on which they'll be used
177 mf_DC[mfi].template setVal<RunOn::Device>(zero,ybx);
178 fab_CF.template setVal<RunOn::Device>(zero,ybx);
179
180 vert_mean_3d(ybx,0,1,v,Hz,DV_avg1,DC,CF,pm,mskv,nnew,N);
181 }
182
183 // Apply physical boundary conditions to u and v
184 (*physbcs[lev])(mf_u,*mf_msku,0,1,mf_u.nGrowVect(),t_old[lev],xvel_bc(),0,*xvel_old[lev]);
185 (*physbcs[lev])(mf_v,*mf_mskv,0,1,mf_v.nGrowVect(),t_old[lev],yvel_bc(),0,*yvel_old[lev]);
186
187#ifdef REMORA_USE_NETCDF
188 // Fill the data which is stored in the boundary data read from netcdf files
189 if (solverChoice.boundary_from_netcdf)
190 {
191 fill_from_bdyfiles(lev, mf_u,*mf_msku,t_old[lev],xvel_bc(),BdyVars::u,0,0,*xvel_old[lev],dt_lev);
192 fill_from_bdyfiles(lev, mf_v,*mf_mskv,t_old[lev],yvel_bc(),BdyVars::v,0,0,*yvel_old[lev],dt_lev);
193 }
194
195 if (solverChoice.do_rivers) {
196 river_source_transport->update_interpolated_to_time(model_time(t_old[lev]));
197 }
198#endif
199 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
200 {
201 Array4<Real > const& u = mf_u.array(mfi);
202 Array4<Real > const& v = mf_v.array(mfi);
203
204 Array4<Real > const& DC = mf_DC.array(mfi);
205
206 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
207
208 Array4<Real const> const& DU_avg1 = mf_DU_avg1->const_array(mfi);
209 Array4<Real const> const& DV_avg1 = mf_DV_avg1->const_array(mfi);
210
211 Array4<Real const> const& DU_avg2 = mf_DU_avg2->const_array(mfi);
212 Array4<Real const> const& DV_avg2 = mf_DV_avg2->const_array(mfi);
213
214 Array4<Real> const& ubar = mf_ubar->array(mfi);
215 Array4<Real> const& vbar = mf_vbar->array(mfi);
216
217 Array4<Real> const& Huon = mf_Huon->array(mfi);
218 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
219
220 Array4<Real const> const& pm = mf_pm->const_array(mfi);
221 Array4<Real const> const& pn = mf_pn->const_array(mfi);
222 Array4<Real const> const& msku = mf_msku->const_array(mfi);
223 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
224
225 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
226
227 FArrayBox fab_FC(gbx2,1,amrex::The_Async_Arena());
228 auto FC = fab_FC.array();
229
230#ifdef REMORA_USE_NETCDF
231 if (solverChoice.do_rivers) {
232 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
233 Array4<Real const> const& z_w = mf_z_w->const_array(mfi);
234 Array4<int const> const& river_pos = vec_river_position[lev]->const_array(mfi);
235 Array4<Real const> const& river_transport = river_source_transport->fab_interp->array();
236 int* river_direction_d = river_direction.data();
237 ParallelFor(gbx1, [=] AMREX_GPU_DEVICE(int i, int j, int k)
238 {
239 int iriver = river_pos(i,j,0);
240 if (iriver >= 0) {
241 if (river_direction_d[iriver] == 0) {
242 Real on_u = two / (pn(i,j,0)+pn(i-1,j,0));
243 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)));
244 u(i,j,k) = cff * river_transport(iriver,0,k);
245 } else {
246 Real om_v = two / (pm(i,j,0)+pm(i,j-1,0));
247 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)));
248 v(i,j,k) = cff * river_transport(iriver,0,k);
249 }
250 }
251 });
252 }
253#endif
254
255#if 0
256 // Reset to zero on the box on which they'll be used
257 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(xbx,IntVect(0,0,1)));
258 fab_CF.template setVal<RunOn::Device>(zero,grow(xbx,IntVect(0,0,1)));
259
260 update_massflux_3d(xbx,1,0,u,ubar,Huon,Hz,pn,DU_avg1,DU_avg2,DC,FC,nnew);
261
262 // Reset to zero on the box on which they'll be used
263 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(ybx,IntVect(0,0,1)));
264 fab_CF.template setVal<RunOn::Device>(zero,grow(ybx,IntVect(0,0,1)));
265
266 update_massflux_3d(ybx,0,1,v,vbar,Hvom,Hz,pm,DV_avg1,DV_avg2,DC,FC,nnew);
267
268#else
269
270 // Reset to zero on the box on which they'll be used
271 fab_FC.template setVal<RunOn::Device>(zero,gbx2);
272 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(gbx2,IntVect(0,0,1)));
273 update_massflux_3d(lev,gbx2,1,0,u,ubar,Huon,Hz,pn,DU_avg1,DU_avg2,DC,FC,msku,nnew);
274
275 // Reset to zero on the box on which they'll be used
276 fab_FC.template setVal<RunOn::Device>(zero,gbx2);
277 mf_DC[mfi].template setVal<RunOn::Device>(zero,grow(gbx2,IntVect(0,0,1)));
278 update_massflux_3d(lev,gbx2,0,1,v,vbar,Hvom,Hz,pm,DV_avg1,DV_avg2,DC,FC,mskv,nnew);
279
280#endif
281 }
282
283 // WE BELIEVE THESE VALUES SHOULD ALREADY BE FILLED
284 // mf_Huon->FillBoundary(geom[lev].periodicity());
285 // mf_Hvom->FillBoundary(geom[lev].periodicity());
286
287 // ************************************************************************
288 // This should fill both temp and salt with temp/salt currently in cons_old
289 // ************************************************************************
290
291 MultiFab mf_W(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW+1,NGROW+1,0));
292 mf_W.setVal(zero);
293 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
294 {
295 Array4<Real> const& Huon = mf_Huon->array(mfi);
296 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
297
298 Array4<Real const> const& z_w = mf_z_w->const_array(mfi);
299 Array4<Real const> const& h = mf_h->const_array(mfi);
300
301 Box bx = mfi.tilebox();
302 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
303 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
304 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
305
306 Box tbxp1 = bx;
307 Box tbxp11 = bx;
308 Box tbxp2 = bx;
309 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
310 tbxp2.grow(IntVect(NGROW,NGROW,0));
311 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
312
313 FArrayBox fab_FC(surroundingNodes(gbx2,2),1,amrex::The_Async_Arena());
314 FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
315 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
316
317 auto W = mf_W.array(mfi);
318
319 //
320 //------------------------------------------------------------------------
321 // Vertically integrate horizontal mass flux divergence.
322 //------------------------------------------------------------------------
323 //
324 //Should really use gbx3uneven
325 //TODO: go over these boxes and compare to other spots where we do the same thing
326 Box gbx1D = gbx1;
327 gbx1D.makeSlab(2,0);
328
329 compute_W(gbx1D, W, Huon, Hvom, z_w, h, N);
330 }
331
332 const int nstp = (iic) % 2;
333 nnew = 1-nstp;
334 if (solverChoice.vert_mixing_type == VertMixingType::GLS) {
335 gls_corrector(lev, vec_gls[lev].get(), vec_tke[lev].get(), mf_W, vec_Akv[lev].get(),
336 vec_Akt[lev].get(),vec_Akk[lev].get(), vec_Akp[lev].get(), vec_mskr[lev].get(),
337 vec_msku[lev].get(), vec_mskv[lev].get(),
338 nstp, nnew, N, dt_lev);
339 }
340 advance_biology(lev, *cons_old[lev], mf_cons, N, dt_lev);
341 nnew = 0;
342
343 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
344 {
345 Array4<Real> const& Hz = mf_Hz->array(mfi);
346
347 Array4<Real> const& Huon = mf_Huon->array(mfi);
348 Array4<Real> const& Hvom = mf_Hvom->array(mfi);
349
350 Array4<Real const> const& pm = mf_pm->const_array(mfi);
351 Array4<Real const> const& pn = mf_pn->const_array(mfi);
352 Array4<Real const> const& mskr = mf_mskr->const_array(mfi);
353 Array4<Real const> const& msku = mf_msku->const_array(mfi);
354 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
355
356 Box bx = mfi.tilebox();
357 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
358 Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
359
360 Box tbxp1 = bx;
361 Box tbxp11 = bx;
362 Box tbxp2 = bx;
363 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
364 tbxp2.grow(IntVect(NGROW,NGROW,0));
365 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
366
367 FArrayBox fab_FC(surroundingNodes(gbx2,2),1,amrex::The_Async_Arena());
368 FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
369 FArrayBox fab_CF(gbx21,1,amrex::The_Async_Arena());
370
371 auto FC = fab_FC.array();
372 auto W = mf_W.array(mfi);
373
374 // Advective tracer fluxes for the coarse-fine correction, all components at once.
375 FArrayBox fab_fx_reg, fab_fy_reg, fab_fz_reg;
376 if (refluxing) {
377 fab_fx_reg.resize(surroundingNodes(bx,0), ncons, The_Async_Arena());
378 fab_fy_reg.resize(surroundingNodes(bx,1), ncons, The_Async_Arena());
379 fab_fz_reg.resize(surroundingNodes(bx,2), ncons, The_Async_Arena());
380 fab_fx_reg.template setVal<RunOn::Device>(zero);
381 fab_fy_reg.template setVal<RunOn::Device>(zero);
382 // No refinement in the vertical, so nothing crosses a z interface.
383 fab_fz_reg.template setVal<RunOn::Device>(zero);
384 }
385 //
386 //-----------------------------------------------------------------------
387 // rhs_t_3d
388 //-----------------------------------------------------------------------
389 //
390 for (int i_comp=0; i_comp < ncons; i_comp++)
391 {
392#ifdef REMORA_USE_NETCDF
393 FArrayBox* fab_river_source;
394 if (solverChoice.do_rivers_cons[i_comp]) {
395 fab_river_source = river_source_cons[i_comp]->fab_interp;
396 }
397 const Array4<const int>& river_pos = (solverChoice.do_rivers) ? vec_river_position[lev]->const_array(mfi) : Array4<const int>();
398 const Array4<const Real>& river_source = (solverChoice.do_rivers_cons[i_comp]) ? fab_river_source->const_array() : Array4<const Real>();
399#else
400 const Array4<const int >& river_pos = Array4<const int>();
401 const Array4<const Real>& river_source = Array4<const Real>();
402#endif
403 Array4<Real> const& sstore = mf_sstore->array(mfi, i_comp);
404 rhs_t_3d(lev,bx, mf_cons.array(mfi,i_comp), sstore, Huon, Hvom,
405 Hz, pn, pm, W, FC, mskr, msku, mskv, river_pos, river_source, nrhs, nnew, N,dt_lev,
406 refluxing ? fab_fx_reg.array(i_comp) : Array4<Real>(),
407 refluxing ? fab_fy_reg.array(i_comp) : Array4<Real>());
408 }
409
410 if (refluxing) {
411 std::array<FArrayBox const*, AMREX_SPACEDIM> flux
412 { &fab_fx_reg, &fab_fy_reg, &fab_fz_reg };
413 const Real* dx = geom[lev].CellSize();
414 if (lev < finest_level) {
415 getAdvFluxReg(lev+1)->CrseAdd(mfi, flux, dx, dt_lev, RunOn::Device);
416 }
417 if (lev > 0) {
418 getAdvFluxReg(lev)->FineAdd(mfi, flux, dx, dt_lev, RunOn::Device);
419 }
420 // The flux FABs are temporaries scoped to this tile, so don't move on until the
421 // register has finished reading them. ERF does the same at its own CrseAdd.
422 Gpu::streamSynchronize();
423 }
424 } // mfi
425
426 FillPatch(lev, t_old[lev], mf_cons, cons_new, BCVars::cons_bc, BdyVars::t,0,true,false,0,0,dt_lev,*cons_old[lev]);
427
428 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
429 {
430 Array4<Real> const& AK = mf_AK.array(mfi);
431 Array4<Real> const& DC = mf_DC.array(mfi);
432
433 Array4<Real> const& Hzk = mf_Hzk.array(mfi);
434 Array4<Real const> const& Hz = mf_Hz->const_array(mfi);
435
436 Box bx = mfi.tilebox();
437
438 // Copy the tilebox
439 Box tbxp1 = bx;
440 Box tbxp11 = bx;
441 Box tbxp2 = bx;
442 Box tbxp21 = bx;
443 //make only gbx be grown to match multifabs
444 tbxp21.grow(IntVect(NGROW,NGROW,NGROW-1));
445 tbxp2.grow(IntVect(NGROW,NGROW,0));
446 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
447 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
448
449 FArrayBox fab_FC(tbxp2,1,amrex::The_Async_Arena());
450 FArrayBox fab_BC(tbxp2,1,amrex::The_Async_Arena());
451 FArrayBox fab_CF(tbxp21,1,amrex::The_Async_Arena());
452 FArrayBox fab_W(tbxp2,1,amrex::The_Async_Arena());
453
454 auto FC = fab_FC.array();
455 auto BC = fab_BC.array();
456 auto CF = fab_CF.array();
457
458 for (int i_comp=0; i_comp < ncons; i_comp++) {
459 vert_visc_3d(bx,0,0,mf_cons.array(mfi,i_comp),Hz,Hzk,
460 AK,mf_Akt->array(mfi,akt_comp(i_comp)),BC,DC,FC,CF,nnew,N,dt_lev);
461 }
462 } // MFiter
463 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]);
464
465#ifdef REMORA_USE_NETCDF
466 // Nudge every tracer that was given a climatology toward it. temp and salt are just
467 // components 0 and 1 of this loop, so their behavior is unchanged.
468 for (int icomp = 0; icomp < ncons; ++icomp) {
469 if (!solverChoice.do_cons_clim_nudg[icomp]) { continue; }
470
471 cons_clim_data_from_file[icomp]->update_interpolated_to_time(model_time(t_old[lev]), lev, cons_new[lev], geom, ref_ratio);
472
473 for ( MFIter mfi(mf_cons, TilingIfNotGPU()); mfi.isValid(); ++mfi )
474 {
475 Box bx = mfi.growntilebox(IntVect(1,1,0));
476 Array4<const Real> const& cons_nudg_coeff = vec_nudg_coeff[BdyVars::cons(icomp)][lev]->const_array(mfi);
477 Array4<const Real> const& cons_clim = cons_clim_data_from_file[icomp]->get_interpolated_mf(lev)->const_array(mfi);
478 Array4< Real> const& cons = cons_new[lev]->array(mfi, icomp);
479 Array4<const Real> const& Hz = vec_Hz[lev]->const_array(mfi);
480 Array4<const Real> const& pm = vec_pm[lev]->const_array(mfi);
481 Array4<const Real> const& pn = vec_pn[lev]->const_array(mfi);
482
483 apply_clim_nudg(bx, 0, 0, cons, cons, cons_clim, cons_nudg_coeff, Hz, pm, pn, dt_lev);
484 }
485 }
486#endif
487}
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:1839
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:1485
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:1873
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:1486
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:1754
REMORA::getAdvFluxReg
amrex::YAFluxRegister * getAdvFluxReg(int lev)
flux register between lev and lev-1
Definition REMORA.H:1760
REMORA::istep
amrex::Vector< int > istep
which step?
Definition REMORA.H:1733
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:2017
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:1671
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:1889
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:1673
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:1657
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:1660
REMORA::river_source_transport
std::unique_ptr< NCTimeSeriesRiver > river_source_transport
Data container for momentum transport in rivers.
Definition REMORA.H:1662
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::compute_W
static void compute_W(const amrex::Box &bx2d, const amrex::Array4< amrex::Real > &W, const amrex::Array4< amrex::Real const > &Huon, const amrex::Array4< amrex::Real const > &Hvom, const amrex::Array4< amrex::Real const > &z_w, const amrex::Array4< amrex::Real const > &h, const int N)
Vertical volume flux W from the horizontal mass flux divergence.
Definition REMORA_compute_W.cpp:19
REMORA::t_old
amrex::Vector< amrex::Real > t_old
old time at each level, in seconds since start_time
Definition REMORA.H:1746
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:2210
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:72
BdyVars::u
static constexpr int u
Definition REMORA_IndexDefines.H:70
BdyVars::v
static constexpr int v
Definition REMORA_IndexDefines.H:71
BdyVars::cons
int cons(int icomp) noexcept
Definition REMORA_IndexDefines.H:75
SolverChoice::do_rivers_cons
amrex::Vector< int > do_rivers_cons
Definition REMORA_DataStruct.H:936
SolverChoice::vert_mixing_type
VertMixingType vert_mixing_type
Definition REMORA_DataStruct.H:976
SolverChoice::do_cons_clim_nudg
amrex::Vector< int > do_cons_clim_nudg
Definition REMORA_DataStruct.H:1126
SolverChoice::coupling_type
CouplingType coupling_type
Definition REMORA_DataStruct.H:948