File solver.hxx

Defines

BOUT_NO_USING_NAMESPACE_BOUTGLOBALS

Typedefs

using rhsfunc = int (*)(BoutReal)

RHS function pointer.

using PhysicsPrecon = int (*)(BoutReal t, BoutReal gamma, BoutReal delta)

User-supplied preconditioner function.

using Jacobian = int (*)(BoutReal t)

User-supplied Jacobian function.

using TimestepMonitorFunc = int (*)(Solver *solver, BoutReal simtime, BoutReal lastdt)

Solution monitor, called each timestep.

using SolverType = std::string

Enums

enum SOLVER_VAR_OP

Values:

LOAD_VARS

LOAD_DERIVS

SET_ID

SAVE_VARS

SAVE_DERIVS

enum MonitorPosition

A type to set where in the list monitors are added.

Values:

BACK

FRONT

Variables

constexpr auto SOLVERCVODE = "cvode"

constexpr auto SOLVERPVODE = "pvode"

constexpr auto SOLVERIDA = "ida"

constexpr auto SOLVERPETSC = "petsc"

constexpr auto SOLVERSLEPC = "slepc"

constexpr auto SOLVERKARNIADAKIS = "karniadakis"

constexpr auto SOLVERRK4 = "rk4"

constexpr auto SOLVEREULER = "euler"

constexpr auto SOLVERRK3SSP = "rk3ssp"

constexpr auto SOLVERPOWER = "power"

constexpr auto SOLVERARKODE = "arkode"

constexpr auto SOLVERIMEXBDF2 = "imexbdf2"

constexpr auto SOLVERSNES = "snes"

constexpr auto SOLVERRKGENERIC = "rkgeneric"

class Solver

#include <solver.hxx>

Interface to integrators, mainly for time integration

Creation

Solver is a base class and can’t be created directly:

Solver *solver = Solver(); // Error

Instead, use the create() static function:

Solver *solver = Solver::create(); // ok

By default this will use the options in the “solver” section of the options, equivalent to:

Options *opts = Options::getRoot()->getSection("solver");
Solver *solver = Solver::create(opts);

To use a different set of options, for example if there are multiple solvers, use a different option section:

Options *opts = Options::getRoot()->getSection("anothersolver");
Solver *anothersolver = Solver::create(opts);

Problem specification

The equations to be solved are specified in a PhysicsModel object

class MyProblem : public PhysicsModel {
  protected:
    // This function called once at beginning
    int init(bool restarting) {
      SOLVE_FOR(f);   // Specify variables to solve
      // Set f to initial value if needed
      // otherwise read from input options
      return 0;
    }

    // This function called to evaluate time derivatives
    int rhs(BoutReal t) {
       ddt(f) = 1.0; // Calculate time derivatives
       return 0;
    }
  private:
    Field3D f; // A variable to evolve
}

The init() and rhs() functions must be defined, but there are other functions which can be defined. See PhysicsModel documentation for details.

Create an object, then add to the solver:

MyProblem *prob = MyProblem();
solver->setModel(prob);

Running simulation

To run a calculation

solver->solve();

This will use NOUT and TIMESTEP in the solver options (specified during creation). If these are not present then the global options will be used.

To specify NOUT and TIMESTEP, pass the values to solve:

solver->solve(NOUT, TIMESTEP);

Subclassed by ArkodeSolver, CvodeSolver, EulerSolver, IdaSolver, IMEXBDF2, KarniadakisSolver, PetscSolver, PowerSolver, PvodeSolver, RK3SSP, RK4Solver, RKGenericSolver, SlepcSolver, SNESSolver, SplitRK

Public Functions

Solver(Options *opts = nullptr)

virtual ~Solver()

void setModel(PhysicsModel *model)

Specify physics model to solve. Currently only one model can be evolved by a Solver.

virtual void setRHS(rhsfunc f)

Set the RHS function.

void setPrecon(PhysicsPrecon f)

Specify a preconditioner (optional)

virtual void setJacobian(Jacobian j)

Specify a Jacobian (optional)

void setSplitOperator(rhsfunc fC, rhsfunc fD)

Split operator solves.

void addMonitor(Monitor *monitor, MonitorPosition pos = MonitorPosition::FRONT)

Add a monitor to be called regularly

The frequency at which the monitor is called is set by Monitor::timestep. By default this is every output timestep. When a new Monitor with a smaller timestep is added, the solver attemps to adjust its internal timestep to match, as well as adjusting the timesteps of the current set of Monitors to be multiples of the new timestep. If this is not possible, addMonitor will throw an exception.

