13 BL_PROFILE(
"REMORA::setup_step()");
26 int nvars = S_old.nComp();
36 const BoxArray& ba = S_old.boxArray();
37 const DistributionMapping& dm = S_old.DistributionMap();
43 MultiFab source(ba,dm,nvars,1);
44 source.setVal(0.0_rt);
51 MultiFab mf_AK(ba,dm,1,IntVect(
NGROW,
NGROW,0));
54 MultiFab mf_logdrg_tmp(ba,dm,1,IntVect(
NGROW,
NGROW,0));
56 MultiFab* mf_z_r =
vec_z_r[lev].get();
57 MultiFab* mf_z_w =
vec_z_w[lev].get();
59 MultiFab* mf_pm =
vec_pm[lev].get();
60 MultiFab* mf_pn =
vec_pn[lev].get();
61 MultiFab* mf_fcor =
vec_fcor[lev].get();
63 MultiFab* mf_gls =
vec_gls[lev].get();
64 MultiFab* mf_tke =
vec_tke[lev].get();
68 MultiFab mf_rho(ba,dm,1,IntVect(
NGROW,
NGROW,0));
69 std::unique_ptr<MultiFab>& mf_rhoS =
vec_rhoS[lev];
70 std::unique_ptr<MultiFab>& mf_rhoA =
vec_rhoA[lev];
71 std::unique_ptr<MultiFab>& mf_bvf =
vec_bvf[lev];
72 std::unique_ptr<MultiFab>& mf_ru =
vec_ru[lev];
73 std::unique_ptr<MultiFab>& mf_rv =
vec_rv[lev];
74 std::unique_ptr<MultiFab>& mf_rufrc =
vec_rufrc[lev];
75 std::unique_ptr<MultiFab>& mf_rvfrc =
vec_rvfrc[lev];
76 std::unique_ptr<MultiFab>& mf_sustr =
vec_sustr[lev];
77 std::unique_ptr<MultiFab>& mf_svstr =
vec_svstr[lev];
78 std::unique_ptr<MultiFab>& mf_rdrag =
vec_rdrag[lev];
79 std::unique_ptr<MultiFab>& mf_rdrag2 =
vec_rdrag2[lev];
80 std::unique_ptr<MultiFab>& mf_ZoBot =
vec_ZoBot[lev];
81 std::unique_ptr<MultiFab>& mf_bustr =
vec_bustr[lev];
82 std::unique_ptr<MultiFab>& mf_bvstr =
vec_bvstr[lev];
84 std::unique_ptr<MultiFab>& mf_mskr =
vec_mskr[lev];
85 std::unique_ptr<MultiFab>& mf_msku =
vec_msku[lev];
86 std::unique_ptr<MultiFab>& mf_mskv =
vec_mskv[lev];
87 std::unique_ptr<MultiFab>& mf_mskp =
vec_mskp[lev];
89 MultiFab mf_rw(ba,dm,1,IntVect(
NGROW,
NGROW,0));
91 std::unique_ptr<MultiFab>& mf_visc2_p =
vec_visc2_p[lev];
92 std::unique_ptr<MultiFab>& mf_visc2_r =
vec_visc2_r[lev];
96 mf_rho.setVal(0.e34_rt,IntVect(AMREX_D_DECL(
NGROW-1,
NGROW-1,0)));
97 mf_rhoS->setVal(0.e34_rt,IntVect(AMREX_D_DECL(
NGROW-1,
NGROW-1,0)));
98 mf_rhoA->setVal(0.e34_rt,IntVect(AMREX_D_DECL(
NGROW-1,
NGROW-1,0)));
106 mf_rw.setVal(0.0_rt);
110 int iic =
istep[lev];
113 MultiFab::Copy(S_new,S_old,0,0,S_new.nComp(),S_new.nGrowVect());
114 MultiFab::Copy(U_new,U_old,0,0,U_new.nComp(),U_new.nGrowVect());
115 MultiFab::Copy(V_new,V_old,0,0,V_new.nComp(),V_new.nGrowVect());
116 MultiFab::Copy(W_new,W_old,0,0,W_new.nComp(),W_new.nGrowVect());
127 auto N = Geom(lev).Domain().size()[2]-1;
129 const auto dlo = amrex::lbound(Geom(lev).Domain());
130 const auto dhi = amrex::ubound(Geom(lev).Domain());
132 for ( MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi )
135 Array4<Real const>
const& Hz =
vec_Hz[lev]->const_array(mfi);
136 Array4<Real >
const& Huon =
vec_Huon[lev]->array(mfi);
137 Array4<Real >
const& Hvom =
vec_Hvom[lev]->array(mfi);
139 Array4<Real const>
const& z_w = mf_z_w->const_array(mfi);
140 Array4<Real const>
const& z_r = mf_z_r->const_array(mfi);
141 Array4<Real const>
const& uold = U_old.const_array(mfi);
142 Array4<Real const>
const& vold = V_old.const_array(mfi);
143 Array4<Real >
const& rho = mf_rho.array(mfi);
144 Array4<Real >
const& rhoA = mf_rhoA->array(mfi);
145 Array4<Real >
const& rhoS = mf_rhoS->array(mfi);
146 Array4<Real >
const& bvf = mf_bvf->array(mfi);
150 Array4<Real const>
const& pm = mf_pm->const_array(mfi);
151 Array4<Real const>
const& pn = mf_pn->const_array(mfi);
153 Array4<Real const>
const& mskr = mf_mskr->const_array(mfi);
155 Box bx = mfi.tilebox();
156 Box gbx1 = mfi.growntilebox(IntVect(
NGROW-1,
NGROW-1,0));
157 Box gbx2 = mfi.growntilebox(IntVect(
NGROW,
NGROW,0));
158 Box ugbx2 = mfi.grownnodaltilebox(0,IntVect(
NGROW,
NGROW,0));
159 Box vgbx2 = mfi.grownnodaltilebox(1,IntVect(
NGROW,
NGROW,0));
168 FArrayBox fab_FC(gbx2,1,amrex::The_Async_Arena());
169 FArrayBox fab_FX(gbx2,1,amrex::The_Async_Arena());
170 FArrayBox fab_FE(gbx2,1,amrex::The_Async_Arena());
171 FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
172 FArrayBox fab_CF(gbx2,1,amrex::The_Async_Arena());
179 ParallelFor(ugbx2, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
181 Real on_u = 2.0_rt / (pn(i-1,j,0)+pn(i,j,0));
182 Huon(i,j,k)=0.5_rt*(Hz(i,j,k)+Hz(i-1,j,k))*uold(i,j,k)* on_u;
185 ParallelFor(vgbx2, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
187 Real om_v= 2.0_rt / (pm(i,j-1,0)+pm(i,j,0));
188 Hvom(i,j,k)=0.5_rt*(Hz(i,j,k)+Hz(i,j-1,k))*vold(i,j,k)* om_v;
191 Array4<Real const>
const& state_old = S_old.const_array(mfi);
192 rho_eos(gbx2,state_old,rho,rhoA,rhoS,bvf,alpha,beta,Hz,z_w,z_r,h,mskr,N);
203 vec_evap[lev]->FillBoundary(geom[lev].periodicity());
207 for ( MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi )
209 Array4<Real >
const& stflx =
vec_stflx[lev]->array(mfi);
210 Array4<Real >
const& btflx =
vec_btflx[lev]->array(mfi);
211 Array4<Real const>
const& stflux =
vec_stflux[lev]->array(mfi);
212 Array4<Real const>
const& btflux =
vec_btflux[lev]->array(mfi);
213 Box gbx2 = mfi.growntilebox(IntVect(
NGROW,
NGROW,0));
216 ParallelFor(gbx2D, [=] AMREX_GPU_DEVICE (
int i,
int j,
int ) {
223 for ( MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi )
225 Array4<Real >
const& stflx =
vec_stflx[lev]->array(mfi);
226 Array4<Real >
const& btflx =
vec_btflx[lev]->array(mfi);
227 Array4<Real const>
const& stflux =
vec_stflux[lev]->const_array(mfi);
228 Array4<Real const>
const& salt_old = S_old.const_array(mfi,
Salt_comp);
229 Box gbx2 = mfi.growntilebox(IntVect(
NGROW,
NGROW,0));
232 ParallelFor(gbx2D, [=] AMREX_GPU_DEVICE (
int i,
int j,
int ) {
242 for ( MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi )
244 Array4<Real >
const& bustr = mf_bustr->array(mfi);
245 Array4<Real >
const& bvstr = mf_bvstr->array(mfi);
246 Array4<Real const>
const& uold = U_old.const_array(mfi);
247 Array4<Real const>
const& vold = V_old.const_array(mfi);
248 Array4<Real >
const& logdrg_tmp = mf_logdrg_tmp.array(mfi);
249 Array4<Real const>
const& z_r = mf_z_r->const_array(mfi);
250 Array4<Real const>
const& z_w = mf_z_w->const_array(mfi);
252 Box gbx2 = mfi.growntilebox(IntVect(
NGROW,
NGROW,0));
255 Box ubx1 = mfi.grownnodaltilebox(0,IntVect(
NGROW-1,
NGROW-1,0));
258 Box vbx1 = mfi.grownnodaltilebox(1,IntVect(
NGROW-1,
NGROW-1,0));
263 Array4<Real const>
const& rdrag = mf_rdrag->const_array(mfi);
264 ParallelFor(ubx1D, [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
266 bustr(i,j,0) = 0.5_rt * (rdrag(i-1,j,0)+rdrag(i,j,0))*(uold(i,j,0));
268 ParallelFor(vbx1D, [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
270 bvstr(i,j,0) = 0.5_rt * (rdrag(i,j-1,0)+rdrag(i,j,0))*(vold(i,j,0));
273 Array4<Real const>
const& rdrag2 = mf_rdrag2->const_array(mfi);
274 ParallelFor(ubx1D, [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
276 Real avg_v = 0.25_rt * (vold(i,j,0) + vold(i,j+1,0) + vold(i-1,j,0) + vold(i-1,j+1,0));
277 Real vel_mag = std::sqrt(uold(i,j,0)*uold(i,j,0) + avg_v * avg_v);
278 bustr(i,j,0) = 0.5_rt * (rdrag2(i-1,j,0) + rdrag2(i,j,0)) * uold(i,j,0) * vel_mag;
280 ParallelFor(vbx1D, [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
282 Real avg_u = 0.25_rt * (uold(i,j,0) + uold(i+1,j,0) + uold(i,j-1,0) + uold(i+1,j-1,0));
283 Real vel_mag = std::sqrt(avg_u * avg_u + vold(i,j,0) * vold(i,j,0));
284 bvstr(i,j,0) = 0.5_rt * (rdrag2(i,j-1,0) + rdrag2(i,j,0)) * vold(i,j,0) * vel_mag;
287 Array4<Real const>
const& ZoBot = mf_ZoBot->const_array(mfi);
288 ParallelFor(gbx2D, [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
290 Real logz = 1.0_rt / (std::log((z_r(i,j,0) - z_w(i,j,0)) / ZoBot(i,j,0)));
292 logdrg_tmp(i,j,0) = std::min(Cdb_max,std::max(Cdb_min,cff));
294 ParallelFor(ubx1D, [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
296 Real avg_v = 0.25_rt * (vold(i,j,0) + vold(i,j+1,0) + vold(i-1,j,0) + vold(i-1,j+1,0));
297 Real vel_mag = std::sqrt(uold(i,j,0)*uold(i,j,0) + avg_v * avg_v);
298 bustr(i,j,0) = 0.5_rt * (logdrg_tmp(i-1,j,0)+logdrg_tmp(i,j,0)) * uold(i,j,0) * vel_mag;
300 ParallelFor(vbx1D, [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
302 Real avg_u = 0.25_rt * (uold(i,j,0) + uold(i+1,j,0) + uold(i,j-1,0) + uold(i+1,j-1,0));
303 Real vel_mag = std::sqrt(avg_u * avg_u + vold(i,j,0) * vold(i,j,0));
304 bvstr(i,j,0) = 0.5_rt * (logdrg_tmp(i,j-1,0) + logdrg_tmp(i,j,0)) * vold(i,j,0) * vel_mag;
318 MultiFab mf_W(convert(ba,IntVect(0,0,1)),dm,1,IntVect(
NGROW+1,
NGROW+1,0));
324 prestep(lev, U_old, V_old, U_new, V_new,
325 mf_ru.get(), mf_rv.get(),
327 mf_DC, mf_z_r, mf_z_w, mf_h, mf_pm, mf_pn,
328 mf_sustr.get(), mf_svstr.get(), mf_bustr.get(), mf_bvstr.get(),
329 mf_msku.get(), mf_mskv.get(),
330 iic, ntfirst, nnew, nstp, nrhs, N, dt_lev);
334 mf_W.FillBoundary(geom[lev].periodicity());
337#ifdef REMORA_USE_NETCDF
345 for ( MFIter mfi(S_old, TilingIfNotGPU()); mfi.isValid(); ++mfi )
347 Array4<Real const>
const& Hz =
vec_Hz[lev]->const_array(mfi);
348 Array4<Real>
const& Huon =
vec_Huon[lev]->array(mfi);
349 Array4<Real>
const& Hvom =
vec_Hvom[lev]->array(mfi);
350 Array4<Real>
const& z_r = (mf_z_r)->array(mfi);
351 Array4<Real>
const& z_w = (mf_z_w)->array(mfi);
352 Array4<Real const>
const& uold = U_old.const_array(mfi);
353 Array4<Real const>
const& vold = V_old.const_array(mfi);
354 Array4<Real>
const& u = U_new.array(mfi);
355 Array4<Real>
const& v = V_new.array(mfi);
356 Array4<Real>
const& rho = (mf_rho).array(mfi);
357 Array4<Real>
const& ru = (mf_ru)->array(mfi);
358 Array4<Real>
const& rv = (mf_rv)->array(mfi);
359 Array4<Real>
const& rufrc = (mf_rufrc)->array(mfi);
360 Array4<Real>
const& rvfrc = (mf_rvfrc)->array(mfi);
361 Array4<Real>
const& W = (mf_W).array(mfi);
362 Array4<Real>
const& sustr = (mf_sustr)->array(mfi);
363 Array4<Real>
const& svstr = (mf_svstr)->array(mfi);
364 Array4<Real>
const& bustr = (mf_bustr)->array(mfi);
365 Array4<Real>
const& bvstr = (mf_bvstr)->array(mfi);
366 Array4<Real>
const& visc2_p = (mf_visc2_p)->array(mfi);
367 Array4<Real>
const& visc2_r = (mf_visc2_r)->array(mfi);
369 Array4<Real const>
const& pm = mf_pm->const_array(mfi);
370 Array4<Real const>
const& pn = mf_pn->const_array(mfi);
371 Array4<Real const>
const& fcor = mf_fcor->const_array(mfi);
373 Array4<Real const>
const& msku = mf_msku->const_array(mfi);
374 Array4<Real const>
const& mskv = mf_mskv->const_array(mfi);
375 Array4<Real const>
const& mskp = mf_mskp->const_array(mfi);
377 Box bx = mfi.tilebox();
381 Box xbx = mfi.nodaltilebox(0);
382 Box ybx = mfi.nodaltilebox(1);
383 Box xbx_adj = mfi.nodaltilebox(0);
384 Box ybx_adj = mfi.nodaltilebox(1);
386 auto xbx_lo = lbound(xbx_adj);
387 auto xbx_hi = ubound(xbx_adj);
388 auto ybx_lo = lbound(ybx_adj);
389 auto ybx_hi = ubound(ybx_adj);
391 if (xbx_lo.x == dlo.x) {
392 xbx_adj.growLo(0,-1);
393 }
else if (xbx_hi.x == dhi.x) {
394 xbx_adj.growHi(0,-1);
397 if (ybx_lo.y == dlo.y) {
398 ybx_adj.growLo(1,-1);
399 }
else if (ybx_hi.y == dhi.y) {
400 ybx_adj.growHi(1,-1);
403 Box gbx1 = mfi.growntilebox(IntVect(
NGROW-1,
NGROW-1,0));
404 Box gbx2 = mfi.growntilebox(IntVect(
NGROW,
NGROW,0));
406 Box utbx = mfi.nodaltilebox(0);
407 Box vtbx = mfi.nodaltilebox(1);
420 tbxp1D.makeSlab(2,0);
422 tbxp2D.makeSlab(2,0);
424 FArrayBox fab_FC(surroundingNodes(tbxp2,2),1,amrex::The_Async_Arena());
425 FArrayBox fab_FX(gbx2,1,amrex::The_Async_Arena());
426 FArrayBox fab_FE(gbx2,1,amrex::The_Async_Arena());
427 FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
428 FArrayBox fab_CF(gbx2,1,amrex::The_Async_Arena());
430 FArrayBox fab_fomn(tbxp2D,1,amrex::The_Async_Arena());
432 auto FC=fab_FC.array();
434 auto fomn=fab_fomn.array();
437 ParallelFor(tbxp2D, [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
439 fomn(i,j,0) = fcor(i,j,0)*(1.0_rt/(pm(i,j,0)*pn(i,j,0)));
443 ParallelFor(gbx2, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
448 prsgrd(tbxp1,gbx1,utbx,vtbx,ru,rv,pn,pm,rho,FC,Hz,z_r,z_w,msku,mskv,nrhs,N);
452 Array4<Real>
const& s_arr = S_new.array(mfi);
453 Array4<Real>
const& s_arr_rhs = S_old.array(mfi);
454 Array4<Real>
const& diff2_arr =
vec_diff2[lev]->array(mfi);
456 t3dmix(bx, s_arr, s_arr_rhs, diff2_arr, Hz, pm, pn, msku, mskv, dt_lev, ncomp);
459 t3dmix(bx, S_new.array(mfi,
Scalar_comp), S_old.array(mfi,
Scalar_comp), diff2_arr_scalar, Hz, pm, pn, msku, mskv, dt_lev, 1);
468 coriolis(xbx, ybx, uold, vold, ru, rv, Hz, fomn, nrhs, nrhs);
471#ifdef REMORA_USE_NETCDF
478 apply_clim_nudg(xbx_adj, 1, 0, ru, uold, uclim, u_nudg_coeff, Hz, pm, pn);
479 apply_clim_nudg(ybx_adj, 0, 1, rv, vold, vclim, v_nudg_coeff, Hz, pm, pn);
486 rhs_uv_3d(xbx, ybx, uold, vold, ru, rv, rufrc, rvfrc,
487 sustr, svstr, bustr, bvstr, Huon, Hvom,
488 pm, pn, W, FC, nrhs, N);
492 uv3dmix(xbx, ybx, u, v, uold, vold, rufrc, rvfrc, visc2_p, visc2_r, Hz, pm, pn, mskp, nrhs, nnew, dt_lev);
496 int nnew = (iic +1)% 2;
499 gls_prestep(lev, mf_gls, mf_tke, mf_W, mf_msku.get(), mf_mskv.get(),
500 nstp, nnew, iic, ntfirst, N, dt_lev);
507 FillPatch(lev, time, *
vec_sstore[lev], GetVecOfPtrs(
vec_sstore),
BCVars::cons_bc,
BdyVars::t,0,
true,
true,0,0,dt_lev,*
cons_old[lev]);
510 vec_Huon[lev]->FillBoundary(geom[lev].periodicity());
511 vec_Hvom[lev]->FillBoundary(geom[lev].periodicity());