Setting up a problem

First review the note on Unit systems in criptic.

To set up a problem in criptic, the first step is to create a sub-directory of Src/Prob that will hold your problem setup. Several examples of problem setups are provided in the repository, for example Prob/ThickTarget, which implements the thick target example problem described in the criptic method paper. All the tests described in The criptic test suite provide examples of problem setups as well. It is probably useful to refer to these examples when looking through this guide.

To set up a problem, the user must provide a file Prob.cpp in the problem directory that contains an implementation of the initProb function in the criptic namespace. The signature for this function is:

Prob *initProb(const criptic::ParmParser& pp,
               const criptic::Geometry& geom,
               criptic::RngThread& rng)

The routine takes as arguments the contents of the parameter file (see Parameter files), which are stored in the ParmParser object, the problem geometry, which is stored in the Geometry object, and the random number generator, which is stored in the RngThread object.

The purpose of initProb is to return a point to an object of class Prob, which provides routines to specify:

The Src/Prob directory contains a default implementation of the Prob class that sets defaults for each of the above steps, but in general users will want to create their own derived classes that will override the defaults for at least some of these in order to set up the particlar problem they wish to simulate. Examples are provided in all of the Prob/Test sub-directories. For example, a user who wishes to simulate a problem with no sources present and using the default CR propagation model, but wishes to specify the background gas configuration and the initial CR packet distribution, could define a derived class that overrides the initGas (Initializing the background gas state) and initPackets (Initializing cosmic ray packets) methods of Prob, and return that from the initProb function. Full details on the methods used to control various aspects of problem setup are provided below. The parent Prob class stores protected references to the ParmParser, Geometry, and RngThread objects, so these are accessible to the child classes.

In addition, the problem directory can optionally contain a Definitions.H file that overrides some code defaults. See Compile-time definitions.

Initializing the background gas state

In order to specify the background gas state (where “gas” here includes not just the gas itself, but also the magnetic and radiation fields), the Prob class provides the initGas routine. The call signature for this routine is:

gas::Gas* Prob::initGas();

The routine must return a pointer to an object of class Gas. In most realistic calculations users probably wish to implement their own descriptions of the background gas, and instructions for doing so are provided in Describing the background gas. However, some standard options are provided.

The default implementation of initGas, which is used if the user does not override this default in the child class, returns on object of type UniformGas, which describes a uniform gas distribution whose properties (e.g., density, magnetic field) are set through the input deck. The precise implementation is:

virtual gas::Gas *Prob::initGas() {
  return new gas::UniformGas(pp);
}

For a UniformGas object, the following parameters are read from the input deck (all a required unless a default value is specified):

  • gas.denstiy: total mass density of the gas

  • gas.ionDensity: mass density of ionized component of the gas

  • gas.magField: the magnetic field (set as a three-component vector)

  • gas.velocity: gas velocity (set as a three-component vector)

  • gas.dx: a characteristic length scale for the gas; defaults to 0. This parameter exists because for some types of CR propagation model the code requires a characteristic length scale to estimate how frequently it is required to re-calculate things like the CR pressure, and a uniform gas otherwise has not characteristic size scales; it has no other effect.

  • gas.comp: a shortcut to set the gas composition; setting this to the keywords atomic, molecular, or ionized will set a generic Milky Way-like atomic, molecular, or ionized chemical composition; defaults to atomic if not otherwise set.

  • gas.xH0, gas.xHp, gas.xHe0, gas.xHep, gas.xHep2, gas.xe, gas.Z: as an alternative to setting gas.comp, users can manually specify the abundnaces of various components – see (The GasData class) for definitions of these abundances; these parameters are optional, and if specified override the abundances set by gas.comp

  • gas.TBB: a vector of values describing the temperatures of any dilute blackbody radiation fields present

  • gas.WBB: a vector of values describing the dilution factors of any dilute blackbody radiation fields present; must have the same number of entries as gas.TBB

Initializing cosmic ray propagation

In order to specify the cosmic ray propagation model, which gives the diffusion coefficients and streaming velocity as a function of cosmic ray properties, gas properties, etc., the Prob class provides the method:

propagation::Propagation *Prob::initProp();

The function must return a pointer to an object of class Propagation. As with the initGas routine (Initializing the background gas state), for many realistic applications users will probably want to implement their own propagation models, and details of how to do so are provided in Describing the cosmic ray propagation model. However, some pre-defined options are also provided.

The default implementation of this method, provided in the base Prob class and used if not overridden in a derived class, returns a pointer to an instance of PropPowerlaw, a class that described a powerlaw propagation model. The precise implementation is:

virtual propagation::Propagation *initProp() {
   return new propagation::ProbPowerlaw(pp);
}

