MRISplitIntegrator< T > Class Template Reference

Split integrator for MRI simulations handling slow and fast timescales. More...

#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)
	Initialize integrator storage. More...
	~MRISplitIntegrator ()=default
	MRISplitIntegrator (MRISplitIntegrator &&) noexcept=default
MRISplitIntegrator &	operator= (MRISplitIntegrator &&other) noexcept=default
	MRISplitIntegrator (const MRISplitIntegrator &other)=delete
MRISplitIntegrator &	operator= (const MRISplitIntegrator &other)=delete
void	setNcompCons (int _ncomp_cons)
	Set the number of conservative components. More...
void	setAnelastic (int _anelastic)
	Set whether to use the anelastic integrator. More...
void	setAnelasticType (AnelasticType _anelastic_type=AnelasticType::RK2)
	Set which two-stage scheme the anelastic integrator uses. More...
void	setNoSubstepping (int _no_substepping)
	Set whether acoustic substepping is disabled. More...
void	setForceFirstStageSingleSubstep (int _force_stage1_single_substep)
	Force the first RK stage to perform only a single substep. More...
void	set_slow_rhs_pre (std::function< void(T &, T &, T &, const double, const double, const double, const int)> F)
	Set the pre-substepping slow RHS function. More...
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)
	Set the acoustic substepping function. More...
void	set_slow_fast_timestep_ratio (const int timestep_ratio=1)
	Set the ratio of slow to fast timestep sizes. More...
int	get_slow_fast_timestep_ratio ()
	Get the current slow-to-fast timestep ratio. More...
void	set_no_substep (std::function< void(T &, T &, T &, const double, const double, int)> F)
	Set the function to be called when acoustic substepping is disabled. More...
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)
	Advance the state from time to time + time_step. More...
void	map_data (std::function< void(T &)> Map)
	Apply a mapping function to all internal stored data. More...

Private Member Functions
void	initialize_data (const T &S_data)
	Allocate internal storage for integrator variables. More...

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 = 0.0
	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 = -1
	Should we not do acoustic substepping (-1 until setNoSubstepping is called) More...
int	anelastic = -1
	Should we use the anelastic integrator (-1 until setAnelastic is called) More...
AnelasticType	anelastic_type = AnelasticType::RK2
	Which two-stage scheme the anelastic integrator uses. More...
int	ncomp_cons = -1
	How many components in the cell-centered MultiFab (-1 until setNcompCons is called) More...
int	force_stage1_single_substep = -1
	Do we follow the recommendation to only perform a single substep in the first RK stage (-1 until setForceFirstStageSingleSubstep is called) 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
T *	S_sum = nullptr
T *	F_slow = nullptr

Detailed Description

template<class T>
class MRISplitIntegrator< T >

Split integrator for MRI simulations handling slow and fast timescales.

Template Parameters

T	State type.

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

    {
        initialize_data(S_data);
    }

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

Advance the state from time to time + time_step.

Parameters

[in,out]	S_old	Current state.
[in,out]	S_new	State to be updated.
[in]	time	Current simulation time.
[in]	time_step	Timestep size.

Returns

