ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_TI_no_substep_fun.H
Go to the documentation of this file.
1 /**
2  * Wrapper for advancing the solution with the slow RHS in the absence of acoustic substepping
3  */
4  auto no_substep_fun = [&,bogus_large_value_d=bogus_large_value](Vector<MultiFab>& S_sum,
5  Vector<MultiFab>& S_old,
6  Vector<MultiFab>& F_slow,
7  const double time_for_fp, const double slow_dt_d,
8  const int nrk)
9  {
10  BL_PROFILE("no_substep_fun");
11  amrex::ignore_unused(nrk);
12  int n_data = IntVars::NumTypes;
13 
14  amrex::Real slow_dt = static_cast<Real>(slow_dt_d);
15 
16  const auto& dxInv = fine_geom.InvCellSizeArray();
17 
18  const amrex::GpuArray<int, IntVars::NumTypes> scomp_fast = {0,0,0,0};
19  const amrex::GpuArray<int, IntVars::NumTypes> ncomp_fast = {2,1,1,1};
20 
21  if (verbose) amrex::Print() << " No-substepping time integration at level " << level
22  << std::setprecision(timeprecision)
23  << " to " << time_for_fp
24  << " with dt = " << slow_dt << std::endl;
25 
26  // Update S_sum = S_stage only for the fast variables
27 #ifdef _OPENMP
28 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
29 #endif
30  {
31  //
32  // NOTE: we must not tile if the terrain is moving because the update of zmom at
33  // k = 0 below calls WFromOmega, which reads the x- and y-momenta at
34  // (i+1,j) and (i,j+1). Those faces are written on nodaltilebox(0) and
35  // nodaltilebox(1), which exclude a tile's upper faces, so with tiling
36  // they would still hold the values from the previous stage -- and, under
37  // OpenMP, would be written by another thread at the same time.
38  // (The substepping version in ERF_Substep_MT.cpp leaves tiling off for
39  // this same reason.)
40  //
41  const bool do_tiling = (solverChoice.terrain_type == TerrainType::MovingFittedMesh)
42  ? false : TilingIfNotGPU();
43 
44  for ( MFIter mfi(S_sum[IntVars::cons],do_tiling); mfi.isValid(); ++mfi)
45  {
46  const Box bx = mfi.tilebox();
47  Box tbx = mfi.nodaltilebox(0);
48  Box tby = mfi.nodaltilebox(1);
49  Box tbz = mfi.nodaltilebox(2);
50 
51  Vector<Array4<Real> > ssum_h(n_data);
52  Vector<Array4<Real> > sold_h(n_data);
53  Vector<Array4<Real> > fslow_h(n_data);
54 
55  for (int i = 0; i < n_data; ++i) {
56  ssum_h[i] = S_sum[i].array(mfi);
57  sold_h[i] = S_old[i].array(mfi);
58  fslow_h[i] = F_slow[i].array(mfi);
59  }
60 
61  Gpu::AsyncVector<Array4<Real> > sold_d(n_data);
62  Gpu::AsyncVector<Array4<Real> > ssum_d(n_data);
63  Gpu::AsyncVector<Array4<Real> > fslow_d(n_data);
64 
65  Gpu::copy(Gpu::hostToDevice, sold_h.begin(), sold_h.end(), sold_d.begin());
66  Gpu::copy(Gpu::hostToDevice, ssum_h.begin(), ssum_h.end(), ssum_d.begin());
67  Gpu::copy(Gpu::hostToDevice, fslow_h.begin(), fslow_h.end(), fslow_d.begin());
68 
69  Array4<Real>* sold = sold_d.dataPtr();
70  Array4<Real>* ssum = ssum_d.dataPtr();
71  Array4<Real>* fslow = fslow_d.dataPtr();
72 
73  // Moving terrain
74  if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
75  {
76  const Array4<const Real>& dJ_old = detJ_cc[level]->const_array(mfi);
77  const Array4<const Real>& dJ_new = detJ_cc_new[level]->const_array(mfi);
78  const Array4<const Real>& dJ_stg = detJ_cc_src[level]->const_array(mfi);
79 
80  const Array4<const Real>& z_nd_old = z_phys_nd[level]->const_array(mfi);
81  const Array4<const Real>& z_nd_new = z_phys_nd_new[level]->const_array(mfi);
82 
83  const Array4<const Real>& mf_u = mapfac[level][MapFacType::u_x]->const_array(mfi);
84  const Array4<const Real>& mf_v = mapfac[level][MapFacType::v_y]->const_array(mfi);
85 
86  const Array4<Real >& z_t_arr = z_t_rk[level]->array(mfi);
87 
88  // We have already scaled the slow source term to have the extra factor of dJ
89  ParallelFor(bx, ncomp_fast[IntVars::cons],
90  [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) {
91  const int n = scomp_fast[IntVars::cons] + nn;
92  ssum[IntVars::cons](i,j,k,n) = dJ_old(i,j,k) * sold[IntVars::cons](i,j,k,n)
93  + slow_dt * dJ_stg(i,j,k) * fslow[IntVars::cons](i,j,k,n);
94  ssum[IntVars::cons](i,j,k,n) /= dJ_new(i,j,k);
95  });
96 
97  // We have already scaled the slow source term to have the extra factor of dJ
98  ParallelFor(tbx, tby,
99  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
100  Real h_zeta_old = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd_old);
101  Real h_zeta_new = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd_new);
102  ssum[IntVars::xmom](i,j,k) = ( h_zeta_old * sold[IntVars::xmom](i,j,k)
103  + slow_dt * fslow[IntVars::xmom](i,j,k) ) / h_zeta_new;
104  },
105  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
106  Real h_zeta_old = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd_old);
107  Real h_zeta_new = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd_new);
108  ssum[IntVars::ymom](i,j,k) = ( h_zeta_old * sold[IntVars::ymom](i,j,k)
109  + slow_dt * fslow[IntVars::ymom](i,j,k) ) / h_zeta_new;
110  });
111  ParallelFor(tbz,
112  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
113  if (k == 0) {
114  // Here we take advantage of the fact that moving terrain has a slip wall
115  // so we can just use the new value at (i,j,0).
116  Real rho_on_face = ssum[IntVars::cons](i,j,k,Rho_comp);
117  ssum[IntVars::zmom](i,j,k) = WFromOmega(i,j,k,rho_on_face*z_t_arr(i,j,k),
118  ssum[IntVars::xmom], ssum[IntVars::ymom],
119  mf_u, mf_v, z_nd_new,dxInv);
120  } else {
121  Real dJ_old_kface = myhalf * (dJ_old(i,j,k) + dJ_old(i,j,k-1));
122  Real dJ_new_kface = myhalf * (dJ_new(i,j,k) + dJ_new(i,j,k-1));
123  ssum[IntVars::zmom](i,j,k) = ( dJ_old_kface * sold[IntVars::zmom](i,j,k)
124  + slow_dt * fslow[IntVars::zmom](i,j,k) ) / dJ_new_kface;
125  }
126  });
127 
128  } else { // Fixed or no terrain
129  const Array4<const Real>& dJ_old = detJ_cc[level]->const_array(mfi);
130  ParallelFor(bx, ncomp_fast[IntVars::cons],
131  [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) {
132  const int n = scomp_fast[IntVars::cons] + nn;
133  if (dJ_old(i,j,k) > zero) {
134  ssum[IntVars::cons](i,j,k,n) = sold[IntVars::cons](i,j,k,n) + slow_dt *
135  ( fslow[IntVars::cons](i,j,k,n) );
136  } else {
137  ssum[IntVars::cons](i,j,k,n) = sold[IntVars::cons](i,j,k,n);
138  }
139  });
140 
141  // Commenting out the update is a HACK while developing the EB capability
142  if (solverChoice.terrain_type == TerrainType::EB)
143  {
144  const Array4<const Real>& vfrac_u = (get_eb(level).get_u_const_factory())->getVolFrac().const_array(mfi);
145  const Array4<const Real>& vfrac_v = (get_eb(level).get_v_const_factory())->getVolFrac().const_array(mfi);
146  const Array4<const Real>& vfrac_w = (get_eb(level).get_w_const_factory())->getVolFrac().const_array(mfi);
147 
148  ParallelFor(tbx, tby, tbz,
149  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
150  if (vfrac_u(i,j,k) > zero) {
151  // ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k);
152  ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k)
153  + slow_dt * fslow[IntVars::xmom](i,j,k);
154  } else {
155  ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k);
156  }
157  },
158  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
159  if (vfrac_v(i,j,k) > zero) {
160  // ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k);
161  ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k)
162  + slow_dt * fslow[IntVars::ymom](i,j,k);
163  } else {
164  ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k);
165  }
166  },
167  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
168  if (vfrac_w(i,j,k) > zero) {
169  // ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k);
170  ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k)
171  + slow_dt * fslow[IntVars::zmom](i,j,k);
172  } else {
173  ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k);
174  }
175  });
176 
177  } else {
178  ParallelFor(tbx, tby, tbz,
179  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
180  ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k)
181  + slow_dt * fslow[IntVars::xmom](i,j,k);
182  },
183  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
184  ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k)
185  + slow_dt * fslow[IntVars::ymom](i,j,k);
186  },
187  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
188  ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k)
189  + slow_dt * fslow[IntVars::zmom](i,j,k);
190  });
191  } // no EB
192  } // not moving terrain
193  } // mfi
194 
195  if (solverChoice.terrain_type == TerrainType::EB)
196  {
197  BL_PROFILE("no_substep_fun_redistribute");
198  // Redistribute cons states (cell-centered)
199 
200  Vector<MultiFab> dUdt(IntVars::NumTypes);
201  dUdt[IntVars::cons].define(ba, dm, F_slow[IntVars::cons].nComp(), F_slow[IntVars::cons].nGrow(), MFInfo(), EBFactory(level));
202  for (int i = 1; i <= AMREX_SPACEDIM; ++i) {
203  dUdt[i].define(F_slow[i].boxArray(), F_slow[i].DistributionMap(), F_slow[i].nComp(), F_slow[i].nGrow(), MFInfo());
204  }
205  for (int i = 0; i <= AMREX_SPACEDIM; ++i) {
206  dUdt[i].setVal(0, 0, ncomp_fast[i], dUdt[i].nGrow());
207  MultiFab::Copy(dUdt[i], F_slow[i], 0, 0, F_slow[i].nComp(), 0);
208  dUdt[i].setDomainBndry(bogus_large_value_d, 0, ncomp_fast[i], fine_geom);
209  dUdt[i].FillBoundary(fine_geom.periodicity());
210  }
211 
212  BCRec const* bc_ptr_d = domain_bcs_type_d.data();
213  IntVect ngv_pert = S_old[IntVars::cons].nGrowVect();
214 
215  MultiFab S_pert_cons(ba, dm, S_old[IntVars::cons].nComp(), ngv_pert, MFInfo(), EBFactory(level));
216  MultiFab rt0(ba, dm, 1, ngv_pert, MFInfo(), EBFactory(level));
217 
218  // Copy S_old into S_pert
219  MultiFab::Copy(S_pert_cons, S_old[IntVars::cons],0,0,S_pert_cons.nComp(), ngv_pert);
220 
221  // Copy theta0 into rt0
222  MultiFab::Copy(rt0,base_state[level],BaseState::th0_comp,0,1,ngv_pert);
223 
224  // Multiply theta0 by rho0 in rt0
225  MultiFab::Multiply(rt0,base_state[level],BaseState::r0_comp,0,1,ngv_pert);
226 
227  MultiFab::Subtract(S_pert_cons, base_state[level], BaseState::r0_comp,0,1,ngv_pert);
228  MultiFab::Subtract(S_pert_cons, rt0 , 0,1,1,ngv_pert);
229 
230  // Update F_slow for perturbational fast quantities (rho - rho0) and (rho_theta - rho_theta_0)
231  redistribute_term(ncomp_fast[IntVars::cons], fine_geom, F_slow[IntVars::cons], dUdt[IntVars::cons],
232  S_pert_cons, EBFactory(level), bc_ptr_d, slow_dt);
233 
234  // Update F_slow for momenta
235  redistribute_term(ncomp_fast[IntVars::xmom], fine_geom, F_slow[IntVars::xmom], dUdt[IntVars::xmom],
236  S_old[IntVars::xmom], *(get_eb(level).get_u_const_factory()), bc_ptr_d, slow_dt, IntVars::xmom);
237  redistribute_term(ncomp_fast[IntVars::ymom], fine_geom, F_slow[IntVars::ymom], dUdt[IntVars::ymom],
238  S_old[IntVars::ymom], *(get_eb(level).get_v_const_factory()), bc_ptr_d, slow_dt, IntVars::ymom);
239  redistribute_term(ncomp_fast[IntVars::zmom], fine_geom, F_slow[IntVars::zmom], dUdt[IntVars::zmom],
240  S_old[IntVars::zmom], *(get_eb(level).get_w_const_factory()), bc_ptr_d, slow_dt, IntVars::zmom);
241 
242 
243  // Update state using the updated F_slow.
244  for ( MFIter mfi(S_sum[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
245  {
246  const Box bx = mfi.tilebox();
247  Box tbx = mfi.nodaltilebox(0);
248  Box tby = mfi.nodaltilebox(1);
249  Box tbz = mfi.nodaltilebox(2);
250 
251  Vector<Array4<Real> > ssum_h(n_data);
252  Vector<Array4<Real> > sold_h(n_data);
253  Vector<Array4<Real> > fslow_h(n_data);
254 
255  for (int i = 0; i < n_data; ++i) {
256  ssum_h[i] = S_sum[i].array(mfi);
257  sold_h[i] = S_old[i].array(mfi);
258  fslow_h[i] = F_slow[i].array(mfi);
259  }
260 
261  Gpu::AsyncVector<Array4<Real> > sold_d(n_data);
262  Gpu::AsyncVector<Array4<Real> > ssum_d(n_data);
263  Gpu::AsyncVector<Array4<Real> > fslow_d(n_data);
264 
265  Gpu::copy(Gpu::hostToDevice, sold_h.begin(), sold_h.end(), sold_d.begin());
266  Gpu::copy(Gpu::hostToDevice, ssum_h.begin(), ssum_h.end(), ssum_d.begin());
267  Gpu::copy(Gpu::hostToDevice, fslow_h.begin(), fslow_h.end(), fslow_d.begin());
268 
269  Array4<Real>* sold = sold_d.dataPtr();
270  Array4<Real>* ssum = ssum_d.dataPtr();
271  Array4<Real>* fslow = fslow_d.dataPtr();
272 
273  // const Array4<const Real>& vfrac_c = detJ_cc[level]->const_array(mfi);
274  const Array4<const Real>& vfrac_c = (get_eb(level).get_const_factory())->getVolFrac().const_array(mfi);
275 
276  ParallelFor(bx, ncomp_fast[IntVars::cons], [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn)
277  {
278  const int n = scomp_fast[IntVars::cons] + nn;
279  if (vfrac_c(i,j,k) > zero) {
280  ssum[IntVars::cons](i,j,k,n) = sold[IntVars::cons](i,j,k,n) + slow_dt *
281  ( fslow[IntVars::cons](i,j,k,n) );
282  }
283  });
284 
285  const Array4<const Real>& vfrac_u = (get_eb(level).get_u_const_factory())->getVolFrac().const_array(mfi);
286  const Array4<const Real>& vfrac_v = (get_eb(level).get_v_const_factory())->getVolFrac().const_array(mfi);
287  const Array4<const Real>& vfrac_w = (get_eb(level).get_w_const_factory())->getVolFrac().const_array(mfi);
288 
289  ParallelFor(tbx, tby, tbz,
290  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
291  if (vfrac_u(i,j,k) > zero) {
292  ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k)
293  + slow_dt * fslow[IntVars::xmom](i,j,k);
294  }
295  },
296  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
297  if (vfrac_v(i,j,k) > zero) {
298  ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k)
299  + slow_dt * fslow[IntVars::ymom](i,j,k);
300  }
301  },
302  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
303  if (vfrac_w(i,j,k) > zero) {
304  ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k)
305  + slow_dt * fslow[IntVars::zmom](i,j,k);
306 
307  }
308  });
309  } // MFIter
310  } // EB
311  } // omp
312 
313  // Even if we update all the conserved variables we don't need
314  // to fillpatch the slow ones every acoustic substep
315  apply_bcs(S_sum, time_for_fp, S_sum[IntVars::cons].nGrow(), S_sum[IntVars::xmom].nGrow(),
316  fast_only=true, vel_and_mom_synced=false);
317 
318  if (solverChoice.anelastic[level]) {
319  bool have_tb = (thin_xforce[0] || thin_yforce[0] || thin_zforce[0]);
320  if (have_tb) {
321  project_velocity_tb(level, slow_dt, S_sum);
322  } else {
323  project_momenta(level, time_for_fp, slow_dt, S_sum);
324  }
325  }
326  };
bogus_large_value
constexpr amrex::Real bogus_large_value
Definition: ERF_Constants.H:26
zero
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
myhalf
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
v_y
@ v_y
Definition: ERF_DataStruct.H:29
u_x
@ u_x
Definition: ERF_DataStruct.H:28
redistribute_term
void redistribute_term(int ncomp, const Geometry &geom, MultiFab &result, MultiFab &result_tmp, MultiFab const &state, EBFactType const &ebfact, BCRec const *bc, double local_dt_d, int const igrid)
Apply EB state redistribution to result_tmp and write the redistributed result.
Definition: ERF_EBRedistribute.cpp:21
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
ParallelFor
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
no_substep_fun
auto no_substep_fun
Definition: ERF_TI_no_substep_fun.H:4
apply_bcs
auto apply_bcs
Definition: ERF_TI_utils.H:34
Compute_h_zeta_AtIface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtIface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:258
Compute_h_zeta_AtJface
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_h_zeta_AtJface(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:328
WFromOmega
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real WFromOmega(int &i, int &j, int &k, amrex::Real omega, const amrex::Array4< const amrex::Real > &u_arr, const amrex::Array4< const amrex::Real > &v_arr, const amrex::Array4< const amrex::Real > &mf_u, const amrex::Array4< const amrex::Real > &mf_v, const amrex::Array4< const amrex::Real > &z_nd, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dxInv)
Definition: ERF_TerrainMetrics.H:845
BaseState::th0_comp
@ th0_comp
Definition: ERF_IndexDefines.H:79
BaseState::r0_comp
@ r0_comp
Definition: ERF_IndexDefines.H:76
IntVars::NumTypes
@ NumTypes
Definition: ERF_IndexDefines.H:236
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:234
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:232
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:235
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:233
MicVar_WDM6::nn
@ nn
Definition: ERF_WDM6.H:31