For the default PropPowerlaw model, the propagation coefficients take the following powerlaw forms:

\[\begin{split}K_\parallel & = & K_{\parallel,0} (p / m_p c)^{\alpha_\parallel} \\ K_\perp & = & K_{\perp,0} (p / m_p c)^{\alpha_\perp} \\ w & = & w_0 (p / m_p c)^{\alpha_w} \\ K_{\mu\mu} & = & K_{\mu\mu,0} (p / m_p c)^{\alpha_{\mu\mu}},\end{split}\]

where \(K_\parallel\) is the spatial diffusion coefficient parallel to the magnetic field, \(K_\perp\) is the spatial diffusion coefficient perpendicular to the magnetic field, \(w\) is the streaming speed, \(K_{\mu\mu}\) is the pitch angle diffusion coefficient, and \(p\) is the cosmic ray momentum. The coefficients subscripted by 0 and the indices \(\alpha\) are determined by parameters read from the input deck as follows:

  • cr.kPar0: the value of \(K_{\parallel,0}\)

  • cr.kParIdx: the value of \(\alpha_{\parallel}\)

  • cr.kPerp0: the value of \(K_{\perp,0}\)

  • cr.kPerpIdx: the value of \(\alpha_{\perp}\)

  • cr.kMu0: the value of \(K_{\mu\mu,0}\) (only used if pitch angle evolution is turned on; see Simulating pitch angle evolution)

  • cr.kMuIdx: the value of \(\alpha_{\mu\mu}\) (only used if pitch angle evolution is turned on; see Simulating pitch angle evolution)

  • cr.vStr0: a value that sets \(w_0\); however, the way that it is interpreted depends on other parameters (see below)

  • cr.vStrIdx: the value of \(\alpha_\mathrm{str}\)

  • cr.vAStream: if this is set to a non-zero integer value, then cr.vStr0 is assumed to be in units of the local ion Alfvén speed, i.e., cr.vStr0  3.5 will be interpreted to mean that \(w_0 = 3.5 v_{A,i}\), where \(v_{A,i}\) is the ion Alfvén speed; if this is set to zero, the value of cr.vStr0 is assumed to give the absolute value of \(w_0\) in physical units (cm/s or m/s depending on whether the code is using CGS or MKS units); defaults to 1

  • cr.varStreamDir: if this is set to a non-zero integer value, cosmic rays are assumed to stream in the direction opposite the local CR pressure gradient, i.e., \(\mathrm{sgn}(w) = -\mathrm{sgn}(\nabla P_\mathrm{CR} \cdot \mathbf{B})\); otherwise CRs stream in a direction specified by the sign of cr.vStr0, with positive values corresponding to streaming parallel to the magnetic field and negative values to streaming anti-parallel; defaults to 1

All parameters are required unless a default value is specified.

Initializing cosmic ray packets

A user can specify a population of CR packets that is present at the start of a simulation using the initPackets method, which has the call signature:

void Prob::initPackets(std::vector<RealVec>& x,
                       std::vector<CRPacket>& packets);

