$$Date$$

Abstract

The Method of Lines is a way of separating the time integration from the rest of an evolution scheme. This thorn is intended to take care of all the bookwork and provide some basic time integration methods, allowing for easy coupling of diﬀerent thorns.

The Method of Lines (MoL) converts a (system of) partial diﬀerential equation(s) into an ordinary diﬀerential equation containing some spatial diﬀerential operator. As an example, consider writing the hyperbolic system of PDE’s

$${\partial}_{t}q+{A}^{i}\left(q\right){\partial}_{i}B\left(q\right)=s\left(q\right)$$ | (1) |

in the alternative form

$${\partial}_{t}q=L\left(q\right),$$ | (2) |

which (assuming a given discretization of space) is an ODE.

Given this separation of the time and space discretizations, well known stable ODE integrators such as Runge-Kutta can be used to do the time integration. This is more modular (allowing for simple extensions to higher order methods), more stable (as instabilities can now only arise from the spatial discretization or the equations themselves) and also avoids the problems of retaining high orders of convergence when coupling diﬀerent physical models.

MoL can be used for hyperbolic, parabolic and even elliptic problems (although I deﬁnitely don’t recommend the latter). As it currently stands it is set up for systems of equations in the ﬁrst order type form of equation (2). If you want to implement a multilevel scheme such as leapfrog it is not obvious to me that MoL is the thing to use. However if you have lots of thorns that you want to interact, for example ADM_BSSN and a hydro code plus maybe EM or a scalar ﬁeld, and they can easily be written in this sort of form, then you probably want to use MoL.

This thorn is meant to provide a simple interface that will implement the MoL inside Cactus as transparently as possible. It will initially implement only the optimal Runge-Kutta time integrators (which are TVD up to RK3, so suitable for hydro) up to fourth order and iterated Crank Nicholson. All of the interaction with the MoL thorn should occur directly through the scheduler. For example, all synchronization steps should now be possible at the schedule level. This is essential for interacting cleanly with diﬀerent drivers, especially to make Mesh Reﬁnement work.

For more information on the Method of Lines the most comprehensive references are the works of Jonathan Thornburg [1, 2] - see especially section 7.3 of the thesis. From the CFD viewpoint the review of ENO methods by Shu, [3], has some information. For relativistic ﬂuids the paper of Neilsen and Choptuik [4] is also quite good.

For those who used the old version of MoL, this version is unfortunately slightly more eﬀort to use. That is, for most methods you’ll now have to set 4 parameters instead of just one.

If you already have a thorn that uses the method of lines, then there are four main parameters that are relevant to change the integration method. The keyword MoL_ODE_Method chooses between the diﬀerent methods. Currently supported are RK2, RK3, ICN, ICN-Avg and Generic. These are second order Runge-Kutta, third order Runge-Kutta, Iterative Crank Nicholson, Iterative Crank Nicholson with averaging, and the generic Shu-Osher type Runge-Kutta methods. To switch between the diﬀerent types of generic methods there is also the keyword Generic_Type which is currently restricted to RK for the standard TVD Runge-Kutta methods (ﬁrst to fourth order) and ICN for the implementation of the Iterative Crank Nicholson method in generic form.

Full descriptions of the currently implemented methods are given in section 4.

The parameter MoL_Intermediate_Steps controls the number of intermediate steps for the ODE solver. For the generic Runge-Kutta solvers it controls the order of accuracy of the method. For the ICN methods this parameter controls the number of iterations taken, which does not check for stability. This parameter defaults to 3.

The parameter MoL_Num_Scratch_Levels controls the amount of scratch space used. If this is insuﬃcient for the method selected there will be an error at parameter checking time. This parameter defaults to 0, as no scratch space is required for the eﬃcient ICN and Runge-Kutta 2 and 3 solvers. For the generic solvers this must be at least MoL_Intermediate_Steps - 1.

Another parameter is MoL_Memory_Always_On which switches on memory for the scratch space always if true and only during evolution if false. This defaults to true for speed reasons; the memory gains are likely to be limited unless you’re doing something very memory intensive at initialization or analysis.

There is also a parameter MoL_NaN_Check that will check your RHS grid functions for NaNs using the NaNChecker thorn from CactusUtils. This will make certain that you ﬁnd the exact grid function computing the ﬁrst NaN; of course, this may not be the real source of your problem.

The parameter disable_prolongation only does anything if you are using mesh reﬁnement, and in particular Carpet. With mesh reﬁnement it may be necessary to disable prolongation in intermediate steps of MoL. This occurs when evolving systems containing second spatial derivatives. This is done by default in MoL. If your system is purely ﬁrst order in space and time you may wish to set this to "no".

Ideally, initial data thorns should always set initial data at all time levels. However, sometimes initial data thorns fail to do this. In this case you can do one of three things:

- Fix the initial data thorn. This is the best solution.
- If you’re using Carpet, it has some facilities to take forward/backward time steps to initialize multiple time levels. See the Carpet parameters init_each_timelevel and init_3_timelevels for details.
- Finally, if you set (the MoL parameter) initial_data_is_crap, MoL will copy the current time
level of all variables it knows about (more precisely, using the terminology of section 2.2, all evolved,
save-and-restore, and constrained variables which have multiple time levels) to all the past time
levels. Note that this copies the same data to each past time level; this will be wrong
if your spacetime is time-dependent!
If enabled, the copy happens in the CCTK_POSTINITIAL schedule bin. By default this happens before the MoL_PostStep schedule group; the parameter copy_ID_after_MoL_PostStep can be used to change this to after MoL_PostStep.

To port an existing thorn using the method of lines, or to write a new thorn using it, should hopefully be relatively simple. As an example, within the MoL arrangement is WaveMoL which duplicates the WaveToy thorn given by CactusWave in a form suitable for use by MoL. In this section, “the thorn” will mean the user thorn doing the physics.

We start with some terminology. Grid functions are split into four cateogories.

- The ﬁrst are those that are evolved using a MoL form. That is, a right hand side is calculated and the variable updated using it. These we call evolved variables.
- The second category are those variables that are set by a thorn at every intermediate step of the evolution, usually to respect the constraints. Examples of these include the primitive variables in a hydrodynamics code. Another example would be the gauge variables if these were set by constraints at every intermediate step (which is slightly artiﬁcial; the usual example would be the use of maximal slicing, which is only applied once every $N$ complete steps). These are known as constrained variables.
- The third category are those variables that a thorn depends on but does not set or evolve. An example would include the metric terms considered from a thorn evolving matter. Due to the way that MoL deals with these, they are known as Save and Restore variables.
- The ﬁnal category are those variables that do not interact with MoL. These would include temporary variables for analysis or setting up the initial data. These can safely be ignored.

As a generic rule of thumb, variables for which you have a time evolution equation are evolved (obviously), variables which your thorn sets but does not evolve are constrained, and any other variables which your thorn reads during evolution is a Save and Restore variable.

MoL needs to know every GF that falls in one of the ﬁrst three groups. If a GF is evolved by one thorn but is a constrained variable in another (for example, the metric in full GR Hydro) then each thorn should register the function as they see it. For example, the hydro thorn will register the metric as a Save and Restore variable and the spacetime thorn will register the metric as an evolved variable. The diﬀerent variable categories are given the priority evolved, constrained, Save and Restore. So if a variable is registered as belonging in two diﬀerent categories, it is always considered by MoL to belong to the category with the highest priority.

MoL needs to know the total number of GFs in each category at parameter time. To do this, your thorn needs to use some accumulator parameters from MoL. As an example, here are the paramaters from WaveMoL:

shares: MethodOfLines

USES CCTK_INT MoL_Num_Evolved_Vars

USES CCTK_INT MoL_Num_Constrained_Vars

USES CCTK_INT MoL_Num_SaveAndRestore_Vars

restricted:

CCTK_INT WaveMoL_MaxNumEvolvedVars \

"The maximum number of evolved variables used by WaveMoL" \

ACCUMULATOR-BASE=MethodofLines::MoL_Num_Evolved_Vars

{

5:5 :: "Just 5: phi and the four derivatives"

} 5

CCTK_INT WaveMoL_MaxNumConstrainedVars \

"The maximum number of constrained variables used by WaveMoL" \

ACCUMULATOR-BASE=MethodofLines::MoL_Num_Constrained_Vars

{

0:1 :: "A small range, depending on testing or not"

} 1

CCTK_INT WaveMoL_MaxNumSandRVars \

"The maximum number of save and restore variables used by WaveMoL" \

ACCUMULATOR-BASE=MethodofLines::MoL_Num_SaveAndRestore_Vars

{

0:1 :: "A small range, depending on testing or not"

} 1

USES CCTK_INT MoL_Num_Evolved_Vars

USES CCTK_INT MoL_Num_Constrained_Vars