Adding new Monitors after the Solver has been initialised is only possible if their timestep is a multiple of the Solver’s timestep. Smaller timesteps will throw an exception.

void removeMonitor(Monitor *monitor)

Remove a monitor function previously added.

void addTimestepMonitor(TimestepMonitorFunc monitor)

Add a monitor function to be called every timestep.

void removeTimestepMonitor(TimestepMonitorFunc monitor)

Remove a previously added timestep monitor.

void add(Field2D &v, const std::string &name)

Add a variable to be solved. This must be done in the initialisation stage, before the simulation starts.

void add(Field3D &v, const std::string &name)

void add(Vector2D &v, const std::string &name)

void add(Vector3D &v, const std::string &name)

virtual bool constraints()

Returns true if constraints available.

void constraint(Field2D &v, Field2D &C_v, std::string name)

Add constraint functions (optional). These link a variable v to a control parameter C_v such that v is adjusted to keep C_v = 0.

void constraint(Field3D &v, Field3D &C_v, std::string name)

void constraint(Vector2D &v, Vector2D &C_v, std::string name)

void constraint(Vector3D &v, Vector3D &C_v, std::string name)

virtual void setMaxTimestep(BoutReal dt)

Set a maximum internal timestep (only for explicit schemes)

virtual BoutReal getCurrentTimestep()

Return the current internal timestep.

int solve(int nout = -1, BoutReal dt = 0.0)

Start the solver. By default solve() uses options to determine the number of steps and the output timestep. If nout and dt are specified here then the options are not used

Parameters

• nout: Number of output timesteps
• dt: The time between outputs

int init(int nout, BoutReal tstep)

Initialise the solver NOTE: nout and tstep should be passed to run, not init. Needed because of how the PETSc TS code works

virtual int run() = 0

Run the solver, calling monitors nout times, at intervals of tstep. This function is called by solve(), and is specific to each solver type

This should probably be protected, since it shouldn’t be called by users.

virtual void resetInternalFields()

Should wipe out internal field vector and reset from current field object data.

virtual int n2Dvars() const

Number of 2D variables. Vectors count as 3.

virtual int n3Dvars() const

Number of 3D variables. Vectors count as 3.

int resetRHSCounter()

Get and reset the number of calls to the RHS function.

int resetRHSCounter_e()

Same but for explicit timestep counter - for IMEX.

int resetRHSCounter_i()

Same but fur implicit timestep counter - for IMEX.

bool splitOperator()

Test if this solver supports split operators (e.g. implicit/explicit)

void outputVars(Datafile &outputfile, bool save_repeat = true)

Add evolving variables to output (dump) file or restart file

Parameters

• outputfile: The file to add variable to
• save_repeat: If true, add variables with time dimension

Public Members

bool canReset = {false}

Public Static Functions

Solver *create(Options *opts = nullptr)

Create a Solver object. This uses the “type” option in the given Option section to determine which solver type to create.

Solver *create(const SolverType &type, Options *opts = nullptr)

Create a Solver object, specifying the type.

static void setArgs(int &c, char **&v)

Pass the command-line arguments. This static function is called by BoutInitialise, and puts references into protected variables. These may then be used by Solvers to control behavior

Public Static Attributes

constexpr MonitorPosition BACK = MonitorPosition::BACK

constexpr MonitorPosition FRONT = MonitorPosition::FRONT

Protected Functions

int getLocalN()

Calculate the number of evolving variables on this processor.

template<class T>
bool contains(const std::vector<VarStr<T>> &vars, const std::string &name)

Does vars contain a field with name?

int run_rhs(BoutReal t)

Run the user’s RHS function.

int run_convective(BoutReal t)

Calculate only the convective parts.

NOTE: This calls add_mms_sources.

int run_diffusive(BoutReal t, bool linear = true)

Calculate only the diffusive parts.

int call_monitors(BoutReal simtime, int iter, int NOUT)

Calls all monitor functions

There are two important things to note about how iter is passed along to each monitor:

• The solvers all start their iteration numbering from zero, so the initial state is calculated at iter = -1
• Secondly, iter is passed along to each monitor relative to that monitor’s period

In practice, this means that each monitor is called like:

monitor->call(solver, simulation_time,
              ((iter + 1) / monitor->period) - 1,
              NOUT / monitor->period);

e.g. for a monitor with period 10, passing iter = 9 will result in it being called with a value of (9 + 1)/10 - 1 == 0

int call_timestep_monitors(BoutReal simtime, BoutReal lastdt)

bool have_user_precon()

Do we have a user preconditioner?

int run_precon(BoutReal t, BoutReal gamma, BoutReal delta)

void load_vars(BoutReal *udata)

