ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
MRISplitIntegrator< T > Class Template Reference

#include <ERF_MRI.H>

Collaboration diagram for MRISplitIntegrator< T >:

Public Member Functions

 MRISplitIntegrator ()=default
 
 MRISplitIntegrator (const T &S_data)
 
void initialize (const T &S_data)
 
 ~MRISplitIntegrator ()=default
 
 MRISplitIntegrator (MRISplitIntegrator &&) noexcept=default
 
MRISplitIntegratoroperator= (MRISplitIntegrator &&other) noexcept=default
 
 MRISplitIntegrator (const MRISplitIntegrator &other)=delete
 
MRISplitIntegratoroperator= (const MRISplitIntegrator &other)=delete
 
void setNcompCons (int _ncomp_cons)
 
void setAnelastic (int _anelastic)
 
void setNoSubstepping (int _no_substepping)
 
void setForceFirstStageSingleSubstep (int _force_stage1_single_substep)
 
void set_slow_rhs_pre (std::function< void(T &, T &, T &, const double, const double, const double, const int)> F)
 
void set_slow_rhs_post (std::function< void(T &, T &, T &, T &, const double, const double, const double, const int)> F)
 
void set_acoustic_substepping (std::function< void(int, int, int, T &, const T &, T &, T &, const double, const double, const amrex::Real, const double, const double)> F)
 
void set_slow_fast_timestep_ratio (const int timestep_ratio=1)
 
int get_slow_fast_timestep_ratio ()
 
void set_no_substep (std::function< void(T &, T &, T &, const double, const double, int)> F)
 
std::function< void(T &, const T &, const double, const double)> get_rhs ()
 
double advance (T &S_old, T &S_new, double time, const double time_step)
 
void map_data (std::function< void(T &)> Map)
 

Private Member Functions

void initialize_data (const T &S_data)
 

Private Attributes

std::function< void(T &, const T &, const double, const double)> rhs
 rhs is the right-hand-side function the integrator will use. More...
 
std::function< void(T &, T &, T &, const double, const double, const double, const int)> slow_rhs_pre
 
std::function< void(T &, T &, T &, T &, const double, const double, const double, const int)> slow_rhs_post
 
std::function< void(int, int, int, T &, const T &, T &, T &, const double, const double, const amrex::Real, const double, const double)> acoustic_substepping
 
double timestep
 Integrator timestep size (Real) More...
 
int slow_fast_timestep_ratio = 0
 The ratio of slow timestep size / fast timestep size (int) More...
 
int no_substepping
 Should we not do acoustic substepping. More...
 
int anelastic
 Should we use the anelastic integrator. More...
 
int ncomp_cons
 How many components in the cell-centered MultiFab. More...
 
int force_stage1_single_substep
 Do we follow the recommendation to only perform a single substep in the first RK stage. More...
 
std::function< void(T &, T &, T &, const double, const double, int)> no_substep
 The no_substep function is called when we have no acoustic substepping. More...
 
amrex::Vector< std::unique_ptr< T > > T_store
 
TS_sum
 
TF_slow
 

Constructor & Destructor Documentation

◆ MRISplitIntegrator() [1/4]

template<class T >
MRISplitIntegrator< T >::MRISplitIntegrator ( )
default

◆ MRISplitIntegrator() [2/4]

template<class T >
MRISplitIntegrator< T >::MRISplitIntegrator ( const T S_data)
inline
87  {
88  initialize_data(S_data);
89  }
void initialize_data(const T &S_data)
Definition: ERF_MRI.H:69
Here is the call graph for this function:

◆ ~MRISplitIntegrator()

template<class T >
MRISplitIntegrator< T >::~MRISplitIntegrator ( )
default

◆ MRISplitIntegrator() [3/4]

template<class T >
MRISplitIntegrator< T >::MRISplitIntegrator ( MRISplitIntegrator< T > &&  )
defaultnoexcept

◆ MRISplitIntegrator() [4/4]

template<class T >
MRISplitIntegrator< T >::MRISplitIntegrator ( const MRISplitIntegrator< T > &  other)
delete

Member Function Documentation

◆ advance()