The actual timestep taken.

    {
    BL_PROFILE_REGION("MRI_advance");
        using namespace amrex;
 
        // *******************************************************************************
        // !no_substepping: we only update the fast variables every fast timestep, then update
        //                  the slow variables after the acoustic sub-stepping.  This has
        //                  2 calls to slow_rhs so that we can update the slow variables
        //                  with the velocity field after the acoustic substepping using
        //                  the time-averaged velocity from the substepping
        // no_substepping:  we don't do any acoustic subcycling so we only make one call per RK
        //                  stage to slow_rhs
        // *******************************************************************************
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(no_substepping >= 0 && anelastic >= 0 &&
                                         force_stage1_single_substep >= 0 && ncomp_cons > 0,
            "MRISplitIntegrator::advance: setNoSubstepping, setAnelastic, "
            "setForceFirstStageSingleSubstep and setNcompCons must be called before advance");
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(S_sum != nullptr && F_slow != nullptr,
            "MRISplitIntegrator::advance: integrator storage must be initialized before advance");
 
        timestep = time_step;
 
        const int substep_ratio = get_slow_fast_timestep_ratio();
 
        if (!no_substepping) {
            AMREX_ALWAYS_ASSERT(substep_ratio > 1 && substep_ratio % 2 == 0);
        }
 
        // Assume before advance() that S_old is valid data at the current time ("time" argument)
        // And that if data is a MultiFab, both S_old and S_new contain ghost cells for evaluating a stencil based RHS
        // We need this from S_old. This is convenient for S_new to have so we can use it
        // as scratch space for stage values without creating a new scratch MultiFab with ghost cells.
 
        // NOTE: In the following, we use S_new to hold S*, S**, and finally, S^(n+1) at the new time
        // DEFINITIONS:
        // S_old  = S^n
        // S_sum  = S(t)
        // F_slow = F(S_stage)
 
        int n_data = IntVars::NumTypes;
 
        /**********************************************/
        /* RK3 Integration with Acoustic Sub-stepping */
        /**********************************************/
        Vector<int> num_vars = {ncomp_cons, 1, 1, 1};
        for (int i(0); i<n_data; ++i)
        {
            // Copy old -> new
            MultiFab::Copy(S_new[i],S_old[i],0,0,num_vars[i],S_old[i].nGrowVect());
        }
 
        // Timestep taken by the fast integrator
        double dtau;
 
        // How many timesteps taken by the fast integrator
        int nsubsteps;
 
        // This is the final time of the full timestep (also the 3rd RK stage)
        // Real new_time = time + timestep;
 
        double time_stage = time;
        double old_time_stage;
 
        if (!anelastic) {
          // RK3 for compressible integrator
          for (int nrk = 0; nrk < 3; nrk++)
          {
            // Capture the time we got to in the previous RK step
            old_time_stage = time_stage;
 
            if (nrk == 0) {
                if (force_stage1_single_substep || no_substepping) {
                    nsubsteps = 1;
                    dtau      = timestep / three;
                } else {
                    // Clamp to 1: substep_ratio is only required to be even, so a ratio of 2
                    // would otherwise give zero substeps and leave S_sum stale in this stage
                    nsubsteps = std::max(substep_ratio/3, 1);
                    dtau      = (timestep / three) / static_cast<double>(nsubsteps);
                }
                time_stage = time + timestep / three;
            }
            if (nrk == 1) {
                if (no_substepping) {
                    nsubsteps = 1;
                    dtau      = myhalf * timestep;
                } else {
                    nsubsteps = std::max(substep_ratio/2, 1);
                    dtau      = (myhalf * timestep) / static_cast<double>(nsubsteps);
                }
                time_stage = time + timestep / two;
            }
            if (nrk == 2) {
                if (no_substepping) {
                    nsubsteps = 1;
                    dtau      = timestep;
                } else {
                    nsubsteps = substep_ratio;
                    dtau      = timestep / static_cast<double>(nsubsteps);
                }
                time_stage = time + timestep;
            }
 
            // step 1 starts with S_stage = S^n  and we always start substepping at the old time
            // step 2 starts with S_stage = S^*  and we always start substepping at the old time
            // step 3 starts with S_stage = S^** and we always start substepping at the old time
 
            slow_rhs_pre(*F_slow, S_old, S_new, time, old_time_stage, time_stage, nrk);
 
            amrex::Real inv_fac = one / static_cast<amrex::Real>(nsubsteps);
 
            // ****************************************************
            // Acoustic substepping
            // ****************************************************
            if (!no_substepping)
            {
                // *******************************************************************************
                // Update the fast variables
                // *******************************************************************************
                for (int ks = 0; ks < nsubsteps; ++ks)
                {
                    acoustic_substepping(ks, nsubsteps, nrk, *F_slow, S_old, S_new, *S_sum, dtau, timestep, inv_fac,
                                         time + ks*dtau, time + (ks+1) * dtau);
 
                } // ks
 
            } else {
                no_substep(*S_sum, S_old, *F_slow, time + nsubsteps*dtau, nsubsteps*dtau, nrk);
            }
 
            // ****************************************************
            // Evaluate F_slow(S_stage) only for the slow variables
            // Note that we are using the current stage versions (in S_new) of the slow variables
            //      (because we didn't update the slow variables in the substepping)
            //      and the momenta slow_rhs_pre set from the same stage, plus the average of
            //      the acoustic substep updates when substepping, to advect them
            //      (the mass fluxes that advance rho over the stage)
            // ****************************************************
            slow_rhs_post(*F_slow, S_old, S_new, *S_sum, time, old_time_stage, time_stage, nrk);
          } // nrk
 
        } else {
            // Two-stage integrator for the anelastic equations.
            //
            // SSP-RK2 (Heun): both stages advance a full timestep from S^n, and the
            //     second stage averages the slow source with the one from the first,
            //     which erf_slow_rhs_{pre,post} do when (anelastic_type == RK2 && nrk == 1).
            //
            // MidPoint: the first stage advances only a half timestep,
            //     S^*   = S^n + (dt/2) F(S^n)
            //     S^n+1 = S^n +  dt    F(S^*)
            //     and the second stage takes no average. The vertical diffusion
            //     is second order in time with a tridiagonal solve in the first
            //     stage only and none in the second.
            const bool l_midpoint = (anelastic_type == AnelasticType::MidPoint);
 
            for (int nrk = 0; nrk < 2; nrk++)
            {
                // Capture the time we got to in the previous RK step
                old_time_stage = time_stage;
 
                // Set the timestep for this stage -- both stages advance from the old time
                nsubsteps  = 1;
                dtau       = (l_midpoint && (nrk == 0)) ? myhalf * timestep : timestep;
                time_stage = time + nsubsteps * dtau;
 
                slow_rhs_pre(*F_slow, S_old, S_new, time, old_time_stage, time_stage, nrk);
 
                no_substep(*S_sum, S_old, *F_slow, time + nsubsteps*dtau, nsubsteps*dtau, nrk);
 
                // ****************************************************
                // Evaluate F_slow(S_stage) only for the slow variables
                // Note that we are using the current stage versions (in S_new) of the slow variables
                //      and the momenta slow_rhs_pre took from the same stage to advect them;
                //      S_sum holds the projected state at the end of the stage
                // ****************************************************
                slow_rhs_post(*F_slow, S_old, S_new, *S_sum, time, old_time_stage, time_stage, nrk);
            } // nrk
        }
 
        // Return timestep
        return timestep;
    }

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

    {
        return rhs;
    }