void load_derivs(BoutReal *udata)

void save_vars(BoutReal *udata)

void save_derivs(BoutReal *dudata)

void set_id(BoutReal *udata)

Field3D globalIndex(int localStart)

Returns a Field3D containing the global indices.

auto getMonitors() const

Get the list of monitors.

Protected Attributes

Options *options = {nullptr}

Settings to use during initialisation (set by constructor)

int NPES = {1}

Number of processors.

int MYPE = {0}

This processor’s index.

std::vector<VarStr<Field2D>> f2d

Vectors of variables to evolve.

std::vector<VarStr<Field3D>> f3d

std::vector<VarStr<Vector2D>> v2d

std::vector<VarStr<Vector3D>> v3d

bool has_constraints = {false}

Can this solver handle constraints? Set to true if so.

bool initialised = {false}

Has init been called yet?

BoutReal simtime = {0.0}

Current simulation time.

int iteration = {0}

Current iteration (output time-step) number.

bool monitor_timestep = {false}

Should timesteps be monitored?

BoutReal max_dt = {-1.0}

Maximum internal timestep.

Protected Static Attributes

int *pargc = nullptr

Number of command-line arguments.

char ***pargv = nullptr

Command-line arguments.

Private Functions

void add_mms_sources(BoutReal t)

void calculate_mms_error(BoutReal t)

void pre_rhs(BoutReal t)

Should be run before user RHS is called.

void post_rhs(BoutReal t)

Should be run after user RHS is called.

void loop_vars_op(Ind2D i2d, BoutReal *udata, int &p, SOLVER_VAR_OP op, bool bndry)

Loading data from BOUT++ to/from solver.

Perform an operation at a given Ind2D (jx,jy) location, moving data between BOUT++ and CVODE.

void loop_vars(BoutReal *udata, SOLVER_VAR_OP op)

Loop over variables and domain. Used for all data operations for consistency.

bool varAdded(const std::string &name)

Check if a variable has already been added.

BoutReal adjustMonitorPeriods(Monitor *monitor)

(Possibly) adjust the periods of monitor, and the monitors timesteps, returning the new Solver timestep

void finaliseMonitorPeriods(int &NOUT, BoutReal &output_timestep)

Fix all the monitor periods based on output_timestep, as well as adjusting NOUT and output_timestep to be consistent

Private Members

int rhs_ncalls = {0}

Number of calls to the RHS function.

int rhs_ncalls_e = {0}

Number of calls to the explicit (convective) RHS function.

int rhs_ncalls_i = {0}

Number of calls to the implicit (diffusive) RHS function.

int default_monitor_period = {1}

Default sampling rate at which to call monitors - same as output to screen.

BoutReal internal_timestep = {-1}

timestep - shouldn’t be changed after init is called.

PhysicsModel *model = {nullptr}

Physics model being evolved.

rhsfunc phys_run = {nullptr}

The user’s RHS function.

PhysicsPrecon prefunc = {nullptr}

The user’s preconditioner function.

bool split_operator = {false}

Is the physics model using separate convective (explicit) and diffusive (implicit) RHS functions?

rhsfunc phys_conv = {nullptr}

Convective part (if split operator)

rhsfunc phys_diff = {nullptr}

Diffusive part (if split operator)

bool is_nonsplit_model_diffusive = {true}

Should non-split physics models be treated as diffusive?

bool mms = {false}

Enable sources and solutions for Method of Manufactured Solutions.

bool mms_initialise = {false}

Initialise variables to the manufactured solution.

std::list<Monitor *> monitors

List of monitor functions.

std::list<TimestepMonitorFunc> timestep_monitors

List of timestep monitor functions.

Friends

template<class T>
bool operator==(const VarStr<T> &var, const std::string &name)

Does var represent field name?

template<class T>
int local_N_sum(int value, const VarStr<T> &f)

Helper function for getLocalN: return the number of points to evolve in f, plus the accumulator value

If f.evolve_bndry, includes the boundary (NB: not guard!) points

FIXME: This could be a lambda local to getLocalN with an auto argument in C++14

template<class T>
struct VarStr

#include <solver.hxx>

A structure to hold an evolving variable.

Public Members

bool constraint = {false}

T *var = {nullptr}

Does F_var represent a constraint?

T *F_var = {nullptr}

The evolving variable.

std::unique_ptr<T> MMS_err = {nullptr}

The time derivative or constraint on var.

CELL_LOC location = {CELL_DEFAULT}

Error for MMS.

bool covariant = {false}

For fields and vector components.

bool evolve_bndry = {false}

For vectors.

std::string name

Are the boundary regions being evolved?