template<class T >
double MRISplitIntegrator< T >::advance ( T S_old,
T S_new,
double  time,
const double  time_step 
)
inline
173  {
174  BL_PROFILE_REGION("MRI_advance");
175  using namespace amrex;
176 
177  // *******************************************************************************
178  // !no_substepping: we only update the fast variables every fast timestep, then update
179  // the slow variables after the acoustic sub-stepping. This has
180  // two calls to slow_rhs so that we can update the slow variables
181  // with the velocity field after the acoustic substepping using
182  // the time-averaged velocity from the substepping
183  // no_substepping: we don't do any acoustic subcycling so we only make one call per RK
184  // stage to slow_rhs
185  // *******************************************************************************
186  timestep = time_step;
187 
188  const int substep_ratio = get_slow_fast_timestep_ratio();
189 
190  if (!no_substepping) {
191  AMREX_ALWAYS_ASSERT(substep_ratio > 1 && substep_ratio % 2 == 0);
192  }
193 
194  // Assume before advance() that S_old is valid data at the current time ("time" argument)
195  // And that if data is a MultiFab, both S_old and S_new contain ghost cells for evaluating a stencil based RHS
196  // We need this from S_old. This is convenient for S_new to have so we can use it
197  // as scratch space for stage values without creating a new scratch MultiFab with ghost cells.
198 
199  // NOTE: In the following, we use S_new to hold S*, S**, and finally, S^(n+1) at the new time
200  // DEFINITIONS:
201  // S_old = S^n
202  // S_sum = S(t)
203  // F_slow = F(S_stage)
204 
205  int n_data = IntVars::NumTypes;
206 
207  /**********************************************/
208  /* RK3 Integration with Acoustic Sub-stepping */
209  /**********************************************/
210  Vector<int> num_vars = {ncomp_cons, 1, 1, 1};
211  for (int i(0); i<n_data; ++i)
212  {
213  // Copy old -> new
214  MultiFab::Copy(S_new[i],S_old[i],0,0,num_vars[i],S_old[i].nGrowVect());
215  }
216 
217  // Timestep taken by the fast integrator
218  double dtau;
219 
220  // How many timesteps taken by the fast integrator
221  int nsubsteps;
222 
223  // This is the final time of the full timestep (also the 3rd RK stage)
224  // Real new_time = time + timestep;
225 
226  double time_stage = time;
227  double old_time_stage;
228 
229  if (!anelastic) {
230  // RK3 for compressible integrator
231  for (int nrk = 0; nrk < 3; nrk++)
232  {
233  // Capture the time we got to in the previous RK step
234  old_time_stage = time_stage;
235 
236  if (nrk == 0) {
238  nsubsteps = 1;
239  dtau = timestep / three;
240  } else {
241  // Clamp to 1: substep_ratio is only required to be even, so a ratio of 2
242  // would otherwise give zero substeps and leave S_sum stale in this stage
243  nsubsteps = std::max(substep_ratio/3, 1);
244  dtau = (timestep / three) / static_cast<double>(nsubsteps);
245  }
246  time_stage = time + timestep / three;
247  }
248  if (nrk == 1) {
249  if (no_substepping) {
250  nsubsteps = 1;
251  dtau = myhalf * timestep;
252  } else {
253  nsubsteps = substep_ratio/2;
254  dtau = (myhalf * timestep) / static_cast<double>(nsubsteps);
255  }
256  time_stage = time + timestep / two;
257  }
258  if (nrk == 2) {
259  if (no_substepping) {
260  nsubsteps = 1;
261  dtau = timestep;
262  } else {
263  nsubsteps = substep_ratio;
264  dtau = timestep / static_cast<double>(nsubsteps);
265  }
266  time_stage = time + timestep;
267  }
268 
269  // step 1 starts with S_stage = S^n and we always start substepping at the old time
270  // step 2 starts with S_stage = S^* and we always start substepping at the old time
271  // step 3 starts with S_stage = S^** and we always start substepping at the old time
272 
273  slow_rhs_pre(*F_slow, S_old, S_new, time, old_time_stage, time_stage, nrk);
274 
275  amrex::Real inv_fac = one / static_cast<amrex::Real>(nsubsteps);
276 
277  // ****************************************************
278  // Acoustic substepping
279  // ****************************************************
280  if (!no_substepping)
281  {
282  // *******************************************************************************
283  // Update the fast variables
284  // *******************************************************************************
285  for (int ks = 0; ks < nsubsteps; ++ks)
286  {
287  acoustic_substepping(ks, nsubsteps, nrk, *F_slow, S_old, S_new, *S_sum, dtau, timestep, inv_fac,
288  time + ks*dtau, time + (ks+1) * dtau);
289 
290  } // ks
291 
292  } else {
293  no_substep(*S_sum, S_old, *F_slow, time + nsubsteps*dtau, nsubsteps*dtau, nrk);
294  }
295 
296  // ****************************************************
297  // Evaluate F_slow(S_stage) only for the slow variables
298  // Note that we are using the current stage versions (in S_new) of the slow variables
299  // (because we didn't update the slow variables in the substepping)
300  // but we are using the "new" versions (in S_sum) of the velocities
301  // (because we did update the fast variables in the substepping)
302  // ****************************************************
303  slow_rhs_post(*F_slow, S_old, S_new, *S_sum, time, old_time_stage, time_stage, nrk);
304  } // nrk
305 
306  } else {
307  // RK2 for anelastic integrator
308  for (int nrk = 0; nrk < 2; nrk++)
309  {
310  // Capture the time we got to in the previous RK step
311  old_time_stage = time_stage;
312 
313  if (nrk == 0) { nsubsteps = 1; dtau = timestep; time_stage = time + timestep; }
314  if (nrk == 1) { nsubsteps = 1; dtau = timestep; time_stage = time + timestep; }
315 
316  slow_rhs_pre(*F_slow, S_old, S_new, time, old_time_stage, time_stage, nrk);
317 
318  no_substep(*S_sum, S_old, *F_slow, time + nsubsteps*dtau, nsubsteps*dtau, nrk);
319 
320  // ****************************************************
321  // Evaluate F_slow(S_stage) only for the slow variables
322  // Note that we are using the current stage versions (in S_new) of the slow variables
323  // (because we didn't update the slow variables in the substepping)
324  // but we are using the "new" versions (in S_sum) of the velocities
325  // (because we did update the fast variables in the substepping)
326  // ****************************************************
327  slow_rhs_post(*F_slow, S_old, S_new, *S_sum, time, old_time_stage, time_stage, nrk);
328  } // nrk
329  }
330 
331  // Return timestep
332  return timestep;
333  }
constexpr amrex::Real three
Definition: ERF_Constants.H:11
constexpr amrex::Real two
Definition: ERF_Constants.H:10
constexpr amrex::Real one
Definition: ERF_Constants.H:9
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
amrex::Real Real
Definition: ERF_ShocInterface.H:19
T * F_slow
Definition: ERF_MRI.H:67
int anelastic
Should we use the anelastic integrator.
Definition: ERF_MRI.H:47
int force_stage1_single_substep
Do we follow the recommendation to only perform a single substep in the first RK stage.
Definition: ERF_MRI.H:57
int ncomp_cons
How many components in the cell-centered MultiFab.
Definition: ERF_MRI.H:52
double timestep
Integrator timestep size (Real)
Definition: ERF_MRI.H:32
std::function< void(T &, T &, T &, const double, const double, const double, const int)> slow_rhs_pre
Definition: ERF_MRI.H:23
std::function< void(T &, T &, T &, T &, const double, const double, const double, const int)> slow_rhs_post
Definition: ERF_MRI.H:24
std::function< void(T &, T &, T &, const double, const double, int)> no_substep
The no_substep function is called when we have no acoustic substepping.
Definition: ERF_MRI.H:62
int get_slow_fast_timestep_ratio()
Definition: ERF_MRI.H:157
T * S_sum
Definition: ERF_MRI.H:66
std::function< void(int, int, int, T &, const T &, T &, T &, const double, const double, const amrex::Real, const double, const double)> acoustic_substepping
Definition: ERF_MRI.H:27
int no_substepping
Should we not do acoustic substepping.
Definition: ERF_MRI.H:42
@ NumTypes
Definition: ERF_IndexDefines.H:198
Definition: ERF_ConsoleIO.cpp:15