The two arguments are two initially-empty vectors, x and packets, into which cosmic rays can be placed: x is a vector of cosmic ray positions (of type RealVec), and packets is a vector of cosmic ray properties (of type CRPacket. Users can alter these objects by inserting into them the positions and properties of any cosmic ray packets that are present at the start of the simulation. The default implementation of initPackets in the Prob class, which is used unless overridden by a user-specific problem setup class, simply leaves these vectors empty so that no packets are present initially.

For each cosmic ray packet added in this routine, the user must set seven fields of the CRPacket class:

  • type: the type of particle the packet represents; allowed values are partTypes::proton, partTypes::electron, and partTypes::positron; see partTypes for details.

  • src: the unique ID of the cosmic ray source responsible for producing this packet; should generally be set to nullIdx for packets present in the initial conditions, since these do not come from any source

  • tInj: time at which the packet was injected; should generally be set to 0 for packets present in the initial conditions

  • p: the momentum of the packet, in code units where momenta are measured in units of \(m_p c\). Users who wish to describe packets in terms of kinetic energy instead of momentum may find it convenient to use the set_T method provided as part of the CRPacket class.

  • w: the weight of the packet, meaning the number of individual cosmic ray particles that it represents

  • wInj: the weight of the packet at the time it was injected; should generally be set equal to w for packets present in the initial conditions

  • gr: the grammage the packet has traversed; should generally be set to 0 for packets present in the initial conditions

A final caution is that, in an MPI calculation, the initPackets function will be called by each MPI rank, which can result in multiple copies of each initial packet being injected. To prevent this, the packet creation code can be wrapped in an if block of the form:

if (MPIUtil::IOProc) {

   (... code to create packets ...)

}

This will prevent the creation code from being executed on more than one MPI rank. See MPIUtil for details.

Initializing cosmic ray sources

A user can specify a population of CR sources that is present at the start of a simulation using the initSources method, which has the call signature:

void Prob::initSources(std::vector<RealVec>& x,
                       std::vector<CRSource>& sources);

Usage is identical to initPackets (Initializing cosmic ray packets), except that the routine provides a vector of CRSource rather than CRPacket. As with initPackets, the default implementation is that no sources are added. If the user overrides this behavior and adds initial sources, for each such source hte user must initialize some of the fields of the CRSource class. These are:

  • type: the type of particle the source injects; allowed values are partTypes::proton, partTypes::electron, and partTypes::positron; see partTypes for details. Note that each source can inject only a single type of particle; if a given source should inject multiple tpes of particle, represent this as two sources whose positions and properties are identical except for type.

  • p0: minimum possible momentum that this source injects, in code units where momenta are measured in units of \(m_p c\). Users who wish to describe sources in terms of kinetic energy instead of momentum may find it convenient to use the set_T0 method provided as part of the CRSource class.

  • p1: same as p0, but giving the upper limit to the allowed momentum. Note that p1 can be exactly equal to p0, in which case the source is treated as injecting monoenergetic cosmic rays with momentum exactly equal to p0.

  • q: index on the powerlaw momentum distribution for the source; cosmic rays injected by the source have a momentum distribution \(dn/dp \propto p^q\).

  • k: the rate coefficient for cosmic ray injection, with units of inverse momentum times time; the source injects cosmic rays at a rate per unit time per unit momentum given by \(dn/dp = k (p/m_p c)^q\). Instead of setting this field manually, it is often convenient to call the setLum method of CRSource, which sets k in order to produce a specified total cosmic ray luminosity integrated over all momenta, with units of energy per unit time.

  • v: the velocity of the source, in units of cm s-1 (or m s-1 if the code was compiled in SI mode – see Compile-time definitions); defaults to zero

  • a: the acceleration of the source, in units of cm s-2 (or m s-2 if the code was compiled in SI mode – see Compile-time definitions); defaults to zero

Problem-specific work

In addition to initializing and specifying the gas and cosmic ray propagation models, the Prob class provides routines to carry out arbitrary user-specified work at various points during a simulation. Users can perform work at the following times:

  • On problem initialization, via the userSetup method.

  • At the end of each time step, via the userWork method.

  • Immediately after reading a checkpoint, via the userRead method.

  • Immediately after writing a checkpoint, via the userWrite method.

  • At the end of the simulation before writing the final checkpoint, via the userFinalize method.

The default implementations of each of these methods provided in the Prob clas are empty, but users can override this behavior in their derived problem setup classes.

The call signatures are:

virtual void Prob::userSetup();
virtual void Prob::userWork(const Real t,
                            Real& dt,
                            gas::Gas& gasBG,
                            propagation::Propagation& prop,
                            CRTree& tree);
virtual void Prob::userRead(const std::string& filename,
                            const int step,
                            Real& t,
                            Real& dt,
                            CRTree& tree);
virtual void Prob::userWrite(const int step,
                             const Real t,
                             const Real dt,
                             const int chkNum,
                             const std::string& baseChkName,
                             const CRTree& tree);
virtual void Prob::userFinalize(const Real t,
                                gas::Gas& gasBG,
                                propagation::Propagation& prop,
                                CRTree& tree);

See Prob for full documentation of the arguments provided to each of these methods. Note that, with the exception of userSetup, these methods provide the user with full access to the CRTree class, which contains the lists of all packets and sources; users can access these using the getPacketPos, getPacketData, getSourcePos, and getSourceData methods of CRTree. This access in turn allows users to perform arbitary calculations on these data, and, for those methods where the tree is not marked const, to implement arbitrary changes to the properties of packets and sources. For an example of how this can be applied, the userWork function in the Prob/Test/LoopDiff test problem is used to update the positions, velocities, and accelerations of sources that move over time.

Compile-time definitions

In addition to the required Prob.cpp file, a problem setup directory can also contain an optional file Definitions.H that overrides some default compile-time definitions. Default values are set in Src/Definitions.H. These definitions specify, among other things, the unit system that criptic uses (CGS/Gaussian units versus SI units), whether to use single or double precision for real numbers, parameters that control performance of the kd-tree, etc. See Definitions.H for a full content listing. Users can add a #define statement for any of these quantities to the Definitions.H file in their problem directory, and these values will take precedence over the corresponding default values in Src/Definitions.H. However, for most simulations it is probably best simply to use the defaults.