Changes

Udo Ziegler · 4ccb63fc
--- a/3-NIRVANA-user-guide/3.2-User-interfaces.md
+++ b/3-NIRVANA-user-guide/3.2-User-interfaces.md
@@ -820,3 +820,1301 @@ for the type of coordinate system.
        *Note: In a multi-species gas the ionization fraction is not a
        freely selectable parameter but computed selfconsistently from
        the ionisation structure.*
+### Defining initial conditions
+The user interface `configUser.c` serves to define IC for a problem.
+Defining IC means assigning values to primary physical variables,
+{𝜚, **m**, *e*, **B**, *C*<sub>*c*</sub>, *n*<sub>*s*</sub>}, on the
+mesh. However, only problem-relevant variables have to be defined.
+Primary variables which are irrelevant for the problem need not to be
+assigned. For example, in simulations with isothermal EOS the total
+energy density *e* is irrelevant since no energy equation is solved.
+Assignment of the magnetic field in HD simulations (when `_C.mf`=0) is
+redundant and even not allowed because arrays for the magnetic field
+components are not allocated in the HD case. Likewise, arrays for
+*C*<sub>*c*</sub> and *n*<sub>*s*</sub> are not declared if
+`_C.tracer`=0 respective `_C.species`=0. Furthermore, derived variables
+like the temperature must not be assigned here. If the TESTFIELDS
+infrastructure is used testfield fluctuation variables,
+**b**<sub>*t*</sub>, *t* = 0, *N*<sub>*t*</sub> − 1, are to be
+considered primary.
+The mesh is represented by the master mesh pointer `gm` which is the
+only argument passed to function `configUser()`. The problem-relevant
+variables have to be assigned for each superblock `g` in `gm`. Recall
+the explanations about the \[mesh data
+structure\](3.1-Code-basics\#mesh-data-structure).
+There are two types of primary variables: cell-averaged variables and
+face-averaged variables. Cell-averaged variables are
+𝜚, *m*<sub>*x*</sub>, *m*<sub>*y*</sub>, *m*<sub>*z*</sub>, *Φ*, *n*<sub>*s*</sub>, *C*<sub>*c*</sub>
+and it are represented by cell-centroid coordinates. Face-averaged
+variables are the magnetic field components
+*B*<sub>*x*</sub>, *B*<sub>*y*</sub>, *B*<sub>*z*</sub> and it are
+represented by face-centroid cell coordinates.
+A superblock `g` contains ghost cells added to its active region.
+However, in the definition of IC ghost cells need not to be assigned by
+the user. It is sufficient to assign active cells, i.e. those which are
+integrated. Active cells of a superblock start at lower grid index
+(`g->ixs`,`g->iys`,`g->izs`). The upper index of active cells depends on
+the type of variable and is given in the following table together with
+the representative location of the various variables in a cell:
+| variable                                                                                              | upper index                           | cell location                         |
+|:------------------------------------------------------------------------------------------------------|:--------------------------------------|:--------------------------------------|
+| 𝜚, *m*<sub>*x*</sub>, *m*<sub>*y*</sub>, *m*<sub>*z*</sub>, *Φ*, *n*<sub>*s*</sub>, *C*<sub>*c*</sub> | (`g->ixe`,`g->iye`,`g->ize`)          | (`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`) |
+| *B*<sub>*x*</sub>                                                                                     | (`g->ixe+1,g->iye,g->ize`)            | (`g->x[ix]`,`g->yc[iy]`,`g->zc[iz]`)  |
+| *B*<sub>*y*</sub>                                                                                     | (`g->ixe`,`g->iye`+1,`g->ize`)        | (`g->xf[ix]`,`g->y[iy]`,`g->zc[iz]`)  |
+| *B*<sub>*z*</sub>                                                                                     | (`g->ixe`,`g->iye`,`g->ize`+`_C.idz`) | (`g->xf[ix]`,`g->yc[iy]`,`g->z[iz]`)  |
+*Note 1: The staggered, face-averaged vector components of the magnetic
+field have an upper index increased by one in the respective coordinate
+direction.*
+*Note 2: The code parameter `_C.idz` indicates whether a problem is 2D
+(`_C.idz`=0; axissymmetry in cylindrical/spherical geometry) or 3D
+(`_C.idz`=1).*
+Looping over the active grid cells of a superblock looks like this:
+      /* ASSIGNMENT OF CELL-AVERAGED VARIABLES */
+      for(iz=g->izs; iz<=g->ize; iz++)
+        for(iy=g->iys; iy<=g->iye; iy++)
+          for(ix=g->ixs; ix<=g->ixe; ix++)
+          { 
+            g->rho[iz][iy][ix]=...
+            ...
+          }
+      /* ASSIGNMENT OF FACE-AVERAGED MAGNETIC FIELD COMPONENTS */
+      /* X-COMPONENT */
+      for(iz=g->izs; iz<=g->ize; iz++)
+        for(iy=g->iys; iy<=g->iye; iy++)
+          for(ix=g->ixs; ix<=g->ixe+1; ix++) /* UPPER X-INDEX ! */
+            g->bx[iz][iy][ix]=...
+      /* Y-COMPONENT */
+      for(iz=g->izs; iz<=g->ize; iz++)
+        for(iy=g->iys; iy<=g->iye+1; iy++) /* UPPER Y-INDEX ! */
+          for(ix=g->ixs; ix<=g->ixe; ix++)
+            g->by[iz][iy][ix]=...
+      /* Z-COMPONENT */
+      for(iz=g->izs; iz<=g->ize+_C.idz; iz++) /* UPPER Z-INDEX ! */
+        for(iy=g->iys; iy<=g->iye; iy++)
+          for(ix=g->ixs; ix<=g->ixe; ix++)
+            g->bz[iz][iy][ix]=...
+Since tracer variables, *C*<sub>*c*</sub>, are cell-averaged quantities
+like the mass density it can be assigned at the same place. The tracer
+array index, `ic`, runs from 0 to `_C.tracer`-1 where the code parameter
+`_C.tracer` gives the number of tracers *N*<sub>*c*</sub>. Tracer
+variables are dimensionless and are usually defined in the range
+\[0, 1\].
+Similary, species number densities, *n*<sub>*s*</sub>, are cell-averaged
+quantities. The species array index, `is`, runs from 0 to `_C.species`-1
+where code parameter `_C.species` gives the number of species
+*N*<sub>*s*</sub>. The species number densities must be consistently
+defined along with the mass density expressed by the constraint
+∑<sub>*s*</sub>(*μ*<sub>*s*</sub> + *μ*<sub>*e*</sub>*q*<sub>*s*</sub>)*n*<sub>*s*</sub> = 𝜚/*m*<sub>*u*</sub>,  *s* = 0, *N*<sub>*s*</sub> − 1
+where *μ*<sub>*s*</sub>, *q*<sub>*s*</sub> and *μ*<sub>*e*</sub> are the
+molecular weight and charge of the s-th specie and the electron weight,
+respectively. In principle, number densities can be assigned
+individually on this premise. Alternatively, and much simpler, number
+densities can be computed in compliance with the user-specified mass
+fractions in the parameter file `NCCM.par` by calling the function
+      CHinitSpecies(gm);
+in the *function body* of `configUser.c` after the mass density has been
+assigned. Consistency between the species number densities and mass
+density is thereby automatically guaranteed.
+A frequently encounterd problem in the definition of IC for the magnetic
+field is to find a divergence-free formulation on the mesh consistent
+with the Constrained-Transport integral approach. This may be a
+non-trivial task for non-uniform fields, in particular on an initially
+refined mesh. There are two clean approaches to achieve a
+divergence-free setup:
+\(1\) computation of cell-face-averaged magnetic field components via
+exact integration of the analytical expressions. The so discretized
+field is per se cell-wise divergence-free. For a grid cell
+(`ix`,`iy`,`iz`) on superblock `g` this means
+$$\\mathtt{g-&gt;bx\[iz\]\[iy\]\[ix\]}=\\frac{1}{\\delta\\,\\!\\mathcal{A}\_x}\\int
+\\limits\_\\mathtt{g-&gt;y\[iy\]}^\\mathtt{g-&gt;y\[iy+1\]}
+\\int\\limits\_\\mathtt{g-&gt;z\[iz\]}^\\mathtt{g-&gt;z\[iz+1\]} B\_x(\\mathtt{g-&gt;x\[ix\]},y,z)h\_yh\_zdydz$$
+$$\\mathtt{g-&gt;by\[iz\]\[iy\]\[ix\]}=\\frac{1}{\\delta\\,\\!\\mathcal{A}\_y}\\int
+\\limits\_\\mathtt{g-&gt;x\[ix\]}^\\mathtt{g-&gt;x\[ix+1\]}
+\\int\\limits\_\\mathtt{g-&gt;z\[iz\]}^\\mathtt{g-&gt;z\[iz+1\]} B\_y(x,\\mathtt{g-&gt;y\[iy\]},z)h\_zdxdz$$
+$$\\mathtt{g-&gt;bz\[iz\]\[iy\]\[ix\]}=\\frac{1}{\\delta\\,\\!\\mathcal{A}\_z}\\int
+\\limits\_\\mathtt{g-&gt;x\[ix\]}^\\mathtt{g-&gt;x\[ix+1\]}
+\\int\\limits\_\\mathtt{g-&gt;y\[iy\]}^\\mathtt{g-&gt;y\[iy+1\]} B\_z(x,y,\\mathtt{g-&gt;z\[iz\]})h\_ydxdy$$
+The numerical expressions for the cell face contents are:
+| cell face            | expression                               |
+|:---------------------|:-----------------------------------------|
+| *δ*  𝒜<sub>*x*</sub> | `g->hyh[ix]*g->hyh[ix]*g->dvy[iy]*g->dz` |
+| *δ*  𝒜<sub>*y*</sub> | `g->dax[ix]*g->hzyh[iy]*g->dz`           |
+| *δ*  𝒜<sub>*z*</sub> | `g->dax[ix]*g->dy`                       |
+\(2\) use of the magnetic vector potential **A**, if known, with
+**B** = ∇ × **A**. When discretized in integral form this gives for the
+face-averaged magnetic field components
+`g->bx[iz][iy][ix]` = \[*h*<sub>*y*</sub>*Δ*<sub>`i`*y*</sub>(*h*<sub>*z**y*</sub>*Â*<sub>*z*</sub>)−*h*<sub>*y*</sub>*Δ*<sub>`i`*z*</sub>*Â*<sub>*y*</sub>\]/*δ*  𝒜<sub>*x*</sub>
+`g->by[iz][iy][ix]` = \[*Δ*<sub>`i`*z*</sub>*Â*<sub>*x*</sub>−*h*<sub>*z**y*</sub>*Δ*<sub>`i`*x*</sub>(*h*<sub>*y*</sub>*Â*<sub>*z*</sub>)\]/*δ*  𝒜<sub>*y*</sub>
+`g->bz[iz][iy][ix]` = \[*Δ*<sub>`i`*x*</sub>(*h*<sub>*y*</sub>*Â*<sub>*y*</sub>)−*Δ*<sub>`i`*y*</sub>*Â*<sub>*x*</sub>\]/*δ*  𝒜<sub>*z*</sub>
+where (*Â*<sub>*x*</sub>, *Â*<sub>*y*</sub>, *Â*<sub>*z*</sub>) denote
+the path integrals
+$$\\hat{A}\_x(\\mathtt{ix,iy,iz})=\\int\\limits\_\\mathtt{g-&gt;x\[ix\]}^\\mathtt{g-&gt;x\[ix+1\]}
+A\_x(x,\\mathtt{g-&gt;y\[iy\]},\\mathtt{g-&gt;z\[iz\]})dx$$
+$$\\hat{A}\_y(\\mathtt{ix,iy,iz})=\\int\\limits\_\\mathtt{g-&gt;y\[iy\]}^\\mathtt{g-&gt;y\[iy+1\]}
+A\_y(\\mathtt{g-&gt;x\[ix\]},y,\\mathtt{g-&gt;z\[iz\]})dy$$
+$$\\hat{A}\_z(\\mathtt{ix,iy,iz})=\\int\\limits\_\\mathtt{g-&gt;z\[iz\]}^\\mathtt{g-&gt;z\[iz+1\]}
+A\_z(\\mathtt{g-&gt;x\[ix\]},\\mathtt{g-&gt;y\[iy\]},z)dz$$
+over corresponding cell edges. The difference operator
+*Δ*<sub>`i`*y*</sub>, for instance, is given by
+*Δ*<sub>`i`*y*</sub>*Â*<sub>*x*</sub> = *Â*<sub>*x*</sub>(`i``x``,` `i``y` `+` `1``,` `i``z`) − *Â*<sub>*x*</sub>(`i``x``,` `i``y``,` `i``z`)
+and similarly for others. If the *Â*’s are *unambiguously* computed (or
+approximated) on coinciding cell edges in the grid hierarchy, the
+magnetic field has a consistent, cell-wise divergence-free mesh
+representation.
+*Note 3: The code can not check for the consistency of user-defined
+magnetic field configurations. If inconsistencies exist initially,
+especially on refined meshes – even not obvious – it may lead to
+magnetic monopole generation during run time.*
+If the TESTFIELDS infrastructure is enabled and IC are needed for
+testfield fluctuation variables it must be handled in the same way as
+the magnetic field. Moreover, the same issues concerning a
+divergence-free mesh representation appear as for the magnetic field.
+The array index for the testfield fluctuation variables runs from 0 to
+`_C.testfields`-1 where the code parameter `_C.testfields` defines the
+number of testfields *N*<sub>*t*</sub>.
+The setup of IC in MPI parallel simulations is more complicated than in
+serial simulations because of mesh partitioning. The mesh master pointer
+`gm` passed to the IC user interface `configUser.c` only represents the
+portion of mesh distributed to the calling MPI thread. The user has to
+take into account here the spatial locations of superblocks when
+assigning physical data to the mesh.
+Many examples of how to define IC can be found in the testproblems
+collection located in subdirectories under `/nirvana/testproblems/`.
+Here, the IC for the Orszag-Tang problem (example 1) and a shock-cloud
+interaction problem (example 2) are presented as examples. Both are MHD
+problems.
+**Example 1** (taken from `/nirvana/testproblems/MHD/problem2`; cf. )
+IC for the Orszag-Tang problem simulated in a doubly-periodic square
+domain of length *L* (given by `_C.up[0]-_C.lo[0]`):
+𝜚 = 𝜚<sub>0</sub>
+**m** = 𝜚<sub>0</sub>(−sin(2*π*  *y*/*L*),sin(2*π*  *x*/*L*),0)
+**B** = *B*<sub>0</sub>(−sin(2*π*  *y*/*L*),sin(4*π*  *x*/*L*),0)
+*e* = *p*<sub>0</sub>/(*γ* − 1) + **m**<sup>2</sup>/(2𝜚) + **B**<sup>2</sup>/(2*μ*)
+where 𝜚<sub>0</sub> = 2.77778, *p*<sub>0</sub> = 5/3 and
+*B*<sub>0</sub> = 1. The adiabatic index *γ* = 5/3 and the magnetic
+permeability *μ* = 1.
+    void configUser(GRD **gm)
+    /**
+     * Definition of ICs.
+     *
+     * @param gm  master mesh pointer
+     *
+     * @author  Udo Ziegler, AIP
+     *
+     * @note  Testproblem: MHD/problem2: Orszag-Tang vortex
+     *                     (cf. Ziegler 04, JCP 196, 393)
+     */
+    {
+      GRD *g;
+      double x0,rho0,e0,b0,pih;
+      int ix,iy,iz;
+    /* MY IMPLEMENTATION */
+    _C.permeability=1.;  /* OVERWRITE MAGNETIC PERMEABILITY TO 1 */ 
+    pih=0.5/_C.permeability;
+    x0=2.*PI/(_C.up[0]-_C.lo[0]);   
+    rho0=2.77778;
+    e0=5./3./(_C.gamma-1.);
+    b0=1.;
+    for(g=gm[0]; (g); g=g->next) /* LOOP OVER BASE LEVEL ON PARTITION */
+    {
+      for(iz=0; iz<=g->nz; iz++) 
+        for(iy=0; iy<=g->ny; iy++)
+          for(ix=0; ix<=g->nx+1; ix++)
+            g->bx[iz][iy][ix]=-b0*sin(x0*g->yc[iy]);
+      for(iz=0; iz<=g->nz; iz++) 
+        for(iy=0; iy<=g->ny+1; iy++)
+          for(ix=0; ix<=g->nx; ix++)
+            g->by[iz][iy][ix]=b0*sin(2.*x0*g->xc[ix]);
+      for(iz=0; iz<=g->nz+_C.idz; iz++)
+        for(iy=0; iy<=g->ny; iy++)
+          for(ix=0; ix<=g->nx; ix++)
+            g->bz[iz][iy][ix]=0.;
+      for(iz=0; iz<=g->nz; iz++)
+        for(iy=0; iy<=g->ny; iy++)
+          for(ix=0; ix<=g->nx; ix++) 
+          {
+            g->rho[iz][iy][ix]=rho0;
+            g->mx[iz][iy][ix]=-rho0*sin(x0*g->yc[iy]);
+            g->my[iz][iy][ix]=rho0*sin(x0*g->xc[ix]);
+            g->mz[iz][iy][ix]=0.;
+            g->e[iz][iy][ix]=e0+0.5/rho0
+              *SPQR(g->mx[iz][iy][ix],g->my[iz][iy][ix],g->mz[iz][iy][ix])
+              +SPQR(g->bx[iz][iy][ix],g->by[iz][iy][ix],g->bz[iz][iy][ix])*pih;
+          }
+    }
+    if(MASTER) fprintf(_FLOG,"..INFO(configUser): MHD/PROBLEM2 INITIALIZED.\n");
+    return;
+    }
+Example 1 uses the macro `SPQR(a,b,c)` which shortcuts the algebraic
+expression *a*<sup>2</sup> + *b*<sup>2</sup> + *c*<sup>2</sup>.
+**Example 2** (taken from `/nirvana/testproblems/MHD/problem17`; cf. )
+IC for the shock-cloud interaction problem simulated in a Cartesian box
+given by (*x*, *y*, *z*) ∈ \[ − 1/2, 1/2\]<sup>3</sup>:
+$$(\\varrho,p,v\_x,v\_y,v\_z,B\_x,B\_y,B\_z)=\\left\\{\\begin{array}{ll}
+(3.86859,167.345,0,0,0,0,2.1826182,-2.1826182) & x&lt;0.1\\\\
+(1,1,-11.2536,0,0,0,0.56418958,0.56418958) & x\\ge 0.1
+\\end{array}\\right.$$
+At **x** = (0.3, 0, 0) a spherical clump with radius 0.15 and density of
+10 is embedded and co-moving with its surrounding flow under the
+assumption of pressure equilibrium. The adiabatic index *γ* = 5/3 and
+magnetic permeability *μ* = 1.
+    void configUser(GRD **gm)
+    /**
+     * Definition of ICs.
+     *
+     * @param gm  master mesh pointer
+     *
+     * @author  Udo Ziegler, AIP
+     *
+     * @note  Testproblem: MHD/problem17: Shock-cloud interaction
+     *                     (cf. Ziegler 05, CPC 170, 153)
+     */
+    {
+      GRD *g;
+      double g1,x,y,z,dx,dy,x0,r2,vxR,eL,eR,byL,byR,bx,by,bz,pih;
+      int il,ix,iy,iz;
+    _C.permeability=1.;
+    _C.time_max=0.05;
+    pih=0.5/_C.permeability;
+    g1=_C.gamma-1.;
+    x0=0.1;
+    r2=0.1*0.1;
+    vxR=-11.2536;
+    eL=167.345/g1;
+    eR=1./g1;
+    byL=2.1826182;
+    byR=0.56418958; 
+    for(il=_C.level; il>=0; il--)
+      for(g=gm[il]; (g); g=g->next) 
+      {
+        if(_C.mf) /* MAGNETIC FIELD */
+        {
+          for(iz=0; iz<=g->nz; iz++) 
+          {
+            for(iy=0; iy<=g->ny; iy++)
+              for(ix=0; ix<=g->nx+1; ix++)
+                g->bx[iz][iy][ix]=0.;
+            for(iy=0; iy<=g->ny+1; iy++)
+              for(ix=0; ix<=g->nx; ix++)
+                g->by[iz][iy][ix]=(g->x[ix+1]<x0 ? byL 
+                  : (g->x[ix]>x0 ? byR
+                                 : ((g->x[ix+1]-x0)*byR+(x0-g->x[ix])*byL)*g->dxi));
+          }
+          for(iz=0; iz<=g->nz+_C.idz; iz++)
+            for(iy=0; iy<=g->ny; iy++)
+              for(ix=0; ix<=g->nx; ix++)
+                g->bz[iz][iy][ix]=(g->x[ix+1]<x0 ? -byL 
+                  : (g->x[ix]>x0 ? byR
+                                 : ((g->x[ix+1]-x0)*byR-(x0-g->x[ix])*byL)*g->dxi));
+        }
+        /* HD VARIABLES */
+        for(iz=0; iz<=g->nz; iz++) 
+        {
+          z=g->zc[iz];  
+          for(iy=0; iy<=g->ny; iy++) 
+          {
+            y=g->yc[iy];  
+            dy=(y-g->y[iy])*g->dyi;
+            for(ix=0; ix<=g->nx; ix++) 
+            {
+              x=g->xc[ix];
+              g->mx[iz][iy][ix]=g->my[iz][iy][ix]=g->mz[iz][iy][ix]=0.;
+              if(x<x0) 
+              {
+                g->rho[iz][iy][ix]=3.86859;
+                g->e[iz][iy][ix]=eL;
+              }
+              else 
+              {
+                g->rho[iz][iy][ix]=(x-0.3)*(x-0.3)+y*y+z*z<r2 ? 10. : 1.;
+                g->mx[iz][iy][ix]=g->rho[iz][iy][ix]*vxR;
+                g->e[iz][iy][ix]=eR+0.5*vxR*vxR*g->rho[iz][iy][ix];
+              }
+              if(_C.mf) 
+              {
+                bx=g->bx[iz][iy][ix]+(g->bx[iz][iy][ix+1]
+                   -g->bx[iz][iy][ix])*(g->xc[ix]-g->x[ix])*g->dxi;
+                by=g->by[iz][iy][ix]+(g->by[iz][iy+1][ix]
+                   -g->by[iz][iy][ix])*dy;
+                bz=0.5*(g->bz[iz][iy][ix]+g->bz[iz+_C.idz][iy][ix]);
+                g->e[iz][iy][ix]+=pih*SPQR(bx,by,bz);
+              }
+            }
+          } 
+        }
+      }
+    if(MASTER) fprintf(_FLOG,"..INFO(configUser): MHD/PROBLEM17 INITIALIZED.\n");
+    return;
+    }
+### Defining boundary conditions
+A superblock `g` includes up to 3 ghost cells per side per coordinate
+direction beyond its active region. The number of ghost cells depends on
+the physics and geometry of the problem and is stored in the code
+parameter `_C.ngc`. 3 ghost cells are needed for MHD on
+cylindrical/spherical grids and in simulations with selfgravity, 2 ghost
+cells otherwise.
+#### Standard BC
+There are a number of predefined BC types (standard BC) which can be
+selected in the parameter interface `nirvana.par` under category
+BOUNDARY CONDITION. There, standard types are specified by letters
+{I,O,D,M,A,R,C,F,P}.
+#### User-defined BC
+There is further the option that a user can define its own BC.
+User-defined BC are enabled by specifying letter U in section
+`BOUNDARY CONDITION` in `nirvana.par`. For the actual definition of the
+BC the following user interfaces must be coded, if problem-relevant:
+| BC interface  | BC for                                                                     |
+|:--------------|:---------------------------------------------------------------------------|
+| `bcrhoUser.c` | mass density 𝜚                                                             |
+| `bcmUser.c`   | momentum density **m**                                                     |
+| `bceUser.c`   | total energy density *e*                                                   |
+| `bcetUser.c`  | thermal energy density *ε* (only when dual energy mode)                    |
+| `bcbUser.c`   | magnetic field **B**                                                       |
+| `bcnXUser.c`  | species number densities *n*<sub>*s*</sub>, *s* = 0, *N*<sub>*s*</sub> − 1 |
+| `bcCUser.c`   | tracer variables *C*<sub>*c*</sub>, *c* = 0, *N*<sub>*c*</sub> − 1         |
+| `bctfUser.c`  | testfield fluctuations **b**<sub>*t*</sub>, *t* = 0, *N*<sub>*t*</sub> − 1 |
+A combination of user-defined BC and standard BC at different domain
+boundaries is possible. Examples can be found in the testproblems
+collection located in subdirectories under `/nirvana/testproblems/`,
+e.g., the double Mach reflection problem defined in
+`/nirvana/testproblems/MHD/problem1`.
+At the *x*(*y*,*z*)-lower domain boundary of a superblock `g` boundary
+values have to be assigned from `ix`=0 to `g->ixs`-1 (`iy`=0 to
+`g->iys`-1, `iz`=0 to `g->izs`-1) and full range for the cross indices.
+Recall that `g->ixs` (`g->iys`, `g->izs`) counts the first active cell
+on `g` in *x*(*y*,*z*)-direction. At the *x*(*y*,*z*)-upper domain
+boundary the `ix`(`iy`,`iz`)-index range of ghost cells depends on the
+physical variable as per
+| variable                                                                                                   | ix-range (x-upper)   | iy-range (y-upper)   | iz-range (z-upper)        |
+|:-----------------------------------------------------------------------------------------------------------|:---------------------|:---------------------|:--------------------------|
+| 𝜚, *m*<sub>*x*</sub>, *m*<sub>*y*</sub>, *m*<sub>*z*</sub>, *e*, *ε*, *n*<sub>*s*</sub>, *C*<sub>*c*</sub> | `g->ixe`+1,`g->nx`   | `g->iye`+1,`g->ny`   | `g->ize`+`_C.idz`,`g->nz` |
+| *B*<sub>*x*</sub>,*b*<sub>*x*<sub>*t*</sub></sub>                                                          | `g->ixe`+2,`g->nx`+1 | `g->iye`+1,`g->ny`   | `g->ize`+`_C.idz`,`g->nz` |
+| *B*<sub>*y*</sub>,*b*<sub>*y*<sub>*t*</sub></sub>                                                          | `g->ixe`+1,`g->nx`   | `g->iye`+2,`g->ny`+1 | `g->ize`+`_C.idz`,`g->nz` |
+| *B*<sub>*z*</sub>,*b*<sub>*z*<sub>*t*</sub></sub>                                                          | `g->ixe`+1,`g->nx`   | `g->iye`+1,`g->ny`   | `g->ize`+2,`g->nz`+1      |
+and cross indices again run through their full range. For instance, the
+following code fragment is the loop for assigning the *x*-upper BC for
+*B*<sub>*x*</sub>:
+      /* UPPER X-FACE */
+      for(ix=g->ixe+2; ix<=g->nx+1; ix++) 
+        for(iz=0; iz<=g->nz; iz++) 
+          for(iy=0; iy<=g->ny; iy++) 
+          {
+            g->bx[iz][iy][ix]=...
+          }
+#### BC for the gravitational potential
+BC for the gravitational potential *Φ* cannot be explicitely set.
+Instead, the following rule is adopted which transforms speficied types
+into BC types for *Φ*:
+| selected BC type |  ⇒  | BC type for *Φ* |
+|-----------------:|:---:|:----------------|
+|                P |     | P               |
+|            M,A,R |     | von-Neumann     |
+|            I,O,D |     | Dirichlet       |
+| user-defined (U) |     | Dirichlet       |
+In case of von-Neumann conditions the gradient of potential vanishes
+normal to the corresponding domain boundary, i.e.,
+**e**<sub>*n*</sub> ⋅ ∇*Φ* = 0 where **e**<sub>*n*</sub> is the unit
+normal vector of the domain boundary. In case of Dirichlet conditions
+boundary values for *Φ* are obtained by a multipole expansion for the
+given mass distribution. More details and formulae can be found in the
+physics guide
+([Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/A-PhysicsGuide),
+[PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/PhysicsGuide.pdf)).
+*Note: The selfgravity solver supports only triple-periodic BC or
+non-periodic BC. The combination of periodic BC at one domain boundary
+with non-periodic BC at another domain boundary is not supported.*
+### User-defined coefficients for dissipative processes
+The modules `viscosityCoeffUser.c`, `conductionCoeffUser.c` and
+`diffusionCoeffUser.c` serve as templates for a user-defined coefficient
+of fluid viscosity, thermal conduction and Ohmic diffusion,
+respectively.
+#### Fluid viscosity
+Fluid viscosity enters the momentum equation and energy equation via the
+viscous stress tensor *τ* given by
+*τ* = *ν*\[∇**v**+(∇**v**)<sup>⊤</sup>−2/3(∇⋅**v**)*I*\]
+where *ν* \[`k``g` ⋅ `m`<sup> − 1</sup> ⋅ `s`<sup> − 1</sup>\] is the
+*dynamic* viscosity coefficient, **v** the fluid velocity and *I* the
+identity operator.
+*Note: The *dynamic* viscosity coefficient *ν* to be defined here should
+not be confused with the kinetic coefficient,*
+*ν*<sub>*k**i**n**e**t**i**c*</sub>
+\[`m`<sup>2</sup> ⋅ `s`<sup> − 1</sup>\], *related to the dynamic
+coefficient by* *ν* = 𝜚*ν*<sub>*k**i**n**e**t**i**c*</sub>.
+In the call of `viscosityCoeffUser()` the function arguments are the
+superblock pointer `g` and the array pointer `vis` of type `double***`
+representing the viscosity coefficient *ν*:
+      viscosityCoeffUser(g,vis);
+The user must assign the `vis`-array for the full index range of `g`,
+whence a loop like
+      /* MY VISCOSITY COEFFICIENT */
+      for(iz=0; iz<=g->nz; iz++)
+        for(iy=0; iy<=g->ny; iy++) 
+          for(ix=0; ix<=g->nx; ix++) 
+          {
+            vis[iz][iy][ix]=...
+          }
+must be added in `viscosityCoeffUser.c`. The coefficient should be
+evaluated at cell-centroid coordinates
+(`g->xc[ix],g->yc[iy],g->zc[iz]`). An example definition for a viscosity
+coefficient can be found in testproblem
+`/nirvana/testproblems/VISC/problem1`.
+User-defined viscosity is enabled by appropriate choice in the parameter
+interface `nirvana.par` under the category `PHYSICS SPECIFICATIONS`
+(code parameter: `_C.viscosity`).
+#### Thermal conduction
+Thermal conduction enters the energy equation through a heat flux
+**F**<sub>*C*</sub>. Generally, the presence of a magnetic field
+introduces anisotropic effects with different transport properties along
+and across the magnetic field. **F**<sub>*C*</sub> is described by
+**F**<sub>*C*</sub> =  − *κ*<sub>∥</sub>(∇*T* ⋅ **B̂**)**B̂** − *κ*<sub>⊥</sub>(∇*T*−(∇*T*⋅**B̂**)**B̂**)
+where *κ*<sub>∥</sub> and *κ*<sub>⊥</sub> is the thermal conduction
+coefficient parallel and perpendicular to the magnetic field,
+respectively, and **B̂** = **B**/\|**B**\| is the unit vector in the
+direction of the magnetic field. *κ*<sub>∥</sub> and *κ*<sub>⊥</sub> are
+measured in units
+`J` ⋅ `K`<sup> − 1</sup> ⋅ `m`<sup> − 1</sup> ⋅ `s`<sup> − 1</sup>.
+In the call of `conductionCoeffUser()` the function arguments are the
+superblock pointer `g` and the array pointers `cond`, `cond_perp` of
+type `double***` representing the conduction coefficients
+*κ*<sub>∥</sub>, *κ*<sub>⊥</sub>:
+      conductionCoeffUser(g,cond,cond_perp);
+The user must assign the `cond,cond_perp`-arrays for the full index
+range of `g`, whence a loop like
+      /* MY CONDUCTION COEFFICIENT */
+      for(iz=0; iz<=g->nz; iz++)
+        for(iy=0; iy<=g->ny; iy++) 
+          for(ix=0; ix<=g->nx; ix++) 
+          {
+            cond[iz][iy][ix]=...
+            cond_perp[iz][iy][ix]=...
+          }
+must be added to `conductionCoeffUser.c`. The coefficients should be
+evaluated at cell-centroid coordinates
+(`g->xc[ix],g->yc[iy],g->zc[iz]`). In HD simulations or MHD simulations
+with forced isotropy (when `COND_FORCE_ISO`=`YES` in `nirvanaUser.h`)
+only the `cond`-array has to be assigned with `cond` representing the
+conduction coefficient *κ* in the isotropic heat flux
+**F**<sub>*C*</sub> =  − *κ*∇*T*.
+User-defined thermal conduction is enabled by appropriate choice in the
+parameter interface `nirvana.par` under the category
+`PHYSICS SPECIFICATIONS` (code parameter: `_C.conduction`).
+#### Ohmic diffusion
+Ohmic diffusion enters the induction equation and energy equation as a
+field contribution given by
+**E**<sub>*D*</sub> =  − *η*<sub>*D*</sub>∇ × **B**
+where *η*<sub>*D*</sub> \[`m`<sup>2</sup> ⋅ `s`<sup> − 1</sup>\] is the
+diffusion coefficient.
+In the call of `diffusionCoeffUser()` function arguments are the
+superblock pointer `g` and the array pointer `diff` of type `double***`
+representing the diffusion coefficient *η*<sub>*D*</sub>:
+      diffusionCoeffUser(g,diff);
+The user must assign the `diff`-array for the full index range of `g`,
+whence a loop like
+      /* MY DIFFUSION COEFFICIENT */
+      for(iz=0; iz<=g->nz; iz++)
+        for(iy=0; iy<=g->ny; iy++) 
+          for(ix=0; ix<=g->nx; ix++) 
+          {
+            diff[iz][iy][ix]=...
+          }
+must be added in `diffusionCoeffUser.c`. The coefficient should be
+evaluated at cell-centroid coordinates
+(`g->xc[ix],g->yc[iy],g->zc[iz]`).
+User-defined Ohmic diffusion is enabled by appropriate choice in the
+parameter interface `nirvana.par` under the category
+`PHYSICS SPECIFICATIONS` (code parameter: `_C.diffusion`).
+### User-defined coefficient for ambipolar diffusion
+The module `APdiffusionCoeffUser.c` serves as template for a
+user-defined ambipolar diffusion coefficient. Ambipolar diffusion enters
+the induction equation and energy equation as a field contribution given
+by
+**E**<sub>*A**D*</sub> = *η*<sub>*A**D*</sub>/*μ*\[(∇×**B**)×**B**\] × **B**
+where *η*<sub>*A**D*</sub>
+\[`V` ⋅ `m` ⋅ `A`<sup> − 1</sup> ⋅ `T`<sup> − 2</sup>\] denotes the
+ambipolar diffusion coefficient.
+*Note: The prefactor* *η*<sub>*A**D*</sub>/*μ* *has units*
+`m`<sup>2</sup> ⋅ `s`<sup> − 1</sup> ⋅ `T`<sup> − 2</sup>.
+In the call of `APdiffusionCoeffUser()` function arguments are the
+superblock pointer `g` and the array pointer `APdiff` of type
+`double***` representing the ambipolar diffusion coefficient
+*η*<sub>*A**D*</sub>:
+      APdiffusionCoeffUser(g,APdiff);
+The user must assign the `APdiff`-array for the full index range of `g`,
+whence a loop like
+      /* MY AMBIPOLAR DIFFUSION COEFFICIENT */
+      for(iz=0; iz<=nz; iz++)
+        for(iy=0; iy<=ny; iy++) 
+          for(ix=0; ix<=nx; ix++) 
+          {
+            APdiff[iz][iy][ix]=...
+          }
+must be added in `APdiffusionCoeffUser.c`. The coefficient should be
+evaluated at cell-centroid coordinates
+(`g->xc[ix],g->yc[iy],g->zc[iz]`). An example definition for an
+ambipolar diffusion coefficient can be found in testproblem
+`/nirvana/testproblems/APDIFF/problem1`.
+User-defined ambipolar diffusion is enabled by appropriate choice in the
+parameter interface `nirvana.par` under the category
+`PHYSICS SPECIFICATIONS` (code parameter: `_C.APdiffusion`).
+### User-defined body force
+The module `forceUser.c` serves as template for coding a user-defined
+*specific* body force **f** (force per mass in units
+*N* ⋅ *k**g*<sup> − 1</sup>). The body force then enters the momentum
+equation and energy equation as source term.
+In the call of `forceUser()` the function arguments are the superblock
+pointer `g` and the array pointers `fx`,`fy`,`fz` of type `double***`
+representing the force components:
+      forceUser(g,fx,fy,fz);
+The user must assign the `fx`,`fy`,`fz`-arrays for active grid cells of
+`g`, whence a loop like
+      /* MY BODY FORCE */
+      for(iz=g->izs; iz<=g->ize; iz++)
+        for(iy=g->iys; iy<=g->iye; iy++) 
+          for(ix=g->ixs; ix<=g->ixe; ix++) 
+          {
+            fx[iz][iy][ix]=...
+            fy[iz][iy][ix]=...
+            fz[iz][iy][ix]=...
+          }
+must be added in `forceUser.c`. Strictly speaking, the body force is to
+be understood as a cell-averaged quantity. Therefore, **f** should be
+evaluated at cell-centroid coordinates (`g->xc[ix],g->yc[iy],g->zc[iz]`)
+and the point values assigned to `fx`,`fy`,`fz`. An example definition
+for a body force can be found in testproblem
+`/nirvana/testproblems/MHD/problem21`.
+The body force is enabled by specification in the parameter interface
+`nirvana.par` under the category `PHYSICS SPECIFICATIONS` (code
+parameter: `_C.force`).
+### User-defined cooling/heating function
+The modules `sourceCoolingUser.c` and `sourceHeatingUser.c` serve as
+templates for coding a user-defined cooling function,
+*L*<sub>*c**o**o**l*</sub> ≤ 0, and heating function,
+*L*<sub>*h**e**a**t*</sub> ≥ 0, respectively. Both functions are allowed
+to depend on temperature *T* and density 𝜚. The net heatloss, i.e. the
+sum
+*L*<sub>*c**o**o**l*</sub>(*T*, 𝜚) + *L*<sub>*h**e**a**t*</sub>(*T*, 𝜚)
+of both functions, enters as a source term in the energy equation.
+*L*<sub>*c**o**o**l*</sub> and *L*<sub>*h**e**a**t*</sub> are measured
+in units *J* ⋅ *s*<sup> − 1</sup> ⋅ *m*<sup> − 3</sup>.
+In the call of `sourceCoolingUser()` (`sourceHeatingUser()`) the
+function arguments are the temperature value `T`, density value `rho`,
+the pointer `deriv` to the derivatives flag and the 2-element vector
+`dfc` (`dfh`):
+      f=sourceCoolingUser(T,rho,deriv,dfc);
+      f=sourceHeatingUser(T,rho,deriv,dfh);
+The user must define the return value
+`f` = *L*<sub>*c**o**o**l*</sub>(`T``,` `r``h``o`)
+(*L*<sub>*h**e**a**t*</sub>(`T``,` `r``h``o`)) in `sourceCoolingUser.c`
+(`sourceHeatingUser.c`). The `deriv`-flag is thought to indicate the
+calling function whether a user provides the derivatives of
+*L*<sub>*c**o**o**l*</sub>(*T*, 𝜚) and
+*L*<sub>*h**e**a**t*</sub>(*T*, 𝜚) with respect to *T* and 𝜚 himself.
+The `deriv`-flag is preset to `NO` when entering the code function. If
+the user sets
+      *deriv=YES; 
+the derivatives must be stored by the user in the 2-element vector `dfc`
+(`dfh`) with `dfc[0]` (`dfh[0]`) storing the derivative
+∂*L*<sub>*c**o**o**l*</sub>/∂*T* (∂*L*<sub>*h**e**a**t*</sub>/∂*T*) and
+`dfc[1]` (`dfh[1]`) storing the derivative
+∂*L*<sub>*c**o**o**l*</sub>/∂𝜚 (∂*L*<sub>*h**e**a**t*</sub>/∂𝜚).
+Otherwise derivatives are automatically computed by a finite difference
+approximation.
+*Note: Knowledge of the derivatives is a prerequisite in the numerical
+procedure, the Field-length-based mesh refinement criterion and timestep
+computation.*
+An example definition for a cooling/heating function can be found in
+testproblem `/nirvana/testproblems/HEATLOSS/problem2`.
+User-defined cooling/heating is enabled by appropriate choice in the
+parameter interface `nirvana.par` under the category
+`PHYSICS SPECIFICATIONS` (code parameter: `_C.heatloss`).
+### User-defined equation of state
+#### Analytic EOS
+NIRVANA makes possible a simple user-defined implementation of an
+analytic EOS. This requires the defintion of the macros
+| macro             | thermodynamic quantity |
+|:------------------|:-----------------------|
+| `PUSR(rho,eth)`   | pressure               |
+| `CS2USR(rho,eth)` | square of sound speed  |
+| `TUSR(rho,eth)`   | temperature            |
+| `ETUSR(rho,eth)`  | thermal energy density |
+as a function of mass density 𝜚 (`rho`) and thermal energy density *ε*
+(`eth`). The macros are declared in the header file `nirvanaUser.h`. The
+following code fragment (taken from testproblem
+`/nirvana/testproblems/GRAVITY/problem4`) illustrates the definition of
+a barotropic EOS of the form
+*p*(𝜚) = *a*<sup>2</sup>𝜚\[1+(𝜚/𝜚<sub>0</sub>)<sup>4/3</sup>\]<sup>1/2</sup>
+with *a* = 200 and 𝜚<sub>0</sub> = 10<sup> − 10</sup>.
+    #define PUSR(rho,eth)   (40000.*(rho)*sqrt(1.+pow((rho)/1.e-10,4./3.))+0*(eth))
+    #define CS2USR(rho,eth) (40000.*sqrt(1.+pow((rho)/1.e-10,4./3.)) \
+                            +2./3.*40000.*pow((rho)/1.e-10,4./3.) \
+                            /sqrt(1.+pow((rho)/1.e-10,4./3.))+0*(eth)) 
+    #define TUSR(rho,eth)   (MOK*_C.mean_molecular_weight*PUSR(rho,eth)/rho)
+    #define ETUSR(rho,eth)  (PUSR(rho,eth)/(_C.gamma-1.))
+Since *p* solely depends on 𝜚 the energy density parameter `eth` is
+irrelevant is this case. The `ETUSR` definition makes use of the prior
+definition of `PUSR` and assumes the usual caloric EOS with ratio of
+specific heats *γ* (code variable: `_C.gamma`). No energy equation is
+solved in this testproblem.
+*Note 1: If the energy equation is solved together with an analytic EOS
+the macro `ETUSR` must read*
+    #define ETUSR(rho,eth)  (eth)
+for obvious reasons of consistency.
+*Note 2: The macros are used in function `TDeos()` in module `utilTD.c`
+and are needed in the numerical procedure.*
+An analytic EOS is enabled by appropriate choice in the parameter
+interface `nirvana.par` under category `PHYSICS SPECIFICATIONS` (code
+paramter: `_C.eos`).
+#### Tabulated EOS
+NIRVANA also permits the implementation of a tabulated EOS. It requires
+the creation of look-up tables for the logarithm of pressure, log *p*,
+and logarithm of temperature, log *T*, both as functions of the
+logarithm of mass density 𝜚 and the logarithm of thermal energy density
+*ε*. The look-up tables for a given (log 𝜚, log *ε*)-domain are
+generated by cubic interpolation from user-specified sampling data. The
+sampling data {(log *p*)<sub>*i**j*</sub>}
+({(log *T*)<sub>*i**j*</sub>}) on the grid
+{(log 𝜚)<sub>*i*</sub>} × {(log *ε*)<sub>*j*</sub>} has to be defined by
+the user in the interface `eosTabPressureUser.c`
+(`eosTabTemperatureUser.c`). The discretizations in log 𝜚-space and
+log *ε*-space, $\\{(\\log\\varrho)\_i; i=0,{\\tt n1}-1\\}$ and
+$\\{(\\log\\varepsilon)\_j; j=0,{\\tt n2}-1\\}$, must be equidistant
+where `n1` and `n2` are the number of sampling points.
+Actually, in `eosTabPressureUser()` (`eosTabTemperatureUser()`) the user
+must return a structure of type `UTAB` (declared in header `tabular.h`)
+which collects all relevant sampling data. Denoting by `ut` an instance
+of `UTAB` the structure elements are accessed by `ut.V`. The following
+structure elements `V` must be defined by the user:
+| `V`             | meaning                                                                                |
+|:----------------|:---------------------------------------------------------------------------------------|
+| `n1`            | number of sampling points in log 𝜚                                                     |
+| `n2`            | number of sampling points in log *ε*                                                   |
+| `u1_lo`,`u1_up` | lower,upper value of log 𝜚-range                                                       |
+| `u2_lo`,`u2_up` | lower,upper value of log *ε*-range                                                     |
+| `dt2[i][j]`     | 2D array for sampling data {(log *p*)<sub>*i**j*</sub>} ({(log *T*)<sub>*i**j*</sub>}) |
+The subsequent example taken from testproblem
+`/nirvana/testproblem/MHD/problem1B` provides sampling data for
+tabulation of the caloric EOS *p* = (*γ* − 1)*ε*.
+    UTAB eosTabPressureUser(void)
+    /**
+     * Template for the definition of user-defined sampling data for ln(pressure)
+     * on 2D [ln(density),ln(thermal_energy_density)]-space. The sample data is 
+     * used for creation of the 2D look-up table _TABLP in a tabulated EOS.
+     *
+     * @retval ut  UTAB structure for creating the 2D look-up table _TABLP
+     *
+     * @author  Udo Ziegler, AIP
+     *
+     * @note  The UTAB structure is declared in tabular.h and contains the 
+     *        following elements which has to be assigned here by the user:
+     *                        ut.n1 = number of equidistant points for ln(density);
+     *                        ut.n2 = number of equidistant points for 
+     *                                ln(thermal_energy_density);
+     *          [ut.u1_lo,ut.u1_up] = ln(density) variable range;
+     *          [ut.u2_lo,ut.u2_up] = ln(thermal_energy_density) variable range;
+     *                       ut.dt2 = 2D array for ln(pressure) sampling data.
+     *        Example implementation: tabulated ideal EOS for testproblem 
+     *                                MHD/problem1B.
+     */
+    {
+      UTAB ut;
+      double *let,dlet,lc;
+      int i1,i2,n1,n2;
+    /* MY IMPLEMENTATION: IDEAL EOS FOR PROBLEM MHD/PROBLEM1B ------------------- */
+    ut.n1=n1=100; /* #POINTS IN LN(DENSITY)-DIRECTION */
+    ut.n2=n2=100; /* #POINTS IN LN(THERMAL_ENERGY_DENSITY)-DIRECTION */
+    /* LN(DENSITY)-RANGE */
+    ut.u1_lo=-4.61;
+    ut.u1_up=4.61;
+    /* LN(THERMAL_ENERGY_DENSITY)-RANGE */
+    ut.u2_lo=-9.22;
+    ut.u2_up=9.22;
+    /* LN(THERMAL_ENERGY_DENSITY) VECTOR */
+    let=Array(n2-1);
+    dlet=(ut.u2_up-ut.u2_lo)/(n2-1);
+    let[0]=ut.u2_lo;
+    for(i2=1; i2<n2; i2++) 
+      let[i2]=let[i2-1]+dlet;
+    /* LN(PRESSURE) SAMPLING DATA FOR ln(p)=ln(_C.gamma-1)+ln(et) */
+    ut.dt2=Array2(n1-1,n2-1);
+    lc=log(_C.gamma-1.);
+    for(i1=0; i1<n1; i1++)
+      for(i2=0; i2<n2; i2++)
+        ut.dt2[i1][i2]=lc+let[i2];
+    free(let);
+    return ut;
+    }
+First, the number of sampling points `ut.n1` in (log 𝜚)-direction and
+`ut.n2` in (log *ε*)-direction are set followed by the ranges
+\[`ut.u1_lo`,`ut.u1_up`\] for log 𝜚 and \[`ut.u2_lo`,`ut.u2_up`\] for
+log *ε*, respectively. Then, the 2D array `ut.dt2` is allocated and
+sampling data {(log *p*)<sub>*i**j*</sub>} is assigned. The returned
+`ut` is used by the calling function to finally create the public
+look-up table `_TABLP`. A corresponding user definition for log *T* in
+`eosTabTemperatureUser.c` provides the input for the public look-up
+table `_TABLT`.
+*Note: `_TABLP` and `_TABLT` are actually pointers of struct type `TAB`
+(declared in the header `tabular.h`). It are used in the function
+`TDeos()` in module `utilTD.c` and in the numerical procedure.*
+A tabulated EOS is enabled by appropriate choice in the parameter
+interface `nirvana.par` under category `PHYSICS SPECIFICATIONS` (code
+paramter: `_C.eos`).
+### User-defined initial/restricted mesh refinement
+For certain problems it may be advantageous to start a simulation with a
+pre-refined mesh in some parts of the computational domain or to retrict
+mesh refinement to some parts of it during runtime. Such situations can
+be realize with help of the user interfaces `initDomainUser.c` and
+`checkDomainUser.c`.
+#### Initial mesh refinement
+The module `initDomainUser.c` serves as template to initialize a
+user-defined refined mesh at simulation start. `initDomainUser()` is
+called by the mesh refinement algorithm. Function arguments are the
+coordinates of a testpoint (`x`,`y`,`z`):
+      level=initDomainUser(x,y,z);
+The user must return the mesh refinement level (`level`) with which the
+vicinity of the given testpoint should be resolved with the limitation
+\|`l``e``v``e``l`\| ≤ `_``C``.``i``m``r` where code paramter `_C.imr` is
+the maximum allowed initial mesh refinement level specified by the user
+in the parameter interface `nirvana.par` under the category
+`MESH REFINEMENT`. In practice, the user defines spatial domains in
+`initDomainUser.c` via mathematical relations, checks whether the
+testpoint is contained and then returns the requested refinement level.
+The following code fragment illustrates the procedure for a spherical
+interface at radius 0.1 relative to the center
+(*x*, *y*, *z*) = (0.3, 0, 0) to be resolved with refinement level 4:
+    int initDomainUser(double x, double y, double z)
+    /**
+     * Template for the specification of user-defined spatial domains for the 
+     * purpose of initial mesh refinement.
+     *
+     * @param x,y,z  spatial coords of testpoint
+     *
+     * @retval level  requested refinement level:{|level|<=_C.imr}
+     *
+     * @author  Udo Ziegler, AIP
+     *
+     * @note  Domains must be defined via mathematical relations. The function
+     *        then checks whether a testpoint (x,y,z) lies inside a specified 
+     *        domain. If this is the case the user must return the requested 
+     *        refinement level for this domain. Negative values (level<0) 
+     *        indicate that the domain should be statically refined (i.e. there
+     *        is no derefinement during runtime). For positive values the refined
+     *        domain is subject to subsequent derefinement according to the regular 
+     *        criteria. Multiple domains are allowed.
+     */
+    {
+      int level=0;
+      double R;
+    /* SPHERICAL INTERFACE CHECK */
+    R=sqrt((x-0.3)*(x-0.3)+y*y+z*z);
+    if(R>0.025 && R<0.175) level=1;
+    if(R>0.05 && R<0.15) level=2;
+    if(R>0.075 && R<0.125) level=4;
+    return level;
+    }
+The coding here realizes a buffer zone approach by placing three hollow
+globes of decreasing thickness around the interface. This ensures that
+the interface is safely resolved with the requested degree of
+refinement.
+*Note: The use of buffer zones is strongly recommended in the setup of
+highly refined, localized regions since the refinement algorithm does
+not allow jumps in refinement more than one level.*
+Further examples can be found in the testproblem
+`/nirvana/testproblems/MHD/problem16` and, for a more complicated
+situation, in the testproblem `/nirvana/testproblems/GRAVITY/problem1`.
+User-specified refined domains are normally subject to subsequent
+derefinement according to the standard refinement criteria. However, if
+the returned refinement level is given a negative sign
+(`l``e``v``e``l` &lt; 0) it tells the calling function that the
+user-specified domain should be regarded as statically refined. This
+implies that there is no subsequent derefinement of this domain during
+runtime.
+#### Restricted mesh refinement
+The module `checkDomainUser.c` serves as template for a user-defined
+control mechanism which allows to restrict mesh refinement during
+runtime. `checkDomainUser()` is called by the mesh refinement algorithm.
+Function arguments are the coordinates of a testpoint (`x`,`y`,`z`):
+      level=checkDomainUser(x,y,z);
+The user must return the refinement control parameter `level`
+(initialized to 0) for the vicinity of the testpoint (`x`,`y`,`z`). A
+return value `l``e``v``e``l` = 0 has no effect on mesh refinement. A
+return value `l``e``v``e``l` =  − 1 prohibits mesh refinement. A return
+value `l``e``v``e``l` &gt; 0 forces mesh refinement up to refinement
+level `level` with the limitation `l``e``v``e``l` ≤ `_``C``.``a``m``r`
+where code parameter `_C.amr` is the maximum allowed mesh refinement
+level specified by the user in the parameter interface `nirvana.par`
+under category `MESH REFINEMENT`. Like in the `initDomainUser.c`
+interface the user can define spatial domains via mathematical
+relations, checks whether the testpoint is contained and then returns
+the requested refinement control parameter.
+An example can be found in the testproblem
+`/nirvana/testproblems/GRAVITY/problem3`.
+### Specification of NCCM parameters
+The parameter file `NCCM.par` serves as user interface to the
+multi-species framework/NCCM. `NCCM.par` is grouped into the category
+sections
+-   `SPECIES SPECIFICATION` – selection of species.
+-   `REACTIONS SPECIFICATION` – selection of chemical reactions.
+-   `THERMAL PROCESSES SPECIFICATION` – selection of microphysical
+    thermal processes.
+A section starts with the category name preceded by the ’`>`’-character.
+:
+    >SPECIES SPECIFICATION ---------------------------------------------------------
+    > LABEL      H__D__He_C__N__O__Ne_Mg_Si_Fe   CHARGE  MASSFRAC
+    + H          1  0  0  0  0  0  0  0  0  0    0       7.000e-01
+    + H+         1  0  0  0  0  0  0  0  0  0    1       1.000e-04
+    + He         0  0  1  0  0  0  0  0  0  0    0       2.800e-01
+    + He+        0  0  1  0  0  0  0  0  0  0    1       0.000e+00
+    + He++       0  0  1  0  0  0  0  0  0  0    2       0.000e+00
+            .          .
+            .          .
+            .          .
+The following species are currently contained in the multi-species
+framework:
+-   chemical elements H, D, He, C, N, O, Ne, Mg, Si, Fe and their ions
+-   negatively charged elements H<sup>−</sup>, C<sup>−</sup>,
+    O<sup>−</sup>
+-   molecules H<sub>2</sub>, H<sub>2</sub>+, H<sub>3</sub><sup>+</sup>,
+    HD, O<sub>2</sub>, C<sub>2</sub>, O<sub>2</sub><sup>+</sup>, CH,
+    CH<sup>+</sup>, CH<sub>2</sub>, CH<sub>2</sub><sup>+</sup>,
+    CH<sub>3</sub><sup>+</sup>, OH, OH<sup>+</sup>, CO, CO<sup>+</sup>,
+    H<sub>2</sub>O, H<sub>2</sub>O<sup>+</sup>,
+    H<sub>3</sub>O<sup>+</sup>, HCO<sup>+</sup> and HOC<sup>+</sup>
+Each specie is defined by a unique label (`LABEL`), its chemical
+signature (H\_\_D\_\_He\_C\_\_N\_\_O\_\_Ne\_Mg\_Si\_Fe), i.e., the
+chemical elements a specie consists of, and its charge (`CHARGE`) in
+units of the elementary charge e. These basic species properties should
+*NOT* be changed here by the user. A specie is selected (unselected) by
+setting a `+` (`-`) sign in front of the specie label. Optionally, the
+user can assign each specie its initial mass fraction (`MASSFRAC`) which
+considerably simplifies the definition of IC for the species number
+densities in the user interface `configUser.c` (cf.
+\[here\](3.2-User-interfaces\#user-defined-intial-conditions)). The
+species properties are collected in a substructure `_C.X[is]` of the
+public structure `_C` where `is` denotes the species index.
+*Note 1: If the sum of `MASSFRAC` over all species does not amount to
+unity the numbers are proportionally rescaled. This changes their
+values.*
+*Note 2: A specie electron with label e- is automatically generated by
+the code.*
+`REACTIONS SPECIFICATION`:
+    >REACTIONS SPECIFICATION -------------------------------------------------------
+    > ID   STCOEFF REACTANT [& STCOEFF REACTANT] > STCOEFF PRODUCT [& STCOEFF PRODUCT]
+    + 1    1 H    & 1 e-   > 1 H+   & 2 e-
+    + 2    1 H+   & 1 e-   > 1 H
+    + 3    1 He   & 1 e-   > 1 He+  & 2 e-
+    + 4    1 He+  & 1 e-   > 1 He
+    + 5    1 He+  & 1 e-   > 1 He++ & 2 e-
+    + 6    1 He++ & 1 e-   > 1 He+
+            .          .
+            .          .
+            .          .
+The chemical reactions in section `REACTIONS SPECIFICATION` are grouped
+into the classes:
+-   subnetworks for the ionization structure of each element
+-   subnetworks for H<sub>2</sub> formation and HD formation
+-   C-bearing and O-bearing chemical cycles
+-   cosmic ray reactions
+-   dust-assisted reactions
+Each reaction is characterized by an ID (a running number) followed by
+its definition according to the pattern
+           STCOEFF REACTANT [& STCOEFF REACTANT] > STCOEFF PRODUCT [& STCOEFF PRODUCT]
+with stoichiometric coefficient STCOEFF, reactants REACTANT and products
+PRODUCT. A reaction is selected (unselected) by setting a `+ (-)` sign
+in front of the reaction ID.
+*Note: The code parser checks each reaction for charge conservation and
+mass conservation as well as the existence of the labels of involved
+species. This rection check can be disabled by setting the macro
+`CHECK_REACTIONS` to `NO` in `nirvanaUser.h`.*
+`THERMAL PROCESSES SPECIFICATION`:
+    >THERMAL PROCESSES SPECIFICATION -----------------------------------------------
+    >ON/OFF 
+     Y   >BREMSSTRAHLUNG: H+/He+/He++
+     Y   >H: CI+RR+CE COOLING
+     N   >D: CI+RR+CE COOLING
+     Y   >He: CI+RR+DR+CE COOLING
+     Y   >C: CE COOLING
+     N   >N: CE COOLING
+     Y   >O: CE COOLING
+     N   >Ne: CE COOLING
+     Y   >Mg: CE COOLING
+     Y   >Si: CE COOLING
+     Y   >Fe: CE COOLING
+     Y   >H2: CE+CD+CIE COOLING / H2/H2+/H- FORMATION HEATING
+     N   >HD: CE COOLING
+     Y   >CR HEATING
+     N   >CO: CE COOLING
+     N   >H2O: CE COOLING
+     N   >OH: CE COOLING
+The NCCM currently contains the following microphysical thermal
+processes:
+-   Bremsstrahlung for the H, He and He ions
+-   chemical- and atomic line cooling of H, D and He
+-   fine-structure line cooling of the metals C, N, O, Ne, Mg, Si and Fe
+-   chemical-, rotovibrational line- and collision-induced cooling of
+    H<sub>2</sub>
+-   rotovibrational line cooling of HD
+-   cosmic ray heating
+-   rotovibrational line cooling of CO
+-   rotovibrational line cooling of H<sub>2</sub>O
+-   rotational line cooling of OH
+A thermal process is activated (deactivated) by specifying the value Y
+(N). Every text in a line following the ’`>`’-character is ignored by
+the parser and serves for comments.
+### User-controllable macros
+Besides the parameters collected in the user interfaces `nirvana.par`
+and `NCCM.par` there are a number of further parameters which can be
+controlled by the user. In practice, however, the user rarely has to
+change them. These parameters are defined in the header file
+`nirvanaUser.h` as macros and it are grouped into the categories:
+-   NUMERICS-RELATED MACROS
+-   PHYSICS-RELATED MACROS
+-   NCCM-RELATED MACROS
+-   MACROS RELATED TO ANALYTICAL EOS
+The following table compiles a complete list of macros including a short
+description of their meaning and, in some cases, their permissable
+values (given in brackets).
+**NUMERICS-RELATED MACROS**:
+-   `SPACE_ORDER` (2): spatial order of the numerical scheme
+-   `LIM(a)` ({MM(a),MC(a),VL(a)}): type of slope limiter in the
+    reconstruction procedure (MM=MinMod, MC=Monotonized-Centered,
+    VL=VanLeer)
+-   `LIM_MULTI_DIM` ({YES,NO}): option to use the multi-d limiter.
+-   `PSI` (\[0.5,1\]): parameter in multi-d limiter
+-   `MAXLEVEL` ( &lt; 128): maximum refinement level
+-   `MG_ITMAX`: maximum number of iterations allowed in the multigrid
+    solver
+-   `MG_TYPE` (MULT): used type of multigrid solver
+-   `MG_TOL` (typical: 10<sup> − 6</sup>): error tolerance in the
+    multigrid solver
+-   `RKL_COURANT_EXPL` ( &lt; 0.5): explicit Courant number used in the
+    RKL solver
+-   `RKL_MAX_COURANT` (typical:  &lt; 1000): maximum Courant-like number
+    allowed in the RKL solver
+-   `RKL_DT_LIM`: RKL-related timestep limit as fraction of the
+    dynamical timestep
+-   `HEATLOSS_TOL` (typical: 10<sup> − 5</sup>): relative error
+    tolerance in heatloss solver
+-   `HEATLOSS_ATOL`: absolute error tolerance in heatloss solver
+-   `REACTIONS_TOL` (typical: 10<sup> − 5</sup>): relative error
+    tolerance in NCCM solver
+-   `REACTIONS_ATOL_X`: absolute error tolerance for number densities in
+    NCCM solver
+-   `REACTIONS_ATOL_T`: absolute error tolerance for the temperature in
+    NCCM solver
+**PHYSICS-RELATED MACROS**
+-   `CENTRIFUGAL_FORCE` ({YES,NO}): option to enable/disable the
+    centrifugal force term
+-   `COND_FORCE_ISO` ({YES,NO}): option to force isotropy in thermal
+    conduction
+**NCCM-RELATED MACROS**
+-   `CR_IH`: cosmic ray ionisation rate for hydrogen
+-   `TDUST`: dust temperature
+-   `ADUST`: dust abundance relative to Milky Way
+-   `ISRF_G` (Milky Way: 1.13): Habing strength of interstellar
+    radiation field
+-   `DBNCCM` ({YES,NO}): NCCM database availability flag
+-   `CHECK_REACTIONS` ({YES,NO}): option to check chemical reactions
+    with regard to mass conservation, charge conservation and existence
+    of species labels
+**MACROS RELATED TO ANALYTICAL EOS**
+-   `PUSR(rho,eth)`: analytic expression for the pressure *p* as a
+    function of 𝜚 (`rho`) and *ε* (`eth`)
+-   `CS2USR(rho,eth)`: analytic expression for the square of sound speed
+    *c*<sub>*s*</sub><sup>2</sup>(𝜚, *ε*)
+-   `TUSR(rho,eth)`: analytic expression for the temperature *T*(𝜚, *ε*)
+-   `ETUSR(rho,eth)`: analytic expression of the thermal energy density
+    *ε*(𝜚) as a function of 𝜚 in cases without energy equation.
+    Otherwise identity.