REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_advance_2d.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2
3using namespace amrex;
4
5namespace {
6/** \brief Copy the freshly filled leapfrog component into the other two, in the coarse-fine
7 * ghost band only.
8 *
9 * FillPatch writes only component knew. The solver updates every component inside the fine
10 * grid but nothing writes the coarse-fine ghost band, so its other components stay a full
11 * parent step stale -- and the 2D solver reads zeta(krhs) and ubar(kstp) from exactly those
12 * cells. ROMS's put_refine2d sets the contact points it is about to read, so make every
13 * component agree with the parent state just filled.
14 *
15 * The band is found by marking the valid region and calling FillBoundary: ghosts another box
16 * on this level covers pick up a 1, so the remaining zeros are the coarse-fine and
17 * domain-boundary ghosts. Domain ghosts are excluded so the physical boundary conditions,
18 * applied after the FillPatch, are not overwritten.
19 */
20void fill_ghost_kcomps (MultiFab& mf, int knew, const Geometry& geom)
21{
22 MultiFab valid(mf.boxArray(), mf.DistributionMap(), 1, mf.nGrowVect());
23 valid.setVal(Real(0.0));
24 valid.setVal(Real(1.0), 0, 1, 0);
25 valid.FillBoundary(geom.periodicity());
26
27 const Box dom = amrex::convert(geom.Domain(), mf.boxArray().ixType());
28
29 const bool per_x = geom.isPeriodic(0);
30 const bool per_y = geom.isPeriodic(1);
31
32 for (MFIter mfi(mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
33 Box gbx = mfi.growntilebox();
34 const auto& a = mf.array(mfi);
35 const auto& v = valid.const_array(mfi);
36 const auto dlo = amrex::lbound(dom);
37 const auto dhi = amrex::ubound(dom);
38 ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
39 {
40 if (v(i,j,k) > Real(0.5)) { return; }
41 // Skip ghosts outside the domain only where the domain really ends. A ghost beyond
42 // a periodic edge is not a physical-boundary cell: if the level wraps onto itself
43 // the FillBoundary above already marked it valid, and if it does not -- a patch
44 // covering part of the periodic width -- it is an ordinary coarse-fine ghost that
45 // FillPatch just filled from the parent, needing the same synchronisation as any
46 // other.
47 if (!per_x && (i < dlo.x || i > dhi.x)) { return; }
48 if (!per_y && (j < dlo.y || j > dhi.y)) { return; }
49 const Real val = a(i,j,k,knew);
50 a(i,j,k,0) = val;
51 a(i,j,k,1) = val;
52 a(i,j,k,2) = val;
53 });
54 }
55}
56} // namespace
57/** Nonlinear shallow-water rpimitive equations predictor (Leap-frog) and
58 * corrector (Adams-Moulton) time-stepping engine. Corresponds to Nonlinear/step2d_LF_AM3.h
59 * in ROMS.
60 *
61 * @param[in ] lev level of refinement (coarsest level is 0)
62 * @param[in ] mf_rhoS density perturbation
63 * @param[in ] mf_rhoA vertically-averaged density
64 * @param[inout] mf_ru2d RHS contributions to 2D u-momentum
65 * @param[inout] mf_rv2d RHS contribtuions to 2D v-momentum
66 * @param[inout] mf_rufrc before first predictor, vertical integral of 3D RHS for uvel, converted to forcing terms
67 * @param[inout] mf_rvfrc before first predictor, vertical integral of 3D RHS for vvel, converted to forcing term
68 * @param[inout] mf_Zt_avg1 average of sea surface height over all fast steps
69 * @param[inout] mf_DU_avg1 time-averaged u-flux for 2D equations
70 * @param[inout] mf_DU_avg2 time-averaged u-flux for 3D equation coupling
71 * @param[inout] mf_DV_avg1 time-averaged v-flux for 2D equations
72 * @param[inout] mf_DV_avg2 time-averaged v-flux for 3D equation coupling
73 * @param[inout] mf_rubar RHS of vertically integrated u-momentum
74 * @param[inout] mf_rvbar RHS of vertically integrated v-momentum
75 * @param[inout] mf_rzeta RHS of sea surface height
76 * @param[inout] mf_ubar vertically integrated u-momentum
77 * @param[inout] mf_vbar vertically integrated v-momentum
78 * @param[inout] mf_zeta Sea-surface height
79 * @param[in ] mf_h Bathymetry
80 * @param[in ] mf_pm 1 / dx
81 * @param[in ] mf_pn 1 / dy
82 * @param[in ] mf_fcor Coriolis factor
83 * @param[inout] mf_visc2_p Harmonic viscosity at psi points
84 * @param[inout] mf_visc2_r Harmoic viscosity at rho points
85 * @param[in ] mf_mskr Land-sea mask at rho-points
86 * @param[in ] mf_msku Land-sea mask at u-points
87 * @param[in ] mf_mskv Land-sea mask at v-points
88 * @param[in ] mf_mskp Land-sea mask at psi-points
89 * @param[in ] dtfast_lev Length of current barotropic step
90 * @param[in ] predictor_2d_step Is this a predictor step?
91 * @param[in ] first_2d_step Is this the first barotropic step?
92 * @param[in ] my_iif Which barotropic predictor-corrector pair?
93 * @param[inout] next_indx1 Cached index for
94 */
95
96void
97REMORA::advance_2d (int lev,
98 MultiFab const* mf_rhoS,
99 MultiFab const* mf_rhoA,
100 MultiFab * mf_ru2d,
101 MultiFab * mf_rv2d,
102 MultiFab * mf_rufrc,
103 MultiFab * mf_rvfrc,
104 MultiFab * mf_Zt_avg1,
105 std::unique_ptr<MultiFab>& mf_DU_avg1,
106 std::unique_ptr<MultiFab>& mf_DU_avg2,
107 std::unique_ptr<MultiFab>& mf_DV_avg1,
108 std::unique_ptr<MultiFab>& mf_DV_avg2,
109 std::unique_ptr<MultiFab>& mf_rubar,
110 std::unique_ptr<MultiFab>& mf_rvbar,
111 std::unique_ptr<MultiFab>& mf_rzeta,
112 std::unique_ptr<MultiFab>& mf_ubar,
113 std::unique_ptr<MultiFab>& mf_vbar,
114 MultiFab * mf_zeta,
115 MultiFab const* mf_h,
116 MultiFab const* mf_pm,
117 MultiFab const* mf_pn,
118 MultiFab const* mf_fcor,
119 MultiFab const* mf_visc2_p,
120 MultiFab const* mf_visc2_r,
121 MultiFab const* mf_mskr,
122 MultiFab const* mf_msku,
123 MultiFab const* mf_mskv,
124 MultiFab const* mf_mskp,
125 Real dtfast_lev,
126 bool predictor_2d_step,
127 bool first_2d_step, int my_iif,
128 int & next_indx1)
129{
130 BL_PROFILE("REMORA::advance2d()");
131 int iic = istep[lev];
132 const int nnew = 0;
133 const int nstp = 0;
134 int ntfirst = 0;
135
136 int knew = 3;
137 int krhs = (my_iif + iic) % 2 + 1;
138 int kstp = my_iif <=1 ? iic % 2 + 1 : (iic % 2 + my_iif % 2 + 1) % 2 + 1;
139 int indx1 = krhs;
140 if (predictor_2d_step) {
141 next_indx1 = 3 - indx1;
142 } else {
143 knew = next_indx1;
144 kstp = 3 - knew;
145 krhs = 3;
146 //If it's not the auxiliary time step, set indx1 to next_indx1
147 // NOTE: should this ever not execute?
148 // Include indx1 updates for diagnostic purposes?
149 // if (my_iif<nfast+1)
150 // indx1=next_indx1;
151 }
152 int ptsk = 3-kstp;
153 knew-=1;
154 krhs-=1;
155 kstp-=1;
156 // Include indx1 updates for diagnostic purposes?
157 //indx1-=1;
158 ptsk-=1;
159 auto ba = mf_h->boxArray();
160 auto dm = mf_h->DistributionMap();
161
162 MultiFab mf_DUon(convert(ba,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0));
163 MultiFab mf_DVom(convert(ba,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0));
164
165 int ncomp = 0;
166 int fomn_comp = ncomp++;
167 int Drhs_comp = ncomp++;
168 int Dnew_comp = ncomp++;
169 int zwrk_comp = ncomp++;
170 int gzeta_comp = ncomp++;
171 int gzeta2_comp = ncomp++;
172 int gzetaSA_comp = ncomp++;
173 int Dstp_comp = ncomp++;
174 int rhs_ubar_comp = ncomp++;
175 int rhs_vbar_comp = ncomp++;
176 int rhs_zeta_comp = ncomp++;
177 int zeta_new_comp = ncomp++;
178
179 MultiFab mf(ba,dm,ncomp,IntVect(NGROW+1,NGROW+1,0));
180
181 for ( MFIter mfi(*mf_rhoS, TilingIfNotGPU()); mfi.isValid(); ++mfi )
182 {
183 Array4<Real > const& ubar = mf_ubar->array(mfi);
184 Array4<Real > const& vbar = mf_vbar->array(mfi);
185 Array4<Real > const& zeta = mf_zeta->array(mfi);
186 Array4<Real const> const& h = mf_h->const_array(mfi);
187
188 Array4<Real const> const& pm = mf_pm->const_array(mfi);
189 Array4<Real const> const& pn = mf_pn->const_array(mfi);
190
191 Box bx = mfi.tilebox();
192 Box gbx = mfi.growntilebox();
193 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
194 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
195 Box xgbx2 = mfi.grownnodaltilebox(0, IntVect(NGROW,NGROW,0));
196 Box ygbx2 = mfi.grownnodaltilebox(1, IntVect(NGROW,NGROW,0));
197
198 Box tbxp1 = bx;
199 Box tbxp2 = bx;
200 Box tbxp3 = bx;
201 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
202 tbxp2.grow(IntVect(NGROW,NGROW,0));
203 tbxp3.grow(IntVect(NGROW+1,NGROW+1,0));
204
205 Box bxD = bx ; bxD.makeSlab(2,0);
206 Box gbxD = gbx ; gbxD.makeSlab(2,0);
207 Box gbx1D = gbx1; gbx1D.makeSlab(2,0);
208 Box gbx2D = gbx2; gbx2D.makeSlab(2,0);
209
210 Box tbxp2D = tbxp2;
211 tbxp2D.makeSlab(2,0);
212
213 // step2d work arrays
214 FArrayBox fab_Drhs(makeSlab(tbxp3,2,0),1,The_Async_Arena());
215 auto Drhs=fab_Drhs.array();
216
217 auto DUon = mf_DUon.array(mfi);
218 auto DVom = mf_DVom.array(mfi);
219
220 ParallelFor(makeSlab(tbxp3,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
221 {
222 Drhs(i,j,0)=zeta(i,j,0,krhs)+h(i,j,0);
223 });
224
225 ParallelFor(makeSlab(xgbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
226 {
227 Real on_u = two / (pn(i,j,0)+pn(i-1,j,0));
228 Real cff1= Real(0.5) * on_u *(Drhs(i,j,0)+Drhs(i-1,j,0));
229 DUon(i,j,0)=ubar(i,j,0,krhs)*cff1;
230 });
231
232 ParallelFor(makeSlab(ygbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
233 {
234 Real om_v = two / (pm(i,j,0)+pm(i,j-1,0));
235 Real cff1= Real(0.5) * om_v * (Drhs(i,j,0)+Drhs(i,j-1,0));
236 DVom(i,j,0)=vbar(i,j,0,krhs)*cff1;
237 });
238 }
239
240 // These are needed to pass the tests with bathymetry but I don't quite see why
241 if (do_substep && cf_impose_flux) {
242 set_2d_cf_flux(lev, t_old[lev], mf_DUon, mf_DVom);
243 }
244
245 mf_DUon.FillBoundary(geom[lev].periodicity());
246 mf_DVom.FillBoundary(geom[lev].periodicity());
247
248#ifdef REMORA_USE_NETCDF
249 if (solverChoice.do_m2_clim_nudg) {
250 ubar_clim_data_from_file->update_interpolated_to_time(model_time(t_new[lev]), lev, vec_ubar[lev].get(), geom, ref_ratio);
251 vbar_clim_data_from_file->update_interpolated_to_time(model_time(t_new[lev]), lev, vec_vbar[lev].get(), geom, ref_ratio);
252 }
253#endif
254
255 for ( MFIter mfi(*mf_rhoS, TilingIfNotGPU()); mfi.isValid(); ++mfi )
256 {
257 Array4<Real const> const& rhoS = mf_rhoS->const_array(mfi);
258 Array4<Real const> const& rhoA = mf_rhoA->const_array(mfi);
259 Array4<Real const> const& h = mf_h->const_array(mfi);
260
261 Array4<Real > const& rufrc = mf_rufrc->array(mfi);
262 Array4<Real > const& rvfrc = mf_rvfrc->array(mfi);
263 Array4<Real > const& Zt_avg1 = mf_Zt_avg1->array(mfi);
264 Array4<Real > const& ubar = mf_ubar->array(mfi);
265 Array4<Real > const& vbar = mf_vbar->array(mfi);
266 Array4<Real > const& zeta = mf_zeta->array(mfi);
267 Array4<Real > const& DU_avg1 = (mf_DU_avg1)->array(mfi);
268 Array4<Real > const& DU_avg2 = (mf_DU_avg2)->array(mfi);
269 Array4<Real > const& DV_avg1 = (mf_DV_avg1)->array(mfi);
270 Array4<Real > const& DV_avg2 = (mf_DV_avg2)->array(mfi);
271 Array4<Real > const& ru2d = (mf_ru2d)->array(mfi);
272 Array4<Real > const& rv2d = (mf_rv2d)->array(mfi);
273 Array4<Real > const& rubar = (mf_rubar)->array(mfi);
274 Array4<Real > const& rvbar = (mf_rvbar)->array(mfi);
275 Array4<Real > const& rzeta = (mf_rzeta)->array(mfi);
276 Array4<Real const> const& visc2_p = mf_visc2_p->const_array(mfi);
277 Array4<Real const> const& visc2_r = mf_visc2_r->const_array(mfi);
278
279 Array4<Real const> const& pm = mf_pm->const_array(mfi);
280 Array4<Real const> const& pn = mf_pn->const_array(mfi);
281 Array4<Real const> const& fcor = mf_fcor->const_array(mfi);
282
283 Array4<Real const> const& mskr = mf_mskr->const_array(mfi);
284 Array4<Real const> const& msku = mf_msku->const_array(mfi);
285 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
286 Array4<Real const> const& mskp = mf_mskp->const_array(mfi);
287
288 Box bx = mfi.tilebox();
289 Box gbx = mfi.growntilebox();
290 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
291 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
292 Box gbx3 = mfi.growntilebox(IntVect(NGROW+1,NGROW+1,0));
293 Box xgbx2 = mfi.grownnodaltilebox(0, IntVect(NGROW,NGROW,0));
294 Box ygbx2 = mfi.grownnodaltilebox(1, IntVect(NGROW,NGROW,0));
295
296 Box xbxD = mfi.nodaltilebox(0);
297 xbxD.makeSlab(2,0);
298
299 Box ybxD = mfi.nodaltilebox(1);
300 ybxD.makeSlab(2,0);
301
302 Box tbxp1 = bx; tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
303 Box tbxp2 = bx; tbxp2.grow(IntVect(NGROW,NGROW,0));
304 Box tbxp3 = bx; tbxp3.grow(IntVect(NGROW+1,NGROW+1,0));
305
306 Box bxD = bx; bxD.makeSlab(2,0);
307 Box gbxD = gbx; gbxD.makeSlab(2,0);
308 Box gbx1D = gbx1; gbx1D.makeSlab(2,0);
309 Box gbx2D = gbx2; gbx2D.makeSlab(2,0);
310
311 Box tbxp2D = tbxp2;
312 tbxp2D.makeSlab(2,0);
313
314 auto fomn = mf.array(mfi,fomn_comp);
315 auto Drhs = mf.array(mfi,Drhs_comp);
316 auto Drhs_const = mf.const_array(mfi,Drhs_comp);
317 auto Dnew = mf.array(mfi,Dnew_comp);
318 auto zwrk = mf.array(mfi,zwrk_comp);
319 auto gzeta = mf.array(mfi,gzeta_comp);
320 auto gzeta2 = mf.array(mfi,gzeta2_comp);
321 auto gzetaSA = mf.array(mfi,gzetaSA_comp);
322 auto Dstp = mf.array(mfi,Dstp_comp);
323 auto rhs_ubar = mf.array(mfi,rhs_ubar_comp);
324 auto rhs_vbar = mf.array(mfi,rhs_vbar_comp);
325 auto rhs_zeta = mf.array(mfi,rhs_zeta_comp);
326 auto zeta_new = mf.array(mfi,zeta_new_comp);
327
328 FArrayBox & fab_DUon=mf_DUon[mfi];
329 FArrayBox & fab_DVom=mf_DVom[mfi];
330 auto DUon=fab_DUon.array();
331 auto DVom=fab_DVom.array();
332
333 auto weight1 = vec_weight1.dataPtr();
334 auto weight2 = vec_weight2.dataPtr();
335
336 //From ana_grid.h and metrics.F
337 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int)
338 {
339 rhs_ubar(i,j,0)=zero;
340 });
341
342 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int)
343 {
344 rhs_vbar(i,j,0)=zero;
345 });
346
347 if (solverChoice.use_coriolis) {
348 ParallelFor(tbxp2D, [=] AMREX_GPU_DEVICE (int i, int j, int )
349 {
350 fomn(i,j,0) = fcor(i,j,0)*(one/(pm(i,j,0)*pn(i,j,0)));
351 });
352 }
353
354 ParallelFor(makeSlab(tbxp3,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
355 {
356 Drhs(i,j,0)=zeta(i,j,0,krhs)+h(i,j,0);
357 });
358
359 if(predictor_2d_step)
360 {
361 if(first_2d_step) {
362 Real cff2=(Real(-1.0)/Real(12.0))*weight2[my_iif+1];
363 ParallelFor(makeSlab(gbx3,2,0),
364 [=] AMREX_GPU_DEVICE (int i, int j, int)
365 {
366 Zt_avg1(i,j,0)=zero;
367 });
368 ParallelFor(makeSlab(xgbx2,2,0),
369 [=] AMREX_GPU_DEVICE (int i, int j, int)
370 {
371 DU_avg1(i,j,0)=zero;
372 DU_avg2(i,j,0)=cff2*DUon(i,j,0);
373 });
374 ParallelFor(makeSlab(ygbx2,2,0),
375 [=] AMREX_GPU_DEVICE (int i, int j, int)
376 {
377 DV_avg1(i,j,0)=zero;
378 DV_avg2(i,j,0)=cff2*DVom(i,j,0);
379 });
380 }
381 else {
382 Real cff1_wt1 = weight1[my_iif-1];
383 Real cff2_wt1 = (Real(8.0)/Real(12.0))*weight2[my_iif]-
384 (one/Real(12.0))*weight2[my_iif+1];
385
386 ParallelFor(makeSlab(gbx3,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
387 {
388 Zt_avg1(i,j,0) += cff1_wt1*zeta(i,j,0,krhs);
389 });
390
391 ParallelFor(makeSlab(xgbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
392 {
393 DU_avg1(i,j,0) += cff1_wt1*DUon(i,j,0);
394 DU_avg2(i,j,0) += cff2_wt1*DUon(i,j,0);
395 });
396
397 ParallelFor(makeSlab(ygbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
398 {
399 DV_avg1(i,j,0) += cff1_wt1*DVom(i,j,0);
400 DV_avg2(i,j,0) += cff2_wt1*DVom(i,j,0);
401 });
402 }
403 }
404 else {
405 Real cff2_wt2;
406
407 if (first_2d_step) {
408 cff2_wt2=weight2[my_iif];
409 } else {
410 cff2_wt2=Real(5.0)/Real(12.0)*weight2[my_iif];
411 }
412
413 ParallelFor(makeSlab(xgbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
414 {
415 DU_avg2(i,j,0)=DU_avg2(i,j,0)+cff2_wt2*DUon(i,j,0);
416 });
417
418 ParallelFor(makeSlab(ygbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
419 {
420 DV_avg2(i,j,0)=DV_avg2(i,j,0)+cff2_wt2*DVom(i,j,0);
421 });
422 }
423 //
424 // Do not perform the actual time stepping during the auxiliary
425 // (nfast(ng)+1) time step. Jump to next box
426 //
427
428 if (my_iif>=nfast) {
429 continue; }
430 //Load new free-surface values into shared array at both predictor
431 //and corrector steps
432 //
433 //=======================================================================
434 // Time step free-surface equation.
435 //=======================================================================
436 //
437 // During the first time-step, the predictor step is Forward-Euler
438 // and the corrector step is Backward-Euler. Otherwise, the predictor
439 // step is Leap-frog and the corrector step is Adams-Moulton.
440 //
441
442 Real fac=Real(1000.0)/solverChoice.rho0;
443
444 if (my_iif==0) {
445 Real cff1=dtfast_lev;
446
447 ParallelFor(makeSlab(tbxp1,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
448 {
449 rhs_zeta(i,j,0) = (DUon(i,j,0)-DUon(i+1,j,0))+
450 (DVom(i,j,0)-DVom(i,j+1,0));
451 zeta_new(i,j,0) = (zeta(i,j,0,kstp)+ pm(i,j,0)*pn(i,j,0)*cff1*rhs_zeta(i,j,0)) * mskr(i,j,0);
452 Dnew(i,j,0) = zeta_new(i,j,0)+h(i,j,0);
453
454 //Pressure gradient terms:
455 zwrk(i,j,0)=Real(0.5)*(zeta(i,j,0,kstp)+zeta_new(i,j,0));
456 gzeta(i,j,0)=(fac+rhoS(i,j,0))*zwrk(i,j,0);
457 gzeta2(i,j,0)=gzeta(i,j,0)*zwrk(i,j,0);
458 gzetaSA(i,j,0)=zwrk(i,j,0)*(rhoS(i,j,0)-rhoA(i,j,0));
459 });
460
461 } else if (predictor_2d_step) {
462
463 Real cff1=two * dtfast_lev;
464 Real cff4=Real(4.0) / Real(25.0);
465 Real cff5=one - two*cff4;
466
467 ParallelFor(makeSlab(tbxp1,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
468 {
469 rhs_zeta(i,j,0)=(DUon(i,j,0)-DUon(i+1,j,0))+
470 (DVom(i,j,0)-DVom(i,j+1,0));
471 zeta_new(i,j,0)=(zeta(i,j,0,kstp)+
472 pm(i,j,0)*pn(i,j,0)*cff1*rhs_zeta(i,j,0)) * mskr(i,j,0);
473 Dnew(i,j,0)=zeta_new(i,j,0)+h(i,j,0);
474 //Pressure gradient terms
475 zwrk(i,j,0)=cff5*zeta(i,j,0,krhs)+
476 cff4*(zeta(i,j,0,kstp)+zeta_new(i,j,0));
477 gzeta(i,j,0)=(fac+rhoS(i,j,0))*zwrk(i,j,0);
478 gzeta2(i,j,0)=gzeta(i,j,0)*zwrk(i,j,0);
479 gzetaSA(i,j,0)=zwrk(i,j,0)*(rhoS(i,j,0)-rhoA(i,j,0));
480 });
481
482 } else if (!predictor_2d_step) { //AKA if(corrector_2d_step)
483
484 Real cff1=dtfast_lev * Real(5.0)/Real(12.0);
485 Real cff2=dtfast_lev * Real(8.0)/Real(12.0);
486 Real cff3=dtfast_lev * one/Real(12.0);
487 Real cff4=two/Real(5.0);
488 Real cff5=one-cff4;
489
490 ParallelFor(makeSlab(tbxp1,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
491 {
492 Real cff=cff1*((DUon(i,j,0)-DUon(i+1,j,0))+
493 (DVom(i,j,0)-DVom(i,j+1,0)));
494 zeta_new(i,j,0)=zeta(i,j,0,kstp)+
495 pm(i,j,0)*pn(i,j,0)*(cff+
496 cff2*rzeta(i,j,0,kstp)-
497 cff3*rzeta(i,j,0,ptsk));
498 zeta_new(i,j,0) *= mskr(i,j,0);
499 Dnew(i,j,0)=zeta_new(i,j,0)+h(i,j,0);
500 //Pressure gradient terms
501 zwrk(i,j,0)=cff5*zeta_new(i,j,0)+cff4*zeta(i,j,0,krhs);
502 gzeta(i,j,0)=(fac+rhoS(i,j,0))*zwrk(i,j,0);
503 gzeta2(i,j,0)=gzeta(i,j,0)*zwrk(i,j,0);
504 gzetaSA(i,j,0)=zwrk(i,j,0)*(rhoS(i,j,0)-rhoA(i,j,0));
505 });
506 }
507
508 //
509 // Load new free-surface values into shared array at both predictor
510 // and corrector steps.
511 //
512 //// zeta(knew) only valid at zeta_new, i.e. tbxp1
513 ParallelFor(makeSlab(gbx1,2,0),
514 [=] AMREX_GPU_DEVICE (int i, int j, int )
515 {
516 zeta(i,j,0,knew) = zeta_new(i,j,0);
517 });
518
519 //
520 // If predictor step, load right-side-term into shared array.
521 //
522 if (predictor_2d_step) {
523 ParallelFor(makeSlab(gbx1,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
524 {
525 rzeta(i,j,0,krhs)=rhs_zeta(i,j,0);
526 });
527 }
528
529 //
530 //=======================================================================
531 // Compute right-hand-side for the 2D momentum equations.
532 //=======================================================================
533 //
534/*
535!
536!-----------------------------------------------------------------------
537! Compute pressure gradient terms.
538!-----------------------------------------------------------------------
539!
540*/
541 Real cff1 = Real(0.5) * solverChoice.g;
542 Real cff2 = one / Real(3.0);
543 ParallelFor(xbxD,
544 [=] AMREX_GPU_DEVICE (int i, int j, int )
545 {
546 Real on_u = two / (pn(i,j,0)+pn(i-1,j,0));
547 rhs_ubar(i,j,0)=cff1 * on_u *
548 (( h(i-1,j,0) + h(i,j,0))*
549 (gzeta(i-1,j,0) - gzeta(i,j,0))+
550 ( h(i-1,j,0) - h(i,j,0))*
551 ( gzetaSA(i-1,j,0) + gzetaSA(i,j,0)+
552 cff2*( rhoA(i-1,j,0) - rhoA(i,j,0))*
553 ( zwrk(i-1,j,0) - zwrk(i,j,0)))+
554 (gzeta2(i-1,j,0)- gzeta2(i ,j,0)));
555 });
556
557 ParallelFor(ybxD,
558 [=] AMREX_GPU_DEVICE (int i, int j, int )
559 {
560 Real om_v = two / (pm(i,j,0)+pm(i,j-1,0));
561 rhs_vbar(i,j,0) = cff1*om_v *
562 (( h(i,j-1,0) + h(i,j,0))*
563 (gzeta(i,j-1,0) - gzeta(i,j,0))+
564 ( h(i,j-1,0) - h(i,j,0))*
565 (gzetaSA(i,j-1,0)+ gzetaSA(i,j ,0)+
566 cff2*(rhoA(i,j-1,0)- rhoA(i,j ,0))*
567 (zwrk(i,j-1,0)- zwrk(i,j ,0)))+
568 (gzeta2(i,j-1,0)- gzeta2(i,j ,0)));
569 });
570
571 // Advection terms for 2d ubar, vbar added to rhs_ubar and rhs_vbar
572 //
573 //-----------------------------------------------------------------------
574 // rhs_uv_2d
575 //-----------------------------------------------------------------------
576 //
577 Array4<Real const> const& ubar_const = mf_ubar->const_array(mfi);
578 Array4<Real const> const& vbar_const = mf_vbar->const_array(mfi);
579
580 rhs_uv_2d(lev,xbxD, ybxD, ubar_const, vbar_const, rhs_ubar, rhs_vbar, DUon, DVom, krhs);
581
582 //-----------------------------------------------------------------------
583 // Add Coriolis forcing
584 //-----------------------------------------------------------------------
585 if (solverChoice.use_coriolis) {
586 // Coriolis terms for 2d ubar, vbar added to rhs_ubar and rhs_vbar
587 //
588 //-----------------------------------------------------------------------
589 // coriolis
590 //-----------------------------------------------------------------------
591 //
592 coriolis(xbxD, ybxD, ubar_const, vbar_const, rhs_ubar, rhs_vbar, Drhs, fomn, krhs, 0);
593 }
594
595 if (solverChoice.use_curvilinear_grid) {
596 Array4<Real const> const& dndx = vec_dndx[lev]->const_array(mfi);
597 Array4<Real const> const& dmde = vec_dmde[lev]->const_array(mfi);
598 curvilinear(bxD, xbxD, ybxD, ubar_const, vbar_const, rhs_ubar, rhs_vbar, Drhs, dndx, dmde, krhs, 0);
599 }
600
601 //-----------------------------------------------------------------------
602 //Add in horizontal harmonic viscosity.
603 // Consider generalizing or copying uv3dmix, where Drhs is used instead of Hz and u=>ubar v=>vbar, drop dt terms
604 //-----------------------------------------------------------------------
605 uv3dmix(xbxD, ybxD, ubar, vbar, ubar, vbar, rhs_ubar, rhs_vbar,
606 visc2_p, visc2_r, Drhs_const,
607 pm, pn, mskp, krhs, nnew, zero);
608
609#ifdef REMORA_USE_NETCDF
610 if (solverChoice.do_m2_clim_nudg) {
611 Array4<Real > const& ubar_krhs = mf_ubar->array(mfi, krhs);
612 Array4<Real > const& vbar_krhs = mf_vbar->array(mfi, krhs);
613 Array4<const Real> const& ubar_clim = ubar_clim_data_from_file->get_interpolated_mf(lev)->const_array(mfi);
614 Array4<const Real> const& vbar_clim = vbar_clim_data_from_file->get_interpolated_mf(lev)->const_array(mfi);
615 Array4<const Real> const& ubar_nudg_coeff = vec_nudg_coeff[bdy_ubar()][lev]->const_array(mfi);
616 Array4<const Real> const& vbar_nudg_coeff = vec_nudg_coeff[bdy_vbar()][lev]->const_array(mfi);
617 // Boxes are like this to match the ROMS loops i=IstrU..Iend, j=JstrV..Jend, which
618 // exclude the u-face on the west/east domain edges and the v-face on south/north
619 Box xbxD_adj = clim_nudg_momentum_box(mfi.nodaltilebox(0), 0,
620 Geom(lev).Domain(), Geom(lev).isPeriodic(0));
621 xbxD_adj.makeSlab(2,0);
622 Box ybxD_adj = clim_nudg_momentum_box(mfi.nodaltilebox(1), 1,
623 Geom(lev).Domain(), Geom(lev).isPeriodic(1));
624 ybxD_adj.makeSlab(2,0);
625 apply_clim_nudg(xbxD_adj, 1, 0, rhs_ubar, ubar_krhs, ubar_clim, ubar_nudg_coeff, Drhs_const, pm, pn);
626 apply_clim_nudg(ybxD_adj, 0, 1, rhs_vbar, vbar_krhs, vbar_clim, vbar_nudg_coeff, Drhs_const, pm, pn);
627 }
628#endif
629
630 //-----------------------------------------------------------------------
631 // Coupling from 3d to 2d
632 //-----------------------------------------------------------------------
633 if (first_2d_step&&predictor_2d_step)
634 {
635 if (iic==ntfirst) {
636 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
637 {
638 rufrc(i,j,0) -= rhs_ubar(i,j,0);
639 rhs_ubar(i,j,0) += rufrc(i,j,0);
640 ru2d(i,j,0,nstp) = rufrc(i,j,0);
641 });
642
643 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
644 {
645 rvfrc(i,j,0) -= rhs_vbar(i,j,0);
646 rhs_vbar(i,j,0) += rvfrc(i,j,0);
647 rv2d(i,j,0,nstp) = rvfrc(i,j,0);
648 });
649
650 } else if (iic==(ntfirst+1)) {
651
652 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
653 {
654 rufrc(i,j,0)=rufrc(i,j,0)-rhs_ubar(i,j,0);
655 rhs_ubar(i,j,0)=rhs_ubar(i,j,0)+Real(1.5)*rufrc(i,j,0)-Real(0.5)*ru2d(i,j,0,0);
656 ru2d(i,j,0,1)=rufrc(i,j,0);
657 Real r_swap= ru2d(i,j,0,1);
658 ru2d(i,j,0,1) = ru2d(i,j,0,0);
659 ru2d(i,j,0,0) = r_swap;
660 });
661
662 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
663 {
664 rvfrc(i,j,0)=rvfrc(i,j,0)-rhs_vbar(i,j,0);
665 rhs_vbar(i,j,0)=rhs_vbar(i,j,0)+Real(1.5)*rvfrc(i,j,0)-Real(0.5)*rv2d(i,j,0,0);
666 rv2d(i,j,0,1)=rvfrc(i,j,0);
667 Real r_swap= rv2d(i,j,0,1);
668 rv2d(i,j,0,1) = rv2d(i,j,0,0);
669 rv2d(i,j,0,0) = r_swap;
670 });
671
672 } else {
673 cff1=Real(23.0)/Real(12.0);
674 cff2=Real(16.0)/Real(12.0);
675 Real cff3= Real(5.0)/Real(12.0);
676
677 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
678 {
679 rufrc(i,j,0)=rufrc(i,j,0)-rhs_ubar(i,j,0);
680 rhs_ubar(i,j,0)=rhs_ubar(i,j,0)+
681 cff1*rufrc(i,j,0)-
682 cff2*ru2d(i,j,0,0)+
683 cff3*ru2d(i,j,0,1);
684 ru2d(i,j,0,1)=rufrc(i,j,0);
685 Real r_swap= ru2d(i,j,0,1);
686 ru2d(i,j,0,1) = ru2d(i,j,0,0);
687 ru2d(i,j,0,0) = r_swap;
688 });
689
690 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
691 {
692 rvfrc(i,j,0)=rvfrc(i,j,0)-rhs_vbar(i,j,0);
693 rhs_vbar(i,j,0)=rhs_vbar(i,j,0)+
694 cff1*rvfrc(i,j,0)-
695 cff2*rv2d(i,j,0,0)+
696 cff3*rv2d(i,j,0,1);
697 rv2d(i,j,0,1)=rvfrc(i,j,0);
698
699 Real r_swap= rv2d(i,j,0,1);
700 rv2d(i,j,0,1) = rv2d(i,j,0,0);
701 rv2d(i,j,0,0) = r_swap;
702 });
703 }
704 } else {
705 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
706 {
707 rhs_ubar(i,j,0) += rufrc(i,j,0);
708 });
709
710 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
711 {
712 rhs_vbar(i,j,0) += rvfrc(i,j,0);
713 });
714 }
715
716 //
717 //=======================================================================
718 // Time step 2D momentum equations.
719 //=======================================================================
720 //
721 // Compute total water column depth.
722 //
723 ParallelFor(makeSlab(tbxp3,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
724 {
725 Dstp(i,j,0)=zeta(i,j,0,kstp)+h(i,j,0);
726 });
727
728 //
729 // During the first time-step, the predictor step is Forward-Euler
730 // and the corrector step is Backward-Euler. Otherwise, the predictor
731 // step is Leap-frog and the corrector step is Adams-Moulton.
732 //
733 if (my_iif==0) {
734 cff1=Real(0.5)*dtfast_lev;
735 ParallelFor(xbxD,
736 [=] AMREX_GPU_DEVICE (int i, int j, int )
737 {
738 Real cff=(pm(i,j,0)+pm(i-1,j,0))*(pn(i,j,0)+pn(i-1,j,0));
739 Real Dnew_avg =one/(Dnew(i,j,0)+Dnew(i-1,j,0));
740 ubar(i,j,0,knew)=(ubar(i,j,0,kstp)*
741 (Dstp(i,j,0)+Dstp(i-1,j,0))+
742 cff*cff1*rhs_ubar(i,j,0))*Dnew_avg * msku(i,j,0);
743 });
744 ParallelFor(ybxD,
745 [=] AMREX_GPU_DEVICE (int i, int j, int )
746 {
747 Real cff=(pm(i,j,0)+pm(i,j-1,0))*(pn(i,j,0)+pn(i,j-1,0));
748 Real Dnew_avg=one/(Dnew(i,j,0)+Dnew(i,j-1,0));
749 vbar(i,j,0,knew)=(vbar(i,j,0,kstp)*
750 (Dstp(i,j,0)+Dstp(i,j-1,0))+
751 cff*cff1*rhs_vbar(i,j,0))*Dnew_avg * mskv(i,j,0);
752 });
753
754 } else if (predictor_2d_step) {
755
756 cff1=dtfast_lev;
757 ParallelFor(xbxD,
758 [=] AMREX_GPU_DEVICE (int i, int j, int )
759 {
760 Real cff=(pm(i,j,0)+pm(i-1,j,0))*(pn(i,j,0)+pn(i-1,j,0));
761 Real Dnew_avg=one/(Dnew(i,j,0)+Dnew(i-1,j,0));
762 ubar(i,j,0,knew)=(ubar(i,j,0,kstp)*
763 (Dstp(i,j,0)+Dstp(i-1,j,0))+
764 cff*cff1*rhs_ubar(i,j,0))*Dnew_avg * msku(i,j,0);
765 });
766 ParallelFor(ybxD,
767 [=] AMREX_GPU_DEVICE (int i, int j, int )
768 {
769 Real cff=(pm(i,j,0)+pm(i,j-1,0))*(pn(i,j,0)+pn(i,j-1,0));
770 Real Dnew_avg=one/(Dnew(i,j,0)+Dnew(i,j-1,0));
771 vbar(i,j,0,knew)=(vbar(i,j,0,kstp)*
772 (Dstp(i,j,0)+Dstp(i,j-1,0))+
773 cff*cff1*rhs_vbar(i,j,0))*Dnew_avg * mskv(i,j,0);
774 });
775
776 } else if ((!predictor_2d_step)) {
777
778 cff1=Real(0.5)*dtfast_lev*Real(5.0)/Real(12.0);
779 cff2=Real(0.5)*dtfast_lev*Real(8.0)/Real(12.0);
780 Real cff3=Real(0.5)*dtfast_lev*one/Real(12.0);
781 ParallelFor(xbxD,
782 [=] AMREX_GPU_DEVICE (int i, int j, int )
783 {
784 Real cff=(pm(i,j,0)+pm(i-1,j,0))*(pn(i,j,0)+pn(i-1,j,0));
785 Real Dnew_avg=one/(Dnew(i,j,0)+Dnew(i-1,j,0));
786 ubar(i,j,0,knew)=(ubar(i,j,0,kstp)*
787 (Dstp(i,j,0)+Dstp(i-1,j,0))+
788 cff*(cff1*rhs_ubar(i,j,0)+
789 cff2*rubar(i,j,0,kstp)-
790 cff3*rubar(i,j,0,ptsk)))*Dnew_avg * msku(i,j,0);
791 });
792 ParallelFor(ybxD,
793 [=] AMREX_GPU_DEVICE (int i, int j, int )
794 {
795 Real cff=(pm(i,j,0)+pm(i,j-1,0))*(pn(i,j,0)+pn(i,j-1,0));
796 Real Dnew_avg=one/(Dnew(i,j,0)+Dnew(i,j-1,0));
797 vbar(i,j,0,knew)=(vbar(i,j,0,kstp)*
798 (Dstp(i,j,0)+Dstp(i,j-1,0))+
799 cff*(cff1*rhs_vbar(i,j,0)+
800 cff2*rvbar(i,j,0,kstp)-
801 cff3*rvbar(i,j,0,ptsk)))*Dnew_avg * mskv(i,j,0);
802 });
803 }
804
805 //store rhs_ubar and rhs_vbar to save later
806 //
807 // If predictor step, load right-side-term into shared arrays for
808 // future use during the subsequent corrector step.
809 //
810
811 if (predictor_2d_step) {
812 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
813 {
814 rubar(i,j,0,krhs)=rhs_ubar(i,j,0);
815 });
816 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
817 {
818 rvbar(i,j,0,krhs)=rhs_vbar(i,j,0);
819 });
820 }
821 }
822
823 // Don't do the FillPatch or rivers at the last truncated predictor step.
824 // We may need to move the zeta FillPatch further up
825 if (my_iif<nfast) {
826 int know;
827 Real dt2d;
828 if (my_iif==0) {
829 know = krhs;
830 dt2d = dtfast_lev;
831 } else if (predictor_2d_step) {
832 know = krhs;
833 dt2d = two * dtfast_lev;
834 } else {
835 know = kstp;
836 dt2d = dtfast_lev;
837 }
838
839 MultiFab ubar_know(*vec_ubar[lev], make_alias, know, 1);
840 MultiFab vbar_know(*vec_vbar[lev], make_alias, know, 1);
841 MultiFab zeta_know(*vec_zeta[lev], make_alias, know, 1);
842
843 FillPatch(lev, t_old[lev], *vec_ubar[lev], GetVecOfPtrs(vec_ubar), ubar_bc(), bdy_ubar(),
844 knew, false,true, 0,know, dt2d, ubar_know);
845 FillPatch(lev, t_old[lev], *vec_vbar[lev], GetVecOfPtrs(vec_vbar), vbar_bc(), bdy_vbar(),
846 knew, false,true, 0,know, dt2d, vbar_know);
847 // With the snapshot pair the coarse contribution is interpolated onto t_old[lev],
848 // which timeStep guarantees lies inside the parent's step; without it FillPatch
849 // falls back to the parent frozen at the end of that step. See roll_2d_snapshot.
850 if (cf_time_interp_zeta && lev > 0 && int(vec_zeta_crse_old.size()) >= lev
851 && vec_zeta_crse_old[lev-1] && vec_zeta_crse_new[lev-1]) {
852 FillPatch(lev, t_old[lev], *vec_zeta[lev], GetVecOfPtrs(vec_zeta), zeta_bc(), bdy_zeta(),
853 knew, false,false, 0,know, dt2d, zeta_know,
854 GetVecOfPtrs(vec_zeta_crse_old), GetVecOfPtrs(vec_zeta_crse_new));
855 } else {
856 FillPatch(lev, t_old[lev], *vec_zeta[lev], GetVecOfPtrs(vec_zeta), zeta_bc(), bdy_zeta(),
857 knew, false,false, 0,know, dt2d, zeta_know);
858 }
859
860 // See fill_ghost_kcomps: without this the components the solver reads in the ghost
861 // band lag the one FillPatch just wrote by a whole parent step.
862 if (lev > 0 && cf_fill_all_kcomp) {
863 fill_ghost_kcomps(*vec_zeta[lev], knew, geom[lev]);
864 fill_ghost_kcomps(*vec_ubar[lev], knew, geom[lev]);
865 fill_ghost_kcomps(*vec_vbar[lev], knew, geom[lev]);
866 }
867
868 // Replace the interface faces the FillPatchers just set from the parent's ubar with
869 // the parent's mass flux, which conserves mass. Must follow the FillPatch.
870 //
871 // cf_set_2d_bcs = 2 imposes it only on the first fast step, where ROMS calls
872 // put_refine2d (main3d.F, before the barotropic loop; inside it only composite grids
873 // and the NESTING_DEBUG check run). Every fast step instead re-applies the parent's
874 // step-averaged transport twenty times per step rather than letting the barotropic
875 // solver carry the interface.
876 const bool set_2d_now = (cf_set_2d_bcs == 1) || (cf_set_2d_bcs == 2 && my_iif == 0);
877 if (do_substep && set_2d_now) {
878 set_2d_cf_bcs(lev, t_old[lev], know, knew);
879 // set_2d_cf_bcs writes each box's own interface faces. Where a box boundary meets
880 // the interface, the neighbouring box's ghost copy of that face still holds the
881 // FillPatch value, and the shared face between the two boxes would be advanced
882 // from different neighbours -- the answer then depends on the box layout.
883 vec_ubar[lev]->FillBoundary(knew, 1, geom[lev].periodicity());
884 vec_vbar[lev]->FillBoundary(knew, 1, geom[lev].periodicity());
885 }
886
887#ifdef REMORA_USE_NETCDF
888 if (solverChoice.do_rivers) {
889 river_source_transportbar->update_interpolated_to_time(model_time(t_old[lev]));
890 int* river_direction_d = river_direction.data();
891 for ( MFIter mfi(*mf_rhoS, TilingIfNotGPU()); mfi.isValid(); ++mfi )
892 {
893 Array4<const int > const& river_pos = vec_river_position[lev]->const_array(mfi);
894 Array4<const Real> const& river_transportbar = river_source_transportbar->fab_interp->array();
895 Array4<Real > const& ubar = mf_ubar->array(mfi);
896 Array4<Real > const& vbar = mf_vbar->array(mfi);
897 Array4<Real const> const& zeta = mf_zeta->const_array(mfi);
898 Array4<Real const> const& h = mf_h->const_array(mfi);
899 Array4<Real const> const& pm = mf_pm->const_array(mfi);
900 Array4<Real const> const& pn = mf_pn->const_array(mfi);
901
902 Box gbx1D = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
903 gbx1D.makeSlab(2,0);
904
905 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
906 {
907 int iriver = river_pos(i,j,0);
908 if (iriver >= 0) {
909 if (river_direction_d[iriver] == 0) {
910 Real on_u = two / (pn(i,j,0)+pn(i-1,j,0));
911 Real cff = one / (on_u * Real(0.5) * (zeta(i-1,j,0,knew) + h(i-1,j,0) +
912 zeta(i,j,0,knew) + h(i,j,0)));
913 ubar(i,j,0,knew) = river_transportbar(iriver,0,0) * cff;
914 } else {
915 Real om_v = two / (pm(i,j,0)+pm(i,j-1,0));
916 Real cff = one / (om_v * Real(0.5) * (zeta(i,j-1,0,knew) + h(i,j-1,0) +
917 zeta(i,j,0,knew) + h(i,j,0)));
918 vbar(i,j,0,knew) = river_transportbar(iriver,0,0) * cff;
919 }
920 }
921 });
922 }
923 }
924 FillPatchNoBC(lev, t_old[lev], *vec_ubar[lev], GetVecOfPtrs(vec_ubar), bdy_ubar(),
925 knew, false,false);
926 FillPatchNoBC(lev, t_old[lev], *vec_vbar[lev], GetVecOfPtrs(vec_vbar), bdy_vbar(),
927 knew, false,false);
928#endif
929 }
930}
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
NGROW
#define NGROW
Definition REMORA_IndexDefines.H:30
setVal
mf_h setVal(geomdata.ProbHi(2))
if
if(traditional)
Definition REMORA_InitAnalyticBathymetry_DogboneAnalytic.H:11
cff1
Real cff1
Definition REMORA_InitAnalyticBiology_BioToy.H:2
cff2
Real cff2
Definition REMORA_InitAnalyticBiology_BioToy.H:3
mf_fcor
Real mf_fcor
Definition REMORA_InitAnalyticSMFlux_ChannelTest.H:12
REMORA::nfast
int nfast
Number of fast steps to take.
Definition REMORA.H:1828
REMORA::vec_zeta_crse_new
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_zeta_crse_new
see vec_zeta_crse_old (2D)
Definition REMORA.H:566
REMORA::vec_zeta_crse_old
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_zeta_crse_old
free surface at the start and end of this level's step, every leapfrog component holding the same val...
Definition REMORA.H:564
REMORA::do_substep
int do_substep
Whether to substep fine levels in time.
Definition REMORA.H:1831
REMORA::set_2d_cf_bcs
void set_2d_cf_bcs(int lev, amrex::Real time, int know, int knew)
set this level's barotropic contact points from its parent, mass-conservingly
Definition REMORA_set_2d_cf_bcs.cpp:150
REMORA::clim_nudg_momentum_box
static amrex::Box clim_nudg_momentum_box(const amrex::Box &nodal_bx, int dir, const amrex::Box &domain, bool is_periodic)
Shrink an x- or y-nodal momentum tilebox to the faces ROMS nudges.
Definition REMORA_apply_clim_nudg.cpp:52
REMORA::river_source_transportbar
std::unique_ptr< NCTimeSeriesRiver > river_source_transportbar
Data container for vertically integrated momentum transport in rivers.
Definition REMORA.H:1620
REMORA::zeta_bc
int zeta_bc() const noexcept
Definition REMORA.H:1453
REMORA::vec_weight2
amrex::Vector< amrex::Real > vec_weight2
Weights for calculating avg2 in 2D advance.
Definition REMORA.H:663
REMORA::vec_dmde
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_dmde
d(1/m)/d(eta)
Definition REMORA.H:644
REMORA::bdy_zeta
int bdy_zeta() const noexcept
Definition REMORA.H:1466
REMORA::istep
amrex::Vector< int > istep
which step?
Definition REMORA.H:1689
REMORA::uv3dmix
void uv3dmix(const amrex::Box &xbx, const amrex::Box &ybx, const amrex::Array4< amrex::Real > &u, const amrex::Array4< amrex::Real > &v, const amrex::Array4< amrex::Real const > &uold, const amrex::Array4< amrex::Real const > &vold, const amrex::Array4< amrex::Real > &rufrc, const amrex::Array4< amrex::Real > &rvfrc, const amrex::Array4< amrex::Real const > &visc2_p, const amrex::Array4< amrex::Real const > &visc2_r, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &pn, const amrex::Array4< amrex::Real const > &mskp, int nrhs, int nnew, const amrex::Real dt_lev)
Harmonic viscosity.
Definition REMORA_uv3dmix.cpp:25
REMORA::rhs_uv_2d
void rhs_uv_2d(int lev, const amrex::Box &xbx, const amrex::Box &ybx, const amrex::Array4< amrex::Real const > &uold, const amrex::Array4< amrex::Real const > &vold, const amrex::Array4< amrex::Real > &ru, const amrex::Array4< amrex::Real > &rv, const amrex::Array4< amrex::Real const > &Duon, const amrex::Array4< amrex::Real const > &Dvom, const int nrhs)
RHS terms for 2D momentum.
Definition REMORA_rhs_uv_2d.cpp:18
REMORA::ubar_clim_data_from_file
std::unique_ptr< NCTimeSeries > ubar_clim_data_from_file
Data container for ubar climatology data read from file.
Definition REMORA.H:1603
REMORA::bdy_vbar
int bdy_vbar() const noexcept
Definition REMORA.H:1465
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::t_new
amrex::Vector< amrex::Real > t_new
new time at each level, in seconds since start_time
Definition REMORA.H:1700
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_weight1
amrex::Vector< amrex::Real > vec_weight1
Weights for calculating avg1 in 2D advance.
Definition REMORA.H:661
REMORA::bdy_ubar
int bdy_ubar() const noexcept
Definition REMORA.H:1464
REMORA::cf_impose_flux
int cf_impose_flux
impose the parent's barotropic mass flux on DUon/DVom at the coarse-fine interface,...
Definition REMORA.H:1842
REMORA::FillPatchNoBC
void FillPatchNoBC(int lev, amrex::Real time, amrex::MultiFab &mf_to_be_filled, amrex::Vector< amrex::MultiFab * > const &mfs, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=true, const bool fill_set=false, 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 without applying boun...
Definition REMORA_FillPatch.cpp:211
REMORA::vec_zeta
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_zeta
free surface height (2D)
Definition REMORA.H:578
REMORA::ubar_bc
int ubar_bc() const noexcept
Definition REMORA.H:1451
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::vec_vbar
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_vbar
barotropic y velocity (2D)
Definition REMORA.H:576
REMORA::advance_2d
void advance_2d(int lev, amrex::MultiFab const *mf_rhoS, amrex::MultiFab const *mf_rhoA, amrex::MultiFab *mf_ru2d, amrex::MultiFab *mf_rv2d, amrex::MultiFab *mf_rufrc, amrex::MultiFab *mf_rvfrc, amrex::MultiFab *mf_Zt_avg1, 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_rubar, std::unique_ptr< amrex::MultiFab > &mf_rvbar, std::unique_ptr< amrex::MultiFab > &mf_rzeta, std::unique_ptr< amrex::MultiFab > &mf_ubar, std::unique_ptr< amrex::MultiFab > &mf_vbar, amrex::MultiFab *mf_zeta, amrex::MultiFab const *mf_h, amrex::MultiFab const *mf_pm, amrex::MultiFab const *mf_pn, amrex::MultiFab const *mf_fcor, amrex::MultiFab const *mf_visc2_p, amrex::MultiFab const *mf_visc2_r, amrex::MultiFab const *mf_mskr, amrex::MultiFab const *mf_msku, amrex::MultiFab const *mf_mskv, amrex::MultiFab const *mf_mskp, amrex::Real dtfast_lev, bool predictor_2d_step, bool first_2d_step, int my_iif, int &next_indx1)
Perform a 2D predictor (predictor_2d_step=True) or corrector (predictor_2d_step=False) step.
Definition REMORA_advance_2d.cpp:97
REMORA::vec_ubar
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ubar
barotropic x velocity (2D)
Definition REMORA.H:574
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::vbar_bc
int vbar_bc() const noexcept
Definition REMORA.H:1452
REMORA::vec_dndx
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_dndx
d(1/n)/d(xi)
Definition REMORA.H:642
REMORA::cf_time_interp_zeta
int cf_time_interp_zeta
interpolate the parent in time onto the child's own sub-time, as put_refine2d does,...
Definition REMORA.H:1881
REMORA::vec_nudg_coeff
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:678
REMORA::cf_set_2d_bcs
int cf_set_2d_bcs
0 = never impose the 2D coarse-fine interface condition, 1 = every fast step, 2 = only the first fast...
Definition REMORA.H:1867
REMORA::curvilinear
void curvilinear(const amrex::Box &bx, const amrex::Box &xbx, const amrex::Box &ybx, const amrex::Array4< amrex::Real const > &uold, const amrex::Array4< amrex::Real const > &vold, const amrex::Array4< amrex::Real > &ru, const amrex::Array4< amrex::Real > &rv, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &dndx, const amrex::Array4< amrex::Real const > &dmde, int nrhs, int nr)
Calculate curvilinear advection terms.
Definition REMORA_curvilinear.cpp:21
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::cf_fill_all_kcomp
int cf_fill_all_kcomp
write every leapfrog record of the coarse-fine ghost band, as put_refine2d does when it sets zeta(:,...
Definition REMORA.H:1876
REMORA::set_2d_cf_flux
void set_2d_cf_flux(int lev, amrex::Real time, amrex::MultiFab &mf_DUon, amrex::MultiFab &mf_DVom)
impose the parent's barotropic mass flux on the coarse-fine interface faces
Definition REMORA_set_2d_cf_bcs.cpp:310
REMORA::coriolis
void coriolis(const amrex::Box &xbx, const amrex::Box &ybx, const amrex::Array4< amrex::Real const > &uold, const amrex::Array4< amrex::Real const > &vold, const amrex::Array4< amrex::Real > &ru, const amrex::Array4< amrex::Real > &rv, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &fomn, int nrhs, int nr)
Calculate Coriolis terms.
Definition REMORA_coriolis.cpp:19
REMORA::vbar_clim_data_from_file
std::unique_ptr< NCTimeSeries > vbar_clim_data_from_file
Data container for vbar climatology data read from file.
Definition REMORA.H:1605