.. highlight:: rest Setting up a problem ==================== First review the note on :doc:`units`. 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 :doc:`tests` 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 :doc:`parameters`), 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: * A description of the background gas state (:ref:`ssec-init-gas`) * A description of how cosmic rays propagate (:ref:`ssec-init-propagation`) * The initial set of packets present at the start of the simulation (:ref:`ssec-init-packets`) * The initial set of sources present at the start of the simulation (:ref:`ssec-init-sources`) * Any problem-specific work to be performed during the calculation (:ref:`ssec-problem-work`) 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`` (:ref:`ssec-init-gas`) and ``initPackets`` (:ref:`ssec-init-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 :ref:`ssec-definitions-H`. .. _ssec-init-gas: 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 :doc:`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 (:ref:`ssec-gasdata`) 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`` .. _ssec-init-propagation: 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 (:ref:`ssec-init-gas`), for many realistic applications users will probably want to implement their own propagation models, and details of how to do so are provided in :doc:`propagation`. 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: .. math:: 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}}, where :math:`K_\parallel` is the spatial diffusion coefficient parallel to the magnetic field, :math:`K_\perp` is the spatial diffusion coefficient perpendicular to the magnetic field, :math:`w` is the streaming speed, :math:`K_{\mu\mu}` is the pitch angle diffusion coefficient, and :math:`p` is the cosmic ray momentum. The coefficients subscripted by 0 and the indices :math:`\alpha` are determined by parameters read from the input deck as follows: * ``cr.kPar0``: the value of :math:`K_{\parallel,0}` * ``cr.kParIdx``: the value of :math:`\alpha_{\parallel}` * ``cr.kPerp0``: the value of :math:`K_{\perp,0}` * ``cr.kPerpIdx``: the value of :math:`\alpha_{\perp}` * ``cr.kMu0``: the value of :math:`K_{\mu\mu,0}` (only used if pitch angle evolution is turned on; see :ref:`sec-pitch-angle`) * ``cr.kMuIdx``: the value of :math:`\alpha_{\mu\mu}` (only used if pitch angle evolution is turned on; see :ref:`sec-pitch-angle`) * ``cr.vStr0``: a value that sets :math:`w_0`; however, the way that it is interpreted depends on other parameters (see below) * ``cr.vStrIdx``: the value of :math:`\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 :math:`w_0 = 3.5 v_{A,i}`, where :math:`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 :math:`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., :math:`\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. .. _ssec-init-packets: 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& x, std::vector& 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 :math:`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. .. _ssec-init-sources: 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& x, std::vector& sources); Usage is identical to ``initPackets`` (:ref:`ssec-init-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 :math:`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 :math:`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 :math:`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\ :sup:`-1` (or m s\ :sup:`-1` if the code was compiled in SI mode -- see :ref:`ssec-definitions-H`); defaults to zero * ``a``: the acceleration of the source, in units of cm s\ :sup:`-2` (or m s\ :sup:`-2` if the code was compiled in SI mode -- see :ref:`ssec-definitions-H`); defaults to zero .. _ssec-problem-work: 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. .. _ssec-definitions-H: 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.