USES CCTK_INT MoL_Num_SaveAndRestore_Vars

restricted:

CCTK_INT WaveMoL_MaxNumEvolvedVars \

"The maximum number of evolved variables used by WaveMoL" \

ACCUMULATOR-BASE=MethodofLines::MoL_Num_Evolved_Vars

{

5:5 :: "Just 5: phi and the four derivatives"

} 5

CCTK_INT WaveMoL_MaxNumConstrainedVars \

"The maximum number of constrained variables used by WaveMoL" \

ACCUMULATOR-BASE=MethodofLines::MoL_Num_Constrained_Vars

{

0:1 :: "A small range, depending on testing or not"

} 1

CCTK_INT WaveMoL_MaxNumSandRVars \

"The maximum number of save and restore variables used by WaveMoL" \

ACCUMULATOR-BASE=MethodofLines::MoL_Num_SaveAndRestore_Vars

{

0:1 :: "A small range, depending on testing or not"

} 1

This should give the maximum number of variables that your thorn will register.

Every thorn should register every grid function that it uses even if you expect it to be registered again by a diﬀerent thorn. For example, a hydro thorn would register the metric variables as Save and Restore, whilst the spacetime evolution thorn would register them as evolved (in ADM) or constrained (in ADM_BSSN), both of which have precedence. To register your GFs with MoL schedule a routine in the bin MoL_Register which just contains the relevant function calls. For an evolved variable the GF corresponding to the update term ($L\left(q\right)$ in equation (2)) should be registered at the same time. The appropriate functions are given in section 5.

These functions are provided using function aliasing. For details on using function aliasing, see sections B10.5 and F2.2.3 of the UsersGuide. For the case of real GFs, you simply add the following lines to your interface.ccl:

##########################################

### PROTOTYPES - DELETE AS APPLICABLE! ###

##########################################

CCTK_INT FUNCTION MoLRegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)

CCTK_INT FUNCTION MoLRegisterEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)

CCTK_INT FUNCTION MoLRegisterConstrained(CCTK_INT ConstrainedIndex)

CCTK_INT FUNCTION MoLRegisterSaveAndRestore(CCTK_INT SandRIndex)

CCTK_INT FUNCTION MoLRegisterEvolvedGroup(CCTK_INT EvolvedIndex, \

CCTK_INT RHSIndex)

CCTK_INT FUNCTION MoLRegisterConstrainedGroup(CCTK_INT ConstrainedIndex)

CCTK_INT FUNCTION MoLRegisterSaveAndRestoreGroup(CCTK_INT SandRIndex)

CCTK_INT FUNCTION MoLChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)

CCTK_INT FUNCTION MoLChangeToEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)

CCTK_INT FUNCTION MoLChangeToConstrained(CCTK_INT ConstrainedIndex)

CCTK_INT FUNCTION MoLChangeToSaveAndRestore(CCTK_INT SandRIndex)

CCTK_INT FUNCTION MoLChangeToNone(CCTK_INT RemoveIndex)

#############################################

### USE STATEMENT - DELETE AS APPLICABLE! ###

#############################################

USES FUNCTION MoLRegisterEvolved

USES FUNCTION MoLRegisterEvolvedSlow

USES FUNCTION MoLRegisterConstrained

USES FUNCTION MoLRegisterSaveAndRestore

USES FUNCTION MoLRegisterEvolvedGroup

USES FUNCTION MoLRegisterConstrainedGroup

USES FUNCTION MoLRegisterSaveAndRestoreGroup

USES FUNCTION MoLChangeToEvolved

USES FUNCTION MoLChangeToConstrained

USES FUNCTION MoLChangeToSaveAndRestore

USES FUNCTION MoLChangeToNone

### PROTOTYPES - DELETE AS APPLICABLE! ###

##########################################

CCTK_INT FUNCTION MoLRegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)

CCTK_INT FUNCTION MoLRegisterEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)

CCTK_INT FUNCTION MoLRegisterConstrained(CCTK_INT ConstrainedIndex)

CCTK_INT FUNCTION MoLRegisterSaveAndRestore(CCTK_INT SandRIndex)

CCTK_INT FUNCTION MoLRegisterEvolvedGroup(CCTK_INT EvolvedIndex, \

CCTK_INT RHSIndex)

CCTK_INT FUNCTION MoLRegisterConstrainedGroup(CCTK_INT ConstrainedIndex)

CCTK_INT FUNCTION MoLRegisterSaveAndRestoreGroup(CCTK_INT SandRIndex)

CCTK_INT FUNCTION MoLChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)

CCTK_INT FUNCTION MoLChangeToEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)

CCTK_INT FUNCTION MoLChangeToConstrained(CCTK_INT ConstrainedIndex)

CCTK_INT FUNCTION MoLChangeToSaveAndRestore(CCTK_INT SandRIndex)

CCTK_INT FUNCTION MoLChangeToNone(CCTK_INT RemoveIndex)

#############################################

### USE STATEMENT - DELETE AS APPLICABLE! ###

#############################################

USES FUNCTION MoLRegisterEvolved

USES FUNCTION MoLRegisterEvolvedSlow

USES FUNCTION MoLRegisterConstrained

USES FUNCTION MoLRegisterSaveAndRestore

USES FUNCTION MoLRegisterEvolvedGroup

USES FUNCTION MoLRegisterConstrainedGroup

USES FUNCTION MoLRegisterSaveAndRestoreGroup

USES FUNCTION MoLChangeToEvolved

USES FUNCTION MoLChangeToConstrained

USES FUNCTION MoLChangeToSaveAndRestore

USES FUNCTION MoLChangeToNone

Note that the list of parameters not complete; see the section on parameters for the use of arrays. However, the list of functions is, and is expanded on in section 5. MoL will check whether a group or variable is a GF or an array and whether it is real.

Having done that, one routine (or group of routines) which we’ll here call Thorn_CalcRHS must be deﬁned. This does all the ﬁnite diﬀerencing that you’d usually do, applied to $q$, and ﬁnds the right hand sides which are stored in $L$. This routine should be scheduled in MoL_CalcRHS. The precise order that these are scheduled should not matter, because no updating of any of the user thorns $q$ will be done until after all the RHSs are calculated. Important note: all the ﬁnite diﬀerencing must be applied to the most recent time level $q$ and not to the previous time level ${q}_{p}$ as you would normally do. Don’t worry about setting up the data before the calculation, as MoL will do that automatically.

Finally, if you have some things that have to be done after each update to an intermediate
level, these should be scheduled in MoL_PostStep. Examples of things that need to go here
include the recalculation of primitive variables for hydro codes, the application of boundary
conditions^{1} ,
the solution of elliptic equations (although this would be a very expensive place to solve them, some sets
of equations might require the updating of some variables by constraints in this fashion). When
applying boundary conditions the cleanest thing to do is to write a routine applying the symmetries
to the appropriate GFs and, when calling it from the scheduler, adding the SYNC statement to
the appropriate groups. An example is given by the routine WaveToyMoL_Boundaries in thorn
WaveMoL.

Points to note. The thorn routine Thorn_CalcRHS does not need to know and in fact should deﬁnitely not know where precisely in the MoL step it is. It just needs to know that it is receiving some intermediate data stored in the GFs $q$ and that it should return the RHS $L\left(q\right)$. All the book-keeping to ensure that it is passed the correct intermediate state at that the GFs contain the correct data at the end of the MoL step will be dealt with by the MoL thorn and the ﬂesh.

When using a multirate scheme the thorn routine Thorn_CalcRHS must check the MoL grid scalars MoL::MoL_SlowStep and MoL::MoL_SlowPostStep which MoL sets to true (non-zero) if it is time for the slow sector to compute its RHS or to apply its poststep routine. MoL::MoL_SlowPostStep is always true outside of the RHS loop, eg. in CCTK_POST_RESTRICT.

If you want to try adding a new evolution method to MoL the simplest way is to use the generic table option to specify it completely in the parameter ﬁle - no coding is required at all.

All the generic methods evolve the equation

$${\partial}_{t}q=L\left(q\right)$$ | (3) |

using the following algorithm for an $N$-step method:

$$\begin{array}{rcll}{q}^{\left(0\right)}& =& {q}^{n},& \text{}\\ {q}^{\left(i\right)}& =& \sum _{k=0}^{i-1}\left({\alpha}_{ik}{q}^{\left(k\right)}\right)+\Delta t{\beta}_{i-1}L\left({q}^{\left(i-1\right)}\right),\phantom{\rule{2em}{0ex}}i=1,\dots ,N,& \text{(4)}\text{}\text{}\\ {q}^{n+1}& =& {q}^{\left(N\right)}.& \text{}\end{array}$$