get_slow_fast_timestep_ratio()

template<class T >

int MRISplitIntegrator< T >::get_slow_fast_timestep_ratio ( )

inline

Get the current slow-to-fast timestep ratio.

Returns

The slow/fast timestep ratio.

    {
        return slow_fast_timestep_ratio;
    }

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

Initialize integrator storage.

Parameters

[in]

S_data

Reference state used to determine storage size.

    {
        initialize_data(S_data);
    }

Here is the call graph for this function:

initialize_data()

template<class T >

void MRISplitIntegrator< T >::initialize_data ( const T &  S_data )

inlineprivate

Allocate internal storage for integrator variables.

Parameters

[in]

S_data

Reference state used to determine storage size.

    {
        // TODO: We can optimize memory by making the cell-centered part of S_sum
        //       have only 2 components, not ncomp_cons components
        const bool include_ghost = true;
        // CreateLike appends to T_store, so drop any storage from a previous
        // initialization before S_sum and F_slow are taken from entries 0 and 1
        T_store.clear();
        amrex::IntegratorOps<T>::CreateLike(T_store, S_data, include_ghost);
        S_sum = T_store[0].get();
        amrex::IntegratorOps<T>::CreateLike(T_store, S_data, include_ghost);
        F_slow = T_store[1].get();
        // initializing to zero
        for (long idx = 0; idx < S_sum->size(); idx++) { (*S_sum)[idx].setVal(0); }
        for (long idx = 0; idx < F_slow->size(); idx++) { (*F_slow)[idx].setVal(0); }
    }

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