Referenced by ERF::advance_dycore().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_rhs()

template<class T >
std::function<void(T&, const T&, const double, const double)> MRISplitIntegrator< T >::get_rhs ( )
inline
168  {
169  return rhs;
170  }
std::function< void(T &, const T &, const double, const double)> rhs
rhs is the right-hand-side function the integrator will use.
Definition: ERF_MRI.H:22

◆ get_slow_fast_timestep_ratio()

template<class T >
int MRISplitIntegrator< T >::get_slow_fast_timestep_ratio ( )
inline
158  {
160  }
int slow_fast_timestep_ratio
The ratio of slow timestep size / fast timestep size (int)
Definition: ERF_MRI.H:37

Referenced by MRISplitIntegrator< T >::advance().

Here is the caller graph for this function:

◆ initialize()

template<class T >
void MRISplitIntegrator< T >::initialize ( const T S_data)
inline
92  {
93  initialize_data(S_data);
94  }
Here is the call graph for this function:

◆ initialize_data()

template<class T >
void MRISplitIntegrator< T >::initialize_data ( const T S_data)
inlineprivate
70  {
71  // TODO: We can optimize memory by making the cell-centered part of S_sum
72  // have only 2 components, not ncomp_cons components
73  const bool include_ghost = true;
74  amrex::IntegratorOps<T>::CreateLike(T_store, S_data, include_ghost);
75  S_sum = T_store[0].get();
76  amrex::IntegratorOps<T>::CreateLike(T_store, S_data, include_ghost);
77  F_slow = T_store[1].get();
78  // initializing to zero
79  for (long idx = 0; idx < S_sum->size(); idx++) { (*S_sum)[idx].setVal(0); }
80  for (long idx = 0; idx < F_slow->size(); idx++) { (*F_slow)[idx].setVal(0); }
81  }
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int idx(int i, int j, int k, int nx, int ny)
Definition: ERF_InitForEnsemble.cpp:365
amrex::Vector< std::unique_ptr< T > > T_store
Definition: ERF_MRI.H:65