This method is completely speciﬁed by $N$ (GenericIntermediateSteps) and the $\alpha $ (GenericAlphaCoeffs) and $\beta $ (GenericBetaCoeffs) arrays. The names in parentheses give the keys in a table that MoL will use. This table is created from the string parameter Generic_Method_Descriptor.

As an example, the standard TVD RK2 method that is implemented both in optimized and generic form is written as

$$\begin{array}{rcll}{q}^{\left(1\right)}& =& {q}^{n}+\Delta tL\left({q}^{n}\right),& \text{(5)}\text{}\text{}\\ {q}^{n+1}& =& \frac{1}{2}\left({q}^{n}+{q}^{\left(1\right)}+\Delta tL\left({q}^{\left(1\right)}\right)\right).& \text{(6)}\text{}\text{}\end{array}$$

To implement this using the generic table options, use

methodoflines::MoL_Intermediate_Steps = 2

methodoflines::MoL_Num_Scratch_Levels = 1

methodoflines::Generic_Method_Descriptor = \

"GenericIntermediateSteps = 2 \

GenericAlphaCoeffs = { 1.0 0.0 0.5 0.5 } \

GenericBetaCoeffs = { 1.0 0.5 }"

methodoflines::MoL_Num_Scratch_Levels = 1

methodoflines::Generic_Method_Descriptor = \

"GenericIntermediateSteps = 2 \

GenericAlphaCoeffs = { 1.0 0.0 0.5 0.5 } \

GenericBetaCoeffs = { 1.0 0.5 }"

The number of steps speciﬁed in the table must be the same as MoL_Intermediate_Steps, and the number of scratch levels should be at least MoL_Intermediate_Steps - 1.

The generic methods are somewhat ineﬃcient for use in production runs, so it is frequently better to write an optimized version once you are happy with the method. To do this you should

- write your code into a new ﬁle, called from the scheduler under the alias MoL_Add,
- make certain that at each intermediate step the correct values of cctk_time and cctk_delta_time are set in SetTime.c for mesh reﬁnement, boundary conditions and so on,
- make certain that you check for the number of intermediate steps in ParamCheck.c.

As a fairly extended example of how to use MoL I’ll outline how ADM_BSSN works in this context. The actual implementation of this is given in the thorn AEIThorns/BSSN_MoL.

As normal the required variables are deﬁned in the interface.ccl ﬁle, together with the associated source terms. For example, the conformal factor and source are deﬁned by

real ADM_BSSN_phi type=GF timelevels=2

{

ADM_BS_phi

} "ADM_BSSN_phi"

real ADM_BSSN_sources type=GF

{

...,

adm_bs_sphi,

...

}

{

ADM_BS_phi

} "ADM_BSSN_phi"

real ADM_BSSN_sources type=GF

{

...,

adm_bs_sphi,

...

}

Also in this ﬁle we write the function aliasing prototypes.

Once the sources are deﬁned the registration with MoL is required, for which the essential ﬁle is MoLRegister.c. In the ADM_BSSN system the standard metric coeﬃcients ${g}_{ij}$ are not evolved, and neither are the standard extrinsic curvature components ${K}_{ij}$. However these are used by ADM_BSSN in a number of places, and are calculated from evolved quantities at the appropriate points. In the MoL terminology these variables are constrained. As the appropriate storage is deﬁned in thorn ADMBase, the actual calls have the form

ierr += MoLRegisterConstrained(CCTK_VarIndex("ADMBase::kxx"));

The actual evolved variables include things such as the conformal factor. This, and the appropriate source term, is deﬁned in thorn ADM_BSSN, and so the call has the form

ierr += MoLRegisterEvolved(CCTK_VarIndex("adm_bssn::ADM_BS_phi"),

CCTK_VarIndex("adm_bssn::adm_bs_sphi"));

CCTK_VarIndex("adm_bssn::adm_bs_sphi"));

As well as the evolved variables, and those constrained variables such as the metric, there are the gauge variables. Precisely what status these have depends on how they are set. If harmonic or 1+log slicing is used then the lapse is evolved:

ierr += MoLRegisterEvolved(CCTK_VarIndex("ADMBase::alp"),

CCTK_VarIndex("adm_bssn::adm_bs_salp"));

CCTK_VarIndex("adm_bssn::adm_bs_salp"));

A matter density $\rho $ might not require such a high order scheme and can be evolved using a multi-rate scheme

ierr += MoLRegisterEvolvedSlow(CCTK_VarIndex("GRHydro::dens"),

CCTK_VarIndex("GRHydro::dendsrhs"));

CCTK_VarIndex("GRHydro::dendsrhs"));

If maximal or static slicing is used then the lapse is a constrained
variable^{2} :

ierr += MoLRegisterConstrained(CCTK_VarIndex("ADMBase::alp"));

Finally, if none of the above apply we assume that the lapse is evolved in some unknown fashion, and so it must be registered as a Save and Restore variable:

ierr += MoLRegisterSaveAndRestore(CCTK_VarIndex("ADMBase::alp"));

However, it is perfectly possible that we may wish to change how we deal with the gauge during the evolution. This is dealt with in the ﬁle PreLoop.F. If the slicing changes then the appropriate routine is called. For example, if we want to use 1+log evolution then we call

call CCTK_VarIndex(lapseindex,"ADMBase::alp")

call CCTK_VarIndex(lapserhsindex,"adm_bssn::adm_bs_salp")

ierr = ierr + MoLChangeToEvolved(lapseindex, lapserhsindex)

call CCTK_VarIndex(lapserhsindex,"adm_bssn::adm_bs_salp")

ierr = ierr + MoLChangeToEvolved(lapseindex, lapserhsindex)

It is not required to tell MoL what the lapse is changing from, or indeed if it is changing at all; MoL will work this out for itself.

Finally there are the routines that we wish to apply after every intermediate step. These are ADM_BSSN_removetrA which enforces various constraints (such as the tracefree conformal extrinsic curvature remaining trace free), ADM_BSSN_Boundaries which applies symmetry boundary conditions as well as various others (such as some of the radiative boundary conditions). Note all the calls to SYNC at this point. We also convert from the updated BSSN variables back to the standard ADM variables in ADM_BSSN_StandardVariables, and also update the time derivative of the lapse in ADM_BSSN_LapseChange.

The default method is Iterative Crank-Nicholson. There are many ways of implementing this. The standard "ICN" and "Generic"/"ICN" methods both implement the following, assuming an $N$ iteration method:

$$\begin{array}{rcll}{q}^{\left(0\right)}& =& {q}^{n},& \text{(7)}\text{}\text{}\\ {q}^{\left(i\right)}& =& {q}^{\left(0\right)}+\frac{\Delta t}{2}L\left({q}^{\left(i-1\right)}\right),\phantom{\rule{1em}{0ex}}i=1,\dots ,N-1,& \text{(8)}\text{}\text{}\\ {q}^{\left(N\right)}& =& {q}^{\left(N-1\right)}+\Delta tL\left({q}^{\left(N-1\right)}\right),& \text{(9)}\text{}\text{}\\ {q}^{n+1}& =& {q}^{\left(N\right)}& \text{(10)}\text{}\text{}\end{array}$$

The “averaging” ICN method "ICN-avg" instead calculates intermediate steps before averaging:

$$\begin{array}{rcll}{q}^{\left(0\right)}& =& {q}^{n},& \text{(11)}\text{}\text{}\\ {\stackrel{\u0303}{q}}^{\left(i\right)}& =& \frac{1}{2}\left({q}^{\left(i\right)}+{q}^{n}\right),\phantom{\rule{1em}{0ex}}i=0,\dots ,N-1& \text{(12)}\text{}\text{}\\ {q}^{\left(i\right)}& =& {q}^{\left(0\right)}+\Delta tL\left({\stackrel{\u0303}{q}}^{\left(N-1\right)}\right),& \text{(13)}\text{}\text{}\\ {q}^{n+1}& =& {q}^{\left(N\right)}& \text{(14)}\text{}\text{}\end{array}$$

The Runge-Kutta methods are those typically used in hydrodynamics by, e.g., Shu and others — see [3] for example. Explicitly the ﬁrst order method is the Euler method:

$$\begin{array}{rcll}{q}^{\left(0\right)}& =& {q}^{n},& \text{(15)}\text{}\text{}\\ {q}^{\left(1\right)}& =& {q}^{\left(0\right)}+\Delta tL\left({\stackrel{\u0303}{q}}^{\left(0\right)}\right),& \text{(16)}\text{}\text{}\\ {q}^{n+1}& =& {q}^{\left(1\right)}.& \text{(17)}\text{}\text{}\end{array}$$

The second order method is:

