Skip to content

Changes

Page history
Update 3.2 User interfaces authored by Udo Ziegler's avatar Udo Ziegler
Hide whitespace changes
Inline Side-by-side
3-NIRVANA-user-guide/3.2-User-interfaces.md
View page @ 4ccb63fc
......@@ -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>*Â*<sub>*z*</sub>)−*h*<sub>*y*</sub>*Δ*<sub>`i`*z*</sub>*Â*<sub>*y*</sub>\]/*δ* 𝒜<sub>*x*</sub>
`g->by[iz][iy][ix]` = \[*Δ*<sub>`i`*z*</sub>*Â*<sub>*x*</sub>−*h*<sub>*z**y*</sub>*Δ*<sub>`i`*x*</sub>(*h*<sub>*y*</sub>*Â*<sub>*z*</sub>)\]/*δ* 𝒜<sub>*y*</sub>
`g->bz[iz][iy][ix]` = \[*Δ*<sub>`i`*x*</sub>(*h*<sub>*y*</sub>*Â*<sub>*y*</sub>)−*Δ*<sub>`i`*y*</sub>*Â*<sub>*x*</sub>\]/*δ* 𝒜<sub>*z*</sub>
where (*Â*<sub>*x*</sub>, *Â*<sub>*y*</sub>, *Â*<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>*Â*<sub>*x*</sub> = *Â*<sub>*x*</sub>(`i``x``,` `i``y` `+` `1``,` `i``z`) − *Â*<sub>*x*</sub>(`i``x``,` `i``y``,` `i``z`)
and similarly for others. If the *Â*’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.