Referenced by MRISplitIntegrator< T >::initialize(), and MRISplitIntegrator< T >::MRISplitIntegrator().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ map_data()

template<class T >
void MRISplitIntegrator< T >::map_data ( std::function< void(T &)>  Map)
inline
336  {
337  for (auto& F : T_store) {
338  Map(*F);
339  }
340  }

◆ operator=() [1/2]

template<class T >
MRISplitIntegrator& MRISplitIntegrator< T >::operator= ( const MRISplitIntegrator< T > &  other)
delete

◆ operator=() [2/2]

template<class T >
MRISplitIntegrator& MRISplitIntegrator< T >::operator= ( MRISplitIntegrator< T > &&  other)
defaultnoexcept

◆ set_acoustic_substepping()

template<class T >
void MRISplitIntegrator< T >::set_acoustic_substepping ( std::function< void(int, int, int, T &, const T &, T &, T &, const double, const double, const amrex::Real, const double, const double)>  F)
inline
148  {
150  }

Referenced by ERF::advance_dycore().

Here is the caller graph for this function:

◆ set_no_substep()

template<class T >
void MRISplitIntegrator< T >::set_no_substep ( std::function< void(T &, T &, T &, const double, const double, int)>  F)
inline
163  {
164  no_substep = F;
165  }

Referenced by ERF::advance_dycore().

Here is the caller graph for this function:

◆ set_slow_fast_timestep_ratio()

template<class T >
void MRISplitIntegrator< T >::set_slow_fast_timestep_ratio ( const int  timestep_ratio = 1)
inline
153  {
154  slow_fast_timestep_ratio = timestep_ratio;
155  }

Referenced by ERF::advance_dycore().

Here is the caller graph for this function:

◆ set_slow_rhs_post()

template<class T >
void MRISplitIntegrator< T >::set_slow_rhs_post ( std::function< void(T &, T &, T &, T &, const double, const double, const double, const int)>  F)
inline
140  {
141  slow_rhs_post = F;
142  }

Referenced by ERF::advance_dycore().

Here is the caller graph for this function:

◆ set_slow_rhs_pre()

template<class T >
void MRISplitIntegrator< T >::set_slow_rhs_pre ( std::function< void(T &, T &, T &, const double, const double, const double, const int)>  F)
inline
136  {
137  slow_rhs_pre = F;
138  }

