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  for ( MFIter mfi(S_sum[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
32  {
33  const Box bx = mfi.tilebox();
34  Box tbx = mfi.nodaltilebox(0);
35  Box tby = mfi.nodaltilebox(1);
36  Box tbz = mfi.nodaltilebox(2);
37 
38  Vector<Array4<Real> > ssum_h(n_data);
39  Vector<Array4<Real> > sold_h(n_data);
40  Vector<Array4<Real> > fslow_h(n_data);
41 
42  for (int i = 0; i < n_data; ++i) {
43  ssum_h[i] = S_sum[i].array(mfi);
44  sold_h[i] = S_old[i].array(mfi);
45  fslow_h[i] = F_slow[i].array(mfi);
46  }
47 
48  Gpu::AsyncVector<Array4<Real> > sold_d(n_data);
49  Gpu::AsyncVector<Array4<Real> > ssum_d(n_data);
50  Gpu::AsyncVector<Array4<Real> > fslow_d(n_data);
51 
52  Gpu::copy(Gpu::hostToDevice, sold_h.begin(), sold_h.end(), sold_d.begin());
53  Gpu::copy(Gpu::hostToDevice, ssum_h.begin(), ssum_h.end(), ssum_d.begin());
54  Gpu::copy(Gpu::hostToDevice, fslow_h.begin(), fslow_h.end(), fslow_d.begin());
55 
56  Array4<Real>* sold = sold_d.dataPtr();
57  Array4<Real>* ssum = ssum_d.dataPtr();
58  Array4<Real>* fslow = fslow_d.dataPtr();
59 
60  // Moving terrain
61  if ( solverChoice.terrain_type == TerrainType::MovingFittedMesh )
62  {
63  const Array4<const Real>& dJ_old = detJ_cc[level]->const_array(mfi);
64  const Array4<const Real>& dJ_new = detJ_cc_new[level]->const_array(mfi);
65  const Array4<const Real>& dJ_stg = detJ_cc_src[level]->const_array(mfi);
66 
67  const Array4<const Real>& z_nd_old = z_phys_nd[level]->const_array(mfi);
68  const Array4<const Real>& z_nd_new = z_phys_nd_new[level]->const_array(mfi);
69 
70  const Array4<const Real>& mf_u = mapfac[level][MapFacType::u_x]->const_array(mfi);
71  const Array4<const Real>& mf_v = mapfac[level][MapFacType::v_y]->const_array(mfi);
72 
73  const Array4<Real >& z_t_arr = z_t_rk[level]->array(mfi);
74 
75  // We have already scaled the slow source term to have the extra factor of dJ
76  ParallelFor(bx, ncomp_fast[IntVars::cons],
77  [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) {
78  const int n = scomp_fast[IntVars::cons] + nn;
79  ssum[IntVars::cons](i,j,k,n) = dJ_old(i,j,k) * sold[IntVars::cons](i,j,k,n)
80  + slow_dt * dJ_stg(i,j,k) * fslow[IntVars::cons](i,j,k,n);
81  ssum[IntVars::cons](i,j,k,n) /= dJ_new(i,j,k);
82  });
83 
84  // We have already scaled the slow source term to have the extra factor of dJ
85  ParallelFor(tbx, tby,
86  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
87  Real h_zeta_old = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd_old);
88  Real h_zeta_new = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd_new);
89  ssum[IntVars::xmom](i,j,k) = ( h_zeta_old * sold[IntVars::xmom](i,j,k)
90  + slow_dt * fslow[IntVars::xmom](i,j,k) ) / h_zeta_new;
91  },
92  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
93  Real h_zeta_old = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd_old);
94  Real h_zeta_new = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd_new);
95  ssum[IntVars::ymom](i,j,k) = ( h_zeta_old * sold[IntVars::ymom](i,j,k)
96  + slow_dt * fslow[IntVars::ymom](i,j,k) ) / h_zeta_new;
97  });
98  ParallelFor(tbz,
99  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
100  if (k == 0) {
101  // Here we take advantage of the fact that moving terrain has a slip wall
102  // so we can just use the new value at (i,j,0).
103  Real rho_on_face = ssum[IntVars::cons](i,j,k,Rho_comp);
104  ssum[IntVars::zmom](i,j,k) = WFromOmega(i,j,k,rho_on_face*z_t_arr(i,j,k),
105  ssum[IntVars::xmom], ssum[IntVars::ymom],
106  mf_u, mf_v, z_nd_new,dxInv);
107  } else {
108  Real dJ_old_kface = myhalf * (dJ_old(i,j,k) + dJ_old(i,j,k-1));
109  Real dJ_new_kface = myhalf * (dJ_new(i,j,k) + dJ_new(i,j,k-1));
110  ssum[IntVars::zmom](i,j,k) = ( dJ_old_kface * sold[IntVars::zmom](i,j,k)
111  + slow_dt * fslow[IntVars::zmom](i,j,k) ) / dJ_new_kface;
112  }
113  });
114 
115  } else { // Fixed or no terrain
116  const Array4<const Real>& dJ_old = detJ_cc[level]->const_array(mfi);
117  ParallelFor(bx, ncomp_fast[IntVars::cons],
118  [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) {
119  const int n = scomp_fast[IntVars::cons] + nn;
120  if (dJ_old(i,j,k) > zero) {
121  ssum[IntVars::cons](i,j,k,n) = sold[IntVars::cons](i,j,k,n) + slow_dt *
122  ( fslow[IntVars::cons](i,j,k,n) );
123  } else {
124  ssum[IntVars::cons](i,j,k,n) = sold[IntVars::cons](i,j,k,n);
125  }
126  });
127 
128  // Commenting out the update is a HACK while developing the EB capability
129  if (solverChoice.terrain_type == TerrainType::EB)
130  {
131  const Array4<const Real>& vfrac_u = (get_eb(level).get_u_const_factory())->getVolFrac().const_array(mfi);
132  const Array4<const Real>& vfrac_v = (get_eb(level).get_v_const_factory())->getVolFrac().const_array(mfi);
133  const Array4<const Real>& vfrac_w = (get_eb(level).get_w_const_factory())->getVolFrac().const_array(mfi);
134 
135  ParallelFor(tbx, tby, tbz,
136  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
137  if (vfrac_u(i,j,k) > zero) {
138  // ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k);
139  ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k)
140  + slow_dt * fslow[IntVars::xmom](i,j,k);
141  } else {
142  ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k);
143  }
144  },
145  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
146  if (vfrac_v(i,j,k) > zero) {
147  // ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k);
148  ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k)
149  + slow_dt * fslow[IntVars::ymom](i,j,k);
150  } else {
151  ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k);
152  }
153  },
154  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
155  if (vfrac_w(i,j,k) > zero) {
156  // ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k);
157  ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k)
158  + slow_dt * fslow[IntVars::zmom](i,j,k);
159  } else {
160  ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k);
161  }
162  });
163 
164  } else {
165  ParallelFor(tbx, tby, tbz,
166  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
167  ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k)
168  + slow_dt * fslow[IntVars::xmom](i,j,k);
169  },
170  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
171  ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k)
172  + slow_dt * fslow[IntVars::ymom](i,j,k);
173  },
174  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
175  ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k)
176  + slow_dt * fslow[IntVars::zmom](i,j,k);
177  });
178  } // no EB
179  } // not moving terrain
180  } // mfi
181 
182  if (solverChoice.terrain_type == TerrainType::EB)
183  {
184  BL_PROFILE("no_substep_fun_redistribute");
185  // Redistribute cons states (cell-centered)
186 
187  Vector<MultiFab> dUdt(IntVars::NumTypes);
188  dUdt[IntVars::cons].define(ba, dm, F_slow[IntVars::cons].nComp(), F_slow[IntVars::cons].nGrow(), MFInfo(), EBFactory(level));
189  for (int i = 1; i <= AMREX_SPACEDIM; ++i) {
190  dUdt[i].define(F_slow[i].boxArray(), F_slow[i].DistributionMap(), F_slow[i].nComp(), F_slow[i].nGrow(), MFInfo());
191  }
192  for (int i = 0; i <= AMREX_SPACEDIM; ++i) {
193  dUdt[i].setVal(0, 0, ncomp_fast[i], dUdt[i].nGrow());
194  MultiFab::Copy(dUdt[i], F_slow[i], 0, 0, F_slow[i].nComp(), 0);
195  dUdt[i].setDomainBndry(bogus_large_value_d, 0, ncomp_fast[i], fine_geom);
196  dUdt[i].FillBoundary(fine_geom.periodicity());
197  }
198 
199  BCRec const* bc_ptr_d = domain_bcs_type_d.data();
200  IntVect ngv_pert = S_old[IntVars::cons].nGrowVect();
201 
202  MultiFab S_pert_cons(ba, dm, S_old[IntVars::cons].nComp(), ngv_pert, MFInfo(), EBFactory(level));
203  MultiFab rt0(ba, dm, 1, ngv_pert, MFInfo(), EBFactory(level));
204 
205  // Copy S_old into S_pert
206  MultiFab::Copy(S_pert_cons, S_old[IntVars::cons],0,0,S_pert_cons.nComp(), ngv_pert);
207 
208  // Copy theta0 into rt0
209  MultiFab::Copy(rt0,base_state[level],BaseState::th0_comp,0,1,ngv_pert);
210 
211  // Multiply theta0 by rho0 in rt0
212  MultiFab::Multiply(rt0,base_state[level],BaseState::r0_comp,0,1,ngv_pert);
213 
214  MultiFab::Subtract(S_pert_cons, base_state[level], BaseState::r0_comp,0,1,ngv_pert);
215  MultiFab::Subtract(S_pert_cons, rt0 , 0,1,1,ngv_pert);
216 
217  // Update F_slow for perturbational fast quantities (rho - rho0) and (rho_theta - rho_theta_0)
218  redistribute_term(ncomp_fast[IntVars::cons], fine_geom, F_slow[IntVars::cons], dUdt[IntVars::cons],
219  S_pert_cons, EBFactory(level), bc_ptr_d, slow_dt);
220 
221  // Update F_slow for momenta
222  redistribute_term(ncomp_fast[IntVars::xmom], fine_geom, F_slow[IntVars::xmom], dUdt[IntVars::xmom],
223  S_old[IntVars::xmom], *(get_eb(level).get_u_const_factory()), bc_ptr_d, slow_dt, IntVars::xmom);
224  redistribute_term(ncomp_fast[IntVars::ymom], fine_geom, F_slow[IntVars::ymom], dUdt[IntVars::ymom],
225  S_old[IntVars::ymom], *(get_eb(level).get_v_const_factory()), bc_ptr_d, slow_dt, IntVars::ymom);
226  redistribute_term(ncomp_fast[IntVars::zmom], fine_geom, F_slow[IntVars::zmom], dUdt[IntVars::zmom],
227  S_old[IntVars::zmom], *(get_eb(level).get_w_const_factory()), bc_ptr_d, slow_dt, IntVars::zmom);
228 
229 
230  // Update state using the updated F_slow.
231  for ( MFIter mfi(S_sum[IntVars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi)
232  {
233  const Box bx = mfi.tilebox();
234  Box tbx = mfi.nodaltilebox(0);
235  Box tby = mfi.nodaltilebox(1);
236  Box tbz = mfi.nodaltilebox(2);
237 
238  Vector<Array4<Real> > ssum_h(n_data);
239  Vector<Array4<Real> > sold_h(n_data);
240  Vector<Array4<Real> > fslow_h(n_data);
241 
242  for (int i = 0; i < n_data; ++i) {
243  ssum_h[i] = S_sum[i].array(mfi);
244  sold_h[i] = S_old[i].array(mfi);
245  fslow_h[i] = F_slow[i].array(mfi);
246  }
247 
248  Gpu::AsyncVector<Array4<Real> > sold_d(n_data);
249  Gpu::AsyncVector<Array4<Real> > ssum_d(n_data);
250  Gpu::AsyncVector<Array4<Real> > fslow_d(n_data);
251 
252  Gpu::copy(Gpu::hostToDevice, sold_h.begin(), sold_h.end(), sold_d.begin());
253  Gpu::copy(Gpu::hostToDevice, ssum_h.begin(), ssum_h.end(), ssum_d.begin());
254  Gpu::copy(Gpu::hostToDevice, fslow_h.begin(), fslow_h.end(), fslow_d.begin());
255 
256  Array4<Real>* sold = sold_d.dataPtr();
257  Array4<Real>* ssum = ssum_d.dataPtr();
258  Array4<Real>* fslow = fslow_d.dataPtr();
259 
260  // const Array4<const Real>& vfrac_c = detJ_cc[level]->const_array(mfi);
261  const Array4<const Real>& vfrac_c = (get_eb(level).get_const_factory())->getVolFrac().const_array(mfi);
262 
263  ParallelFor(bx, ncomp_fast[IntVars::cons], [=] AMREX_GPU_DEVICE (int i, int j, int k, int nn)
264  {
265  const int n = scomp_fast[IntVars::cons] + nn;
266  if (vfrac_c(i,j,k) > zero) {
267  ssum[IntVars::cons](i,j,k,n) = sold[IntVars::cons](i,j,k,n) + slow_dt *
268  ( fslow[IntVars::cons](i,j,k,n) );
269  }
270  });
271 
272  const Array4<const Real>& vfrac_u = (get_eb(level).get_u_const_factory())->getVolFrac().const_array(mfi);
273  const Array4<const Real>& vfrac_v = (get_eb(level).get_v_const_factory())->getVolFrac().const_array(mfi);
274  const Array4<const Real>& vfrac_w = (get_eb(level).get_w_const_factory())->getVolFrac().const_array(mfi);
275 
276  ParallelFor(tbx, tby, tbz,
277  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
278  if (vfrac_u(i,j,k) > zero) {
279  ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k)
280  + slow_dt * fslow[IntVars::xmom](i,j,k);
281  }
282  },
283  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
284  if (vfrac_v(i,j,k) > zero) {
285  ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k)
286  + slow_dt * fslow[IntVars::ymom](i,j,k);
287  }
288  },
289  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
290  if (vfrac_w(i,j,k) > zero) {
291  ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k)
292  + slow_dt * fslow[IntVars::zmom](i,j,k);
293 
294  }
295  });
296  } // MFIter
297  } // EB
298  } // omp
299 
300  // Even if we update all the conserved variables we don't need
301  // to fillpatch the slow ones every acoustic substep
302  apply_bcs(S_sum, time_for_fp, S_sum[IntVars::cons].nGrow(), S_sum[IntVars::xmom].nGrow(),
303  fast_only=true, vel_and_mom_synced=false);
304 
305  if (solverChoice.anelastic[level]) {
306  bool have_tb = (thin_xforce[0] || thin_yforce[0] || thin_zforce[0]);
307  if (have_tb) {
308  project_velocity_tb(level, slow_dt, S_sum);
309  } else {
310  project_momenta(level, time_for_fp, slow_dt, S_sum);
311  }
312  }
313  };
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:25
u_x
@ u_x
Definition: ERF_DataStruct.H:24
redistribute_term
void redistribute_term(int ncomp, const Geometry &geom, MultiFab &result, MultiFab &result_tmp, MultiFab const &state, EBFArrayBoxFactory const &ebfact, BCRec const *bc, double local_dt_d)
Definition: ERF_EBRedistribute.cpp:13
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
dxInv
amrex::GpuArray< Real, AMREX_SPACEDIM > dxInv
Definition: ERF_InitCustomPertVels_ParticleTests.H:17
ParallelFor
ParallelFor(grown_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:104
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:144
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:464
BaseState::th0_comp
@ th0_comp
Definition: ERF_IndexDefines.H:76
BaseState::r0_comp
@ r0_comp
Definition: ERF_IndexDefines.H:73
IntVars::NumTypes
@ NumTypes
Definition: ERF_IndexDefines.H:197
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:195
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:193
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:196
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:194