Mesh ODE solver

The MeshODESolver<N> implements a solver for

\[\frac{\partial \vec{\phi}}{\partial t} = \vec{S},\]

where \(\vec{\phi}\) represents \(N\) unknowns on the mesh, and \(\vec{S}\) is the corresponding source term. The class is templated as

/*!
  @brief Class for solving dy/dt = f on an AMR hierarchy. 
  @details The template parameter is the number of variables in y and f.
*/
template <size_t N = 1>
class MeshODESolver

where N indicates the number of variables stored on the mesh.

MeshODESolver<N> is designed to store N variables in each grid cell, without any cell-to-cell coupling. To instantiate the solver, use the full constructor with reference to AmrMesh:

/*!
  @brief Full constructor. 
  @param[in] a_amr AMR core.
*/
MeshODESolver(const RefCountedPtr<AmrMesh>& a_amr) noexcept;

Tip

Source code for the MeshODESolver<N> resides in Source/MeshODESolver. See https://chombo-discharge.github.io/chombo-discharge/doxygen/html/classMeshODESolver.html for the full C++ API.

Setting \(\vec{\phi}\)

Mesh-based

To set \(\vec{\phi}\), one can fetch the \(N\) mesh components from

/*!
  @brief Get the solution vector (left-hand side of equation).
*/
EBAMRCellData&
getPhi() noexcept;

This will return the data holder that holds the cell centered data \(\vec{\phi}\) which the user can iterate through. See Mesh data for examples.

Analytic function

One can set \(\vec{\phi}(\mathbf{x}) = \vec{f}\left(\mathbf{x}\right)\) through the following functions:

/*!
  @brief Set phi for a specific component. 
  @param[in] a_phiFunc Phi function
  @param[in] a_comp    Component
*/
virtual void
setPhi(const std::function<Real(const RealVect& a_pos)>& a_phiFunc, const size_t a_comp) noexcept;

/*!
  @brief Set phi everywhere. 
  @param[in] a_phiFunc Phi function. 
*/
virtual void
setPhi(const std::function<std::array<Real, N>(const RealVect& a_pos)>& a_phiFunc) noexcept;

These differ only in that the first function sets a specific component, whereas the second version sets all \(N\) components.

Setting \(\vec{S}\)

Mesh-based

For a general method of setting the source term one can fetch \(\vec{S}\) through

/*!
  @brief Get the solution vector (left-hand side of equation).
*/
EBAMRCellData&
getRHS() noexcept;

This returns a reference to \(\vec{S}\) which the user can iterate through and set the value in each cell. See Mesh data for explicit examples.

Spatially dependent

The source term can also be set on a component-by-component basis using

/*!
  @brief Set right-hand side for specified component. 
  @param[in] a_srcFunc Source term function.
  @param[in] a_comp    Component
*/
virtual void
setRHS(const std::function<Real(const RealVect& a_pos)>& a_srcFunc, const size_t a_comp) noexcept;

As a function of \(\vec{\phi}\)

In order to compute the source term \(\vec{S}\) as a function of \(\vec{\phi}\), MeshOdeSolver<N> has a function

/*!
  @brief Alias for right-hand side
*/
using RHSFunction = std::function<std::array<Real, N>(const std::array<Real, N>&, const Real)>;

/*!
  @brief Compute right-hand side from left-hand side. I.e. compute f = f(y,t). 
  @param[in] a_rhsFunction Function for computing the right-hand side. 
*/
virtual void
computeRHS(const RHSFunction& a_rhsFunction) noexcept;

which computes the source term \(\vec{S}\) as a function

\[\vec{S} = \vec{f}\left(\vec{\phi},t\right).\]

An example which sets \(\vec{S} = \vec{\phi}\) is given below

auto f = [](const std::array<Real, N>& phi, const Real t) -> std::array<Real, N> {
   const std::array<Real, N> S = phi;

   return S;
};

solver.computeRHS(f);

Regridding

When regridding the MeshODESolver<N>, one must first ensure that the mesh data on the old mesh is stored before calling the regrid function:

/*!
  @brief Perform pre-regrid operations. 
  @param[in] a_lbase          Coarsest level that changed during regrid. 
  @param[in] a_oldFinestLevel Finest grid level before the regrid operation. 
  @note This copies m_phi onto m_cache
*/
virtual void
preRegrid(const int a_lbase, const int a_oldFinestLevel) noexcept;

This must be done before AmrMesh creates the new grids. This will store \(\vec{\phi}\) on the old mesh. After AmrMesh has generated the new grids, \(\vec{\phi}\) can be interpolated onto the new grids by calling

/*!
  @brief Regrid this solver. 
  @param[in] a_lmin           Coarsest level where grids did not change. 
  @param[in] a_oldFinestLevel Finest AMR level before the regrid. 
  @param[in] a_newFinestLevel Finest AMR level after the regrid. 
  @details This linearly interpolates (potentially with limiters) m_phi to the new grids. 
*/
virtual void
regrid(const int a_lmin, const int a_oldFinestLevel, const int a_newFinestLevel) noexcept;

Users can also choose to turn on/off slope limiters when putting the solution on the new mesh.

Important

The source term \(\vec{S}\) is also allocated on the new mesh, but is not interpolated onto the new grids. It must therefore be set by the user after calling the regrid function.

Input options

Several input options are available for configuring the run-time configuration of MeshODESolver<N>, which are listed below

Listing 28 Input options for the MeshODESolver<N> class. All options are run-time configurable.

# ====================================================================================================
# MeshODESolver class options
# ====================================================================================================
MeshODESolver.verbosity             = -1         ## Verbosity level
MeshODESolver.plt_vars              = phi rhs    ## Stuff that can be plotted. 'phi' and/or 'rhs
MeshODESolver.use_regrid_slopes     = false      ## Use or don't use slopes when regridding

I/O

The user can add \(\vec{\phi}\) and \(\vec{S}\) to output files by specifying these in the input script. These variables are named

MeshODESolver.plt_vars = phi rhs

Only phi and rhs are recognized as valid arguments. If choosing to omit output variables for the solver, one can put e.g. MeshODESolver.plt_vars = -1.

Note

MeshODESolver<N> checkpoint files only contain \(\vec{\phi}\).