Referenced by ERF::advance_dycore().

Here is the caller graph for this function:

◆ setAnelastic()

template<class T >
void MRISplitIntegrator< T >::setAnelastic ( int  _anelastic)
inline
121  {
122  anelastic = _anelastic;
123  }

◆ setForceFirstStageSingleSubstep()

template<class T >
void MRISplitIntegrator< T >::setForceFirstStageSingleSubstep ( int  _force_stage1_single_substep)
inline
131  {
132  force_stage1_single_substep = _force_stage1_single_substep;
133  }

◆ setNcompCons()

template<class T >
void MRISplitIntegrator< T >::setNcompCons ( int  _ncomp_cons)
inline
116  {
117  ncomp_cons = _ncomp_cons;
118  }

◆ setNoSubstepping()

template<class T >
void MRISplitIntegrator< T >::setNoSubstepping ( int  _no_substepping)
inline
126  {
127  no_substepping = _no_substepping;
128  }

Member Data Documentation

◆ acoustic_substepping

template<class T >
std::function<void(int, int, int, T&, const T&, T&, T&, const double, const double, const amrex::Real, const double, const double)> MRISplitIntegrator< T >::acoustic_substepping
private

◆ anelastic

template<class T >
int MRISplitIntegrator< T >::anelastic
private

Should we use the anelastic integrator.

Referenced by MRISplitIntegrator< T >::advance(), and MRISplitIntegrator< T >::setAnelastic().

◆ F_slow

template<class T >
T* MRISplitIntegrator< T >::F_slow
private

◆ force_stage1_single_substep

template<class T >
int MRISplitIntegrator< T >::force_stage1_single_substep
private

Do we follow the recommendation to only perform a single substep in the first RK stage.

Referenced by MRISplitIntegrator< T >::advance(), and MRISplitIntegrator< T >::setForceFirstStageSingleSubstep().

◆ ncomp_cons

template<class T >
int MRISplitIntegrator< T >::ncomp_cons
private

How many components in the cell-centered MultiFab.

Referenced by MRISplitIntegrator< T >::advance(), and MRISplitIntegrator< T >::setNcompCons().

◆ no_substep

template<class T >
std::function<void (T&, T&, T&, const double, const double, int)> MRISplitIntegrator< T >::no_substep
private

The no_substep function is called when we have no acoustic substepping.

Referenced by MRISplitIntegrator< T >::advance(), and MRISplitIntegrator< T >::set_no_substep().

◆ no_substepping

template<class T >
int MRISplitIntegrator< T >::no_substepping
private

Should we not do acoustic substepping.

Referenced by MRISplitIntegrator< T >::advance(), and MRISplitIntegrator< T >::setNoSubstepping().

◆ rhs

template<class T >
std::function<void(T&, const T&, const double, const double)> MRISplitIntegrator< T >::rhs
private

rhs is the right-hand-side function the integrator will use.

Referenced by MRISplitIntegrator< T >::get_rhs().

◆ S_sum

template<class T >
T* MRISplitIntegrator< T >::S_sum
private

◆ slow_fast_timestep_ratio

template<class T >
int MRISplitIntegrator< T >::slow_fast_timestep_ratio = 0
private

The ratio of slow timestep size / fast timestep size (int)

Referenced by MRISplitIntegrator< T >::get_slow_fast_timestep_ratio(), and MRISplitIntegrator< T >::set_slow_fast_timestep_ratio().

◆ slow_rhs_post

template<class T >
std::function<void(T&, T&, T&, T&, const double, const double, const double, const int)> MRISplitIntegrator< T >::slow_rhs_post
private

◆ slow_rhs_pre

template<class T >
std::function<void(T&, T&, T&, const double, const double, const double, const int)> MRISplitIntegrator< T >::slow_rhs_pre
private

◆ T_store

template<class T >
amrex::Vector<std::unique_ptr<T> > MRISplitIntegrator< T >::T_store
private

◆ timestep

template<class T >
double MRISplitIntegrator< T >::timestep
private

Integrator timestep size (Real)

Referenced by MRISplitIntegrator< T >::advance().


The documentation for this class was generated from the following file: