TimeStepper
TimeStepper is essentially the problem solving class chombo-discharge.
The class is abstract and is subclassed by the programmer in order to solve new problems (or change time discretizations).
Typically, TimeStepper owns all numerical solvers, has responsibility for setting up solvers, regridding its internal state, and advancing the equations of motion.
Because TimeStepper and not Driver owns the solvers as class members, it will also partially coordinate I/O by providing data that Driver will add to HDF5 files.
Tip
Here is the TimeStepper C++ API
Basic functions
There are numerous functions that must be implemented and coordinated in order to provide a full-fledged TimeStepper for various problems.
Below, we include the current header file for TimeStepper.
/* chombo-discharge
* Copyright © 2021 SINTEF Energy Research.
* Please refer to Copyright.txt and LICENSE in the chombo-discharge root directory.
*/
/*!
@file CD_TimeStepper.H
@brief Declaration of main (abstract) time stepper class.
@author Robert Marskar
*/
#ifndef CD_TimeStepper_H
#define CD_TimeStepper_H
// Chombo includes
#include <CH_HDF5.H>
// Our includes
#include <CD_ComputationalGeometry.H>
#include <CD_MultiFluidIndexSpace.H>
#include <CD_AmrMesh.H>
#include <CD_NamespaceHeader.H>
/*!
@brief Base class for advancing equations.
@details This class is used by Driver for advancing sets of equations. In short, this class should implement a time-stepping routine for a
set of solvers. One should also implement routines for setting up the solvers, allocating necessary memory, regridding etc.
*/
class TimeStepper
{
public:
/*!
@brief Default constructor (does nothing)
*/
TimeStepper();
/*!
@brief Default destructor (does nothing)
*/
virtual ~TimeStepper();
/*!
@brief Set AmrMesh
@param[in] a_amr AmrMesh
*/
void
setAmr(const RefCountedPtr<AmrMesh>& a_amr);
/*!
@brief Set the computational geometry
@param[in] a_computationalGeometry The computational geometry.
*/
void
setComputationalGeometry(const RefCountedPtr<ComputationalGeometry>& a_computationalGeometry);
/*!
@brief Set up solvers
*/
virtual void
setupSolvers() = 0;
/*!
@brief Allocate data for the time stepper and solvers.
*/
virtual void
allocate() = 0;
/*!
@brief Fill solvers with initial data
*/
virtual void
initialData() = 0;
/*!
@brief Post-initialize operations to be performed at end of setup stage.
*/
virtual void
postInitialize() = 0;
/*!
@brief Post-initialize operations to be performed after filling solvers with data read from checkpoint files.
*/
virtual void
postCheckpointSetup() = 0;
/*!
@brief Register realms to be used for the simulation.
*/
virtual void
registerRealms() = 0;
/*!
@brief Register operators to be used for the simulation
*/
virtual void
registerOperators() = 0;
/*!
@brief Parse runtime options
@details Override this routine if your time stepper can use run-time configuration of the solvers that it advances. This can e.g.
be the CFL condition.
*/
virtual void
parseRuntimeOptions();
#ifdef CH_USE_HDF5
/*!
@brief Read header data from checkpoint file.
@param[inout] a_header HDF5 header.
*/
virtual void
readCheckpointHeader(HDF5HeaderData& a_header);
#endif
#ifdef CH_USE_HDF5
/*!
@brief Write header data to checkpoint file.
@param[inout] a_header HDF5 header.
*/
virtual void
writeCheckpointHeader(HDF5HeaderData& a_header) const;
#endif
#ifdef CH_USE_HDF5
/*!
@brief Write checkpoint data to file
@param[inout] a_handle HDF5 file
@param[in] a_lvl Grid level
@details Implement this routine for checkpointing data for restarts. This will
typically call the solvers' checkpointing routine, but more data can be added.
*/
virtual void
writeCheckpointData(HDF5Handle& a_handle, const int a_lvl) const = 0;
#endif
#ifdef CH_USE_HDF5
/*!
@brief Read checkpoint data from file
@param[inout] a_handle HDF5 file
@param[in] a_lvl Grid level
@details Implement this routine for reading data for restarts. This will typically call the solvers' checkponiting routine.
*/
virtual void
readCheckpointData(HDF5Handle& a_handle, const int a_lvl) = 0;
#endif
/*!
@brief Get the number of plot variables for this time stepper.
@details This is necessary because Driver, not TimeStepper, is responsible for allocating the necessary memory.
@return Returns number of plot variables that will be written during writePlotData
*/
virtual int
getNumberOfPlotVariables() const = 0;
/*!
@brief Get plot variable names
*/
virtual Vector<std::string>
getPlotVariableNames() const = 0;
/*!
@brief Write plot data to output holder.
@param[inout] a_output Output data holder.
@param[inout] a_icomp Starting component in a_output to begin at.
@param[in] a_outputRealm Realm where a_output belongs
@param[in] a_level Grid level
*/
virtual void
writePlotData(LevelData<EBCellFAB>& a_output,
int& a_icomp,
const std::string a_outputRealm,
const int a_level) const = 0;
/*!
@brief An option for calling special functions prior to plotting data.
Called by Driver in the IMMEDIATELY before writing the plot file.
*/
virtual void
prePlot();
/*!
@brief An option for calling special functions prior to plotting data.
Called by Driver in the IMMEDIATELY after writing the plot file.
*/
virtual void
postPlot();
/*!
@brief Get computational loads to be checkpointed.
@details This is used by Driver both for setting up load-balanced restarts AND for plotting the computational loads to a file. This routine is
disjoint from loadBalanceBoxes because this routine is not part of a regrid. This means that we are not operating with temporarily load balanced
grids where the but the final ones.
@note The default implementation uses the box volume as a proxy for the load. You should overwrite this if you load balance your application, and also
make sure that the loads returned from this routine are consistent with what you put in loadBalanceBoxes.
Also note that the return vector has the same ordering as the DisjointBoxLayout's boxes on the input grid level. See the implementation for further details.
@param[in] a_realm Realm
@param[in] a_level Grid level
@return Returns computational loads for each box on grid level a_level.
*/
virtual Vector<long int>
getCheckpointLoads(const std::string a_realm, const int a_level) const;
/*!
@brief Compute a time step to be used by Driver.
*/
virtual Real
computeDt() = 0;
/*!
@brief Advancement method. The implementation of this method should advance all equations of motion
@param[in] a_dt Time step to be used for advancement
@return Returns the time step that was used.
@note The return value does not need to equal a_dt. Adaptive time stepping methods will generally return != a_dt.
*/
virtual Real
advance(const Real a_dt) = 0;
/*!
@brief Synchronzie solver times and time steps
@param[in] a_step Time step
@param[in] a_time Time (in seconds)
@param[in] a_dt Time step that was used.
*/
virtual void
synchronizeSolverTimes(const int a_step, const Real a_time, const Real a_dt) = 0;
/*!
@brief Print a step report.
@details This is called by Driver after time step. The routine can be used to display use information about the simulation progress.
*/
virtual void
printStepReport() = 0;
/*!
@brief Perform pre-regrid operations.
@details This should include all copying all data which should be interpolated to the new grids. It can also include deallocating memory in case the regrid
operation takes a lot of memory.
@param[in] a_lmin The coarsest level that changes
@param[in] a_oldFinestLevel The finest level before the regrid.
*/
virtual void
preRegrid(const int a_lmin, const int a_oldFinestLevel) = 0;
/*!
@brief Time stepper regrid method.
@param[in] a_lmin The coarsest level that changed.
@param[in] a_oldFinestLevel The finest level before the regrid.
@param[in] a_newFinestLevel The finest level after the regrid.
*/
virtual void
regrid(const int a_lmin, const int a_oldFinestLevel, const int a_newFinestLevel) = 0;
/*!
@brief Perform post-regrid operations.
@details This includes all operations to be done AFTER interpolating data to new grids.
*/
virtual void
postRegrid() = 0;
/*!
@brief Fuction which can have Driver do regrids at arbitrary points in the simulation.
@details This is called by Driver at every time step and if it returns true, Driver will regrid, regardless of whether or not it is
on the regular "regrid intervals".
*/
virtual bool
needToRegrid();
/*!
@brief Load balancing query for a specified realm. If this returns true for a_realm, load balancing routines will be called during regrids.
@param[in] a_realm Realm name
*/
virtual bool
loadBalanceThisRealm(const std::string a_realm) const;
/*!
@brief Load balance grid boxes for a specified realm.
@param[out] a_procs MPI ranks owning the various grid boxes.
@param[out] a_boxes Grid boxes on every level (obtain them with a_grids[lvl].boxArray())
@param[in] a_realm Realm identifier
@param[in] a_grids Original grids
@param[in] a_lmin Coarsest grid level that changed
@param[in] a_finestLevel New finest grid level
@details This is only called by Driver if TimeStepper::loadBalanceThisRealm(a_realm) returned true. The default implementation
uses volume-based loads for the grid patches. If the user wants to load balance boxes on a realm, this routine must be overwritten and
he should compute loads for the various patches in a_grids and call LoadBalancing::makeBalance on each level. It is up to the user/programmer
to decide if load balancing should be done independently on each level, or if loads per MPI rank are accumulated across levels.
*/
virtual void
loadBalanceBoxes(Vector<Vector<int>>& a_procs,
Vector<Vector<Box>>& a_boxes,
const std::string a_realm,
const Vector<DisjointBoxLayout>& a_grids,
const int a_lmin,
const int a_finestLevel);
protected:
/*!
@brief Class verbosity
*/
int m_verbosity;
/*!
@brief Time step
*/
int m_timeStep;
/*!
@brief TIme
*/
Real m_time;
/*!
@brief Previous time step size
*/
Real m_dt;
/*!
@brief AmrMesh.
*/
RefCountedPtr<AmrMesh> m_amr;
/*!
@brief Computational geometry.
*/
RefCountedPtr<ComputationalGeometry> m_computationalGeometry;
};
#include <CD_NamespaceFooter.H>
#endif
Setup routines
Here, we consider the various setup routines in TimeStepper.
The routines are used by Driver in the simulation setup step, both for fresh simulation setups as well as restarts.
registerRealms
chombo-discharge permits things to happen on different sets of grids where the the grids themselves cover the same physical region, but where MPI ownership of grids might change between grid sets (see Realm for details).
To register a Realm, users will have TimeStepper register realms in the registerRealms() routine, as follows:
void myTimeStepper::registerRealms(){
m_amr->registerRealm(Realm::Primal);
m_amr->registerRealm("particleRealm");
m_amr->registerRealm("otherParticleRealm");
}
The above code will ensure that chombo-discharge generates three Realm, which can be individually load balanced.
Since at least one realm is required, Driver will always register the realm "Primal".
Fundamentally, there is no limitation to the number of realms that can be allocated.
setupSolvers
setupSolvers is used for instantiating the solvers.
This routine is called prior to creating grids, so it is not possible to allocate mesh data for the solvers inside this routine.
The rationale for this design is that AmrMesh must know relatively early which part of the AMR infrastructure that will be instantiated, so the solvers are created before allocating the grids.
It is still quite possible to parse lots of data into the solvers, e.g., setting input variables.
To provide an example, the code snippet below shows the implementation of this routine for the TimeStepper implementation for the Advection-diffusion model:
void
AdvectionDiffusionStepper::setupSolvers()
{
CH_TIME("AdvectionDiffusionStepper::setupSolvers");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::setupSolvers" << endl;
}
CH_assert(!m_solver.isNull());
// Instantiate the species.
m_species = RefCountedPtr<AdvectionDiffusionSpecies>(
new AdvectionDiffusionSpecies(m_initialData, m_mobile, m_diffusive));
// Prep the solver.
m_solver->setVerbosity(m_verbosity);
m_solver->setSpecies(m_species);
m_solver->parseOptions();
m_solver->setPhase(m_phase);
m_solver->setAmr(m_amr);
m_solver->setComputationalGeometry(m_computationalGeometry);
m_solver->setRealm(m_realm);
if (!m_solver->isMobile() && !m_solver->isDiffusive()) {
MayDay::Error("AdvectionDiffusionStepper::setupSolvers - can't turn off both advection AND diffusion");
}
}
registerOperators
Internally, an instantiation of Realm will provide access to the grids and cut-cell information on that Realm, as well as any operators that the user has seen fit to register.
Various operators are available for, e.g., computing gradients, conservative coarsening, ghost cell interpolation, filling a patch with interpolation data, redistribution, particle-mesh operations, and so on.
Since operators always incur overhead and not all applications require all operators, they must be registered.
If a solver needs an operator for, say, piecewise linear ghost cell interpolation, the solver needs to register that operator through the AmrMesh public interface:
m_amr->registerOperator(s_eb_pwl_interp, m_realm, m_phase);
Once an operator has been registered, Realm will automatically instantiate those operators during initialization or regrid operations.
Run-time error messages are issued if an AMR operator is used, but has not been registered.
More commonly, chombo-discharge solvers will contain a routine that registers the operators that the solver needs.
A valid TimeStepper implementation must register all required operators in the function registerOperators(), but this is usually done by letting the solvers register what they need.
Failure to do so will issue a run-time error.
Solvers will typically allocate a subset of these operators, but for multiphysics code that use both fluid and particles, most of these will be in use.
An example of this is given in Listing 7, which shows how this routine is implemented for the Advection-diffusion model:
registerOperators routine for a simple advection-diffusion problem.void
AdvectionDiffusionStepper::registerOperators()
{
CH_TIME("AdvectionDiffusionStepper::registerOperators");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::registerOperators" << endl;
}
// Let the solver do this -- it knows what it needs.
m_solver->registerOperators();
}
allocate
allocate is used for allocating particle and mesh data that is required during simulations.
This step is done after the grids have been initialized by AmrMesh and during regrids.
Again using Advection-diffusion model as an example,
void
AdvectionDiffusionStepper::allocate()
{
CH_TIME("AdvectionDiffusionStepper::allocate");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::allocate" << endl;
}
m_solver->allocate();
}
In the above snippet, TimeStepper only calls the solver allocation function (solvers generally know how to allocate their own internal data).
For more complex problems this routine will probably allocate additonal data that only lives within TimeStepper (and not the solvers).
initialData
initialData is called by Driver setup routines after the allocate step, and has responsibility of setting up the problem with initial data.
This can occasionally be simple, or for coupled problems it might be highly complex.
For discharge problems, this can involve filling the solvers with initial densities, and solving the Poisson equation for obtaining the electric field.
A simpler example is again given by Advection-diffusion model:
void
AdvectionDiffusionStepper::initialData()
{
CH_TIME("AdvectionDiffusionStepper::initialData");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::initialData" << endl;
}
// Fill the solver with initial data from the species.
m_solver->initialData();
// Set velocity, diffusion coefficient, and boundary conditions.
m_solver->setSource(0.0);
m_solver->setEbFlux(0.0);
if (m_solver->isDiffusive()) {
m_solver->setDiffusionCoefficient(m_diffCo);
}
if (m_solver->isMobile()) {
m_solver->setVelocity(m_velocity);
}
// Set flux functions
auto fluxFunc = [](const RealVect a_pos, const Real a_time) {
return 0.0;
};
// m_solver->setDomainFlux(fluxFunc);
}
The above code defers the initialization of the density in the advection-diffusion-reaction solver to the actual solver implementation. Here, one could equally well have fetched the density directly from the solver and set it to something else by simply iterating through the grid cells (or setting it through other means). In addition, source terms are set to zero, as are boundary fluxes.
postInitialize
postInitialize is a special routine that is called after Driver has filled the solvers with initial data and Driver is done with all the initial regrids.
While most data initialization steps can, however, be done in initialData, the function is put there as an open door to the programmer for performing certain post-initialization functions that do not not need be performed in initialData (which is called once per initial regrid).
For example, the Discharge inception model uses this function to compute several relevant quantities after the electric has been obtained in initialData.
postCheckpointSetup
During simulation restarts, Driver will open an HDF5 file and have TimeStepper fill solvers and its own internal data with data from that file.
postCheckpointSetup is a routine which is called immediately after the solvers have performed this step.
Several gas discharge models use this function to compute the electric field from the potential that was saved in the HDF5 file.
I/O routines
TimeStepper contains I/O routines primarily serves two purposes:
To provide data for HDF5 plot files, used for post-processing analysis.
To read and write data from HDF5 checkpoint files, which are used to restart simulations from a specified time step.
In general, plot and checkpoint data do not contain the same data, and TimeStepper therefore requires that plot and checkpoint files are filled separately.
We discuss these below:
getNumberOfPlotVariables
getNumberOfPlotVariables must return the number of components that will be plotted by TimeStepper.
The reason why this routine exists is the Driver will pre-allocate the necessary memory on each AMR level, and TimeStepper will then copy solver data into this data holder.
Specifically, if TimeStepper will plot a single scalar, it must return a value of one.
If it plots a single vector, it must return a value of SpaceDim.
Below we include the implementation of this routine for the Advection-diffusion model:
getNumberOfPlotVariables routine for a simple advection-diffusion problem.int
AdvectionDiffusionStepper::getNumberOfPlotVariables() const
{
CH_TIME("AdvectionDiffusionStepper::getNumberOfPlotVariables");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::getNumberOfPlotVariables" << endl;
}
// Not plotting anything of our own, so return whatever the solver wants to plot.
return m_solver->getNumberOfPlotVariables();
}
getPlotVariableNames
getPlotVariableNames provides a list of plot variable names for the HDF5 file.
This list must have the same length as the returned value of getNumberOfPlotVariables.
Below we include the implementation of this routine for the Advection-diffusion model:
getPlotVariableNames routine for a simple advection-diffusion problem.Vector<std::string>
AdvectionDiffusionStepper::getPlotVariableNames() const
{
CH_TIME("AdvectionDiffusionStepper::getPlotVariableNames");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::getPlotVariableNames" << endl;
}
return m_solver->getPlotVariableNames();
}
Tip
When writing some vector data \(F\) to the HDF5-file, one can write the variables as x-F, y-F, and z-F, and VisIt visualization will automatically recognize \(F\) as a vector field.
writePlotData
writePlotData will write the plot data to the provided data holder.
The function signature is
/*!
@brief Write plot data to output holder.
@param[inout] a_output Output data holder.
@param[inout] a_icomp Starting component in a_output to begin at.
@param[in] a_outputRealm Realm where a_output belongs
@param[in] a_level Grid level
*/
virtual void
writePlotData(LevelData<EBCellFAB>& a_output,
int& a_icomp,
const std::string a_outputRealm,
In this function, a_output is pre-allocated block of memory that TimeStepper will write its components to (beginning at a_icomp).
Note that if TimeStepper writes \(N\) components, the implementation must increment a_icomp by \(N\).
Usually, solvers will have their own writePlotData routines which lets TimeStepper simply call the solver functions.
An example is given below for the Advection-diffusion model:
void
AdvectionDiffusionStepper::writePlotData(LevelData<EBCellFAB>& a_output,
int& a_icomp,
const std::string a_outputRealm,
const int a_level) const
{
CH_TIME("AdvectionDiffusionStepper::writePlotData");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::writePlotData" << endl;
}
CH_assert(a_level >= 0);
CH_assert(a_level <= m_amr->getFinestLevel());
m_solver->writePlotData(a_output, a_icomp, a_outputRealm, a_level);
}
writeCheckpointData
writeCheckpointData must write necessary data for checkpointing the simulation state.
This data is used when restarting simulations from a checkpoint file.
Note that checkpoint data is written on a level-by-level basis.
The function signature is
/*!
@brief Write checkpoint data to file
@param[inout] a_handle HDF5 file
@param[in] a_lvl Grid level
@details Implement this routine for checkpointing data for restarts. This will
typically call the solvers' checkpointing routine, but more data can be added.
*/
virtual void
writeCheckpointData(HDF5Handle& a_handle, const int a_lvl) const = 0;
Usually, the solvers know themselves what data to put in the checkpoint files and these routines are then pretty simple. Below, we again include an example for the Advection-diffusion model:
writeCheckpointData routine for a simple advection-diffusion problem.void
AdvectionDiffusionStepper::writeCheckpointData(HDF5Handle& a_handle, const int a_lvl) const
{
CH_TIME("AdvectionDiffusionStepper::writeCheckpointData");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::writeCheckpointData" << endl;
}
m_solver->writeCheckpointLevel(a_handle, a_lvl);
}
When implementing the function, it is important to add any data that is required when restarting the simulation. It is also beneficial to not add data that is not required since this will just lead to larger files. An example of this is the FieldSolver class, which only checkpoints the potential and not the electric field, as the latter is simply obtained by taking the gradient.
Tip
Computational particles can also be added to the checkpoint files.
readCheckpointData
readCheckpointData is the function that will read data from an HDF5 checkpoint file and populate TimeStepper with this data.
The data is read on a level-by-level basis, with a function signature
/*!
@brief Read checkpoint data from file
@param[inout] a_handle HDF5 file
@param[in] a_lvl Grid level
@details Implement this routine for reading data for restarts. This will typically call the solvers' checkponiting routine.
*/
virtual void
readCheckpointData(HDF5Handle& a_handle, const int a_lvl) = 0;
Solvers will normally already know what data to read into their data members. E.g., the example for the Advection-diffusion model is
readCheckpointData routine for a simple advection-diffusion problem.void
AdvectionDiffusionStepper::readCheckpointData(HDF5Handle& a_handle, const int a_lvl)
{
CH_TIME("AdvectionDiffusionStepper::readCheckpointData");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::readCheckpointData" << endl;
}
m_solver->readCheckpointLevel(a_handle, a_lvl);
}
Advance routines
computeDt
computeDt is a routine that will compute a trial time step when calling the advance method.
We have chosen to call this a trial time step because
Driver might choose to use a smaller time step in order to write plot files at specific times.
When calling the actual advance method (see below), it is possible to return a different time step than the one computed through
computeDt.
The calculation of a time step can be quite involved, depending on the application being imlemented.
Moreover, many TimeStepper implementations will provide hooks for swapping algorithms, and in this case the time step might be limited differently.
For the Advection-diffusion model the implementation is as follows:
Real
AdvectionDiffusionStepper::computeDt()
{
CH_TIME("AdvectionDiffusionStepper::computeDt");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::computeDt" << endl;
}
// TLDR: If we run explicit advection but implicit diffusion then we are only limited by the advective CFL. Otherwise,
// if diffusion is also explicit we need the advection-diffusion limited time step.
// A weird thing, but sometimes we want to be able to force the CFL so that
// we override run-time configurations of the CFL number. This code does that.
Real cfl = 0.0;
if (m_forceCFL > 0.0) {
cfl = m_forceCFL;
}
else {
cfl = m_cfl;
}
Real dt = std::numeric_limits<Real>::max();
switch (m_integrator) {
case Integrator::Heun: {
dt = cfl * m_solver->computeAdvectionDiffusionDt();
break;
}
case Integrator::IMEX: {
dt = cfl * m_solver->computeAdvectionDt();
break;
}
default: {
MayDay::Error("AdvectionDiffusionStepper::computeDt - logic bust");
break;
}
}
dt = std::max(dt, m_minDt);
dt = std::min(dt, m_maxDt);
return dt;
}
In the code above, TimeStepper supports both fully explicit advection-diffusion advances as well as split-step advances with implicit diffusion.
Depending on how the user chooses to run the code, the time step is therefore computed differently.
At the bottom of the above code, hard limits on the time step are also enforced.
advance
The advance method has responsibility for advancing physics module one time step, and is called by Driver.
The function signature is
/*!
@brief Advancement method. The implementation of this method should advance all equations of motion
@param[in] a_dt Time step to be used for advancement
@return Returns the time step that was used.
@note The return value does not need to equal a_dt. Adaptive time stepping methods will generally return != a_dt.
*/
virtual Real
advance(const Real a_dt) = 0;
As mentioned in the documentation for this method, the function takes a trial time step a_dt which is the physical time step.
This time step is the one computed by computeDt.
It is, however, quite possible to advance the equations of motion over a time that does not equal a_dt, which is generally the case for adaptive time stepping methods.
The implementation of the advance method is usually the most time-consuming part of implementation a new TimeStepper, and the implementation of this routine can become substantially complicated.
For simpler problems this routine is relatively straightforward to implement, however, also in a way that involves AMR and cut-cells.
Real
AdvectionDiffusionStepper::advance(const Real a_dt)
{
CH_TIME("AdvectionDiffusionStepper::advance");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::advance" << endl;
}
// State to be advanced.
EBAMRCellData& state = m_solver->getPhi();
switch (m_integrator) {
case Integrator::Heun: {
const bool conservativeOnly = false;
const bool addEbFlux = true;
const bool addDomainFlux = true;
// Transient storage
EBAMRCellData yp;
EBAMRCellData k1;
EBAMRCellData k2;
m_amr->allocate(yp, m_realm, m_phase, 1);
m_amr->allocate(k1, m_realm, m_phase, 1);
m_amr->allocate(k2, m_realm, m_phase, 1);
// Compute k1 coefficient
m_solver->computeDivJ(k1, state, 0.0, conservativeOnly, addEbFlux, addDomainFlux);
DataOps::copy(yp, state);
DataOps::incr(yp, k1, -a_dt);
// Compute k2 coefficient and final state
m_solver->computeDivJ(k2, yp, 0.0, conservativeOnly, addEbFlux, addDomainFlux);
DataOps::incr(state, k1, -0.5 * a_dt);
DataOps::incr(state, k2, -0.5 * a_dt);
break;
}
m_amr->conservativeAverage(state, m_realm, m_phase);
m_amr->interpGhost(state, m_realm, m_phase);
return a_dt;
}
In the above code we have included the part of the advance routine that executes Heun’s method for advancing the scalar.
This code also includes allocation of temporaries for computing the coefficients, storage for the intermediate state, and enforcement of boundary conditions, all of which include AMR.
At the way out of the routine the solution is coarsened and the ghost cells are updated, and the trial time step (a_dt) is returned.
synchronizeSolverTimes
synchronizeSolverTimes is called after the advance method and is used to update the simulation time for all solvers.
Again, this routine exists because there is often a physical time to be tracked by the solvers (e.g., enforcement of time-dependent voltage applications).
This routine simply ensures that TimeStepper and all solvers that TimeStepper owns see the same physical time, number of steps, and time step sizes.
The implementation for Advection-diffusion model is
synchronizeSolverTimes routine for a simple advection-diffusion problem.void
AdvectionDiffusionStepper::synchronizeSolverTimes(const int a_step, const Real a_time, const Real a_dt)
{
CH_TIME("AdvectionDiffusionStepper::synchronizeSolverTimes");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::synchronizeSolverTimes" << endl;
}
m_timeStep = a_step;
m_time = a_time;
m_dt = a_dt;
m_solver->setTime(a_step, a_time, a_dt);
}
printStepReport
printStepReport is called after the advance method, and provides extra information printed to the pout.* files (see Controlling chombo-discharge).
This function is called by Driver after performing a time step, and can be used to print extra information not covered by Driver, such as how the time step was limited, or other information that is useful for monitoring the behavior of TimeStepper.
For example, the current gas discharge models in chombo-discharge print the maximum electric field and density at each time step.
Note that printStepReport has (or should have!) no side-effects that affect the simulation state.
Regrid routines
The regrid routines in TimeStepper must, in combination, be able to transfer the simulation between old and new grids.
For an explanation to how regridding occurs in chombo-discharge, see Regridding.
In particular, when regrids occur the old grids are eventually destroyed so it is necessary to cache the old-grid simulation states so that we have something to interpolate from whan transfer the state to the new grids.
preRegrid
preRegrid should any necessary pre-regrid operations that are necessary in order to call the regrid.
This will virtually always include caching the old-grid simulation state, both for the solvers and also for internal data in TimeStepper
Solvers usually know how to do this, and in some cases this function can be deceptively simple, as illustrated by the implementation of this function in Advection-diffusion model:
void
AdvectionDiffusionStepper::preRegrid(const int a_lbase, const int a_oldFinestLevel)
{
CH_TIME("AdvectionDiffusionStepper::preRegrid");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::preRegrid" << endl;
}
m_solver->preRegrid(a_lbase, a_oldFinestLevel);
}
Other implementations of TimeStepper may have substantially more complicated preRegrid routines.
regrid
regrid is the function that performs an actual regrid operation.
At the time when regrid is called, the old grids are already destroyed and are only available through the cached data.
Solvers are usually implemented with their own regrid routine, and if the only things that need to be regridded are the solvers, the implementation of this routine can be comparatively simple, as illustrated below for the Advection-diffusion model:
void
AdvectionDiffusionStepper::regrid(const int a_lmin, const int a_oldFinestLevel, const int a_newFinestLevel)
{
CH_TIME("AdvectionDiffusionStepper::regrid");
if (m_verbosity > 5) {
pout() << "AdvectionDiffusionStepper::regrid" << endl;
}
// Regrid CDR solver and set up the flow fields
m_solver->regrid(a_lmin, a_oldFinestLevel, a_newFinestLevel);
m_solver->setSource(0.0);
m_solver->setEbFlux(0.0);
if (m_solver->isDiffusive()) {
m_solver->setDiffusionCoefficient(m_diffCo);
}
if (m_solver->isMobile()) {
m_solver->setVelocity(m_velocity);
}
}
For other TimeStepper implementations this routine can become much more complex.
The Îto-KMC plasma model, for example, will regrid not only solver data but also internal mesh and particle data, recompute conductivities, deposit particles, handle superparticles, and prepare the simulation state for the next time step.
postRegrid
The postRegrid is called after regrid has completed and can be used to perform any post-regrid specific procedures.
This function is not a pure function, and an implementation of this function is therefore not a requirement.
Load balancing routines
The default load-balancing method in chombo-discharge is to distribute grid patches equally among ranks, respecting a space-filling Morton curve on each grid level.
However, during the regrid step, Driver will check if meshes should be load balanced using different heuristics.
This load balancing can be done separately for each Realm, and in this case the MPI ranks will have different patch ownership in different grid sets.
If a realm should be load balanced with a different method than the default load balancing scheme, then TimeStepper can take a DisjointBoxLayout which originally load balanced using the patch volume, and regenerate the patch-to-rank ownership for the grids.
This functionality is implemented through two routines:
loadBalanceThisRealmwhich checks if a specific Realm should be load balanced.loadBalanceBoxeswhich load balances the boxes on the specified Realm.
Note that these functions are not pure functions, and it is perfect fine to use their default implementation, in which case each MPI rank gets approximately the same number of grid patches.
loadBalanceThisRealm
The function signature for this function is
/*!
@brief Load balancing query for a specified realm. If this returns true for a_realm, load balancing routines will be called during regrids.
@param[in] a_realm Realm name
*/
virtual bool
loadBalanceThisRealm(const std::string a_realm) const;
This function must return true if the input Realm (a_realm) should be load balanced.
loadBalanceBoxes
If loadBalanceThisRealm returns true, the following function is responsible for actually regenerating the grids:
/*!
@brief Load balance grid boxes for a specified realm.
@param[out] a_procs MPI ranks owning the various grid boxes.
@param[out] a_boxes Grid boxes on every level (obtain them with a_grids[lvl].boxArray())
@param[in] a_realm Realm identifier
@param[in] a_grids Original grids
@param[in] a_lmin Coarsest grid level that changed
@param[in] a_finestLevel New finest grid level
@details This is only called by Driver if TimeStepper::loadBalanceThisRealm(a_realm) returned true. The default implementation
uses volume-based loads for the grid patches. If the user wants to load balance boxes on a realm, this routine must be overwritten and
he should compute loads for the various patches in a_grids and call LoadBalancing::makeBalance on each level. It is up to the user/programmer
to decide if load balancing should be done independently on each level, or if loads per MPI rank are accumulated across levels.
*/
virtual void
loadBalanceBoxes(Vector<Vector<int>>& a_procs,
Vector<Vector<Box>>& a_boxes,
const std::string a_realm,
const Vector<DisjointBoxLayout>& a_grids,
const int a_lmin,
const int a_finestLevel);
This is called if loadBalanceThisRealm evaluates to true, and in this case the TimeStepper should compute a new set of rank ownership for the input grid boxes.
Observe that loadBalanceBoxes occurs for the entire AMR hierarchy, where the outer vector of a_procs and a_boxes is the grid level, and the inner vectors describe the ownership of each box.
The default implementation of this function ensures that when we load balance a level, we account for the accumulated load on coarser levels (see Default implementation of loadBalanceBoxes.).
/*!
@brief Load balancing query for a specified realm. If this returns true for a_realm, load balancing routines will be called during regrids.
@param[in] a_realm Realm name
*/
virtual bool
loadBalanceThisRealm(const std::string a_realm) const;
In the above, we use the Loads class to hold the computational load for each rank, and on each level we compute the load for each patch to be equal to the number of grid cells in the patch.
This is later load balanced in LoadBalancing::makeBalance, which is a routine that ensures that when we assign boxes on some grid level \(l\), we account for loads already assigned on coarser grid levels.