$$\begin{array}{rcll}{q}^{\left(0\right)}& =& {q}^{n},& \text{(18)}\text{}\text{}\\ {q}^{\left(1\right)}& =& {q}^{\left(0\right)}+\Delta tL\left({q}^{\left(0\right)}\right),& \text{(19)}\text{}\text{}\\ {q}^{\left(2\right)}& =& \frac{1}{2}\left({q}^{\left(0\right)}+{q}^{\left(1\right)}+\Delta tL\left({q}^{\left(1\right)}\right)\right),& \text{(20)}\text{}\text{}\\ {q}^{n+1}& =& {q}^{\left(2\right)}.& \text{(21)}\text{}\text{}\end{array}$$

The third order method is:

$$\begin{array}{rcll}{q}^{\left(0\right)}& =& {q}^{n},& \text{(22)}\text{}\text{}\\ {q}^{\left(1\right)}& =& {q}^{\left(0\right)}+\Delta tL\left({q}^{\left(0\right)}\right),& \text{(23)}\text{}\text{}\\ {q}^{\left(2\right)}& =& \frac{1}{4}\left(3{q}^{\left(0\right)}+{q}^{\left(1\right)}+\Delta tL\left({q}^{\left(1\right)}\right)\right),& \text{(24)}\text{}\text{}\\ {q}^{\left(3\right)}& =& \frac{1}{3}\left({q}^{\left(0\right)}+2{q}^{\left(2\right)}+2\Delta tL\left({q}^{\left(2\right)}\right)\right),& \text{(25)}\text{}\text{}\\ {q}^{n+1}& =& {q}^{\left(3\right)}.& \text{(26)}\text{}\text{}\end{array}$$

The fourth order method, which is not strictly TVD, is:

$$\begin{array}{rcll}{q}^{\left(0\right)}& =& {q}^{n},& \text{(27)}\text{}\text{}\\ {q}^{\left(1\right)}& =& {q}^{\left(0\right)}+\frac{1}{2}\Delta tL\left({q}^{\left(0\right)}\right),& \text{(28)}\text{}\text{}\\ {q}^{\left(2\right)}& =& {q}^{\left(0\right)}+\frac{1}{2}\Delta tL\left({q}^{\left(1\right)}\right),& \text{(29)}\text{}\text{}\\ {q}^{\left(3\right)}& =& {q}^{\left(0\right)}+\Delta tL\left({q}^{\left(2\right)}\right),& \text{(30)}\text{}\text{}\\ {q}^{\left(4\right)}& =& \frac{1}{6}\left(-2{q}^{\left(0\right)}+2{q}^{\left(1\right)}+4{q}^{\left(2\right)}+2{q}^{\left(3\right)}+\Delta tL\left({q}^{\left(3\right)}\right)\right),& \text{(31)}\text{}\text{}\\ {q}^{n+1}& =& {q}^{\left(4\right)}.& \text{(32)}\text{}\text{}\end{array}$$

A scheme for coupling diﬀerent parts of a system of equations

$$\begin{array}{rcll}{\partial}_{t}g& =& F\left(g,q\right),& \text{(33)}\text{}\text{}\\ {\partial}_{t}q& =& G\left(g,q\right),& \text{(34)}\text{}\text{}\end{array}$$

representing eg. spacetime and matter variables, respectively, with diﬀerent RK integrators is given by multirate RK schemes (e.g. [5, 6]). Here, we pursuit the simple Ansatz of performing one RK2 intermediate RHS evaluation for two RK4 intermediate RHS evaluations. That is, the additional RK4 intermediate RHS evaluations simple use the results from the last intermediate RK2 step.

To be more explicit, given the equation

$${\partial}_{t}y=f\left(t,y\right)\phantom{\rule{0.3em}{0ex}},$$ | (35) |

where $f$ corresponds to the RHS possibly including spatial derivatives, we write a generic RK scheme according to

$$\begin{array}{rcll}{y}_{n+1}& =& {y}_{n}+\Delta t\sum _{i=1}^{s}{b}_{i}\phantom{\rule{0.3em}{0ex}}{k}_{i}\phantom{\rule{0.3em}{0ex}},& \text{(36)}\text{}\text{}\\ {k}_{i}& =& f\left({t}_{n}+{c}_{i}\Delta t\phantom{\rule{0.3em}{0ex}},{y}_{n}+\Delta t\sum _{j=1}^{s}{a}_{ij}{k}_{j}\right)\phantom{\rule{0.3em}{0ex}}.& \text{(37)}\text{}\text{}\end{array}$$

The coeﬃcients ${b}_{i}$, ${c}_{i}$, and ${a}_{ij}$ can be written in the standard Butcher notation.

In our multirate scheme, we use two diﬀerent sets of coeﬃcients. One set of coeﬃcients determines the RK4 scheme used for integrating the spacetime variables (33), the other set determines the RK2 scheme for the hydro variables (34). The coeﬃcients for the RK2 scheme are arranged such that RHS evaluations coincide with RK4 RHS evaluations. We list the corresponding multirate Butcher tableau in Table 1.

All the functions listed below return error codes in theory. However at this current point in time they always return 0 (success). Any failure to register or change a GF is assumed fatal and MoL will issue a level 0 warning stopping the code. This may change in future, in which case negative return values will indicate errors.

These are all aliased functions. You can get the functions directly through header ﬁles, but this feature may be phased out. Using function aliasing is the recommended method.

Tells MoL that the given GF is in the evolved category with the associated update GF.

Synopsis

C

