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/** Nonlinear shallow-water rpimitive equations predictor (Leap-frog) and
5 * corrector (Adams-Moulton) time-stepping engine. Corresponds to Nonlinear/step2d_LF_AM3.h
6 * in ROMS.
7 *
8 * @param[in ] lev level of refinement (coarsest level is 0)
9 * @param[in ] mf_rhoS density perturbation
10 * @param[in ] mf_rhoA vertically-averaged density
11 * @param[inout] mf_ru2d RHS contributions to 2D u-momentum
12 * @param[inout] mf_rv2d RHS contribtuions to 2D v-momentum
13 * @param[inout] mf_rufrc before first predictor, vertical integral of 3D RHS for uvel, converted to forcing terms
14 * @param[inout] mf_rvfrc before first predictor, vertical integral of 3D RHS for vvel, converted to forcing term
15 * @param[inout] mf_Zt_avg1 average of sea surface height over all fast steps
16 * @param[inout] mf_DU_avg1 time-averaged u-flux for 2D equations
17 * @param[inout] mf_DU_avg2 time-averaged u-flux for 3D equation coupling
18 * @param[inout] mf_DV_avg1 time-averaged v-flux for 2D equations
19 * @param[inout] mf_DV_avg2 time-averaged v-flux for 3D equation coupling
20 * @param[inout] mf_rubar RHS of vertically integrated u-momentum
21 * @param[inout] mf_rvbar RHS of vertically integrated v-momentum
22 * @param[inout] mf_rzeta RHS of sea surface height
23 * @param[inout] mf_ubar vertically integrated u-momentum
24 * @param[inout] mf_vbar vertically integrated v-momentum
25 * @param[inout] mf_zeta Sea-surface height
26 * @param[in ] mf_h Bathymetry
27 * @param[in ] mf_pm 1 / dx
28 * @param[in ] mf_pn 1 / dy
29 * @param[in ] mf_fcor Coriolis factor
30 * @param[inout] mf_visc2_p Harmonic viscosity at psi points
31 * @param[inout] mf_visc2_r Harmoic viscosity at rho points
32 * @param[in ] mf_mskr Land-sea mask at rho-points
33 * @param[in ] mf_msku Land-sea mask at u-points
34 * @param[in ] mf_mskv Land-sea mask at v-points
35 * @param[in ] mf_mskp Land-sea mask at psi-points
36 * @param[in ] dtfast_lev Length of current barotropic step
37 * @param[in ] predictor_2d_step Is this a predictor step?
38 * @param[in ] first_2d_step Is this the first barotropic step?
39 * @param[in ] my_iif Which barotropic predictor-corrector pair?
40 * @param[inout] next_indx1 Cached index for
41 */
42
43void
45 MultiFab const* mf_rhoS,
46 MultiFab const* mf_rhoA,
47 MultiFab * mf_ru2d,
48 MultiFab * mf_rv2d,
49 MultiFab * mf_rufrc,
50 MultiFab * mf_rvfrc,
51 MultiFab * mf_Zt_avg1,
52 std::unique_ptr<MultiFab>& mf_DU_avg1,
53 std::unique_ptr<MultiFab>& mf_DU_avg2,
54 std::unique_ptr<MultiFab>& mf_DV_avg1,
55 std::unique_ptr<MultiFab>& mf_DV_avg2,
56 std::unique_ptr<MultiFab>& mf_rubar,
57 std::unique_ptr<MultiFab>& mf_rvbar,
58 std::unique_ptr<MultiFab>& mf_rzeta,
59 std::unique_ptr<MultiFab>& mf_ubar,
60 std::unique_ptr<MultiFab>& mf_vbar,
61 MultiFab * mf_zeta,
62 MultiFab const* mf_h,
63 MultiFab const* mf_pm,
64 MultiFab const* mf_pn,
65 MultiFab const* mf_fcor,
66 MultiFab const* mf_visc2_p,
67 MultiFab const* mf_visc2_r,
68 MultiFab const* mf_mskr,
69 MultiFab const* mf_msku,
70 MultiFab const* mf_mskv,
71 MultiFab const* mf_mskp,
72 Real dtfast_lev,
74 bool first_2d_step, int my_iif,
75 int & next_indx1)
76{
77 BL_PROFILE("REMORA::advance2d()");
78 int iic = istep[lev];
79 const int nnew = 0;
80 const int nstp = 0;
81 int ntfirst = 0;
82
83 int knew = 3;
84 int krhs = (my_iif + iic) % 2 + 1;
85 int kstp = my_iif <=1 ? iic % 2 + 1 : (iic % 2 + my_iif % 2 + 1) % 2 + 1;
86 int indx1 = krhs;
88 next_indx1 = 3 - indx1;
89 } else {
91 kstp = 3 - knew;
92 krhs = 3;
93 //If it's not the auxiliary time step, set indx1 to next_indx1
94 // NOTE: should this ever not execute?
95 // Include indx1 updates for diagnostic purposes?
96 // if (my_iif<nfast+1)
97 // indx1=next_indx1;
98 }
99 int ptsk = 3-kstp;
100 knew-=1;
101 krhs-=1;
102 kstp-=1;
103 // Include indx1 updates for diagnostic purposes?
104 //indx1-=1;
105 ptsk-=1;
106 auto ba = mf_h->boxArray();
107 auto dm = mf_h->DistributionMap();
108
109 MultiFab mf_DUon(convert(ba,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0));
110 MultiFab mf_DVom(convert(ba,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0));
111
112 int ncomp = 0;
113 int fomn_comp = ncomp++;
114 int Drhs_comp = ncomp++;
115 int Dnew_comp = ncomp++;
116 int zwrk_comp = ncomp++;
117 int gzeta_comp = ncomp++;
118 int gzeta2_comp = ncomp++;
119 int gzetaSA_comp = ncomp++;
120 int Dstp_comp = ncomp++;
121 int rhs_ubar_comp = ncomp++;
122 int rhs_vbar_comp = ncomp++;
123 int rhs_zeta_comp = ncomp++;
124 int zeta_new_comp = ncomp++;
125
126 MultiFab mf(ba,dm,ncomp,IntVect(NGROW+1,NGROW+1,0));
127
128 for ( MFIter mfi(*mf_rhoS, TilingIfNotGPU()); mfi.isValid(); ++mfi )
129 {
130 Array4<Real > const& ubar = mf_ubar->array(mfi);
131 Array4<Real > const& vbar = mf_vbar->array(mfi);
132 Array4<Real > const& zeta = mf_zeta->array(mfi);
133 Array4<Real const> const& h = mf_h->const_array(mfi);
134
135 Array4<Real const> const& pm = mf_pm->const_array(mfi);
136 Array4<Real const> const& pn = mf_pn->const_array(mfi);
137
138 Box bx = mfi.tilebox();
139 Box gbx = mfi.growntilebox();
140 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
141 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
142 Box xgbx2 = mfi.grownnodaltilebox(0, IntVect(NGROW,NGROW,0));
143 Box ygbx2 = mfi.grownnodaltilebox(1, IntVect(NGROW,NGROW,0));
144
145 Box tbxp1 = bx;
146 Box tbxp2 = bx;
147 Box tbxp3 = bx;
148 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
149 tbxp2.grow(IntVect(NGROW,NGROW,0));
150 tbxp3.grow(IntVect(NGROW+1,NGROW+1,0));
151
152 Box bxD = bx ; bxD.makeSlab(2,0);
153 Box gbxD = gbx ; gbxD.makeSlab(2,0);
154 Box gbx1D = gbx1; gbx1D.makeSlab(2,0);
155 Box gbx2D = gbx2; gbx2D.makeSlab(2,0);
156
157 Box tbxp2D = tbxp2;
158 tbxp2D.makeSlab(2,0);
159
160 // step2d work arrays
161 FArrayBox fab_Drhs(makeSlab(tbxp3,2,0),1,The_Async_Arena());
162 auto Drhs=fab_Drhs.array();
163
164 auto DUon = mf_DUon.array(mfi);
165 auto DVom = mf_DVom.array(mfi);
166
167 ParallelFor(makeSlab(tbxp3,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
168 {
169 Drhs(i,j,0)=zeta(i,j,0,krhs)+h(i,j,0);
170 });
171
172 ParallelFor(makeSlab(xgbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
173 {
174 Real on_u = two / (pn(i,j,0)+pn(i-1,j,0));
175 Real cff1= Real(0.5) * on_u *(Drhs(i,j,0)+Drhs(i-1,j,0));
176 DUon(i,j,0)=ubar(i,j,0,krhs)*cff1;
177 });
178
179 ParallelFor(makeSlab(ygbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
180 {
181 Real om_v = two / (pm(i,j,0)+pm(i,j-1,0));
182 Real cff1= Real(0.5) * om_v * (Drhs(i,j,0)+Drhs(i,j-1,0));
183 DVom(i,j,0)=vbar(i,j,0,krhs)*cff1;
184 });
185 }
186
187 // These are needed to pass the tests with bathymetry but I don't quite see why
188 mf_DUon.FillBoundary(geom[lev].periodicity());
189 mf_DVom.FillBoundary(geom[lev].periodicity());
190
191#ifdef REMORA_USE_NETCDF
193 ubar_clim_data_from_file->update_interpolated_to_time(t_new[lev], lev, vec_ubar[lev].get(), geom, ref_ratio);
194 vbar_clim_data_from_file->update_interpolated_to_time(t_new[lev], lev, vec_vbar[lev].get(), geom, ref_ratio);
195 }
196#endif
197
198 for ( MFIter mfi(*mf_rhoS, TilingIfNotGPU()); mfi.isValid(); ++mfi )
199 {
200 Array4<Real const> const& rhoS = mf_rhoS->const_array(mfi);
201 Array4<Real const> const& rhoA = mf_rhoA->const_array(mfi);
202 Array4<Real const> const& h = mf_h->const_array(mfi);
203
204 Array4<Real > const& rufrc = mf_rufrc->array(mfi);
205 Array4<Real > const& rvfrc = mf_rvfrc->array(mfi);
206 Array4<Real > const& Zt_avg1 = mf_Zt_avg1->array(mfi);
207 Array4<Real > const& ubar = mf_ubar->array(mfi);
208 Array4<Real > const& vbar = mf_vbar->array(mfi);
209 Array4<Real > const& zeta = mf_zeta->array(mfi);
210 Array4<Real > const& DU_avg1 = (mf_DU_avg1)->array(mfi);
211 Array4<Real > const& DU_avg2 = (mf_DU_avg2)->array(mfi);
212 Array4<Real > const& DV_avg1 = (mf_DV_avg1)->array(mfi);
213 Array4<Real > const& DV_avg2 = (mf_DV_avg2)->array(mfi);
214 Array4<Real > const& ru2d = (mf_ru2d)->array(mfi);
215 Array4<Real > const& rv2d = (mf_rv2d)->array(mfi);
216 Array4<Real > const& rubar = (mf_rubar)->array(mfi);
217 Array4<Real > const& rvbar = (mf_rvbar)->array(mfi);
218 Array4<Real > const& rzeta = (mf_rzeta)->array(mfi);
219 Array4<Real const> const& visc2_p = mf_visc2_p->const_array(mfi);
220 Array4<Real const> const& visc2_r = mf_visc2_r->const_array(mfi);
221
222 Array4<Real const> const& pm = mf_pm->const_array(mfi);
223 Array4<Real const> const& pn = mf_pn->const_array(mfi);
224 Array4<Real const> const& fcor = mf_fcor->const_array(mfi);
225
226 Array4<Real const> const& mskr = mf_mskr->const_array(mfi);
227 Array4<Real const> const& msku = mf_msku->const_array(mfi);
228 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
229 Array4<Real const> const& mskp = mf_mskp->const_array(mfi);
230
231 Box bx = mfi.tilebox();
232 Box gbx = mfi.growntilebox();
233 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
234 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
235 Box gbx3 = mfi.growntilebox(IntVect(NGROW+1,NGROW+1,0));
236 Box xgbx2 = mfi.grownnodaltilebox(0, IntVect(NGROW,NGROW,0));
237 Box ygbx2 = mfi.grownnodaltilebox(1, IntVect(NGROW,NGROW,0));
238
239 Box xbxD = mfi.nodaltilebox(0);
240 xbxD.makeSlab(2,0);
241
242 Box ybxD = mfi.nodaltilebox(1);
243 ybxD.makeSlab(2,0);
244
245 Box tbxp1 = bx; tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
246 Box tbxp2 = bx; tbxp2.grow(IntVect(NGROW,NGROW,0));
247 Box tbxp3 = bx; tbxp3.grow(IntVect(NGROW+1,NGROW+1,0));
248
249 Box bxD = bx; bxD.makeSlab(2,0);
250 Box gbxD = gbx; gbxD.makeSlab(2,0);
251 Box gbx1D = gbx1; gbx1D.makeSlab(2,0);
252 Box gbx2D = gbx2; gbx2D.makeSlab(2,0);
253
254 Box tbxp2D = tbxp2;
255 tbxp2D.makeSlab(2,0);
256
257 auto fomn = mf.array(mfi,fomn_comp);
258 auto Drhs = mf.array(mfi,Drhs_comp);
259 auto Drhs_const = mf.const_array(mfi,Drhs_comp);
260 auto Dnew = mf.array(mfi,Dnew_comp);
261 auto zwrk = mf.array(mfi,zwrk_comp);
262 auto gzeta = mf.array(mfi,gzeta_comp);
263 auto gzeta2 = mf.array(mfi,gzeta2_comp);
264 auto gzetaSA = mf.array(mfi,gzetaSA_comp);
265 auto Dstp = mf.array(mfi,Dstp_comp);
266 auto rhs_ubar = mf.array(mfi,rhs_ubar_comp);
267 auto rhs_vbar = mf.array(mfi,rhs_vbar_comp);
268 auto rhs_zeta = mf.array(mfi,rhs_zeta_comp);
269 auto zeta_new = mf.array(mfi,zeta_new_comp);
270
271 FArrayBox & fab_DUon=mf_DUon[mfi];
272 FArrayBox & fab_DVom=mf_DVom[mfi];
273 auto DUon=fab_DUon.array();
274 auto DVom=fab_DVom.array();
275
276 auto weight1 = vec_weight1.dataPtr();
277 auto weight2 = vec_weight2.dataPtr();
278
279 //From ana_grid.h and metrics.F
280 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int)
281 {
282 rhs_ubar(i,j,0)=zero;
283 });
284
285 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int)
286 {
287 rhs_vbar(i,j,0)=zero;
288 });
289
291 ParallelFor(tbxp2D, [=] AMREX_GPU_DEVICE (int i, int j, int )
292 {
293 fomn(i,j,0) = fcor(i,j,0)*(one/(pm(i,j,0)*pn(i,j,0)));
294 });
295 }
296
297 ParallelFor(makeSlab(tbxp3,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
298 {
299 Drhs(i,j,0)=zeta(i,j,0,krhs)+h(i,j,0);
300 });
301
303 {
304 if(first_2d_step) {
305 Real cff2=(Real(-1.0)/Real(12.0))*weight2[my_iif+1];
307 [=] AMREX_GPU_DEVICE (int i, int j, int)
308 {
309 Zt_avg1(i,j,0)=zero;
310 });
312 [=] AMREX_GPU_DEVICE (int i, int j, int)
313 {
314 DU_avg1(i,j,0)=zero;
315 DU_avg2(i,j,0)=cff2*DUon(i,j,0);
316 });
318 [=] AMREX_GPU_DEVICE (int i, int j, int)
319 {
320 DV_avg1(i,j,0)=zero;
321 DV_avg2(i,j,0)=cff2*DVom(i,j,0);
322 });
323 }
324 else {
325 Real cff1_wt1 = weight1[my_iif-1];
326 Real cff2_wt1 = (Real(8.0)/Real(12.0))*weight2[my_iif]-
327 (one/Real(12.0))*weight2[my_iif+1];
328
329 ParallelFor(makeSlab(gbx3,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
330 {
331 Zt_avg1(i,j,0) += cff1_wt1*zeta(i,j,0,krhs);
332 });
333
334 ParallelFor(makeSlab(xgbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
335 {
336 DU_avg1(i,j,0) += cff1_wt1*DUon(i,j,0);
337 DU_avg2(i,j,0) += cff2_wt1*DUon(i,j,0);
338 });
339
340 ParallelFor(makeSlab(ygbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
341 {
342 DV_avg1(i,j,0) += cff1_wt1*DVom(i,j,0);
343 DV_avg2(i,j,0) += cff2_wt1*DVom(i,j,0);
344 });
345 }
346 }
347 else {
348 Real cff2_wt2;
349
350 if (first_2d_step) {
352 } else {
353 cff2_wt2=Real(5.0)/Real(12.0)*weight2[my_iif];
354 }
355
356 ParallelFor(makeSlab(xgbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
357 {
358 DU_avg2(i,j,0)=DU_avg2(i,j,0)+cff2_wt2*DUon(i,j,0);
359 });
360
361 ParallelFor(makeSlab(ygbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
362 {
363 DV_avg2(i,j,0)=DV_avg2(i,j,0)+cff2_wt2*DVom(i,j,0);
364 });
365 }
366 //
367 // Do not perform the actual time stepping during the auxiliary
368 // (nfast(ng)+1) time step. Jump to next box
369 //
370
371 if (my_iif>=nfast) {
372 continue; }
373 //Load new free-surface values into shared array at both predictor
374 //and corrector steps
375 //
376 //=======================================================================
377 // Time step free-surface equation.
378 //=======================================================================
379 //
380 // During the first time-step, the predictor step is Forward-Euler
381 // and the corrector step is Backward-Euler. Otherwise, the predictor
382 // step is Leap-frog and the corrector step is Adams-Moulton.
383 //
384
385 Real fac=Real(1000.0)/solverChoice.rho0;
386
387 if (my_iif==0) {
388 Real cff1=dtfast_lev;
389
390 ParallelFor(makeSlab(tbxp1,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
391 {
392 rhs_zeta(i,j,0) = (DUon(i,j,0)-DUon(i+1,j,0))+
393 (DVom(i,j,0)-DVom(i,j+1,0));
394 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);
395 Dnew(i,j,0) = zeta_new(i,j,0)+h(i,j,0);
396
397 //Pressure gradient terms:
398 zwrk(i,j,0)=Real(0.5)*(zeta(i,j,0,kstp)+zeta_new(i,j,0));
399 gzeta(i,j,0)=(fac+rhoS(i,j,0))*zwrk(i,j,0);
400 gzeta2(i,j,0)=gzeta(i,j,0)*zwrk(i,j,0);
401 gzetaSA(i,j,0)=zwrk(i,j,0)*(rhoS(i,j,0)-rhoA(i,j,0));
402 });
403
404 } else if (predictor_2d_step) {
405
406 Real cff1=two * dtfast_lev;
407 Real cff4=Real(4.0) / Real(25.0);
408 Real cff5=one - two*cff4;
409
410 ParallelFor(makeSlab(tbxp1,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
411 {
412 rhs_zeta(i,j,0)=(DUon(i,j,0)-DUon(i+1,j,0))+
413 (DVom(i,j,0)-DVom(i,j+1,0));
414 zeta_new(i,j,0)=(zeta(i,j,0,kstp)+
415 pm(i,j,0)*pn(i,j,0)*cff1*rhs_zeta(i,j,0)) * mskr(i,j,0);
416 Dnew(i,j,0)=zeta_new(i,j,0)+h(i,j,0);
417 //Pressure gradient terms
418 zwrk(i,j,0)=cff5*zeta(i,j,0,krhs)+
419 cff4*(zeta(i,j,0,kstp)+zeta_new(i,j,0));
420 gzeta(i,j,0)=(fac+rhoS(i,j,0))*zwrk(i,j,0);
421 gzeta2(i,j,0)=gzeta(i,j,0)*zwrk(i,j,0);
422 gzetaSA(i,j,0)=zwrk(i,j,0)*(rhoS(i,j,0)-rhoA(i,j,0));
423 });
424
425 } else if (!predictor_2d_step) { //AKA if(corrector_2d_step)
426
427 Real cff1=dtfast_lev * Real(5.0)/Real(12.0);
428 Real cff2=dtfast_lev * Real(8.0)/Real(12.0);
429 Real cff3=dtfast_lev * one/Real(12.0);
430 Real cff4=two/Real(5.0);
431 Real cff5=one-cff4;
432
433 ParallelFor(makeSlab(tbxp1,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
434 {
435 Real cff=cff1*((DUon(i,j,0)-DUon(i+1,j,0))+
436 (DVom(i,j,0)-DVom(i,j+1,0)));
437 zeta_new(i,j,0)=zeta(i,j,0,kstp)+
438 pm(i,j,0)*pn(i,j,0)*(cff+
439 cff2*rzeta(i,j,0,kstp)-
440 cff3*rzeta(i,j,0,ptsk));
441 zeta_new(i,j,0) *= mskr(i,j,0);
442 Dnew(i,j,0)=zeta_new(i,j,0)+h(i,j,0);
443 //Pressure gradient terms
444 zwrk(i,j,0)=cff5*zeta_new(i,j,0)+cff4*zeta(i,j,0,krhs);
445 gzeta(i,j,0)=(fac+rhoS(i,j,0))*zwrk(i,j,0);
446 gzeta2(i,j,0)=gzeta(i,j,0)*zwrk(i,j,0);
447 gzetaSA(i,j,0)=zwrk(i,j,0)*(rhoS(i,j,0)-rhoA(i,j,0));
448 });
449 }
450
451 //
452 // Load new free-surface values into shared array at both predictor
453 // and corrector steps.
454 //
455 //// zeta(knew) only valid at zeta_new, i.e. tbxp1
457 [=] AMREX_GPU_DEVICE (int i, int j, int )
458 {
459 zeta(i,j,0,knew) = zeta_new(i,j,0);
460 });
461
462 //
463 // If predictor step, load right-side-term into shared array.
464 //
465 if (predictor_2d_step) {
466 ParallelFor(makeSlab(gbx1,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
467 {
468 rzeta(i,j,0,krhs)=rhs_zeta(i,j,0);
469 });
470 }
471
472 //
473 //=======================================================================
474 // Compute right-hand-side for the 2D momentum equations.
475 //=======================================================================
476 //
477/*
478!
479!-----------------------------------------------------------------------
480! Compute pressure gradient terms.
481!-----------------------------------------------------------------------
482!
483*/
484 Real cff1 = Real(0.5) * g;
485 Real cff2 = one / Real(3.0);
487 [=] AMREX_GPU_DEVICE (int i, int j, int )
488 {
489 Real on_u = two / (pn(i,j,0)+pn(i-1,j,0));
490 rhs_ubar(i,j,0)=cff1 * on_u *
491 (( h(i-1,j,0) + h(i,j,0))*
492 (gzeta(i-1,j,0) - gzeta(i,j,0))+
493 ( h(i-1,j,0) - h(i,j,0))*
494 ( gzetaSA(i-1,j,0) + gzetaSA(i,j,0)+
495 cff2*( rhoA(i-1,j,0) - rhoA(i,j,0))*
496 ( zwrk(i-1,j,0) - zwrk(i,j,0)))+
497 (gzeta2(i-1,j,0)- gzeta2(i ,j,0)));
498 });
499
501 [=] AMREX_GPU_DEVICE (int i, int j, int )
502 {
503 Real om_v = two / (pm(i,j,0)+pm(i,j-1,0));
504 rhs_vbar(i,j,0) = cff1*om_v *
505 (( h(i,j-1,0) + h(i,j,0))*
506 (gzeta(i,j-1,0) - gzeta(i,j,0))+
507 ( h(i,j-1,0) - h(i,j,0))*
508 (gzetaSA(i,j-1,0)+ gzetaSA(i,j ,0)+
509 cff2*(rhoA(i,j-1,0)- rhoA(i,j ,0))*
510 (zwrk(i,j-1,0)- zwrk(i,j ,0)))+
511 (gzeta2(i,j-1,0)- gzeta2(i,j ,0)));
512 });
513
514 // Advection terms for 2d ubar, vbar added to rhs_ubar and rhs_vbar
515 //
516 //-----------------------------------------------------------------------
517 // rhs_uv_2d
518 //-----------------------------------------------------------------------
519 //
520 Array4<Real const> const& ubar_const = mf_ubar->const_array(mfi);
521 Array4<Real const> const& vbar_const = mf_vbar->const_array(mfi);
522
524
525 //-----------------------------------------------------------------------
526 // Add Coriolis forcing
527 //-----------------------------------------------------------------------
529 // Coriolis terms for 2d ubar, vbar added to rhs_ubar and rhs_vbar
530 //
531 //-----------------------------------------------------------------------
532 // coriolis
533 //-----------------------------------------------------------------------
534 //
536 }
537
539 Array4<Real const> const& dndx = vec_dndx[lev]->const_array(mfi);
540 Array4<Real const> const& dmde = vec_dmde[lev]->const_array(mfi);
542 }
543
544 //-----------------------------------------------------------------------
545 //Add in horizontal harmonic viscosity.
546 // Consider generalizing or copying uv3dmix, where Drhs is used instead of Hz and u=>ubar v=>vbar, drop dt terms
547 //-----------------------------------------------------------------------
548 uv3dmix(xbxD, ybxD, ubar, vbar, ubar, vbar, rhs_ubar, rhs_vbar,
550 pm, pn, mskp, krhs, nnew, zero);
551
552#ifdef REMORA_USE_NETCDF
554 Array4<Real > const& ubar_krhs = mf_ubar->array(mfi, krhs);
555 Array4<Real > const& vbar_krhs = mf_vbar->array(mfi, krhs);
556 Array4<const Real> const& ubar_clim = ubar_clim_data_from_file->get_interpolated_mf(lev)->const_array(mfi);
557 Array4<const Real> const& vbar_clim = vbar_clim_data_from_file->get_interpolated_mf(lev)->const_array(mfi);
560 // Boxes are like this to match the ROMS loops i=IstrU..Iend, j=JstrV..Jend, which
561 // exclude the u-face on the west/east domain edges and the v-face on south/north
562 Box xbxD_adj = clim_nudg_momentum_box(mfi.nodaltilebox(0), 0,
563 Geom(lev).Domain(), Geom(lev).isPeriodic(0));
564 xbxD_adj.makeSlab(2,0);
565 Box ybxD_adj = clim_nudg_momentum_box(mfi.nodaltilebox(1), 1,
566 Geom(lev).Domain(), Geom(lev).isPeriodic(1));
567 ybxD_adj.makeSlab(2,0);
570 }
571#endif
572
573 //-----------------------------------------------------------------------
574 // Coupling from 3d to 2d
575 //-----------------------------------------------------------------------
577 {
578 if (iic==ntfirst) {
579 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
580 {
581 rufrc(i,j,0) -= rhs_ubar(i,j,0);
582 rhs_ubar(i,j,0) += rufrc(i,j,0);
583 ru2d(i,j,0,nstp) = rufrc(i,j,0);
584 });
585
586 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
587 {
588 rvfrc(i,j,0) -= rhs_vbar(i,j,0);
589 rhs_vbar(i,j,0) += rvfrc(i,j,0);
590 rv2d(i,j,0,nstp) = rvfrc(i,j,0);
591 });
592
593 } else if (iic==(ntfirst+1)) {
594
595 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
596 {
597 rufrc(i,j,0)=rufrc(i,j,0)-rhs_ubar(i,j,0);
598 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);
599 ru2d(i,j,0,1)=rufrc(i,j,0);
600 Real r_swap= ru2d(i,j,0,1);
601 ru2d(i,j,0,1) = ru2d(i,j,0,0);
602 ru2d(i,j,0,0) = r_swap;
603 });
604
605 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
606 {
607 rvfrc(i,j,0)=rvfrc(i,j,0)-rhs_vbar(i,j,0);
608 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);
609 rv2d(i,j,0,1)=rvfrc(i,j,0);
610 Real r_swap= rv2d(i,j,0,1);
611 rv2d(i,j,0,1) = rv2d(i,j,0,0);
612 rv2d(i,j,0,0) = r_swap;
613 });
614
615 } else {
616 cff1=Real(23.0)/Real(12.0);
617 cff2=Real(16.0)/Real(12.0);
618 Real cff3= Real(5.0)/Real(12.0);
619
620 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
621 {
622 rufrc(i,j,0)=rufrc(i,j,0)-rhs_ubar(i,j,0);
623 rhs_ubar(i,j,0)=rhs_ubar(i,j,0)+
624 cff1*rufrc(i,j,0)-
625 cff2*ru2d(i,j,0,0)+
626 cff3*ru2d(i,j,0,1);
627 ru2d(i,j,0,1)=rufrc(i,j,0);
628 Real r_swap= ru2d(i,j,0,1);
629 ru2d(i,j,0,1) = ru2d(i,j,0,0);
630 ru2d(i,j,0,0) = r_swap;
631 });
632
633 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
634 {
635 rvfrc(i,j,0)=rvfrc(i,j,0)-rhs_vbar(i,j,0);
636 rhs_vbar(i,j,0)=rhs_vbar(i,j,0)+
637 cff1*rvfrc(i,j,0)-
638 cff2*rv2d(i,j,0,0)+
639 cff3*rv2d(i,j,0,1);
640 rv2d(i,j,0,1)=rvfrc(i,j,0);
641
642 Real r_swap= rv2d(i,j,0,1);
643 rv2d(i,j,0,1) = rv2d(i,j,0,0);
644 rv2d(i,j,0,0) = r_swap;
645 });
646 }
647 } else {
648 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
649 {
650 rhs_ubar(i,j,0) += rufrc(i,j,0);
651 });
652
653 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
654 {
655 rhs_vbar(i,j,0) += rvfrc(i,j,0);
656 });
657 }
658
659 //
660 //=======================================================================
661 // Time step 2D momentum equations.
662 //=======================================================================
663 //
664 // Compute total water column depth.
665 //
666 ParallelFor(makeSlab(tbxp3,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
667 {
668 Dstp(i,j,0)=zeta(i,j,0,kstp)+h(i,j,0);
669 });
670
671 //
672 // During the first time-step, the predictor step is Forward-Euler
673 // and the corrector step is Backward-Euler. Otherwise, the predictor
674 // step is Leap-frog and the corrector step is Adams-Moulton.
675 //
676 if (my_iif==0) {
677 cff1=Real(0.5)*dtfast_lev;
679 [=] AMREX_GPU_DEVICE (int i, int j, int )
680 {
681 Real cff=(pm(i,j,0)+pm(i-1,j,0))*(pn(i,j,0)+pn(i-1,j,0));
682 Real Dnew_avg =one/(Dnew(i,j,0)+Dnew(i-1,j,0));
683 ubar(i,j,0,knew)=(ubar(i,j,0,kstp)*
684 (Dstp(i,j,0)+Dstp(i-1,j,0))+
685 cff*cff1*rhs_ubar(i,j,0))*Dnew_avg * msku(i,j,0);
686 });
688 [=] AMREX_GPU_DEVICE (int i, int j, int )
689 {
690 Real cff=(pm(i,j,0)+pm(i,j-1,0))*(pn(i,j,0)+pn(i,j-1,0));
691 Real Dnew_avg=one/(Dnew(i,j,0)+Dnew(i,j-1,0));
692 vbar(i,j,0,knew)=(vbar(i,j,0,kstp)*
693 (Dstp(i,j,0)+Dstp(i,j-1,0))+
694 cff*cff1*rhs_vbar(i,j,0))*Dnew_avg * mskv(i,j,0);
695 });
696
697 } else if (predictor_2d_step) {
698
701 [=] AMREX_GPU_DEVICE (int i, int j, int )
702 {
703 Real cff=(pm(i,j,0)+pm(i-1,j,0))*(pn(i,j,0)+pn(i-1,j,0));
704 Real Dnew_avg=one/(Dnew(i,j,0)+Dnew(i-1,j,0));
705 ubar(i,j,0,knew)=(ubar(i,j,0,kstp)*
706 (Dstp(i,j,0)+Dstp(i-1,j,0))+
707 cff*cff1*rhs_ubar(i,j,0))*Dnew_avg * msku(i,j,0);
708 });
710 [=] AMREX_GPU_DEVICE (int i, int j, int )
711 {
712 Real cff=(pm(i,j,0)+pm(i,j-1,0))*(pn(i,j,0)+pn(i,j-1,0));
713 Real Dnew_avg=one/(Dnew(i,j,0)+Dnew(i,j-1,0));
714 vbar(i,j,0,knew)=(vbar(i,j,0,kstp)*
715 (Dstp(i,j,0)+Dstp(i,j-1,0))+
716 cff*cff1*rhs_vbar(i,j,0))*Dnew_avg * mskv(i,j,0);
717 });
718
719 } else if ((!predictor_2d_step)) {
720
721 cff1=Real(0.5)*dtfast_lev*Real(5.0)/Real(12.0);
722 cff2=Real(0.5)*dtfast_lev*Real(8.0)/Real(12.0);
723 Real cff3=Real(0.5)*dtfast_lev*one/Real(12.0);
725 [=] AMREX_GPU_DEVICE (int i, int j, int )
726 {
727 Real cff=(pm(i,j,0)+pm(i-1,j,0))*(pn(i,j,0)+pn(i-1,j,0));
728 Real Dnew_avg=one/(Dnew(i,j,0)+Dnew(i-1,j,0));
729 ubar(i,j,0,knew)=(ubar(i,j,0,kstp)*
730 (Dstp(i,j,0)+Dstp(i-1,j,0))+
731 cff*(cff1*rhs_ubar(i,j,0)+
732 cff2*rubar(i,j,0,kstp)-
733 cff3*rubar(i,j,0,ptsk)))*Dnew_avg * msku(i,j,0);
734 });
736 [=] AMREX_GPU_DEVICE (int i, int j, int )
737 {
738 Real cff=(pm(i,j,0)+pm(i,j-1,0))*(pn(i,j,0)+pn(i,j-1,0));
739 Real Dnew_avg=one/(Dnew(i,j,0)+Dnew(i,j-1,0));
740 vbar(i,j,0,knew)=(vbar(i,j,0,kstp)*
741 (Dstp(i,j,0)+Dstp(i,j-1,0))+
742 cff*(cff1*rhs_vbar(i,j,0)+
743 cff2*rvbar(i,j,0,kstp)-
744 cff3*rvbar(i,j,0,ptsk)))*Dnew_avg * mskv(i,j,0);
745 });
746 }
747
748 //store rhs_ubar and rhs_vbar to save later
749 //
750 // If predictor step, load right-side-term into shared arrays for
751 // future use during the subsequent corrector step.
752 //
753
754 if (predictor_2d_step) {
755 ParallelFor(xbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
756 {
757 rubar(i,j,0,krhs)=rhs_ubar(i,j,0);
758 });
759 ParallelFor(ybxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
760 {
761 rvbar(i,j,0,krhs)=rhs_vbar(i,j,0);
762 });
763 }
764 }
765
766 // Don't do the FillPatch or rivers at the last truncated predictor step.
767 // We may need to move the zeta FillPatch further up
768 if (my_iif<nfast) {
769 int know;
770 Real dt2d;
771 if (my_iif==0) {
772 know = krhs;
774 } else if (predictor_2d_step) {
775 know = krhs;
776 dt2d = two * dtfast_lev;
777 } else {
778 know = kstp;
780 }
781
782 MultiFab ubar_know(*vec_ubar[lev], make_alias, know, 1);
783 MultiFab vbar_know(*vec_vbar[lev], make_alias, know, 1);
784 MultiFab zeta_know(*vec_zeta[lev], make_alias, know, 1);
786 knew, false,true, 0,know, dt2d, ubar_know);
788 knew, false,true, 0,know, dt2d, vbar_know);
790 knew, false,false, 0,know, dt2d, zeta_know);
791
792#ifdef REMORA_USE_NETCDF
794 river_source_transportbar->update_interpolated_to_time(t_old[lev]);
796 for ( MFIter mfi(*mf_rhoS, TilingIfNotGPU()); mfi.isValid(); ++mfi )
797 {
798 Array4<const int > const& river_pos = vec_river_position[lev]->const_array(mfi);
800 Array4<Real > const& ubar = mf_ubar->array(mfi);
801 Array4<Real > const& vbar = mf_vbar->array(mfi);
802 Array4<Real const> const& zeta = mf_zeta->const_array(mfi);
803 Array4<Real const> const& h = mf_h->const_array(mfi);
804 Array4<Real const> const& pm = mf_pm->const_array(mfi);
805 Array4<Real const> const& pn = mf_pn->const_array(mfi);
806
807 Box gbx1D = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
808 gbx1D.makeSlab(2,0);
809
810 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
811 {
812 int iriver = river_pos(i,j,0);
813 if (iriver >= 0) {
814 if (river_direction_d[iriver] == 0) {
815 Real on_u = two / (pn(i,j,0)+pn(i-1,j,0));
816 Real cff = one / (on_u * Real(0.5) * (zeta(i-1,j,0,knew) + h(i-1,j,0) +
817 zeta(i,j,0,knew) + h(i,j,0)));
818 ubar(i,j,0,knew) = river_transportbar(iriver,0,0) * cff;
819 } else {
820 Real om_v = two / (pm(i,j,0)+pm(i,j-1,0));
821 Real cff = one / (om_v * Real(0.5) * (zeta(i,j-1,0,knew) + h(i,j-1,0) +
822 zeta(i,j,0,knew) + h(i,j,0)));
823 vbar(i,j,0,knew) = river_transportbar(iriver,0,0) * cff;
824 }
825 }
826 });
827 }
828 }
830 knew, false,false);
832 knew, false,false);
833#endif
834 }
835}
constexpr amrex::Real two
constexpr amrex::Real one
constexpr amrex::Real zero
constexpr amrex::Real g
#define NGROW
mf_h setVal(geomdata.ProbHi(2))
int nfast
Number of fast steps to take.
Definition REMORA.H:1675
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.
std::unique_ptr< NCTimeSeriesRiver > river_source_transportbar
Data container for vertically integrated momentum transport in rivers.
Definition REMORA.H:1487
int zeta_bc() const noexcept
Definition REMORA.H:1318
amrex::Vector< amrex::Real > vec_weight2
Weights for calculating avg2 in 2D advance.
Definition REMORA.H:628
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_dmde
d(1/m)/d(eta)
Definition REMORA.H:609
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)
Fill a new MultiFab by copying in phi from valid region and filling ghost cells without applying boun...
int bdy_zeta() const noexcept
Definition REMORA.H:1331
amrex::Vector< int > istep
which step?
Definition REMORA.H:1552
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.
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.
std::unique_ptr< NCTimeSeries > ubar_clim_data_from_file
Data container for ubar climatology data read from file.
Definition REMORA.H:1470
int bdy_vbar() const noexcept
Definition REMORA.H:1330
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.
amrex::Vector< amrex::Real > t_new
new time at each level
Definition REMORA.H:1556
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1717
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > vec_river_position
iMultiFab for river positions; contents are indices of rivers
Definition REMORA.H:1494
amrex::Vector< amrex::Real > vec_weight1
Weights for calculating avg1 in 2D advance.
Definition REMORA.H:626
int bdy_ubar() const noexcept
Definition REMORA.H:1329
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_zeta
free surface height (2D)
Definition REMORA.H:551
int ubar_bc() const noexcept
Definition REMORA.H:1316
amrex::Gpu::DeviceVector< int > river_direction
Vector over rivers of river direction: 0: u-face; 1: v-face; 2: w-face.
Definition REMORA.H:1496
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_vbar
barotropic y velocity (2D)
Definition REMORA.H:549
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.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ubar
barotropic x velocity (2D)
Definition REMORA.H:547
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())
Fill a new MultiFab by copying in phi from valid region and filling ghost cells.
int vbar_bc() const noexcept
Definition REMORA.H:1317
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_dndx
d(1/n)/d(xi)
Definition REMORA.H:607
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:643
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.
amrex::Vector< amrex::Real > t_old
old time at each level
Definition REMORA.H:1558
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.
std::unique_ptr< NCTimeSeries > vbar_clim_data_from_file
Data container for vbar climatology data read from file.
Definition REMORA.H:1472