Apply a mapping function to all internal stored data.

Parameters

[in]

Map

Mapping function to apply to each stored state object.

    {
        for (auto& F : T_store) {
            Map(*F);
        }
    }

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

Set the acoustic substepping function.

Parameters

[in]

F

Function to perform the acoustic substepping update.

    {
        acoustic_substepping = F;
    }

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

Set the function to be called when acoustic substepping is disabled.

Parameters

[in]

F

Function to use in place of substepping.

    {
        no_substep = F;
    }

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

Set the ratio of slow to fast timestep sizes.

Parameters

[in]

timestep_ratio

Ratio of slow/fast timesteps.

    {
        slow_fast_timestep_ratio = timestep_ratio;
    }

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

    {
        slow_rhs_post = F;
    }

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

Set the pre-substepping slow RHS function.

Parameters

[in]

F

Function to compute the pre-substepping slow RHS.

    {
        slow_rhs_pre = F;
    }

Referenced by ERF::advance_dycore().

Here is the caller graph for this function:

setAnelastic()

template<class T >

void MRISplitIntegrator< T >::setAnelastic ( int  _anelastic )

inline

Set whether to use the anelastic integrator.

Parameters

[in]

_anelastic

Integer flag (1 for anelastic, 0 for compressible).

    {
        anelastic = _anelastic;
    }

setAnelasticType()

template<class T >

void MRISplitIntegrator< T >::setAnelasticType ( AnelasticType  _anelastic_type = AnelasticType::RK2 )

inline

Set which two-stage scheme the anelastic integrator uses.

Parameters

[in]

_anelastic_type

AnelasticType::RK2 or AnelasticType::MidPoint.

    {
        anelastic_type = _anelastic_type;
    }

setForceFirstStageSingleSubstep()

template<class T >

void MRISplitIntegrator< T >::setForceFirstStageSingleSubstep ( int  _force_stage1_single_substep )

inline

Force the first RK stage to perform only a single substep.

Parameters

[in]

_force_stage1_single_substep

Integer flag to enable this behavior.

    {
        force_stage1_single_substep = _force_stage1_single_substep;
    }

setNcompCons()

template<class T >

void MRISplitIntegrator< T >::setNcompCons ( int  _ncomp_cons )

inline

Set the number of conservative components.

Parameters

[in]

_ncomp_cons

Number of conservative components.

    {
        ncomp_cons = _ncomp_cons;
    }

setNoSubstepping()

template<class T >

void MRISplitIntegrator< T >::setNoSubstepping ( int  _no_substepping )

inline

Set whether acoustic substepping is disabled.

Parameters

[in]

_no_substepping

Integer flag (1 to disable substepping).

    {
        no_substepping = _no_substepping;
    }

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

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

anelastic

template<class T >

int MRISplitIntegrator< T >::anelastic = -1

private

Should we use the anelastic integrator (-1 until setAnelastic is called)

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

anelastic_type

template<class T >

AnelasticType MRISplitIntegrator< T >::anelastic_type = AnelasticType::RK2

private

Which two-stage scheme the anelastic integrator uses.

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

F_slow

template<class T >

T* MRISplitIntegrator< T >::F_slow = nullptr

private

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

force_stage1_single_substep

template<class T >

int MRISplitIntegrator< T >::force_stage1_single_substep = -1

private

Do we follow the recommendation to only perform a single substep in the first RK stage (-1 until setForceFirstStageSingleSubstep is called)

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

ncomp_cons

template<class T >

int MRISplitIntegrator< T >::ncomp_cons = -1

private

How many components in the cell-centered MultiFab (-1 until setNcompCons is called)

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 = -1

private

Should we not do acoustic substepping (-1 until setNoSubstepping is called)

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 = nullptr

private

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

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

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

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

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

T_store

template<class T >

amrex::Vector<std::unique_ptr<T> > MRISplitIntegrator< T >::T_store

private

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

timestep

template<class T >

double MRISplitIntegrator< T >::timestep = 0.0

private

Integrator timestep size (Real)

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

---

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

• Source/TimeIntegration/ERF_MRI.H