CCTK_INT ierr = MoLRegisterEvolved(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

CCTK_INT ierr = MoLRegisterEvolvedSlow(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

CCTK_INT RHSIndex)

CCTK_INT ierr = MoLRegisterEvolvedSlow(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

Fortran

CCTK_INT ierr = MoLRegisterEvolved(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

CCTK_INT ierr = MoLRegisterEvolvedSlow(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

CCTK_INT RHSIndex)

CCTK_INT ierr = MoLRegisterEvolvedSlow(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

Result

Currently if there is an error, MoL will issue a level 0 warning. No sensible return codes exist.
0 success

Parameters

EvolvedIndex Index of the GF to be evolved.

RHSIndex Index of the associated update GF.

Discussion

Should be called in a function scheduled in MoL_Register. Use the Slow variant to register the slow sector of a multirate scheme.

See Also

CCTK_VarIndex() Get the variable index.

MoLRegisterSaveAndRestore() Register Save and Restore variables.

MoLRegisterConstrained() Register constrained variables.

MoLChangeToEvolved() Change a variable at runtime to be evolved.

MoLChangeToEvolvedSlow() Change a variable at runtime to be evolved in the slow sector.

Examples

C

ierr = MoLRegisterEvolved(CCTK_VarIndex("wavetoymol::phi"),

CCTK_VarIndex("wavetoymol::phirhs"));

CCTK_VarIndex("wavetoymol::phirhs"));

Fortran

call CCTK_VarIndex(EvolvedIndex, "wavetoymol::phi")

call CCTK_VarIndex(RHSIndex, "wavetoymol::phirhs")

ierr = MoLRegisterEvolved(EvolvedIndex, RHSIndex)

call CCTK_VarIndex(RHSIndex, "wavetoymol::phirhs")

ierr = MoLRegisterEvolved(EvolvedIndex, RHSIndex)

Tells MoL that the given GF is in the constrained category.

Synopsis

C

CCTK_INT ierr = MoLRegisterConstrained(CCTK_INT ConstrainedIndex)

Fortran

CCTK_INT ierr = MoLRegisterConstrained(CCTK_INT ConstrainedIndex)

Result

Currently if there is an error, MoL will issue a level 0 warning. No sensible return codes exist.
0 success

Parameters

ConstrainedIndex Index of the constrained GF.

Discussion

Should be called in a function scheduled in MoL_Register.

See Also

CCTK_VarIndex() Get the variable index.

MoLRegisterEvolved() Register evolved variables.

MoLRegisterSaveAndRestore() Register Save and Restore variables.

MoLChangeToConstrained() Change a variable at runtime to be constrained.

Examples

C

ierr = MoLRegisterConstrained(CCTK_VarIndex("ADMBase::alp"));

Fortran

call CCTK_VarIndex(ConstrainedIndex, "ADMBase::alp")

ierr = MoLRegisterConstrained(ConstrainedIndex)

ierr = MoLRegisterConstrained(ConstrainedIndex)

Tells MoL that the given GF is in the Save and Restore category.

Synopsis

C

CCTK_INT ierr = MoLRegisterSaveAndRestore(CCTK_INT SandRIndex)

Fortran

CCTK_INT ierr = MoLRegisterSaveAndRestore(CCTK_INT SandRIndex)

Result

Currently if there is an error, MoL will issue a level 0 warning. No sensible return codes exist.
0 success

Parameters

SandRIndex Index of the Save and Restore GF.

Discussion

Should be called in a function scheduled in MoL_Register.

See Also

CCTK_VarIndex() Get the variable index.

MoLRegisterEvolved() Register evolved variables.

MoLRegisterConstrained() Register constrained variables.

MoLChangeToSaveAndRestore() Change a variable at runtime to be Save and Restore.

Examples

C

ierr = MoLRegisterSaveAndRestore(CCTK_VarIndex("ADMBase::alp"));

Fortran

call CCTK_VarIndex(SandRIndex, "ADMBase::alp")

ierr = MoLRegisterSaveAndRestore(SandRIndex)

ierr = MoLRegisterSaveAndRestore(SandRIndex)

Tells MoL that the given group is in the evolved category with the associated update group.

Synopsis

C

CCTK_INT ierr = MoLRegisterEvolvedGroup(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

CCTK_INT ierr = MoLRegisterEvolvedGroupSlow(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

CCTK_INT RHSIndex)

CCTK_INT ierr = MoLRegisterEvolvedGroupSlow(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

Fortran

CCTK_INT ierr = MoLRegisterEvolvedGroup(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

CCTK_INT ierr = MoLRegisterEvolvedGroupSlow(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

CCTK_INT RHSIndex)

CCTK_INT ierr = MoLRegisterEvolvedGroupSlow(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

Result

Parameters

EvolvedIndex Index of the group to be evolved.

RHSIndex Index of the associated update group.

Discussion

Should be called in a function scheduled in MoL_Register. Use the Slow variant to register the slow sector of a multirate scheme.

See Also

CCTK_GroupIndex() Get the group index.

MoLRegisterSaveAndRestoreGroup() Register Save and Restore variables.

MoLRegisterConstrainedGroup() Register constrained variables.

Examples

C

ierr = MoLRegisterEvolvedGroup(CCTK_GroupIndex("wavetoymol::scalarevolvemol"),

CCTK_GroupIndex("wavetoymol::scalarevolvemolrhs"));

CCTK_GroupIndex("wavetoymol::scalarevolvemolrhs"));

Fortran

call CCTK_GroupIndex(EvolvedIndex, "wavetoymol::scalarevolvemol")

call CCTK_GroupIndex(RHSIndex, "wavetoymol::scalarevolvemolrhs")

ierr = MoLRegisterEvolvedGroup(EvolvedIndex, RHSIndex)

call CCTK_GroupIndex(RHSIndex, "wavetoymol::scalarevolvemolrhs")

ierr = MoLRegisterEvolvedGroup(EvolvedIndex, RHSIndex)

Tells MoL that the given group is in the constrained category.

Synopsis

C

CCTK_INT ierr = MoLRegisterConstrainedGroup(CCTK_INT ConstrainedIndex)

Fortran

CCTK_INT ierr = MoLRegisterConstrainedGroup(CCTK_INT ConstrainedIndex)

Result

Parameters

ConstrainedIndex Index of the constrained group.

Discussion

Should be called in a function scheduled in MoL_Register.

See Also

CCTK_GroupIndex() Get the group index.

MoLRegisterEvolvedGroup() Register evolved variables.

MoLRegisterSaveAndRestoreGroup() Register Save and Restore variables.

MoLChangeToConstrained() Change a variable at runtime to be constrained.

Examples

C

ierr = MoLRegisterConstrainedGroup(CCTK_GroupIndex("ADMBase::lapse"));

Fortran

call CCTK_GroupIndex(ConstrainedIndex, "ADMBase::lapse")

ierr = MoLRegisterConstrainedGroup(ConstrainedIndex)

ierr = MoLRegisterConstrainedGroup(ConstrainedIndex)

Tells MoL that the given group is in the Save and Restore category.

Synopsis

C

CCTK_INT ierr = MoLRegisterSaveAndRestoreGroup(CCTK_INT SandRIndex)

Fortran

CCTK_INT ierr = MoLRegisterSaveAndRestoreGroup(CCTK_INT SandRIndex)

Result

Parameters

SandRIndex Index of the save and restore group.

Discussion

Should be called in a function scheduled in MoL_Register.

See Also

CCTK_GroupIndex() Get the group index.

MoLRegisterEvolvedGroup() Register evolved variables.

MoLRegisterConstrainedGroup() Register constrained variables.

Examples

C

ierr = MoLRegisterSaveAndRestoreGroup(CCTK_GroupIndex("ADMBase::shift"));

Fortran

call CCTK_GroupIndex(SandRIndex, "ADMBase::shift")

ierr = MoLRegisterSaveAndRestoreGroup(SandRIndex)

ierr = MoLRegisterSaveAndRestoreGroup(SandRIndex)

Sets a GF to belong to the evolved category, with the associated update GF. Not used for the initial setting.

Synopsis

C

CCTK_INT ierr = MoLChangeToEvolved(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

CCTK_INT ierr = MoLChangeToEvolvedSlow(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

CCTK_INT RHSIndex)

CCTK_INT ierr = MoLChangeToEvolvedSlow(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

Fortran

CCTK_INT ierr = MoLChangeToEvolvedSlow(CCTK_INT EvolvedIndex,

CCTK_INT RHSIndex)

CCTK_INT RHSIndex)

Result

Parameters

EvolvedIndex Index of the evolved GF.

RHSIndex Index of the associated update GF.

Discussion

Should be called in a function scheduled in MoL_PreStep. Note that this function was designed to allow mixed slicings for thorn ADMBase. This set of functions is largely untested and should be used with great care.

See Also

CCTK_VarIndex() Get the variable index.

MoLRegisterEvolved() Register evolved variables.

MoLChangeToSaveAndRestore() Change a variable at runtime to be Save and Restore.

MoLChangeToConstrained() Change a variable at runtime to be constrained.

Examples

C

ierr = MoLChangeToEvolved(CCTK_VarIndex("ADMBase::alp"),

CCTK_VarIndex("adm_bssn::adm_bs_salp"));

CCTK_VarIndex("adm_bssn::adm_bs_salp"));

Fortran

call CCTK_VarIndex(EvolvedIndex, "ADMBase::alp")

call CCTK_VarIndex(RHSIndex,"adm_bssn::adm_bs_salp")

ierr = MoLChangeToEvolved(EvolvedIndex, RHSIndex)

call CCTK_VarIndex(RHSIndex,"adm_bssn::adm_bs_salp")

ierr = MoLChangeToEvolved(EvolvedIndex, RHSIndex)

Sets a GF to belong to the constrained category. Not used for the initial setting.

Synopsis

C

CCTK_INT ierr = MoLChangeToConstrained(CCTK_INT EvolvedIndex)

Fortran

CCTK_INT ierr = MoLChangeToConstrained(CCTK_INT EvolvedIndex)

Result

Parameters

ConstrainedIndex Index of the constrained GF.

Discussion

Should be called in a function scheduled in MoL_PreStep. Note that this function was designed to allow mixed slicings for thorn ADMBase. This set of functions is largely untested and should be used with great care.

See Also

CCTK_VarIndex() Get the variable index.

MoLRegisterConstrained() Register constrained variables.

MoLChangeToSaveAndRestore() Change a variable at runtime to be Save and Restore.

MoLChangeToEvolved() Change a variable at runtime to be evolved.

Examples

C

ierr = MoLChangeToConstrained(CCTK_VarIndex("ADMBase::alp"));

Fortran

call CCTK_VarIndex(EvolvedIndex, "ADMBase::alp")

ierr = MoLChangeToConstrained(EvolvedIndex)

ierr = MoLChangeToConstrained(EvolvedIndex)

Sets a GF to belong to the Save and Restore category. Not used for the initial setting.

Synopsis

C

CCTK_INT ierr = MoLChangeToSaveAndRestore(CCTK_INT SandRIndex)

Fortran

CCTK_INT ierr = MoLChangeToSaveAndRestore(CCTK_INT SandRIndex)

Result

Parameters

SandRIndex Index of the Save and Restore GF.

Discussion

Should be called in a function scheduled in MoL_PreStep. Note that this function was designed to allow mixed slicings for thorn ADMBase. This set of functions is largely untested and should be used with great care.

See Also

CCTK_VarIndex() Get the variable index.

MoLRegisterSaveAndRestore() Register Save and Restore variables.

MoLChangeToEvolved() Change a variable at runtime to be evolved.

MoLChangeToConstrained() Change a variable at runtime to be constrained.

Examples

C

ierr = MoLChangeToSaveAndRestore(CCTK_VarIndex("ADMBase::alp"));

Fortran

call CCTK_VarIndex(SandRIndex, "ADMBase::alp")

ierr = MoLChangeToSaveAndRestore(SandRIndex)

ierr = MoLChangeToSaveAndRestore(SandRIndex)

Sets a GF to belong to the “unknown” category. Not used for the initial setting.

Synopsis

C

CCTK_INT ierr = MoLChangeToNone(CCTK_INT RemoveIndex)

Fortran

CCTK_INT ierr = MoLChangeToNone(CCTK_INT RemoveIndex)

Result

Parameters

RemoveIndex Index of the “unknown” GF.

Discussion

See Also

CCTK_VarIndex() Get the variable index.

MoLChangeToEvolved() Change a variable at runtime to be evolved.

MoLChangeToSaveAndRestore() Change a variable at runtime to be Save and Restore.

MoLChangeToConstrained() Change a variable at runtime to be constrained.

Examples

C

ierr = MoLChangeToNone(CCTK_VarIndex("ADMBase::alp"));

Fortran

call CCTK_VarIndex(RemoveIndex, "ADMBase::alp")

ierr = MoLChangeToNone(RemoveIndex)

ierr = MoLChangeToNone(RemoveIndex)

Queries MoL for the index of the update variable for given GF in the evolved category.

Synopsis

C

CCTK_INT RHSindex = MoLQueryEvolvedRHS(CCTK_INT EvolvedIndex)

Fortran

CCTK_INT RHSindex = MoLQueryEvolvedRHS(CCTK_INT EvolvedIndex)

Result

If the grid function passed does not exists, MoL will issue a level 0 warning. If the grid function is not of an evolved type (fast
or slow sector) $-1$
will be returned. Otherwise the variable index of the update GF is returned.
$>0$
variable index of update GF

Parameters

EvolvedIndex Index of the GF whose update GF is to be returned.

Discussion

Both slow and fast evolved variables can be queried.

See Also

CCTK_VarIndex() Get the variable index.

MoLRegisterEvolved() Register evolved variables.

MoLChangeToEvolvedSlow() Change a variable at runtime to be evolved in the slow sector.

Examples

C

rhsindex = MoLQueryEvolvedRHS(CCTK_VarIndex("wavetoymol::phi"));

Fortran

call CCTK_VarIndex(EvolvedIndex, "wavetoymol::phi")

rhsindex = MoLQueryEvolvedRHS(EvolvedIndex)

rhsindex = MoLQueryEvolvedRHS(EvolvedIndex)

[1] J. Thornburg. Numerical Relativity in Black Hole Spacetimes. Unpublished thesis, University of British Columbia. 1993. Available from http://www.aei.mpg.de/~jthorn/phd/html/phd.html.

[2] J. Thornburg. A 3+1 Computational Scheme for Dynamic Spherically Symmetric Black Hole Spacetimes – II: Time Evolution. Preprint gr-qc/9906022, submitted to Phys. Rev. D.

[3] C. Shu. High Order ENO and WENO Schemes for Computational Fluid Dynamics. In T. J. Barth and H. Deconinck, editors High-Order Methods for Computational Physics. Springer, 1999. A related online version can be found under Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws at http://www.icase.edu/library/reports/rdp/97/97-65RDP.tex.refer.html.

[4] D. W. Neilsen and M. W. Choptuik. Ultrarelativistic ﬂuid dynamics. Class. Quantum Grav., 17: 733–759, 2000.

[5] M. Schlegel, O. Knoth, M. Arnold, and R. Wolke Journal of Computational and Applied Mathematics, 226, 345, 2009.

[6] E. Constantinescu and A. Sandu SIAM J. Sci. Comput., 33, 239, 2007.

ab_initially_reduce_order | Scope: private | BOOLEAN |

Description: Reduce order of accuracy initially so that no past timelevels of initial data are
required
| ||

Default: yes | ||

ab_type | Scope: private | KEYWORD |

Description: If using the the AB method, which sort
| ||

Range | Default: 1 | |

1 | same as forward Euler
| |

2 | second order
| |

3 | third order
| |

4 | fourth order
| |

5 | ﬁfth order
| |

adaptive_stepsize | Scope: private | BOOLEAN |

Description: Choose the time step size adaptively
| ||

Default: no | ||

copy_id_after_mol_poststep | Scope: private | BOOLEAN |

Description: if initial_data_is_crap is true, *when* should we copy the current time level to all
previous time levels: false ==> copy *before* MoL_PostStep (default, matches old behavior) true
==> copy *after* MoL_PostStep (maybe preferable for new code)
| ||

Default: no | ||

disable_prolongation | Scope: private | BOOLEAN |

Description: If Mesh reﬁnement is enabled should we use buﬀer zones in intermediate steps?
| ||

Default: yes | ||

generic_method_descriptor | Scope: private | STRING |

Description: A string used to create a table containing the description of the generic method
| ||

Range | Default: GenericIntermediateSteps = 2 GenericAlphaCoeﬀs = { 1.0 0.0 0.5 0.5 } GenericBetaCoeﬀs = { 1.0 0.5 } | |

.* | Should contain the Alpha and Beta arrays, and the
number of intermediate steps
| |

generic_type | Scope: private | KEYWORD |

Description: If using the generic method, which sort
| ||

Range | Default: RK | |

RK | One of the standard TVD Runge-Kutta methods
| |

ICN | Iterative Crank Nicholson as a generic method
| |

Table | Given from the generic method descriptor parameter
| |

Classic RK3 | Eﬃcient RK3 - classical version
| |

icn_avg_swapped | Scope: private | BOOLEAN |

Description: Use swapped averages in ICN method?
| ||

Default: no | ||

icn_avg_theta | Scope: private | REAL |

Description: theta of averaged ICN method, usually 0.5
| ||

Range | Default: 0.5 | |

0:1 | 0 <= theta <= 1
| |

init_rhs_zero | Scope: private | BOOLEAN |

Description: Initialise the RHS to zero
| ||

Default: yes | ||

initial_data_is_crap | Scope: private | BOOLEAN |

Description: If the initial data routine fails to set up the previous time levels, copy the current
backwards
| ||

Default: no | ||

maximum_absolute_error | Scope: private | REAL |

Description: Maximum allowed absolute error for adaptive stepsize control
| ||

Range | Default: 1.0e-6 | |

0.0:*) | ||

maximum_decrease | Scope: private | REAL |

Description: Maximum stepsize decrease factor
| ||

Range | Default: 10.0 | |

(1.0:*) | should be larger than one
| |

maximum_increase | Scope: private | REAL |

Description: Maximum stepsize increase factor
| ||

Range | Default: 5.0 | |

(1.0:*) | should be larger than one
| |

maximum_relative_error | Scope: private | REAL |

Description: Maximum allowed relative error for adaptive stepsize control
| ||

Range | Default: 1.0e-6 | |

0.0:*) | ||

mol_intermediate_steps | Scope: private | INT |

Description: Number of intermediate steps taken by the ODE method
| ||

Range | Default: 3 | |

1:* | Anything greater than 1
| |

mol_memory_always_on | Scope: private | BOOLEAN |

Description: Do we keep the scratch arrays allocated all the time?
| ||

Default: yes | ||

mol_nan_check | Scope: private | BOOLEAN |

Description: Should the RHS GFs be checked for NaNs?
| ||

Default: no | ||

mol_tiny | Scope: private | REAL |

Description: Eﬀective local machine zero; required by generic solvers
| ||

Range | Default: 1.e-15 | |

0:* | Defaults to 1.e-15
| |

ode_method | Scope: private | KEYWORD |

Description: The ODE method use by MoL to do time integration
| ||

Range | Default: ICN | |

Generic | Generic Shu-Osher Runge-Kutta type
| |

ICN | Iterative Crank Nicholson | |

ICN-avg | Iterative Crank Nicholson with averaging
| |

Euler | Euler | |

RK2 | Eﬃcient RK2
| |

RK2-central | Central RK2 | |

RK3 | Eﬃcient RK3
| |

RK4 | Eﬃcient RK4
| |

RK45 | RK45 (Fehlberg) with error estimation
| |

RK45CK | RK45CK (Cash-Karp) with error estimation
| |

RK65 | RK65 with error estimation
| |

RK87 | RK87 with error estimation
| |

AB | Adams-Bashforth
| |

RK2-MR-2:1 | 2nd order 2:1 multirate RK scheme based on
RK2 due to Schlegel et al 2009. This requires
init_RHS_zero=’no’.
| |

RK4-MR-2:1 | 3rd order 2:1 multirate RK scheme based on
RK43 due to Schlegel et al 2009. This requires
init_RHS_zero=’no’.
| |

RK4-RK2 | RK4 as fast method and RK2 as slow method
| |

rhs_error_weight | Scope: private | REAL |

Description: Weight of the RHS in the relative error calculation
| ||

Range | Default: 1.0 | |

0.0:* | should be between zero and one
| |

run_mol_poststep_in_post_recover_variables | Scope: private | BOOLEAN |

Description: Schedule the PostStep parts after recovery so that symmetries are automatically
done correctly.
| ||

Default: yes | ||

safety_factor | Scope: private | REAL |

Description: Safety factor for stepsize control
| ||

Range | Default: 0.9 | |

(0.0:*) | should be less than one
| |

set_id_boundaries | Scope: private | BOOLEAN |

Description: Should boundaries be overwritten (via synchronization, prolongation, boundary
conditions) by MoL?
| ||

Default: yes | ||

skip_initial_copy | Scope: private | BOOLEAN |

Description: Skip initial copy from previous to current time level
| ||

Default: no | ||

verbose | Scope: private | KEYWORD |

Description: How verbose should MoL be?
| ||

Range | Default: normal | |

none | No output at all (not implemented)
| |

normal | Standard verbosity
| |

register | List the variables registered as well
| |

extreme | Everything you never wanted to know
| |

mol_max_evolved_array_size | Scope: restricted | INT |

Description: The maximum total size of any grid arrays to be evolved
| ||

Range | Default: (none) | |

0:* | Anything non negative. Accumulated by other thorns
| |

mol_num_arrayconstrained_vars | Scope: restricted | INT |

Description: The maximum number of array constrained variables with timelevels that MoL needs
to know about (DPRECATED)
| ||

Range | Default: (none) | |

0:* | Anything non negative. Added to by other thorns.
| |

mol_num_arrayevolved_vars | Scope: restricted | INT |

Description: The maximum number of array variables to be evolved by MoL (DPRECATED)
| ||

Range | Default: (none) | |

0:* | Anything non negative. Added to by other thorns.
| |

mol_num_arraysaveandrestore_vars | Scope: restricted | INT |

Description: The maximum number of array variables to be evolved outside of MoL but that MoL
needs to know about (DPRECATED)
| ||

Range | Default: (none) | |

0:* | Anything non negative. Added to by other thorns.
| |

mol_num_constrained_vars | Scope: restricted | INT |

Description: The maximum number of constrained variables with timelevels that MoL needs to
know about (DPRECATED)
| ||

Range | Default: (none) | |

0:* | Anything non negative. Added to by other thorns.
| |

mol_num_evolved_vars | Scope: restricted | INT |

Description: The maximum number of variables to be evolved by MoL (DPRECATED)
| ||

Range | Default: (none) | |

0:* | Anything non negative. Added to by other thorns.
| |

mol_num_evolved_vars_slow | Scope: restricted | INT |

Description: The maximum number of ’slow’ variables to be evolved by MoL (DPRECATED)
| ||

Range | Default: (none) | |

0:* | Anything non negative. Added to by other thorns.
| |

mol_num_saveandrestore_vars | Scope: restricted | INT |

Description: The maximum number of variables to be evolved outside of MoL but that MoL needs
to know about (DPRECATED)
| ||

Range | Default: (none) | |

0:* | Anything non negative. Added to by other thorns.
| |

mol_num_scratch_levels | Scope: restricted | INT |

Description: Number of scratch levels required by the ODE method
| ||

Range | Default: (none) | |

0:* | Anything non negative
| |

cctk_initial_time | Scope: shared from CACTUS | REAL |

Implements:

methodoﬂines

Group Names | Variable Names | Details | |

rkalphacoeﬃcients | compact | 0 | |

RKAlphaCoeﬃcients | dimensions | 2 | |

distribution | CONSTANT | ||

group type | ARRAY | ||

size | MOL_INTERMEDIATE_STEPS | ||

size | MOL_NUM_SCRATCH_LEVELS+1 | ||

tags | Checkpoint=”no” | ||

timelevels | 1 | ||

variable type | REAL | ||

rkbetacoeﬃcients | compact | 0 | |

RKBetaCoeﬃcients | dimensions | 1 | |

distribution | CONSTANT | ||

group type | ARRAY | ||

size | MOL_INTERMEDIATE_STEPS | ||

tags | Checkpoint=”no” | ||

timelevels | 1 | ||

variable type | REAL | ||

mol_counters | compact | 0 | |

MoL_Intermediate_Step | description | The counter for the time integration method | |

MoL_Stepsize_Bad | dimensions | 0 | |

MoL_SlowStep | distribution | CONSTANT | |

MoL_SlowPostStep | group type | SCALAR | |

tags | Checkpoint=”no” | ||

timelevels | 1 | ||

variable type | INT | ||

mol_original_time | compact | 0 | |

Original_Time | description | The original time and delta time which are reset by MoL during evolution | |

Original_Delta_Time | dimensions | 0 | |

distribution | CONSTANT | ||

group type | SCALAR | ||

tags | Checkpoint=”no” | ||

timelevels | 1 | ||

variable type | REAL | ||

scratchspaceslow | compact | 0 | |

ScratchSpaceSlow | description | The original time and delta time which are reset by MoL during evolution | |

dimensions | 3 | ||

distribution | DEFAULT | ||

group type | GF | ||

tags | Prolongation=”None” Checkpoint=”no” | ||

timelevels | 99 | ||

vararray_size | MoL_Num_Scratch_Levels | ||

variable type | REAL | ||

scratchspace | compact | 0 | |

ScratchSpace | description | The original time and delta time which are reset by MoL during evolution | |

dimensions | 3 | ||

distribution | DEFAULT | ||

group type | GF | ||

tags | Prolongation=”None” Checkpoint=”no” | ||

timelevels | 99 | ||

vararray_size | MoL_Num_Scratch_Levels | ||

variable type | REAL | ||

Group Names | Variable Names | Details | |

sandrscratchspace | compact | 0 | |

SandRScratchSpace | description | The original time and delta time which are reset by MoL during evolution | |

dimensions | 3 | ||

distribution | DEFAULT | ||

group type | GF | ||

tags | Prolongation=”None” Checkpoint=”no” | ||

timelevels | 99 | ||

variable type | REAL | ||

errorestimate | compact | 0 | |

ErrorEstimate | description | The original time and delta time which are reset by MoL during evolution | |

dimensions | 3 | ||

distribution | DEFAULT | ||

group type | GF | ||

tags | Prolongation=”None” Checkpoint=”no” | ||

timelevels | 99 | ||

variable type | REAL | ||

errorscalars | compact | 0 | |

Error | description | Global error estimate | |

Count | dimensions | 0 | |

EstimatedDt | distribution | CONSTANT | |

group type | SCALAR | ||

timelevels | 1 | ||

variable type | REAL | ||

Provides:

MoLRegisterEvolved to

MoLRegisterEvolvedSlow to

MoLRegisterConstrained to

MoLRegisterSaveAndRestore to

MoLRegisterEvolvedGroup to

MoLRegisterEvolvedGroupSlow to

MoLRegisterConstrainedGroup to

MoLRegisterSaveAndRestoreGroup to

MoLChangeToEvolved to

MoLChangeToEvolvedSlow to

MoLChangeToConstrained to

MoLChangeToSaveAndRestore to

MoLChangeToNone to

MoLQueryEvolvedRHS to

MoLNumIntegratorSubsteps to

This section lists all the variables which are assigned storage by thorn CactusNumerical/MoL. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

Always: | Conditional: |

MoL_Counters | RKAlphaCoeﬃcients |

MoL_Original_Time | RKBetaCoeﬃcients |

ErrorScalars | |

ErrorScalars | |

ErrorScalars | |

CCTK_STARTUP (conditional)

mol_startup

startup banner

Language: | c | |

Type: | function | |

CCTK_PARAMCHECK (conditional)

mol_paramcheck

basic parameter checking

Language: | c | |

Type: | function | |

CCTK_INITIAL (conditional)

mol_startloop

initialise the step size control

Language: | c | |

Options: | level | |

Type: | function | |

CCTK_EVOL (conditional)

mol_evolution

a single cactus evolution step using mol

Type: | group | |

While: | mol::mol_stepsize_bad | |

MoL_Evolution (conditional)

mol_startstep

mol internal setup for the evolution step

Type: | group | |

MoL_StartStep (conditional)

mol_setcounter

set the counter for the ode method to loop over

Language: | c | |

Options: | level | |

Type: | function | |

MoL_StartStep (conditional)

mol_settime

ensure the correct time and timestep are used

Language: | c | |

Options: | level | |

Type: | function | |

MoL_StartStep (conditional)

mol_allocatescratchspace

allocate storage for scratch levels

Language: | c | |

Options: | level | |

Type: | function | |

MoL_Evolution (conditional)

mol_prestep

physics thorns can schedule preloop setup routines in here

After: | mol_startstep | |

Before: | mol_step | |

Type: | group | |

MoL_Evolution (conditional)

mol_initialcopy

ensure the data is in the correct timelevel

After: | mol_prestep | |

Before: | mol_step | |

Language: | c | |

Type: | function | |

MoL_Evolution (conditional)

mol_step

the loop over the intermediate steps for the ode integrator

After: | mol_prestep | |

Type: | group | |

While: | mol::mol_intermediate_step | |

MoL_Step (conditional)

mol_icnaverage

averages the time levels for the averaging icn method

Before: | mol_calcrhs | |

Language: | c | |

Type: | function | |

CCTK_WRAGH (conditional)

mol_setupindexarrays

set up the mol bookkeeping index arrays

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_initrhs

initialise the rhs functions

Before: | mol_calcrhs | |

Language: | c | |

Type: | function | |

MoL_Step

mol_calcrhs

physics thorns schedule the calculation of the discrete spatial operator in here

Type: | group | |

MoL_Step

mol_postrhs

modify rhs functions

After: | mol_calcrhs | |

Before: | mol_nancheck | |

mol_add | ||

Type: | group | |

MoL_Step

mol_rhsboundaries

any ’ﬁnal’ modiﬁcations to the rhs functions (boundaries etc.)

After: | mol_postrhs | |

Before: | mol_nancheck | |

mol_add | ||

Type: | group | |

MoL_Step (conditional)

mol_nancheck

check the rhs gfs for nans

After: | mol_calcrhs | |

Before: | mol_add | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_genericrkadd

updates calculated with a generic method

After: | mol_calcrhs | |

Before: | mol_poststep | |

mol_poststepmodify | ||

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_euleradd

updates calculated with the euler method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_rk2add

updates calculated with the eﬃcient runge-kutta 2 method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_rk2centraladd

updates calculated with the central runge-kutta 2 method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_rk3add

updates calculated with the eﬃcient runge-kutta 3 method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

CCTK_WRAGH (conditional)

mol_setuprkcoeﬃcients

initialize the generic runge-kutta coeﬃcients

Language: | c | |

Options: | global | |

Storage: | rkalphacoeﬃcients | |

rkbetacoeﬃcients | ||

Type: | function | |

MoL_Step (conditional)

mol_rk4add

updates calculated with the eﬃcient runge-kutta 4 method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_rk45add

updates calculated with the runge-kutta 45 method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_rk65add

updates calculated with the runge-kutta 65 method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_rk87add

updates calculated with the runge-kutta 87 method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_icnadd

updates calculated with the eﬃcient icn method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_icnadd

updates calculated with the averaging icn method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_abadd

updates calculated with the adams-bashforth

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_rk2_mr_2_1_add

updates calculated with the multirate runge-kutta 2 method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_rk4_mr_2_1_add

updates calculated with the multirate runge-kutta 4 method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

MoL_Step (conditional)

mol_rk4_rk2_add

updates calculated with the multirate rk4/rk2 method

After: | mol_calcrhs | |

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

CCTK_WRAGH (conditional)

mol_setschedulestatus

set the ﬂag so it is ok to register with mol

After: | mol_setupindexarrays | |

Language: | c | |

Options: | global | |

Type: | function | |

MoL_Step (conditional)

mol_poststep

the group for physics thorns to schedule boundary calls etc.

After: | mol_add | |

Type: | group | |

CCTK_POSTINITIAL (conditional)

mol_poststep

ensure that everything is correct after the initial data have been set up

Type: | group | |

CCTK_POSTREGRID (conditional)

mol_poststep

ensure that everything is correct after regridding

Type: | group | |

CCTK_POSTRESTRICTINITIAL (conditional)

mol_poststep

ensure that everything is correct after restriction

Type: | group | |

CCTK_POSTRESTRICT (conditional)

mol_poststep

ensure that everything is correct after restriction

Type: | group | |

CCTK_POST_RECOVER_VARIABLES (conditional)

mol_poststep

ensure that everything is correct after recovery

Type: | group | |

MoL_Step (conditional)

mol_poststepmodify

the group for physics thorns to schedule enforcing constraints

After: | mol_add | |

Before: | mol_poststep | |

Type: | group | |

CCTK_POSTINITIAL (conditional)

mol_poststepmodify

the group for physics thorns to schedule enforcing constraints

Before: | mol_poststep | |

Type: | group | |

CCTK_POSTINITIAL (conditional)

mol_pseudoevolution

calculate pseudo-evolved quantities

After: | mol_poststep | |

Type: | group | |

CCTK_EVOL (conditional)

mol_pseudoevolution

calculate pseudo-evolved quantities

After: | mol_evolution | |

Type: | group | |

CCTK_WRAGH (conditional)

mol_register

the group where physics thorns register variables with mol

After: | mol_setschedulestatus | |

Type: | group | |

CCTK_POSTREGRIDINITIAL (conditional)

mol_pseudoevolutionboundaries

apply boundary conditions to pseudo-evolved quantities

After: | mol_poststep | |

Type: | group | |

CCTK_POSTREGRID (conditional)

mol_pseudoevolutionboundaries

apply boundary conditions to pseudo-evolved quantities

After: | mol_poststep | |

Type: | group | |

CCTK_POSTRESTRICTINITIAL (conditional)

mol_pseudoevolutionboundaries

apply boundary conditions to pseudo-evolved quantities

After: | mol_poststep | |

Type: | group | |

CCTK_POSTRESTRICT (conditional)

mol_pseudoevolutionboundaries

apply boundary conditions to pseudo-evolved quantities

After: | mol_poststep | |

Type: | group | |

MoL_Step (conditional)

mol_decrementcounter

alter the counter number

After: | mol_add | |

Before: | mol_poststep | |

mol_poststepmodify | ||

Language: | c | |

Options: | level | |

Type: | function | |

MoL_Step (conditional)

mol_resettime

if necessary, change the time

After: | mol_decrementcounter | |

Before: | mol_poststep | |

mol_poststepmodify | ||

Language: | c | |

Options: | level | |

Type: | function | |

MoL_Step (conditional)

mol_resetdeltatime

if necessary, change the timestep

After: | mol_poststep | |

mol_poststepmodify | ||

Language: | c | |

Options: | level | |

Type: | function | |

MoL_Evolution (conditional)

mol_restoresandr

restoring the save and restore variables to the original state

After: | mol_reduceadaptiveerror | |

mol_ﬁnishloop | ||

Language: | c | |

Type: | function | |

MoL_Evolution (conditional)

mol_freescratchspace

free storage for scratch levels

After: | mol_restoresandr | |

Language: | c | |

Options: | level | |

Type: | function | |

MoL_Evolution (conditional)

mol_initadaptiveerror

control the step size: initialize error check variables

After: | mol_poststep | |

Language: | c | |

Options: | level | |

Type: | function | |

CCTK_WRAGH (conditional)

mol_reportnumbervariables

report how many of each type of variable there are

After: | mol_register | |

Language: | c | |

Options: | meta | |

Type: | function | |

MoL_Evolution (conditional)

mol_ﬁndadaptiveerror

control the step size: compute error check variables

After: | mol_initadaptiveerror | |

Language: | c | |

Type: | function | |

MoL_Evolution (conditional)

mol_reduceadaptiveerror

control the step size: reduce error check variables

After: | mol_ﬁndadaptiveerror | |

Language: | c | |

Options: | level | |

Type: | function | |

CCTK_POSTSTEP (conditional)

mol_setestimateddt

control the step size: set the new timestep

Language: | c | |

Options: | level | |

Type: | function | |

MoL_Evolution (conditional)

mol_ﬁnishloop

control the step size

After: | mol_step | |

Language: | c | |

Options: | level | |

Type: | function | |

CCTK_TERMINATE

mol_freeindexarrays

free the mol bookkeeping index arrays

Before: | driver_terminate | |

Language: | c | |

Type: | function | |

CCTK_POSTINITIAL (conditional)

mol_ﬁllalllevels

a bad routine. ﬁlls all previous timelevels with data copied from the current.

After: | mol_poststep | |

Language: | c | |

Type: | function | |

CCTK_POSTINITIAL (conditional)

mol_ﬁllalllevels

a bad routine. ﬁlls all previous timelevels with data copied from the current.

Before: | mol_poststep | |

Language: | c | |

Type: | function | |

CCTK_EVOL (conditional)

mol_startloop

initialise the step size control

Before: | mol_evolution | |

Language: | c | |

Options: | level | |

Type: | function | |

Alias Name: | Function Name: |

MoL_ABAdd | MoL_Add |

MoL_EulerAdd | MoL_Add |

MoL_GenericRKAdd | MoL_Add |

MoL_ICNAdd | MoL_Add |

MoL_ICNAverage | MoL_Prepare |

MoL_RK2Add | MoL_Add |

MoL_RK2CentralAdd | MoL_Add |

MoL_RK2_MR_2_1_Add | MoL_Add |

MoL_RK3Add | MoL_Add |

MoL_RK45Add | MoL_Add |

MoL_RK4Add | MoL_Add |

MoL_RK4_MR_2_1_Add | MoL_Add |

MoL_RK4_RK2_Add | MoL_Add |

MoL_RK65Add | MoL_Add |

MoL_RK87Add | MoL_Add |