Skip to content

Changes

Page history
version 4.2 authored by Udo Ziegler's avatar Udo Ziegler
Show whitespace changes
Inline Side-by-side
3-NIRVANA-user-guide/3.2-User-interfaces.md
View page @ ea67580c
![logo_short_small](uploads/1e4ee858ca174663b210647640da0976/logo_short_small.png) `4.1` ![logo](uploads/1e4ee858ca174663b210647640da0976/logo_short_small_circle.png) `4.2` | Udo Ziegler, AIP
CONTENTS: CONTENTS:
[Overview](3.2-User-interfaces#overview) [Overview](#overview)
[Specification of main simulation parameters](3.2-User-interfaces#specification-of-main-simulation-parameters) [Specification of main simulation parameters](#specification-of-main-simulation-parameters)
[Defining initial conditions](3.2-User-interfaces#defining-initial-conditions) [Defining initial conditions](#defining-initial-conditions)
[Defining boundary conditions](3.2-User-interfaces#defining-boundary-conditions) [Defining boundary conditions](#defining-boundary-conditions)
[User-defined coefficients for dissipative processes](3.2-User-interfaces#user-defined-coefficients-for-dissipative-processes) [User-defined coefficients for dissipative processes](#user-defined-coefficients-for-dissipative-processes)
[User-defined coefficient for ambipolar diffusion](3.2-User-interfaces#user-defined-coefficient-for-ambipolar-diffusion) [User-defined coefficient for ambipolar diffusion](#user-defined-coefficient-for-ambipolar-diffusion)
[User-defined body force](3.2-User-interfaces#user-defined-body-force) [User-defined body force](#user-defined-body-force)
[User-defined cooling and heating function](3.2-User-interfaces#user-defined-cooling-and-heating-function) [User-defined cooling- and heating function](#user-defined-cooling-and-heating-function)
[User-defined equation of state](3.2-User-interfaces#user-defined-equation-of-state) [User-defined equation of state](#user-defined-equation-of-state)
[User-defined mesh refinement](3.2-User-interfaces#user-defined-mesh-refinement) [User-defined initial mesh refinement and refinement control](#user-defined-initial-mesh-refinement-and-refinement-control)
[User-defined background magnetic field](3.2-User-interfaces#user-defined-background-magnetic-field) [User-defined background magnetic field](#user-defined-background-magnetic-field)
[Specification of NCCM parameters](3.2-User-interfaces#specification-of-nccm-parameters) [User-defined MHD driver](#user-defined-mhddriver)
[User-controllable macros](3.2-User-interfaces#user-controllable-macros) [Specification of NCCM parameters](#specification-of-nccm-parameters)
[User-controllable macros](#user-controllable-macros)
## Overview
### Overview
User interfaces are certain code modules in order to generate and
customize a simulation problem. There are three different types of user User interfaces are code modules in order to generate and customize a
interfaces in NIRVANA: simulation problem. There are three different types of user interfaces
in NIRVANA:
\(1\) User interfaces which demand to add or modify code. These are
recognized by C files ending on `User.c`: **(I)** User interfaces which demand to add or modify code. These are
recognized by files ending on `User.c`:
- `configUser.c` – definition of IC
- `configUser.c` -- definition of IC
- `bcrhoUser.c`,`bcmUser.c`,`bceUser.c`,`bcetUser.c`, - `bcrhoUser.c`,`bcmUser.c`,`bceUser.c`,`bcetUser.c`,
`bcbUser.c`,`bcnXUser.c`,`bcCUser.c`,`bctfUser.c` – user-defined BC `bcbUser.c`,`bcnXUser.c`,`bcCUser.c`,`bctfUser.c` -- user-defined BC
- `viscosityCoeffUser.c`,`conductionCoeffUser.c`,`diffusionCoeffUser.c` - `viscosityCoeffUser.c`,`conductionCoeffUser.c`,`diffusionCoeffUser.c`
– user-defined dissipation coefficients for fluid viscosity, thermal -- user-defined dissipation coefficients for fluid viscosity,
conduction and magnetic diffusion thermal conduction and Ohmic diffusion
- `APdiffusionCoeffUser.c` – user-defined coefficient for ambipolar - `APdiffusionCoeffUser.c` -- user-defined coefficient for ambipolar
diffusion diffusion
- `forceUser.c` – user-defined body force - `forceUser.c` -- user-defined body force
- `sourceCoolingUser.c`,`sourceHeatingUser.c` – user-defined - `sourceCoolingUser.c`,`sourceHeatingUser.c` -- user-defined cooling-
cooling/heating function and heating function
- `eosTabPressureUser.c`,`eosTabTemperatureUser.c` – user-defined - `eosTabPressureUser.c`,`eosTabTemperatureUser.c` -- user-defined
pressure and temperature for a tabulated equation of state pressure and temperature for a tabulated equation of state
- `initDomainUser.c`,`checkDomainUser.c` – user-defined initial- and - `initDomainUser.c`,`checkDomainUser.c` -- user-defined initial mesh
restricted mesh refinement refinement and refinement control
- `sourceB0User.c` - user-defined background magnetic field in the - `sourceB0User.c` -- user-defined background magnetic field in the
B-splitting formalism B-splitting formalism
Defining interfaces requires some familiarity with NIRVANA and, in - `driverMHDUser.c` -- user-defined driver for MHD flows
Coding those interfaces requires some familiarity with NIRVANA and, in
particular, a basic understanding of its data structure. The defintion particular, a basic understanding of its data structure. The defintion
of IC in `configUser.c` is the minimum necessary for a problem setup. of IC in module `configUser.c` is the least necessary for a problem
Sample `...User.c` files can be found in the various testproblems setup. Some examples of `...User.c` files can be found in the
subdirectories under top-level `/nirvana/testproblems`. testproblems subdirectories below `/nirvana/testproblems`.
\(2\) User parameter files (suffix `.par`) which demand editing in order **(II)** User parameter files (ending on `.par`) which demand editing to
to configure simulation properties: configure simulation properties:
- `nirvana.par` – specification of main simulation parameters - `nirvana.par` -- specification of main simulation parameters
- `NCCM.par` – specification of NCCM parameters - `NCCM.par` -- specification of NCCM parameters
Configuration of the main parameter file `nirvana.par` is the minimum Customization of the parameter file `nirvana.par` is the least necessary
necessary for a problem setup. to configure a problem.
\(3\) User header files ending on `User.h`: **(III)** User header files ending on `User.h`:
- `nirvanaUser.h` – header file containing macros for - `nirvanaUser.h` -- header file containing user-contollable macros
user-controllable, secondary simulation parameters
- `User.h` – freely available header file - `User.h` -- freely available header file
## Specification of main simulation parameters ### Specification of main simulation parameters
NIRVANA requires the parameter file `nirvana.par` as input. In this file `nirvana.par` is the place to specify important simulation parameters.
the user specifies main parameters which characterize the properties of The file is divided into several sections grouping parameters into
a simulation. Within `nirvana.par` parameters are grouped into blocks categories:
with category names
- **SIMULATION I/O** – I/O-related parameter - **SIMULATION I/O** -- I/O-related parameter
- **GEOMETRY** – selection of the coordinate system - **GEOMETRY** -- selection of the coordinate system
- **DOMAIN SETTINGS** – physical/numerical domain specifications - **DOMAIN SETTINGS** -- physical/numerical domain specifications
- **BOUNDARY CONDITIONS** – selection of boundary conditions - **BOUNDARY CONDITIONS** -- selection of boundary conditions
- **MESH REFINEMENT** – AMR-related parameters - **MESH REFINEMENT** -- AMR-related parameters
- **USER-SPECIFIC PARAMETERS** – freely usuable parameter - **USER-SPECIFIC PARAMETERS** -- freely usuable parameter
- **SOLVER SPECIFICATIONS** – selection of solver options - **SOLVER SPECIFICATIONS** -- selection of solver options
- **PHYSICS SPECIFICATIONS** – selection of physics options - **PHYSICS SPECIFICATIONS** -- selection of physics options
A parameter section starts with its category name preceded by the A parameter section starts with its category name preceded by a
`>`-character. The associated parameters follow in lines where a line `>`-character. The parameters follow in lines where a line can contain
can contain more than one parameter. Parameters are separated by white more than one parameter. Parameters are separated by white spaces and it
spaces and it can be numbers, single characters or words. Text following can be numbers, single characters or words. Text following the
the `>`-character is ignored by the parser and serves for comments. For `>`-character at the end of a line is comment and is ignored by the
example, the GEOMETRY section in `nirvana.par` looks like parser.
>GEOMETRY ---------------------------------------------------------------------- All parameters are described in detail in the following. On referencing
CART 0.00e+00 0.00e+00 0.00e+00 >geometry:{CART,CYL,SPH},omega[0-2] purposes parameter lines are marked here with a 2-digit running number.
Parameter names set in parentheses follow after. Actually, line numbers
There is just one parameter line here containing 4 parameters: the first and names are not present in `nirvana.par`. Note that almost all
specifies the type of coordinate system being `CART` for Cartesian and parameters in `nirvana.par` are subsumed in a global structure variable
the remaining three parameters give the vector components (doubles) of named `_C` defined in the header file `ctrl.h`. A certain parameter `V`
the angular velocity vector of a rotating frame of reference with in `_C` is then addressed by `_C.V`, e.g., `_C.data_dir`, e.g.,
respect to the inertial frame. It follows a comment after `>`. `_C.geometry` for the type of coordinate system. For some parameters
their range of validity is given either in the form of a set {} of
A complete list of parameters for each category is presented next. possible values or a numeric range.
Parameter lines are marked with a 2-digit running number for a clear
reference (Lines are actually not numbered in `nirvana.par`). Parameter
explanations follow the pattern ’line\# (parameter names within the line
from left to right)’ plus a description of each parameter appearing in a
line. For some parameters their valid range of values are given in round
brackets after the parameter’s name either in form of a set {} or a
numeric range. For some parameters typically used values are given.
Technically, most of the parameters are collected in a global structure
variable named `_C` defined in the header `nirvana.h`. A certain
variable `V` in `_C` is then addressed by `_C.V`, e.g., `_C.geometry`
for the type of coordinate system.
**SIMULATION I/O**: **SIMULATION I/O**:
...@@ -136,59 +126,52 @@ for the type of coordinate system. ...@@ -136,59 +126,52 @@ for the type of coordinate system.
- `mode` ({NEW,RES,MOD,ANA}): running mode of a simulation. Value - `mode` ({NEW,RES,MOD,ANA}): running mode of a simulation. Value
NEW starts a new simulation. Value RES restarts a simulation NEW starts a new simulation. Value RES restarts a simulation
from a snapshot with name specified by the parameter `fname`. from a snapshot with filename given by the parameter `fname`.
Value MOD is like RES but, in addition, allows a modification of Value MOD is like RES but, in addition, allows modification of
read-in simulation data by calling the user interface function the imported data via coding the user interface function
`modifyConfigUser()` before the simulation continues. Value ANA `modifyConfigUser()`. Value ANA reads data from a snapshot and
reads data from a snapshot followed by a call of the user then calls the user interface function `analysisUser()` before
interface function `analysisUser()` and exits thereafter. The the code exits. This mode is intended for data analysis of a
ANA mode is intended for a one-time user-specific analysis of snapshot with all the code functionality available as in a
the snapshot data with full availability of all code regular run. In mode RES, MOD and ANA all remaining parameters
functionalities. In the modes RES, MOD and ANA all parameters in `nirvana.par` following `fname` are ignored.
following `fname` in `nirvana.par` are ignored.
- `fname`: snapshot filename including the directory path relative
- `fname`: qualified name of the snapshot. The qualified name here to the `NIRVANA_RUN` directory as set in `runNIRVANA` (or
is the path relative to the `NIRVANA_RUN` directory as set in `runNIRVANA_PBS`, `runNIRVANA_SLURM`). The base filename has to
the `runNIRVANA` (or `runNIRVANA_MPI.PBS`) script plus the be specified without container id file extension, i.e., `NIR#`
*base* container filename of the snapshot. The base container or `NIRLAST` instead of `NIR#.#` or `NIRLAST.#` (cf. [Snapshot
filename is `NIR#` or `NIRLAST#` with `#`=`_C.mod` but files](#snapshot-files)).
excluding the dot and post-dot container id (cf.
\[Output data\](3.3-Output-data\#snapshot-files)).
- `02` (`_C.mod_max`, `_C.time_max`) - `02` (`_C.mod_max`, `_C.time_max`)
- `_C.mod_max`: maximum number of timestep cycles a simulation - `_C.mod_max`: maximum number of timesteps a simulation should
should run. If the physical evolution time reaches the value run. If the physical evolution time reaches the value
`_C.time_max` before `_C.mod_max` cycles are done the simulation `_C.time_max` before `_C.mod_max` timesteps are done the
stops ahead of schedule. simulation stops ahead of schedule.
- `_C.time_max`: maximum physical evolution time a simulation - `_C.time_max`: maximum physical evolution time a simulation
should run. If `_C.mod_max` cycles are done before the physical should run. If `_C.mod_max` timesteps are done before the
evolution time `_C.time_max` is reached the simulation stops physical evolution time `_C.time_max` is reached the simulation
ahead of schedule. stops ahead of schedule.
- `03` (`_C.freq_log`, `_C.freq_nir`, `_C.freq_ana`, - `03` (`_C.freq_log`, `_C.freq_nir`, `_C.freq_ana`,
`_C.freq_walltime`) `_C.freq_walltime`)
- `_C.freq_log`: interval in units of timestep cycles at which the - `_C.freq_log`: interval in units of timesteps at which the
NIRVANA log file `nirvana.log` (cf. NIRVANA log file `nirvana.log` and monitoring file `nirvana.mon`
[Output data](3.3-Output-data#log-file)) and monitoring file (cf. [Output data](#Output-data)) is updated.
`nirvana.mon` (cf. [Output data](3.3-Output-data#monitoring-file))
is updated.
- `_C.freq_nir`: interval in units of timestep cycles at which - `_C.freq_nir`: interval in units of timesteps at which NIRVANA
NIRVANA snapshots (cf. snapshots (cf. [Snapshot files](#snapshot-files)) are produced.
[Output data](3.3-Output-data#snapshot-files)) are produced.
- `_C.freq_ana`: interval in units of timestep cycles at which the - `_C.freq_ana`: interval in units of timesteps at which the user
user interface function `analysisUser()` is called for interface function `analysisUser()` is called for user-specific
user-specific analysis tasks (cf. analysis tasks (cf. [Data analysis](#data-analysis)).
[Data analysis](3.5-Data-analysis#runtime-data-analysis)).
- `_C.freq_walltime`: interval in seconds at which the special - `_C.freq_walltime`: interval in seconds at which the special
snapshot `NIRLAST.#` (cf. snapshot `NIRLAST.#` (cf. [Snapshot files](#snapshot-files)) is
[Output data](3.3-Output-data#snapshot-files)) is renewed renewed by overwriting existing older files.
overwriting existing older files.
**GEOMETRY**: **GEOMETRY**:
...@@ -202,9 +185,8 @@ for the type of coordinate system. ...@@ -202,9 +185,8 @@ for the type of coordinate system.
- `_C.omega[0-2]`: components of the angular velocity vector of a - `_C.omega[0-2]`: components of the angular velocity vector of a
rotating frame of reference with respect to the inertial frame rotating frame of reference with respect to the inertial frame
(cf. physics guide (cf. [physics
[Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/A-PhysicsGuide), guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
[PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/PhysicsGuide.pdf)).
In case of cylindrical- or spherical geometry only the first In case of cylindrical- or spherical geometry only the first
component `_C.omega[0]` is accepted as relevant representing component `_C.omega[0]` is accepted as relevant representing
rotation around the geometric axis. rotation around the geometric axis.
...@@ -219,35 +201,34 @@ for the type of coordinate system. ...@@ -219,35 +201,34 @@ for the type of coordinate system.
- `01` (`_C.lo[0]`, `_C.up[0]`, `_C.dim[0]`) - `01` (`_C.lo[0]`, `_C.up[0]`, `_C.dim[0]`)
- `_C.lo[0]`,\_C.up\[0\]: lower,upper *x*-coordinate of the - `_C.lo[0]`,\_C.up\[0\]: lower,upper $x$-coordinate of the
computational domain. computational domain.
- `_C.dim[0]`: number of *base-level* grid points in - `_C.dim[0]`: number of *base-level* grid points in
*x*-direction. `_C.dim[0]` must be an integral factor of 4 and $x$-direction. `_C.dim[0]` must be an integral factor of 4, and
excludes ghost cells (which are automatically added). excludes ghost cells which are automatically added by the code.
- `02` (`_C.lo[1]`, `_C.up[1]`, `_C.dim[1]`) - `02` (`_C.lo[1]`, `_C.up[1]`, `_C.dim[1]`)
- `_C.lo[1]`,`_C.up[1]`: lower,upper *y*-coordinate of the - `_C.lo[1]`,`_C.up[1]`: lower,upper $y$-coordinate of the
computational domain. In case of spherical geometry (*y* ≡ *θ*) computational domain. In case of spherical geometry
`_C.lo[1]`,`_C.up[1]` have to be specified in units of *π*. ($y\equiv \theta$) `_C.lo[1]`,`_C.up[1]` have to be specified in
units of $\pi$.
- `_C.dim[1]`: number of *base-level* grid points in - `_C.dim[1]`: number of *base-level* grid points in
*y*-direction. `_C.dim[1]` must be a multiple factor of 4 and $y$-direction. `_C.dim[1]` must be a multiple factor of 4.
excludes ghost cells (which are automatically added).
- `03` (`_C.lo[2]`, `_C.up[2]`, `_C.dim[2]`) - `03` (`_C.lo[2]`, `_C.up[2]`, `_C.dim[2]`)
- `_C.lo[2]`,`_C.up[2]`: lower,upper *z*-coordinate of the - `_C.lo[2]`,`_C.up[2]`: lower,upper $z$-coordinate of the
computational domain. In case of cylindrical- or spherical computational domain. In case of cylindrical- or spherical
geometry (*z* ≡ *ϕ*) `_C.lo[2]`,`_C.up[2]` have to be specified geometry ($z\equiv \phi$) `_C.lo[2]`,`_C.up[2]` have to be
in units of *π*. specified in units of $\pi$.
- `_C.dim[2]`: number of *base-level* grid points in - `_C.dim[2]`: number of *base-level* grid points in
*z*-direction. `_C.dim[2]` must be a multiple factor of 4 and $z$-direction. `_C.dim[2]` must be a multiple factor of 4. If
excludes ghost cells (which are automatically added). If `_C.dim[2]`=0 the simulation is assumed 2D, i.e., axisymmetric
`_C.dim[2]`=0 the simulation is assumed 2D (e.g. axisymmetric in in case of cylindrical- or spherical coordinates.
case of cylindrical- or spherical coordinates).
- `04` (`_C.partitioning_type`, `_C.bnx`, `_C.bny`, `_C.bnz`) - `04` (`_C.partitioning_type`, `_C.bnx`, `_C.bny`, `_C.bnz`)
...@@ -255,16 +236,15 @@ for the type of coordinate system. ...@@ -255,16 +236,15 @@ for the type of coordinate system.
MPI simulations. The value SFC stands for a Space-Filling-Curve MPI simulations. The value SFC stands for a Space-Filling-Curve
domain decomposition and the value BLOCK for a rectangular domain decomposition and the value BLOCK for a rectangular
domain decomposition. Both decomposition types work for unigrid domain decomposition. Both decomposition types work for unigrid
simulations. For AMR simulations a specified value BLOCK is simulations. For AMR simulations BLOCK is ignored and a
ignored and a SFC-decomposition is automatically applied. SFC-decomposition is automatically used instead.
- `_C.bnx`,`_C.bny`,`_C.bnz`: number of domain subdivisions in - `_C.bnx`,`_C.bny`,`_C.bnz`: number of domain subdivisions in
*x*, *y*, *z*-direction in case \_C.partitioning\_type=BLOCK. $x,y,z$-direction in case \_C.partitioning_type=BLOCK. Numbers
Numbers must be chosen such that the grid dimension of must be chosen such that the grid dimension of subdomains is a
subdomains is a multiple factor of 4 in each coordinate multiple factor of 4 in each coordinate direction. Moreover, the
direction. Moreover, the total number of subdomains must equal total number of subdomains must equal the number of MPI threads,
the number of MPI threads, i.e. `_C.bnx`$\times$`_C.bny`$\times$`_C.bnz`=#threads.
`_C.bnx`×`_C.bny`×`_C.bnz`=number of MPI threads.
**BOUNDARY CONDITIONS**: **BOUNDARY CONDITIONS**:
...@@ -275,22 +255,21 @@ for the type of coordinate system. ...@@ -275,22 +255,21 @@ for the type of coordinate system.
- `_C.bc[0-5]` ({I,O,A,M,R,C,P,D,F,U}): selection of boundary - `_C.bc[0-5]` ({I,O,A,M,R,C,P,D,F,U}): selection of boundary
conditions. The type of boundary condition at a physical domain conditions. The type of boundary condition at a physical domain
boundary is characterized by a single uppercase. Uppercases for boundary is characterized by a single capital letter. It are
the domain boundaries are grouped in a 6-letters word with the grouped in a 6-letters word with the individual letter
individual letters representing, from left to right, the representing, from left to right, the lower-$x$ (`_C.bc[0]`),
lower-*x* (`_C.bc[0]`), upper-*x* (`_C.bc[1]`), lower-*y* upper-$x$ (`_C.bc[1]`), lower-$y$ (`_C.bc[2]`), upper-$y$
(`_C.bc[2]`), upper-*y* (`_C.bc[3]`), lower-*z* (`_C.bc[4]`) and (`_C.bc[3]`), lower-$z$ (`_C.bc[4]`) and upper-$z$ (`_C.bc[5]`)
upper-*z* (`_C.bc[5]`) domain boundary. Possible boundary domain boundary. Possible boundary condition types are
condition types are (*u*<sub>*i*/*o*</sub>: inner/outer-domain ($u_{i/o}$: inner/outer-domain boundary values of a variable
boundary values of a variable *u*): $u$):
- `I`: inflow – fluid is allowed to flow inwards through the - `I`: inflow -- fluid is allowed to flow inwards through the
corresponding domain boundary but not to flow outwards. This corresponding domain boundary but not to flow outwards. This
means for the perpendicular velocity component means For the perpendicular velocity component this means
*u*<sub>*o*</sub>=± sgn*u*<sub>*i*</sub>⋅*u*<sub>*i*</sub> $u_o=\pm\operatorname{sgn}u_i\cdot u_i$ where the
where the ’+’(’-’)-sign stands for a lower (upper) domain '+'('-')-sign is valid for a lower (upper) domain boundary.
boundary. Zero derivative, Zero derivative, $u_o=u_i$, is assumed for other
*u*<sub>*o*</sub>=*u*<sub>*i*</sub>, is assumed for other
cell-centroid quantities except the total energy density. cell-centroid quantities except the total energy density.
The total energy density is set consistently under the The total energy density is set consistently under the
assumption of a zero derivative for the thermal energy assumption of a zero derivative for the thermal energy
...@@ -299,61 +278,56 @@ for the type of coordinate system. ...@@ -299,61 +278,56 @@ for the type of coordinate system.
magnetic field component is obtained by a divergence-free magnetic field component is obtained by a divergence-free
extrapolation. extrapolation.
- `O`: outflow – fluid is allowed to flow outwards through the - `O`: outflow -- fluid is allowed to flow outwards through
corresponding domain boundary but not to flow inwards. This the corresponding domain boundary but not to flow inwards.
means for the perpendicular velocity component For the perpendicular velocity component this means
*u*<sub>*o*</sub>=∓sgn *u*<sub>*i*</sub>⋅*u*<sub>*i*</sub> $u_o=\mp\operatorname{sgn}u_i\cdot u_i$ where the
where the ’+’(’-’)-sign stands for a lower (upper) domain '+'('-')-sign is valid for a lower (upper) domain boundary.
boundary. Other variables are likewise to I. Other variables are set likewise to inflow I.
- `M`: mirror symmetry – reflecting conditions, - `M`: mirror symmetry -- reflecting conditions, $u_o=-u_i$,
*u*<sub>*o*</sub>=−*u*<sub>*i*</sub>, for the for the perpendicular components of velocity and magnetic
perpendicular components of velocity and magnetic field and field. Zero derivative for the other variables.
zero derivative for the other variables.
- `A`: antimirror symmetry – same as M for non-magnetic - `A`: anti-mirror symmetry -- same as M for non-magnetic
variables. The magnetic field is forced to have dipole variables. The magnetic field is forced to have dipole
parity with respect to the domain boundary. This means parity with respect to the domain boundary. This means
*u*<sub>*o*</sub>=−*u*<sub>*i*</sub> for the parallel $u_o=-u_i$ for the parallel magnetic field components and
magnetic field components and $u_o=u_i$ for the perpendicular magnetic field component.
*u*<sub>*o*</sub>=*u*<sub>*i*</sub> for the perpendicular
magnetic field component.
- `R`: reflection-on-axis (only of relevance in cylindrical - `R`: reflection-on-axis (only of relevance in cylindrical
geometry for `_C.bc[2]` and in spherical geometry for geometry at the $y$-lower boundary (at $R=0$) or in
`_C.bc[2]`,`_C.bc[3]`) – reflecting conditions at the spherical geometry at the $y$-lower/upper boundaries (at
geometric axis. Same as M except for the azimutal magnetic $\theta=0,\pi$)) -- reflecting conditions at the geometric
field component for which axis. Same as M except for the azimutal magnetic field
*u*<sub>*o*</sub>=−*u*<sub>*i*</sub>. component for which $u_o=-u_i$.
- `C`: reflection-at-center (only of relevance in spherical - `C`: reflection-at-center (only of relevance in spherical
geometry for `_C.bc[0]`) – reflecting conditions at the geometry at the $x$-lower boundary (at $r=0$)) -- reflecting
coordinate origin. Same as M except for the non-radial conditions at the coordinate center. Same as M except for
magnetic field components for which the non-radial magnetic field components for which
*u*<sub>*o*</sub>=−*u*<sub>*i*</sub>. $u_o=-u_i$.
- `P`: periodicity – periodic conditions for all variables. - `P`: periodicity -- periodic conditions for all variables.
- `D`: default – zero derivative in the non-magnetic variables - `D`: default -- zero derivative in the non-magnetic
and perpendicular-to-boundary condition (pseudo-vacuum variables and perpendicular-to-boundary condition
condition) for the magnetic field. (pseudo-vacuum condition) for the magnetic field.
- `F`: free boundary (for cylindrical/spherical geometry - `F`: free boundary (only of relevance in cylindrical
including the geometric axis) – ’natural’ boundary geometry at the $y$-lower boundary (at $R=0$) or in
conditions at the geometric axis. Boundary values are not spherical geometry at the $y$-lower/upper boundaries (at
set explictely but are implicitly given by *π*-shifted $\theta=0,\pi$) and $x$-lower boundary (at $r=0$)) --
values at the geometric axis. Relevant in cylindrical 'natural' boundary condition at the geometric axis. Boundary
geometry for `_C.bc[2]` and in spherical geometry for values are not set explictely but are implicitly given by
`_C.bc[0]`,`_C.bc[2]`,`_C.bc[3]`. $\pi$-shifted values. Note that when a `F`-type boundary
condition is used at a domain boundary not matching the
*Note 1: If the geometric axis is actually not part of the geometrical axis exactly a pseudo-axis is pretented. Note
domain setting F-type BC mimicks a pseudo-axis.* also that restrictions apply in using F-type boundary
conditions in combination with MPI and AMR.
*Note 2: There are restrictions in using F-type BC
concerning MPI and AMR.* - `U`: user-defined boundary conditions (cf. [Defining
boundary conditions](#defining-boundary-conditions)).
- `U`: user-defined boundary conditions (cf.
[Defining BC](3.2-User-interfaces#defining-boundary-conditions)).
**MESH REFINEMENT**: **MESH REFINEMENT**:
...@@ -366,26 +340,27 @@ for the type of coordinate system. ...@@ -366,26 +340,27 @@ for the type of coordinate system.
- `01` (`_C.imr`, `_C.amr`) - `01` (`_C.imr`, `_C.amr`)
- `_C.imr` (≤`MAXLEVEL`): maximal requested refinement level in a - `_C.imr` ($\le$`MAXLEVEL`): maximum refinement level for an
user-defined initially refined mesh (cf. initially refined mesh by the user (cf. [User-defined initial
[User-defined mesh refinement](3.2-User-interfaces#user-defined-mesh-refinement)). mesh refinement and refinement
control](#user-defined-initial-mesh-refinement-and-refinement-control)).
`_C.imr` cannot be larger than the macro `MAXLEVEL` defined in `_C.imr` cannot be larger than the macro `MAXLEVEL` defined in
the header file `nirvanaUser.h`. the header file `nirvanaUser.h`.
- `_C.amr` (≤`MAXLEVEL`): maximum allowed mesh refinement level in - `_C.amr` ($\le$`MAXLEVEL`): allowed maximum mesh refinement
a AMR simulations. level in a AMR simulations. `_C.amr` cannot be larger than the
macro `MAXLEVEL` defined in the header file `nirvanaUser.h`.
- `02` (`_C.amr_eps[0-4]`) - `02` (`_C.amr_eps[0-4]`)
- `_C.amr_eps[0-4]` (typical range: \[0.1, 0.5\]): thresholds for - `_C.amr_eps[0-4]` (typical range: $[0.1,0.5]$): thresholds for
the mass density (`_C.amr_eps[0]`), momentum density the mass density (`_C.amr_eps[0]`), momentum density
(`_C.amr_eps[1]`), energy density (`_C.amr_eps[2]`), magnetic (`_C.amr_eps[1]`), energy density (`_C.amr_eps[2]`), magnetic
field strength (`_C.amr_eps[3]`) and tracer variables field strength (`_C.amr_eps[3]`) and tracer variables
(`_C.amr_eps[4]`) in the derivatives-based mesh refinement (`_C.amr_eps[4]`) in the derivatives-based mesh refinement
criterion (cf. criterion (cf. [AMR](#adaptive-mesh-refinement)). A zero or
[AMR](3.1-Code-basics#adaptive-mesh-refinement)). A zero or negative value indicates that the respective component is not
negative value means that the respective component is disabled used in the mesh refinement check.
for a mesh refinement check.
- `03` (`_C.rhoref`, `_C.vref`, `_C.eref`, `_C.bref`, `_C.Cref`) - `03` (`_C.rhoref`, `_C.vref`, `_C.eref`, `_C.bref`, `_C.Cref`)
...@@ -415,42 +390,45 @@ for the type of coordinate system. ...@@ -415,42 +390,45 @@ for the type of coordinate system.
- `_C.amr_d1` (\[0,1\]): tuning parameter in the derivatives-based - `_C.amr_d1` (\[0,1\]): tuning parameter in the derivatives-based
mesh refinement criterion controlling the transition between a mesh refinement criterion controlling the transition between a
purely gradient-based (`_C.amr_d1=1`) and a purely purely gradient-based criterion (`_C.amr_d1=1`) and a purely
2nd-derivatives-based (`_C.amr_d1=0`) check. 2nd-derivatives-based criterion (`_C.amr_d1=0`).
- `_C.amr_exp` (\[0,2\]): tuning parameter in the - `_C.amr_exp` (\[0,2\]): tuning parameter in the
derivatives-based mesh refinement criterion controlling the derivatives-based mesh refinement criterion controlling the
power-law dependence of mesh refinement on the grid spacing. A power-law dependence of mesh refinement on the grid spacing. A
value `_C.amr_exp>1` (`_C.amr_exp<1`) would mean that mesh value `_C.amr_exp>1` (`_C.amr_exp<1`) means that mesh refinement
refinement becomes progressively more difficult (easier) with becomes progressively more difficult (easier) with increasing
increasing refinement level compared to the standard linear refinement level compared to the standard linear dependence
dependence (`_C.amr_exp=1`). (`_C.amr_exp=1`) (cf. [AMR](#adaptive-mesh-refinement)).
- `05` (`_C.amr_Jeans`, `_C.amr_exp`, `_C.amr_Field`) - `05` (`_C.amr_Jeans`, `_C.amr_exp`, `_C.amr_Field`)
- `_C.amr_Jeans` (typical value: 0.2): threshold in the - `_C.amr_Jeans` (typical value: $0.2$): threshold in the
Jeans-length-based mesh refinement criterion (cf. Jeans-length-based mesh refinement criterion (cf.
[AMR](3.1-Code-basics#adaptive-mesh-refinement)). The value [AMR](#adaptive-mesh-refinement)). The value `_C.amr_Jeans`
`_C.amr_Jeans` defines the fraction of local Jeans length to be defines the fraction of local Jeans length to be resolved by one
resolved by at least one grid cell. A zero or negative value grid cell. A zero or negative value means that the
means that the Jeans-length-based criterion is disabled. Jeans-length-based criterion is disabled.
- `_C.amr_dJeans` (≥0): tuning parameter for the - `_C.amr_dJeans` ($\ge 0$): tuning parameter for the
Jeans-length-based mesh refinement criterion allowing a Jeans-length-based mesh refinement criterion allowing a
systematic reduction of the Jeans threshold with increasing systematic reduction of the Jeans threshold with increasing
refinement level *l* according to the expression `_C.amr_Jeans` refinement level $l$ according to the expression `_C.amr_Jeans`
-*l*\*`_C.amr_dJeans`, i.e., the local Jeans length becomes -$l$\*`_C.amr_dJeans`, i.e., the local Jeans length becomes
gradually higher resolved with incresing *l*. gradually higher resolved with incresing $l$. `_C.amr_dJeans`
must be positiv.
*Note: Tuning of `_C.amr_dJeans` must be such that the actual **Important**: `_C.amr_dJeans` must be chosen with care such
Jeans threshold does not become too small of even negative.* that the actual Jeans threshold never becomes too small or
negative at some refinement level. Otherwise it will trigger
mesh refinement up to the maximum refinement level allowed.
- `_C.amr_Field` (typical value: 0.2): threshold in the - `_C.amr_Field` (typical value: $0.2$): threshold in the
Field-length-based mesh refinement criterion (cf. Field-length-based mesh refinement criterion (cf.
[AMR](3.1-Code-basics#adaptive-mesh-refinement)). The value [AMR](#adaptive-mesh-refinement)). The value `_C.amr_Field`
`_C.amr_Field` defines the fraction of the local Field length to defines the fraction of the local Field length to be resolved by
be resolved by at least one grid cell. A zero or negative value one grid cell. A zero or negative value means that the
means that the Field-length-based criterion is disabled. Field-length-based criterion is disabled.
**USER-SPECIFIC PARAMETERS**: **USER-SPECIFIC PARAMETERS**:
...@@ -462,10 +440,10 @@ for the type of coordinate system. ...@@ -462,10 +440,10 @@ for the type of coordinate system.
05 0.00000e+00 0.00000e+00 0.00000e+00 >param[12-14] / 05 0.00000e+00 0.00000e+00 0.00000e+00 >param[12-14] /
06 0 0 0 0 0 0 0 0 >flag[0-7] / 06 0 0 0 0 0 0 0 0 >flag[0-7] /
- `01` (`_C.param[0-2]`) - `01` (`_C.param[0-2]`)\
`02` (`_C.param[3-5]`) `02` (`_C.param[3-5]`)\
`03` (`_C.param[6-8]`) `03` (`_C.param[6-8]`)\
`04` (`_C.param[9-11]`) `04` (`_C.param[9-11]`)\
`05` (`_C.param[12-14]`) `05` (`_C.param[12-14]`)
- `_C.param[0-14]`: freely usable parameter of `double` type. - `_C.param[0-14]`: freely usable parameter of `double` type.
...@@ -510,15 +488,15 @@ for the type of coordinate system. ...@@ -510,15 +488,15 @@ for the type of coordinate system.
constrained-transport MHD solver. constrained-transport MHD solver.
- CT: flow-upwinded face-to-edge interpolation of the - CT: flow-upwinded face-to-edge interpolation of the
face-centroid electric field fluxes obtained in the FV flux face-centered electric field fluxes obtained in the FV flux
computation part. computation part.
- CCT: 2D upwinding procedure based on a Hamilton-Jacobi-type - CCT: 2D upwinding procedure based on a Hamilton-Jacobi-type
numerical technique (equiv. to genuinely 2D-HLL in the numerical technique (equiv. to genuinely 2D-HLL in the
current implementation). current implementation).
- `_C.mhd_courant` (typical value: ≤ 0.5): CFL number in the - `_C.mhd_courant` (typical value: $\le 0.5$): CFL number in the
expression for the MHD timestep. MHD timestep.
- `02` (`_C.viscosity_solver`, `_C.viscosity_courant`) - `02` (`_C.viscosity_solver`, `_C.viscosity_courant`)
...@@ -526,15 +504,14 @@ for the type of coordinate system. ...@@ -526,15 +504,14 @@ for the type of coordinate system.
the fluid viscosity. the fluid viscosity.
- STD: Runge-Kutta method as selected for the MHD integrator. - STD: Runge-Kutta method as selected for the MHD integrator.
Coupling to the MHD integrator is unsplit. Coupling to MHD integrator is unsplit.
- RKL: second-order, stabilized Runge-Kutta-Legendre method. - RKL: second-order, stabilized Runge-Kutta-Legendre method.
Coupling to the MHD integrator is via Strang-type splitting. Coupling to MHD integrator is via Strang-type splitting.
- `_C.viscosity_courant`: CFL-analog number in the expression for - `_C.viscosity_courant`: CFL-like number in the fluid viscosity
the fluid viscosity timestep. A typical value in case of STD is timestep. A typical value in case of STD is $<0.5$. Values much
&lt;0.5. Values much larger than 1 are possible in the case of larger than 1 are possible in the case of RKL.
RKL.
- `03` (`_C.diffusion_solver`, `_C.diffusion_courant`) - `03` (`_C.diffusion_solver`, `_C.diffusion_courant`)
...@@ -550,24 +527,26 @@ for the type of coordinate system. ...@@ -550,24 +527,26 @@ for the type of coordinate system.
- `06` (`_C.heatloss_max_change`) - `06` (`_C.heatloss_max_change`)
- `_C.heatloss_max_change` (typical value: ≤0.1): maximal - `_C.heatloss_max_change` (typical value: $\le 0.1$): allowed
allowed relative change in temperature due to the heatloss maximal relative change in the temperature due to the heatloss
source term. source term.
- `07` (`_C.reactions_max_changeX`, `_C.reactions_max_changeT`) - `07` (`_C.reactions_max_changeX`, `_C.reactions_max_changeT`)
- `_C.reactions_max_changeX` (typical value: ≤0.1): maximal - `_C.reactions_max_changeX` (typical value: $\le 0.1$): allowed
relative change of species number densities (or total number maximal relative change of species number densities (or total
density) in the time integration of the chemo-thermal rate number density) in the time integration of the chemo-thermal
equations. rate equations.
*Note: The macro `SXN` in module `solveNCCM.c` provides a linear **Important:** The macro `SXN` in module `solveNCCM.c` is a
weight in applying the parameter `_C.reactions_max_change` to linear weight applying `_C.reactions_max_change` to a combined
individual number densities and to the total number density.* criterion based on individual number densities and the total
number density (`SXN`=1 uses individual number densities;
`SXN`=0 uses the total number density).
- `_C.reactions_max_changeT` (typical value: ≤0.1): maximal - `_C.reactions_max_changeT` (typical value: $\le 0.1$): allowed
relative change of temperature in the time integration of the maximal relative change in the temperature in the time
chemo-thermal rate equations. integration of the chemo-thermal rate equations.
- `08` (`_C.dt0_reduce`) - `08` (`_C.dt0_reduce`)
...@@ -592,137 +571,131 @@ for the type of coordinate system. ...@@ -592,137 +571,131 @@ for the type of coordinate system.
- `01` (`_C.mf`, `_C.permeability_rel`) - `01` (`_C.mf`, `_C.permeability_rel`)
- `_C.mf` ({Y,N}): option for activating magnetic fields. The - `_C.mf` ({Y,N}): option to include magnetic fields (MHD). Y
value Y enables the evolution of magnetic fields. The value N enables MHD. N refers to HD where the Lorentz force is dropped
refers to HD where the Lorentz force is dropped from the from the momentum equation and the induction equation is not
momentum equation and the induction equation is not solved. solved.
- `_C.permeability_rel`: relative magnetic permeability - `_C.permeability_rel`: relative magnetic permeability
*μ*<sub>*r**e**l*</sub> $\mu_{rel}$ ($=\mu/\mu_0,\,\mu_0=4\pi\cdot10^{-7}
( = *μ*/*μ*<sub>0</sub>, *μ*<sub>0</sub> = 4*π*⋅10<sup>−7</sup>V⋅m<sup>−1</sup>⋅A<sup>−1</sup>⋅s<sup>−1</sup>). \mathtt{V}\cdot\mathtt{m}^{-1}\cdot\mathtt{A}^{-1}\cdot\mathtt{s}^{-1}$).
*Note: The Gaussian unit system can be mimicked by choosing a **Important:** The Gaussian unit system can be mimicked by
value choosing a value $\mathtt{\_C.permeability\_rel}=10^7$ so that
`_C.permeability_rel`=10<sup>7</sup> the magnetic permeability is $\mu=4\pi$.
so that the magnetic permeability is *μ*=4*π*.*
- `02` (`_C.viscosity`, `_C.viscosity_coeff`) - `02` (`_C.viscosity`, `_C.viscosity_coeff`)
- `_C.viscosity` ({Y,N,U}): option for activating fluid viscosity. - `_C.viscosity` ({Y,N,U}): option for fluid viscosity. Y enables
The value Y enables the fluid viscosity term using a constant the fluid viscosity term using a constant *dynamic* viscosity
*dynamic* viscosity coefficient as given by the parameter coefficient as given by the parameter `_C.viscosity_coeff`
`_C.viscosity_coeff` described next. The value U enables fluid described below. U enables fluid viscosity with a user-defined
viscosity with a user-defined dynamic viscosity coefficient as dynamic viscosity coefficient as coded in the user interface
coded in the user interface `viscosityCoeffUser.c`. The value N function `viscosityCoeffUser()` (cf. [Fluid
disables fluid viscosity. viscosity](#fluid-viscosity)). N disables fluid viscosity.
- `_C.viscosity_coeff`: standard constant *dynamic* viscosity - `_C.viscosity_coeff`: constant *dynamic* viscosity coefficient.
coefficient.
- `03` (`_C.diffusion`, `_C.diffusion_coeff`) - `03` (`_C.diffusion`, `_C.diffusion_coeff`)
- `_C.diffusion` ({Y,N,U}): option for activating magnetic - `_C.diffusion` ({Y,N,U}): option for Ohmic diffusion. Y enables
diffusion. The value Y enables the magnetic diffusion term using Ohmic diffusion using a constant magnetic diffusion coefficient
a constant magnetic diffusion coefficient given by the parameter given by the parameter `_C.diffusion_coeff` described below. U
`_C.diffusion_coeff` described next. The value U enables enables Ohmic diffusion with a user-defined magnetic diffusion
magnetic diffusion with a user-defined magnetic diffusion coefficient as coded in the user interface function
coefficient as coded in the user interface **diffusionCoeffUser()**. (cf. [Ohmic
`diffusionCoeffUser.c`. The value N disables magnetic diffusion](#ohmic-diffusion)). N disables Ohmic diffusion.
diffusion.
- `_C.diffusion_coeff`: standard constant magnetic diffusion - `_C.diffusion_coeff`: constant magnetic diffusion coefficient.
coefficient.
- `04` (`_C.conduction`, `_C.conduction_coeff`, - `04` (`_C.conduction`, `_C.conduction_coeff`,
`_C.conduction_coeff_perp`, `_C.conduction_coeff_sat`) `_C.conduction_coeff_perp`, `_C.conduction_coeff_sat`)
- `_C.conduction` ({Y,N,U,S}): option for activating thermal - `_C.conduction` ({Y,N,U,S}): option for thermal conduction. In
conduction. In the presence of magnetic fields the thermal the presence of magnetic fields the thermal conduction is
conduction is anisotropic with distinct conduction coefficients anisotropic with distinct conduction coefficients parallel and
parallel and perpendicular to the magnetic field. In HD thermal perpendicular to the magnetic field. In HD thermal conduction is
conduction is assumed isotropic. The value Y enables the thermal assumed isotropic. Y enables thermal conduction using a constant
conduction term using a constant conduction coefficient given by conduction coefficient given by the parameter
the parameter `_C.conduction_coeff` (and `_C.conduction_coeff` (and `_C.conduction_coeff_perp` in the
`_C.conduction_coeff_perp` in the anisotropic case) described anisotropic case) described below. U enables thermal conduction
next. The value U enables thermal conduction with a user-defined with a user-defined conduction coefficient as coded in the user
conduction coefficient as coded in the user interface interface function **conductionCoeffUser()** (cf. [Thermal
`conductionCoeffUser.c`. The value S refers to the standard conduction](#thermal-conduction)). S enables thermal conduction
Spitzer conductivity model (cf. physics guide using the standard Spitzer conductivity model predefined in
[Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/A-PhysicsGuide), NIRVANA (cf. [physics
[PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/PhysicsGuide.pdf)). guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
The value N disables thermal conduction. N disables thermal conduction.
- `_C.conduction_coeff`: standard constant thermal conduction - `_C.conduction_coeff`: constant thermal conduction coefficient
coefficient (in the anisotropic case this is the parallel (in the anisotropic case this is the parallel coefficient).
coefficient).
- `_C.conduction_coeff_perp`: thermal conduction coefficient
- `_C.conduction_coeff_perp`: constant thermal conduction perpendicular to the magnetic field (meaningless in isotropic
coefficient perpendicular to the magnetic field (meaningless in conduction).
isotropic conduction).
- `_C.conduction_coeff_sat` (typical value: $0.3$): $\Psi$ is the
- `_C.conduction_coeff_sat` (typical value: 0.3): parameter *Ψ* in parameter in the saturation heat flux model of [@CM77] (cf.
the saturation heat flux model of \[[CM77](3.2-User-interfaces#references)\] (cf. physics guide [physics
[Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/C-NumericsGuide), guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
[PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/NumericsGuide.pdf)).
- `05` (`_C.APdiffusion`, `_C.APdiffusion_coeff`) - `05` (`_C.APdiffusion`, `_C.APdiffusion_coeff`)
- `_C.APdiffusion` ({Y,N,U}): option for activating ambipolar - `_C.APdiffusion` ({Y,N,U}): option for ambipolar diffusion. Y
diffusion. The value Y enables the ambipolar diffusion term enables ambipolar diffusion using a constant ambipolar diffusion
using a constant ambipolar diffusion coefficient given by the coefficient given by the parameter `_C.APdiffusion_coeff`
parameter `_C.APdiffusion_coeff` described next. The value U described below. U enables ambipolar diffusion with a
enables ambipolar diffusion with a user-defined anbipolar user-defined anbipolar diffusion coefficient as coded in the
diffusion coefficient as coded in the user interface user interface function **APdiffusionCoeffUser()** (cf.
`APdiffusionCoeffUser.c`. The value N disables ambipolar [Ambipolar
diffusion. diffusion](#user-defined-coefficient-for-ambipolar-diffusion)).
N disables ambipolar diffusion.
- `_C.APdiffusion_coeff`: standard constant ambipolar diffusion
- `_C.APdiffusion_coeff`: constant ambipolar diffusion
coefficient. coefficient.
- `06` (`_C.energy`, `_C.energy_dual_sw`) - `06` (`_C.energy`, `_C.energy_dual_sw`)
- `_C.energy` ({Y,N,DUAL}): option for activation of the energy - `_C.energy` ({Y,N,DUAL}): option for the energy equation. Y
equation. The value Y means that the total energy equation is means that a total energy equation is solved in step with
solved in time with the gasdynamics equations. For value N no gasdynamics. For N no energy equation is solved. DUAL enables
energy equation is solved. The value DUAL enables a dual energy the dual energy formalism where both the total energy equation
formalism where both the total energy equation and the thermal and the thermal energy equation is solved simultaneously. This
energy equation is solved simultaneously. This allows a more allows a more robust computation of the pressure in regions of
robust computation of the pressure (in view of possitivity) in high Mach number and/or low plasma-$\beta$. In DUAL mode the
regions of high Mach number and/or low plasma-*β*. Actually, the
thermal energy density (hence, the pressure) is computed from thermal energy density (hence, the pressure) is computed from
the solution of the thermal energy equation if the the solution of the thermal energy equation if the
thermal-to-total energy density ratio falls below a certain thermal-to-total energy density ratio falls below a certain
threshold (next parameter). Otherwise the thermal energy density threshold (next parameter). Otherwise the thermal energy density
is computed from the solution of the total energy equation as is computed from the solution of the total energy equation. Note
usual. that the dual energy formalism violates the energy conservation
principle since syncronization between the total- and thermal
*Note 1: Disabling the energy equation is not consistent with energy density is necessary.
the use of an adiabatic EOS.*
*Note 2: The dual energy formalism violates the energy **Important:** Disabling the energy equation is not compatible
conservation principle since syncronization steps between total with the use of an adiabatic EOS, for instance.
and thermal energy densities are necessary to maintain
consistency.*
- `_C.energy_dual_sw` (\[0, 1\]; typical value: 0.01): threshold - `_C.energy_dual_sw` ($[0,1]$; typical value: $0.01$): threshold
value for the thermal-to-total energy density ratio in the dual value for the thermal-to-total energy density ratio in the dual
energy formalism (when `_C.energy`=DUAL). energy formalism (when `_C.energy`=DUAL).
- `07` (`_C.force`) - `07` (`_C.force`)
- `_C.force` ({Y,N}): option for activating a body force term. The - `_C.force` ({Y,N}): option for an external body force term. Y
value Y enables a body force in the momentum equation. The body enables a body force in the momentum equation. The body force
force has to be defined by the user in the user interface has to be defined by the user in the user interface function
`forceUser.c`. The value N disables body forces. `forceUser()` (cf. [User-defined body
force](#user-defined-body-force)). N disables body forces.
- `08` (`_C.gravity`) - `08` (`_C.gravity`)
- `_C.gravity` ({Y,N}): option for activating selfgravity. The - `_C.gravity` ({Y,N,S}): option for selfgravity. Y enables
value Y enables selfgravity. A Poisson equation is solved for selfgravity. S enables selfgravity using the screening mass
the gravitational potential (cf. physics guide approach to solve for Dirichlet boundary conditions. S works
[Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/C-NumericsGuide), only for cylindrical/spherical geometry with $2\pi$ periodicity
[PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/NumericsGuide.pdf)). in azimuthal direction (cf. [physics
The value N disables selfgravity. guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
N disables selfgravity.
- `09` (`_C.eos`,`_C.gamma`,`_C.polytropic_constant`,`_C.temperature`) - `09` (`_C.eos`,`_C.gamma`,`_C.polytropic_constant`,`_C.temperature`)
...@@ -730,33 +703,33 @@ for the type of coordinate system. ...@@ -730,33 +703,33 @@ for the type of coordinate system.
state. Possible values are: state. Possible values are:
- ADIA: adiabatic EOS. This choice requires specificiation of - ADIA: adiabatic EOS. This choice requires specificiation of
the adiabatic index (parameter: `_C.gamma`). the adiabatic index `_C.gamma`.
- ISO: isothermal EOS. This choice requires specification of - ISO: isothermal EOS. This choice requires specification of
the temperature (parameter: `_C.temperature`). the temperature `_C.temperature`.
- POLY: polytropic EOS. This choice requires specification of - POLY: polytropic EOS. This choice requires specification of
the polytropic exponent (parameter: `_C.gamma`) and the polytropic exponent `_C.gamma` and polytropic constant
polytropic constant (parameter: `_C.polytropic_constant`). `_C.polytropic_constant`.
- USER: user-defined analytic EOS. This choice requires the - USER: user-defined analytic EOS. This choice requires the
specification of special macros (cf. specification of special macros (cf. [User-defined
[User-defined EOS](3.2-User-interfaces#user-defined-equation-of-state)) EOS](#user-defined-equation-of-state))
- TAB: user-defined tabulated EOS. This choice requires the - TAB: user-defined tabulated EOS. This choice requires the
generation of look-up tables (cf. generation of look-up tables (cf. [User-defined
[User-defined EOS](3.2-User-interfaces#user-defined-equation-of-state)). EOS](#user-defined-equation-of-state)).
*Note: The choice of an isothermal EOS or polytropic EOS is not **Important:** The choice of an isothermal EOS or polytropic EOS
compatible with solving for an energy equation.* is not compatible with solving for an energy equation.
- `_C.gamma`: adiabatic index (polytropic exponent) in the - `_C.gamma`: adiabatic index (polytropic exponent) in the
adiabatic (polytropic) EOS. adiabatic (polytropic) EOS. Note that in a multi-species gas the
adiabatic index is not a freely selectable parameter but is
approximated by the expression
*Note: In a multi-species gas the adiabatic index is not a $\gamma=\left[5(n_\mathrm{H}+n_\mathrm{He}+n_\mathrm{e})+7n_{\mathrm{H}_2}\right]
freely selectable parameter but is approximated by the /\left[3(n_\mathrm{H}+n_\mathrm{He}+n_\mathrm{e})+5n_{\mathrm{H}_2}\right]$.
expression*
*γ*=\[5(*n*<sub>*H*</sub>+*n*<sub>*He*</sub>+*n*<sub>*e*</sub>)+7*n*<sub>*H*<sub>2</sub></sub>\]/\[3(*n*<sub>*H*</sub>+*n*<sub>*He*</sub>+*n*<sub>*e*</sub>)+5*n*<sub>*H*<sub>2</sub></sub>\]
- `_C.polytropic_constant`: polytropic constant in the polytropic - `_C.polytropic_constant`: polytropic constant in the polytropic
EOS. EOS.
...@@ -765,18 +738,17 @@ for the type of coordinate system. ...@@ -765,18 +738,17 @@ for the type of coordinate system.
- `10` (`_C.heatloss`) - `10` (`_C.heatloss`)
- `_C.heatloss` ({N,ISM,U}): option for activating an external - `_C.heatloss` ({N,ISM,U}): option for an external heatloss
heatloss term (cooling/heating function) in the energy equation. source term (cooling/heating function) in the energy equation. Y
The value N disables the heatloss term. The value ISM enables enables the heatloss term. N disables it. ISM enables the
the heatloss term applying a predefined interstellar medium heatloss term and applies a predefined interstellar medium
cooling/heating function (cf. physics guide cooling/heating function (cf. [physics
[Wiki](https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/A-PhysicsGuide), guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
[PDF](https://gitlab.aip.de/ziegler/NIRVANA/doc/pdf/PhysicsGuide.pdf)). U enables the heatloss term assuming a user-defined
The value U enables the heatloss term assuming a user-defined
cooling/heating function coded in the interfaces cooling/heating function coded in the interfaces
`sourceCoolingUser.c` for cooling and `sourceHeatingUser.c` **sourceCoolingUser()** for cooling and **sourceHeatingUser()**
for heating (cf. for heating (cf. [User-defined
[User-defined cooling/heating function](3.2-User-interfaces#user-defined-cooling-and-heating-function)). cooling/heating](#user-defined-cooling-and-heating-function)).
- `11` (`_C.tracer`) - `11` (`_C.tracer`)
...@@ -785,108 +757,125 @@ for the type of coordinate system. ...@@ -785,108 +757,125 @@ for the type of coordinate system.
- `12` (`_C.testfields`) - `12` (`_C.testfields`)
- `_C.testfields`: number of testfield equations associated with - `_C.testfields`: number of testfields associated with the code
the code functionality TESTFIELDS. functionality TESTFIELDS.
- `13` (`NN`, `NN`, `_C.energy_reactions`, `_C.mean_molecular_weight`, - `13` ( NN, NN, `_C.energy_reactions`, `_C.mean_molecular_weight`,
`_C.mean_ionization`) `_C.mean_ionization`)
- NN ({Y,N}): option for activating the multi-species framework. - NN ({Y,N}): option for a multi-species framework. Y enables a
The value Y enables the multi-species framework. In the multi-species framework. In the multi-species framework an
multi-species framework an advection equation is solved for each advection equation is solved for each gas specie. Species are
specie defined in the `SPECIES` section of the user interface defined in the `SPECIES` section in the parameter file
`NCCM.par`. The value N refers to a standard single-component `NCCM.par`. N refers to using a canonical single-species gas.
gas.
- NN ({Y,N}): option for a chemical reaction network. Y enables
- NN ({Y,N}): option for activating a chemical reaction network chemical reactions between species if a multi-species framework
between species. The value Y enables chemical reactions if the is enabled. In this case an advection-reaction equation is
multi-species framework is enabled. An advection-reaction solved for each gas specie subject to chemical reactions.
equation is solved for each specie defined in the `SPECIES` Chemical reactions are specified in the `REACTIONS` section in
section of the user interface `NCCM.par` subject to chemical `NCCM.par`.
reactions as specified in the `REACTIONS` section of the same
file. - `_C.energy_reactions` ({Y,N}): option for chemo-thermal
processes implemented in the NCCM. Y enables thermal processes
- `_C.energy_reactions` ({Y,N}): option for activating the thermal if both a multi-species framework and chemical reaction network
processes implemented in the NCCM. The value Y enables thermal is enabled. Thermal processes are selected in the
processes if the multi-species framework and chemical reaction `THERMAL PROCESSES` section of `NCCM.par`. N disables
network are both enabled. Thermal processes have to be selected chemo-thermal processes.
in the `THERMAL PROCESSES` section of `NCCM.par`.
- `_C.mean_molecular_weight`: mean molecular weight. - `_C.mean_molecular_weight`: mean molecular weight.
*Note: In a multi-species gas the mean molecular weight is not a **Important:** In a multi-species framework the mean molecular
freely selectable parameter but computed selfconsistently from weight is not a freely selectable parameter but it is computed
the species composition.* selfconsistently from the gas composition.
- `_C.mean_ionization`: mean ionization fraction. - `_C.mean_ionization`: mean ionization fraction.
*Note: In a multi-species gas the ionization fraction is not a **Important:** In a multi-species framework the ionization
freely selectable parameter but computed selfconsistently from fraction is not a freely selectable parameter but is computed
the ionisation structure.* selfconsistently from the ionisation degree of gas species.
## Defining initial conditions ### Defining initial conditions
The user interface `configUser.c` serves to define IC for a problem. The module `configUser.c` serves as user interface for the definition of
Defining IC means assigning values to primary physical variables, initial conditions (IC). To do this the user must program the function
{𝜚, **m**, *e*, **B**, *C*<sub>*c*</sub>, *n*<sub>*s*</sub>}, on the `configUser()` in the module. Defining IC means assigning values to
mesh. However, only problem-relevant variables have to be defined. primary physical variables on the mesh, that are
Primary variables which are irrelevant for the problem need not to be $\{\varrho,\mathbf{m},e,\mathbf{B},C_\mc,n_\ms\}$ -- gas density,
assigned. For example, in simulations with isothermal EOS the total momentum,total energy density, magnetic field, tracer
energy density *e* is irrelevant since no energy equation is solved. ($c=0,{\tt\_C.tracer}-1$), species number densities
Assignment of the magnetic field in HD simulations (when `_C.mf`=0) is ($s=0,{\tt\_C.species}-1$). The parameter `_C.tracer` (`_C.species`)
redundant and even not allowed because arrays for the magnetic field denotes the number of tracer variables (species). However, those primary
components are not allocated in the HD case. Likewise, arrays for variables which are irrelevant for the problem under study need not to
*C*<sub>*c*</sub> and *n*<sub>*s*</sub> are not declared if be assigned. For example, in simulations with isothermal EOS the total
`_C.tracer`=0 respective `_C.species`=0. Furthermore, derived variables energy density $e$ is irrelevant since no energy equation is solved. In
like the temperature must not be assigned here. If the TESTFIELDS HD simulations, the assignment of the magnetic field is redundant, of
infrastructure is used testfield fluctuation variables, course, and would lead to an error because arrays for the magnetic field
**b**<sub>*t*</sub>, *t*=0,*N*<sub>*t*</sub>−1, are to be components are not allocated in the HD case. Likewise, arrays for tracer
considered primary. and species, $C_\mc$ and $n_\ms$, are not declared unless the parameters
`_C.tracer`$>0$ respective ${\tt\_C.species}>0$. If the TESTFIELDS
*Note: In case magnetic field splitting is enabled infrastructure is used (`_C.testfields`$>0$) testfields fluctuation
(cf. [Functionality overview](3.1-Code-basics#functionality-overview)) variables, $\mathbf{b}_\mt,\,\mt=0,{\tt\_C.testfields}-1$, are also
the variable* **B** *represents the residual magnetic field which considered as primary variables and must be assigned by the user. The
must be initialized instead of the total field. Likewise, *e* represents parameter `_C.testfields` denotes the number of testfields.
the residual total energy density.*
**Important:** In case magnetic field splitting is enabled (cf. [Code
The mesh is represented by the master mesh pointer `gm` which is the features](#code-features)) the variable $\mathbf{B}$ represents the
only argument passed to function `configUser()`. The problem-relevant residual magnetic field, i.e., total field minus background field. In
variables have to be assigned for each superblock `g` in `gm`. Recall addition, $e$ represents the residual total energy density, i.e, has a
in this context the explanations about the [mesh data magnetic energy density contribution solely from the residual magnetic
structure](3.1-Code-basics\#mesh-data-structure). field.
There are two types of primary variables: cell-averaged variables and The function
configUser(gm)
is passed the master mesh pointer `gm` which contains all the mesh and
data structures. In serial simulations `gm` represents the global mesh.
In MPI simulations, `gm` only represents the mesh partition on a thread.
Primary variables have to be assigned for each mesh superblock `g` in
`gm`. Looping over the mesh is looping over all refinement levels,
il=0,`_C.level`, and all its superblocks:
for(il=0; il<=_C.level; il++)
for(g=gm[il]; (g); g=g->next)
{
...
}
Please recall in this context the introduction to NIRVANA's [Mesh data
structure](#mesh-data-structure).
There are two types of mesh variables: cell-averaged variables and
face-averaged variables. Cell-averaged variables are 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> $\varrho,m_x,m_y,m_z,n_\ms,C_\mc$. It are described by cell-centroid
and it are represented by cell-centroid coordinates. Face-averaged coordinates. Face-averaged variables are the magnetic field components
variables are the magnetic field components $B_x,B_y,B_z$ (and testfields). It are described by cell face-centroid
*B*<sub>*x*</sub>, *B*<sub>*y*</sub>, *B*<sub>*z*</sub> and it are coordinates.
represented by face-centroid cell coordinates.
Although a superblock `g` contains ghost cells (taking care of boundary
A superblock `g` contains ghost cells added to its active region. values) it is sufficient to only assign cells of the active region.
However, in the definition of IC ghost cells need not to be assigned by Active cells of a superblock start at lower grid index
the user. It is sufficient to assign active cells, i.e. those which are (`g->ixs`,`g->iys`,`g->izs`). The upper index depends on the type of
integrated. Active cells of a superblock start at lower grid index variable, whether cell-averaged or face-averaged, and is given in the
(`g->ixs`,`g->iys`,`g->izs`). The upper index of active cells depends on following table. The table also compiles the coordinates of the
the type of variable and is given in the following table together with variables in a cell:
the representative location of the various variables in a cell:
| variable | upper index | cell coordinate location |
| 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]`) | | 𝜚, *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>*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>*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]`) | | *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 Note that the face-averaged vector components of the magnetic field (and
field have an upper index increased by one in the respective coordinate likewise the testfields components) have an upper index increased by one
direction.* in the respective coordinate direction compared to the cell-averaged
variables. The code parameter `_C.idz` indicates whether a problem is 2D
*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`=0; axissymmetry in cylindrical/spherical geometry) or 3D
(`_C.idz`=1).* (`_C.idz`=1).
Looping over the active grid cells of a superblock looks like this: Here is an example of looping over the active region of a superblock:
/* ASSIGNMENT OF CELL-AVERAGED VARIABLES */ /* ASSIGNMENT OF CELL-AVERAGED VARIABLES */
for(iz=g->izs; iz<=g->ize; iz++) for(iz=g->izs; iz<=g->ize; iz++)
...@@ -895,6 +884,14 @@ Looping over the active grid cells of a superblock looks like this: ...@@ -895,6 +884,14 @@ Looping over the active grid cells of a superblock looks like this:
{ {
g->rho[iz][iy][ix]=... g->rho[iz][iy][ix]=...
... ...
if(_C.tracer) /* SET TRACER VARIABLES */
for(ic=0; ic<_C.tracer; ic++)
g->C[ic][iz][iy][ix]=...
if(_C.species) /* SET SPECIES */
for(is=0; is<_C.species; is++)
g->nX[is][iz][iy][ix]=...
} }
/* ASSIGNMENT OF FACE-AVERAGED MAGNETIC FIELD COMPONENTS */ /* ASSIGNMENT OF FACE-AVERAGED MAGNETIC FIELD COMPONENTS */
...@@ -916,48 +913,50 @@ Looping over the active grid cells of a superblock looks like this: ...@@ -916,48 +913,50 @@ Looping over the active grid cells of a superblock looks like this:
for(ix=g->ixs; ix<=g->ixe; ix++) for(ix=g->ixs; ix<=g->ixe; ix++)
g->bz[iz][iy][ix]=... g->bz[iz][iy][ix]=...
Since tracer variables, *C*<sub>*c*</sub>, are cell-averaged quantities The tracer array index, `ic`, runs from 0 to `_C.tracer`-1 whereas the
like the mass density it can be assigned at the same place. The tracer species array index, `is`, runs from 0 to `_C.species`-1. Tracer
array index, `ic`, runs from 0 to `_C.tracer`-1 where the code parameter variables are assumed dimensionless and it are usually defined in the
`_C.tracer` gives the number of tracers *N*<sub>*c*</sub>. Tracer range $[0,1]$. The mesh of the superblock is indexed by
variables are dimensionless and are usually defined in the range (`ix`,`iy`,`iz`).
\[0,1\].
Similary, species number densities, *n*<sub>*s*</sub>, are cell-averaged **Important:** The species number densities must be consistently defined
quantities. The species array index, `is`, runs from 0 to `_C.species`-1 along with the mass density because of the conservation condition
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 $$\sum_s(\mu_s+\mu_e q_s)n_s = \varrho/m_u,\quad s=0,{\tt \_C.species}-1$$
where *μ*<sub>*s*</sub>, *q*<sub>*s*</sub> and *μ*<sub>*e*</sub> are the where $\mu_s$, $q_s$ and $\mu_\mathrm{e}$ are the molecular weight,
molecular weight and charge of the s-th specie and the electron weight, charge of the s-th specie and the electron weight, respectively. In
respectively. In principle, number densities can be assigned principle, number densities can be assigned individually in this way.
individually on this premise. Alternatively, and much simpler, number However, it may be simpler to specify the species mass fractions in the
densities can be computed in compliance with the user-specified mass parameter file `NCCM.par` and then to call the function
fractions in the parameter file `NCCM.par` by calling the function
CHinitSpecies(gm); CHinitSpecies(gm);
in the *function body* of `configUser.c` after the mass density has been in the *body* (outside of `g`-loops) of `configUser()` after the mass
assigned. Consistency between the species number densities and mass density has been assigned. It guarantees consistency between the species
density is thereby automatically guaranteed. number densities and mass density.
A frequently encounterd problem in the definition of IC for the magnetic **On defining a divergence-free magnetic field.** A frequently
field is to find a divergence-free formulation on the mesh consistent encounterd problem in the definition of IC for the magnetic field is to
with the Constrained-Transport integral approach. This may be a fulfill the divergence contraint ($\nabla\cdot\mathbf{B}=0$) on the mesh
non-trivial task for non-uniform fields, in particular on an initially in a Finite-Volume sense. This may be a non-trivial task, in particular,
refined mesh. There are two clean approaches to achieve a for non-uniform fields. There are two clean approaches to achieve a
divergence-free setup: divergence-free setup:
\(1\) computation of cell-face-averaged magnetic field components via **(I)** Explicit computation of cell-face-averaged magnetic field
exact integration of the analytical expressions. The so discretized components by exact integration of given analytical expressions. For a
field is per se cell-wise divergence-free. For a grid cell grid cell (,,) on superblock `g` the discretized components would then
(`ix`,`iy`,`iz`) on superblock `g` the discritized components read read
$$\mathtt{g->bx[iz][iy][ix]}=\frac{1}{\d \mathcal{A}_x}\int
![ic_b](uploads/411b3e817fec94afa0a2861772697e71/ic_b.png) \limits_\mathtt{g->y[iy]}^\mathtt{g->y[iy+1]}
The numerical expressions for the cell face contents are: \int\limits_\mathtt{g->z[iz]}^\mathtt{g->z[iz+1]} B_x(\mathtt{g->x[ix]},y,z)h_yh_zdydz$$
$$\mathtt{g->by[iz][iy][ix]}=\frac{1}{\d \mathcal{A}_y}\int
\limits_\mathtt{g->x[ix]}^\mathtt{g->x[ix+1]}
\int\limits_\mathtt{g->z[iz]}^\mathtt{g->z[iz+1]} B_y(x,\mathtt{g->y[iy]},z)h_zdxdz$$
$$\mathtt{g->bz[iz][iy][ix]}=\frac{1}{\d \mathcal{A}_z}\int
\limits_\mathtt{g->x[ix]}^\mathtt{g->x[ix+1]}
\int\limits_\mathtt{g->y[iy]}^\mathtt{g->y[iy+1]} B_z(x,y,\mathtt{g->z[iz]})h_ydxdy$$
where the numerical expressions for the cell faces are:
| cell face | expression | | cell face | expression |
|:---------------------|:-----------------------------------------| |:---------------------|:-----------------------------------------|
...@@ -965,68 +964,83 @@ The numerical expressions for the cell face contents are: ...@@ -965,68 +964,83 @@ The numerical expressions for the cell face contents are:
| *δ*𝒜<sub>*y*</sub> | `g->dax[ix]*g->hzyh[iy]*g->dz` | | *δ*𝒜<sub>*y*</sub> | `g->dax[ix]*g->hzyh[iy]*g->dz` |
| *δ*𝒜<sub>*z*</sub> | `g->dax[ix]*g->dy` | | *δ*𝒜<sub>*z*</sub> | `g->dax[ix]*g->dy` |
\(2\) use of the magnetic vector potential **A**, if known, with This way of approach is in most cases not feasible because the 2D
**B** = ∇ × **A**. When discretized in integral form this gives for the integration over cell faces turns out to be too complicated.
face-averaged magnetic field components
**(II)** Use of the magnetic vector potential, $\mathbf{A}$, if known.
The face-averaged magnetic field components are then obtained from the
integral form of $\mathbf{B}=\nabla\times \mathbf{A}$:
`g->bx[iz][iy][ix]`$=\left[h_y\Delta_\miy (h_{zy}\hat{A}_z)
-h_y\Delta_\miz \hat{A}_y\right]/\d \mathcal{A}_x$
`g->by[iz][iy][ix]`$=\left[\Delta_\miz \hat{A}_x
-h_{zy}\Delta_\mix (h_y\hat{A}_z)\right]/\d \mathcal{A}_y$
`g->bx[iz][iy][ix]` = \[*h*<sub>*y*</sub>*Δ*<sub>`iy`</sub>(*h*<sub>*zy*</sub>*Â*<sub>*z*</sub>)−*h*<sub>*y*</sub>*Δ*<sub>`iz`</sub>*Â*<sub>*y*</sub>\]/*δ*𝒜<sub>*x*</sub> `g->bz[iz][iy][ix]`$=\left[\Delta_\mix (h_y\hat{A}_y)
-\Delta_\miy \hat{A}_x\right]/\d \mathcal{A}_z$
`g->by[iz][iy][ix]` = \[*Δ*<sub>`iz`</sub>*Â*<sub>*x*</sub>−*h*<sub>*zy*</sub>*Δ*<sub>`ix`</sub>(*h*<sub>*y*</sub>*Â*<sub>*z*</sub>)\]/*δ*𝒜<sub>*y*</sub> where the difference operators $\Delta_\mix$, $\Delta_\miy$ and
$\Delta_\miz$ are given by
`g->bz[iz][iy][ix]` = \[*Δ*<sub>`ix`</sub>(*h*<sub>*y*</sub>*Â*<sub>*y*</sub>)−*Δ*<sub>`iy`</sub>*Â*<sub>*x*</sub>\]/*δ*𝒜<sub>*z*</sub> $\Delta_{\mix}X=X(\mathtt{ix+1,iy,iz})-X(\mathtt{ix,iy,iz})$,
where (*Â*<sub>*x*</sub>, *Â*<sub>*y*</sub>, *Â*<sub>*z*</sub>) denote $\Delta_{\miy}X=X(\mathtt{ix,iy+1,iz})-X(\mathtt{ix,iy,iz})$,
the path integrals
![ic_A](uploads/d7670a3203d0f3a1e2a9a759e69e2e77/ic_A.png)
over corresponding cell edges. The difference operator
*Δ*<sub>`iy`</sub>, for instance, is given by
*Δ*<sub>`iy`</sub>*Â*<sub>*x*</sub> = *Â*<sub>*x*</sub>(`ix`,`iy`+1,`iz`) − *Â*<sub>*x*</sub>(`ix`,`iy`,`iz`) $\Delta_{\miz}X=X(\mathtt{ix,iy,iz+1})-X(\mathtt{ix,iy,iz})$
and similarly for others. If the *Â*’s are *unambiguously* computed (or and $(\hat{A}_x,\hat{A}_y,\hat{A}_z)$ are the cell-edge integrals
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 $$\hat{A}_x(\mathtt{ix,iy,iz})=\int\limits_\mathtt{g->x[ix]}^\mathtt{g->x[ix+1]}
magnetic field configurations. If inconsistencies exist initially, A_x(x,\mathtt{g->y[iy]},\mathtt{g->z[iz]})dx$$
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 $$\hat{A}_y(\mathtt{ix,iy,iz})=\int\limits_\mathtt{g->y[iy]}^\mathtt{g->y[iy+1]}
testfield fluctuation variables it must be handled in the same way as A_y(\mathtt{g->x[ix]},y,\mathtt{g->z[iz]})dy$$
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 $$\hat{A}_z(\mathtt{ix,iy,iz})=\int\limits_\mathtt{g->z[iz]}^\mathtt{g->z[iz+1]}
serial simulations because of mesh partitioning. The mesh master pointer A_z(\mathtt{g->x[ix]},\mathtt{g->y[iy]},z)dz$$
`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 The $\hat{A}$'s are usually easier to compute than the face-averaged
collection located in subdirectories under `/nirvana/testproblems/`. magnetic field components directly. Even if the $\hat{A}$'s are
Here, the IC for the Orszag-Tang problem (example 1) and a shock-cloud approximated but otherwise set *unambiguously* on cell edges in the grid
interaction problem (example 2) are presented as examples. Both are MHD hierarchy, the magnetic field derived from the above discrete curl of
problems. vector potential has a divergence-free mesh representation.
**Example 1** (taken from `/nirvana/testproblems/MHD/problem2`; cf. \[[Zie04](3.2-User-interfaces#references)\]) Note that NIRVANA can not check for the consistency of a user-defined
magnetic field configuration. If inconsistencies exist initially,
especially on a refined mesh -- even not obvious -- it may lead to
magnetic monopole generation during run time.
IC for the Orszag-Tang problem simulated in a doubly-periodic square When the TESTFIELDS infrastructure is enabled and IC are needed for
domain of length *L* (given by `_C.up[0]-_C.lo[0]`): testfields fluctuation variables it must be handled in the same way as
the magnetic field. In particular, the same issue concerning a
divergence-free representation on the mesh applies.
𝜚 = 𝜚<sub>0</sub> **Important:** For MPI simulations in case the IC require global mesh
**m** = 𝜚<sub>0</sub>(−sin(2*π* *y*/*L*),sin(2*π* *x*/*L*),0) operations parallel coding using MPI functionality is needed.
**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 Many examples of how to define IC can be found in NIRVANA's collection
*B*<sub>0</sub>=1. The adiabatic index *γ*=5/3 and the magnetic of testproblems. It are located in subdirectories below
permeability *μ*=1. `/nirvana/testproblems/`. Here, the definition of IC are shown examplary
for two problems: the Orszag-Tang vortex problem (example 1) and a
shock-cloud interaction problem (example 2).
**Example 1** (taken from `/nirvana/testproblems/MHD/problem2`; cf.
[@Z04])
IC for the Orszag-Tang vortex problem simulated in a doubly-periodic
square domain of length $L$ (given by `_C.up[0]-_C.lo[0]`):
$\varrho = \varrho_0$
$\mathbf{m}=\varrho_0 \left(-\sin (2\pi\,\!y/L),\sin (2\pi\,\!x/L),0\right)$
$\mathbf{B}=B_0 \left(-\sin (2\pi\,\!y/L),\sin (4\pi\,\!x/L),0\right)$
$e=p_0/(\gamma -1)+\mathbf{m}^2/(2\varrho)+\mathbf{B}^2/(2\mu)$
where $\varrho_0=2.77778$, $p_0=5/3$ and $B_0=1$. The adiabatic index
$\gamma =5/3$ and the magnetic permeability $\mu=1$.
void configUser(GRD **gm) void configUser(GRD **gm)
/** /**
...@@ -1092,22 +1106,30 @@ permeability *μ*=1. ...@@ -1092,22 +1106,30 @@ permeability *μ*=1.
return; return;
} }
Example 1 uses the macro `SPQR(a,b,c)` which shortcuts the algebraic The macro `SPQR(a,b,c)` defined in header file `nirvana.h` shortcuts the
expression *a*<sup>2</sup> + *b*<sup>2</sup> + *c*<sup>2</sup>. algebraic expression $a^2+b^2+c^2$.
**Example 2** (taken from `/nirvana/testproblems/MHD/problem17`; cf. \[[Zie05](3.2-User-interfaces#references)\]) **Example 2** (taken from `/nirvana/testproblems/MHD/problem17`; cf.
[@Z05])
IC for the shock-cloud interaction problem simulated in a Cartesian box
given by (*x*, *y*, *z*) ∈ \[ − 1/2, 1/2\]<sup>3</sup>: IC for the shock-cloud interaction problem simulated in a 3D Cartesian
![ic_example2](uploads/e31ebefeb64b91e0ef74d839ebac4de8/ic_example2.png) box given by $(x,y,z)\in [-1/2,1/2]^3$:
At **x**=(0.3,0,0) a spherical clump with radius 0.15 and density of $$(\varrho,p,v_x,v_y,v_z,B_x,B_y,B_z)=\left\{\begin{array}{ll}
10 is embedded and co-moving with its surrounding flow under the (3.86859,167.345,0,0,0,0,2.1826182,-2.1826182) & x<0.1\\
assumption of pressure equilibrium. The adiabatic index *γ*=5/3 and (1,1,-11.2536,0,0,0,0.56418958,0.56418958) & x\ge 0.1
magnetic permeability *μ*=1. \end{array}\right.$$ At ${\mathbf 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 $\gamma =5/3$ and magnetic permeability $\mu=1$. The
grid is initially refined with 3 refinement levels, `_C.level`=3, in the
vicinity of the shock and cloud (cf. [User-defined initial mesh
refinement and refinement
control](#user-defined-initial-mesh-refinement-and-refinement-control)
how to refine a grid initially)
void configUser(GRD **gm) void configUser(GRD **gm)
/** /**
* Definition of ICs. * Definition of IC.
* *
* @param gm master mesh pointer * @param gm master mesh pointer
* *
...@@ -1205,101 +1227,105 @@ magnetic permeability *μ*=1. ...@@ -1205,101 +1227,105 @@ magnetic permeability *μ*=1.
return; return;
} }
## Defining boundary conditions ### Defining boundary conditions
A superblock `g` includes up to 3 ghost cells per side per coordinate 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 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 the physics and geometry of the problem and is stored in the code
parameter `_C.ngc`. 3 ghost cells are needed for MHD on parameter `_C.ngc`. 3 ghost cells are needed for MHD on
cylindrical/spherical grids and in simulations with selfgravity, 2 ghost cylindrical/spherical grids and in simulations with selfgravity, 2 ghost
cells otherwise. cells otherwise. If a superblock shares a physical domain boundary ghost
cells there have to be assigned with boundary values.
There are 3 categories of boundary conditions (BC): standard,
user-defined and those for the gravitational potential in simulations
with selfgravity.
#### Standard BC #### Standard BC
There are a number of predefined BC types (standard BC) which can be Standard BC are predefined and can be selected in the parameter file
selected in the parameter interface `nirvana.par` under category `nirvana.par` below category BOUNDARY CONDITION. Standard types are
BOUNDARY CONDITION. There, standard types are specified by letters specified by letters out of the set {I,O,D,M,A,R,C,F,P,U}. Each domain
{I,O,D,M,A,R,C,F,P} (cf. [Specification of main simulation parameters](3.2-User-interfaces#specification-of-main-simulation-parameters)). boundary is assigned a letter. The possible BC types are described in
detail in [Specification of main simulation
parameters](#specification-of-main-simulation-parameters). Capital U is
special and indicates that user-defined BC are to be used at that domain
boundary.
#### User-defined BC #### User-defined BC
There is further the option that a user can define its own BC. When U specified BC at the corresponding domain boundary are assumed
User-defined BC are enabled by specifying letter U in section user-defined. In that case the code calls user-programmable functions in
`BOUNDARY CONDITION` in `nirvana.par`. For the actual definition of the the following modules:
BC the following user interfaces must be coded, if problem-relevant:
| module | function | BC for |
| BC interface | BC for | |:--------------|:--------------|:---------------------------------------------------------------------------|
|:--------------|:---------------------------------------------------------------------------| | `bcrhoUser.c` | `bcrhoUser()` | mass density 𝜚 |
| `bcrhoUser.c` | mass density 𝜚 | | `bcmUser.c` | `bcmUser()` | momentum density **m** |
| `bcmUser.c` | momentum density **m** | | `bceUser.c` | `bceUser()` | total energy density *e* |
| `bceUser.c` | total energy density *e* | | `bcetUser.c` | `bcetUser()` | thermal energy density *ε* (only when dual energy mode) |
| `bcetUser.c` | thermal energy density *ε* (only when dual energy mode) | | `bcbUser.c` | `bcbUser()` | magnetic field **B** |
| `bcbUser.c` | magnetic field **B** | | `bcnXUser.c` | `bcnXUser()` | species number densities *n*<sub>*s*</sub>, *s* = 0, *N*<sub>*s*</sub> − 1 |
| `bcnXUser.c` | species number densities *n*<sub>*s*</sub>, *s* = 0, *N*<sub>*s*</sub> − 1 | | `bcCUser.c` | `bcCUser()` | tracer variables *C*<sub>*c*</sub>, *c* = 0, *N*<sub>*c*</sub> − 1 |
| `bcCUser.c` | tracer variables *C*<sub>*c*</sub>, *c* = 0, *N*<sub>*c*</sub> − 1 | | `bctfUser.c` | `bctfUser()` | testfield fluctuations **b**<sub>*t*</sub>, *t* = 0, *N*<sub>*t*</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 A combination of user-defined BC and standard BC at different domain
boundaries is possible. Examples can be found in the testproblems boundaries is possible. Examples can be found in NIRVANA's collection of
collection located in subdirectories under `/nirvana/testproblems/`, testproblems located in subdirectories below `/nirvana/testproblems/`,
e.g., the double Mach reflection problem defined in e.g., the double Mach reflection problem defined in
`/nirvana/testproblems/MHD/problem1`. `/nirvana/testproblems/MHD/problem1`. Note that there is no need to
define BC for variables which are not used by the problem.
**Important:** Although not a primary variable, defining BC for the
thermal energy density in `bcetUser()` is mandatory if the dual energy
formalism is enabled (cf. [Specification of main simulation
parameters](#specification-of-main-simulation-parameters)).
At the *x*(*y*,*z*)-lower domain boundary of a superblock `g` boundary Boundary values have to be assigned in ghost cells running from to
values have to be assigned from `ix`=0 to `g->ixs`-1 (`iy`=0 to `g->ixs`-1 ( to `g->iys`-1, to `g->izs`-1) in $x$($y$,$z$)-direction at
`g->iys`-1, `iz`=0 to `g->izs`-1) and full range for the cross indices. the lower domain boundary. At upper domain boundaries the index range
Recall that `g->ixs` (`g->iys`, `g->izs`) counts the first active cell depends on the type of variable as listed in the following table:
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) | | variable | ix-range | iy-range | iz-range |
|:-----------------------------------------------------------------------------------------------------------|:---------------------|:---------------------|:--------------------------| |:-----------------------------------------------------------------------------------------------------------|:---------------------|:---------------------|:--------------------------|
| 𝜚, *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` | | 𝜚, *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>*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>*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 | | *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
BC for the gravitational potential *Φ* cannot be explicitely set. BC for the gravitational potential $\Phi$ cannot be explicitely set
Instead, the following rule is adopted which transforms speficied types except when specifying U in `nirvana.par`. Otherwise, the following rule
into BC types for *Φ*: is adopted to transform MHD BC types into BC types for $\Phi$:
| selected BC type | ⇒ | BC type for *Φ* | | selected BC type | ⇒ | BC type for *Φ* |
|-----------------:|:---:|:----------------| |------------------|:----|:---------------------------------------|
| P | | P | | P | | P (periodicity) |
| M,A,R | | von-Neumann | | M,A,R | | von-Neumann |
| I,O,D | | Dirichlet | | I,O,D | | Dirichlet |
| user-defined (U) | | Dirichlet | | F | | not supported |
| U | | user-defined *Φ* (function `phiUser()` |
In case of von-Neumann conditions (M,A,R) the gradient of potential
vanishes normal to the corresponding domain boundary, i.e.,
$\mathbf{e}_n\cdot \nabla\Phi=0$ where $\mathbf{e}_n$ is the unit normal
vector of the domain boundary. In case of Dirichlet conditions (I,O,D)
boundary values for $\Phi$ are obtained by a multipole expansion for the
given mass distribution (cf. [physics
guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
In case of user-defined conditions U the selfgravity solver calls the
function `phiUser()` expecting a user-defined gravitational potential.
In case of von-Neumann conditions the gradient of potential vanishes **Important:** A combination of periodic BC with non-periodic BC is not
normal to the corresponding domain boundary, i.e., possible. The solver only supports triple-periodic BC when periodic.
**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 When using the screening mass approach for cylindrical/spherical
non-periodic BC. The combination of periodic BC at one domain boundary problems with $2\pi$-periodicity BC are computed to discretization error
with non-periodic BC at another domain boundary is not supported.* accuracy by the gravity solver itself. In this case above rules do not
apply and any specified BC are ignored.
## User-defined coefficients for dissipative processes ### User-defined coefficients for dissipative processes
The modules `viscosityCoeffUser.c`, `conductionCoeffUser.c` and The modules `viscosityCoeffUser.c`, `conductionCoeffUser.c` and
`diffusionCoeffUser.c` serve as templates for a user-defined coefficient `diffusionCoeffUser.c` serve as templates for a user-defined coefficient
...@@ -1309,28 +1335,30 @@ respectively. ...@@ -1309,28 +1335,30 @@ respectively.
#### Fluid viscosity #### Fluid viscosity
Fluid viscosity enters the momentum equation and energy equation via the Fluid viscosity enters the momentum equation and energy equation via the
viscous stress tensor *τ* given by viscous stress tensor $\tau$ given by
*τ* = *ν*\[∇**v**+(∇**v**)<sup>⊤</sup>−2/3(∇⋅**v**)*I*\] $$\tau=\nu \left[\nabla\mathbf{v}+(\nabla\mathbf{v})^\top
-2/3(\nabla\cdot\mathbf{v})\mathrm{I}\right]$$
where *ν* \[kg⋅m<sup>−1</sup>⋅s<sup>−1</sup>\] is the where $\nu$ in units
*dynamic* viscosity coefficient, **v** the fluid velocity and *I* the \[$\mathtt{kg}\cdot\mathtt{m}^{-1}\cdot\mathtt{s}^{-1}$\] is the
identity operator. *dynamic* viscosity coefficient, $\mathbf{v}$ the fluid velocity and
$\mathrm{I}$ the identity operator.
*Note: The *dynamic* viscosity coefficient *ν* to be defined here should **Important:** The viscosity coefficient $\nu$ here should not be
not be confused with the kinetic coefficient,* confused with the kinetic viscosity coefficient,
*ν*<sub>*kinetic*</sub> $\nu_{\mathrm{kinetic}}$ with unit \[m$^2\cdot$s$^{-1}$\], which is
\[m<sup>2</sup>⋅s<sup>−1</sup>\], *related to the dynamic related to the dynamic coefficient by
coefficient by* *ν* = 𝜚*ν*<sub>*kinetic*</sub>. $\nu=\varrho\nu_{\mathrm{kinetic}}$.
In the call of `viscosityCoeffUser()` the function arguments are the A user-defined coefficient, $\nu$, has to be assigned in module
superblock pointer `g` and the array pointer `vis` of type `double***` `viscosityUser.c` in the function
representing the viscosity coefficient *ν*:
viscosityCoeffUser(g,vis); viscosityCoeffUser(g,vis);
The user must assign the `vis`-array for the full index range of `g`, taking arguments `g`, a superblock pointer (not the mesh master
whence a loop like pointer), and the array pointer `vis` of type `double***` representing
the dynamic viscosity coefficient. Here is an example code fragment
/* MY VISCOSITY COEFFICIENT */ /* MY VISCOSITY COEFFICIENT */
for(iz=0; iz<=g->nz; iz++) for(iz=0; iz<=g->nz; iz++)
...@@ -1340,41 +1368,42 @@ whence a loop like ...@@ -1340,41 +1368,42 @@ whence a loop like
vis[iz][iy][ix]=... vis[iz][iy][ix]=...
} }
must be added in `viscosityCoeffUser.c`. The coefficient should be Spatial indexing covers here the full domain of the superblock including
evaluated at cell-centroid coordinates ghost cells. The coefficient is assumed a cell-centroid quantity, i.e.,
(`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`). An example definition for a viscosity it should be defined at location (`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`).
coefficient can be found in testproblem
`/nirvana/testproblems/VISC/problem1`.
User-defined viscosity is enabled by appropriate choice in the parameter Fluid viscosity with user-defined viscosity coefficient is enabled by
interface `nirvana.par` under the category PHYSICS SPECIFICATIONS the appropriate choice of parameter `_C.viscosity` in file `nirvana.par`
(code parameter: `_C.viscosity`). below category PHYSICS SPECIFICATIONS. An example can be found in
testproblem `/nirvana/testproblems/VISC/problem1`.
#### Thermal conduction #### Thermal conduction
Thermal conduction enters the energy equation through a heat flux Thermal conduction enters the energy equation through a heat flux
**F**<sub>*C*</sub>. Generally, the presence of a magnetic field $\mathbf{F}_{\mathrm{C}}$. Generally, the presence of a magnetic field
introduces anisotropic effects with different transport properties along introduces anisotropic effects with different transport properties along
and across the magnetic field. **F**<sub>*C*</sub> is described by and across the magnetic field. $\mathbf{F}_\mathrm{C}$ is described by
**F**<sub>*C*</sub> = −*κ*<sub>∥</sub>(∇*T*⋅**B̂**)**B̂**−*κ*<sub>⊥</sub>(∇*T*−(∇*T*⋅**B̂**)**B̂**) $$\mathbf{F}_\mathrm{C}=-\kappa_\parallel (\nabla T\cdot\mathbf{\hat{B}})\mathbf{\hat{B}}
-\kappa_\perp \left(\nabla T-(\nabla T\cdot\mathbf{\hat{B}})\mathbf{\hat{B}}\right)$$
where *κ*<sub>∥</sub> and *κ*<sub>⊥</sub> is the thermal conduction where $\kappa_\parallel$ and $\kappa_\perp$ is the thermal conduction
coefficient parallel and perpendicular to the magnetic field, coefficient parallel and perpendicular to the magnetic field,
respectively, and **B̂** = **B**/\|**B**\| is the unit vector in the respectively, and $\mathbf{\hat{B}}=\mathbf{B}/|\mathbf{B}|$ is the unit
direction of the magnetic field. *κ*<sub>∥</sub> and *κ*<sub>⊥</sub> are vector in the direction of the magnetic field. $\kappa_\parallel$ and
measured in units $\kappa_\perp$ are measured in units
J⋅K<sup>−1</sup>⋅m<sup>−1</sup>⋅s<sup>−1</sup>. J$\cdot$K$^{-1}\cdot$m$^{-1}\cdot$s$^{-1}$ (cf. [physics
guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
In the call of `conductionCoeffUser()` the function arguments are the User-defined coefficients, $\kappa_\parallel$ and $\kappa_\perp$, have
superblock pointer `g` and the array pointers `cond`, `cond_perp` of to be assigned in module `conductionCoeffUser.c` in the function
type `double***` representing the conduction coefficients
*κ*<sub>∥</sub>, *κ*<sub>⊥</sub>:
conductionCoeffUser(g,cond,cond_perp); conductionCoeffUser(g,cond,cond_perp);
The user must assign the `cond,cond_perp`-arrays for the full index taking arguments `g`, a superblock pointer (not the mesh master
range of `g`, whence a loop like pointer), and the array pointers `cond` and `cond_perp` of type
`double***` representing the parallel and perpendicular conduction
coefficients, respectively. Here is an example code fragment:
/* MY CONDUCTION COEFFICIENT */ /* MY CONDUCTION COEFFICIENT */
for(iz=0; iz<=g->nz; iz++) for(iz=0; iz<=g->nz; iz++)
...@@ -1386,36 +1415,37 @@ range of `g`, whence a loop like ...@@ -1386,36 +1415,37 @@ range of `g`, whence a loop like
cond_perp[iz][iy][ix]=... cond_perp[iz][iy][ix]=...
} }
must be added to `conductionCoeffUser.c`. The coefficients should be Spatial indexing covers here the full domain of the superblock including
evaluated at cell-centroid coordinates ghost cells. Both coefficients are assumed cell-centroid quantities,
(`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`). In HD simulations or MHD simulations i.e., it should be defined at location
with forced isotropy (when `COND_FORCE_ISO`=`YES` in `nirvanaUser.h`) (`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`). In HD simulations and in MHD
only the `cond`-array has to be assigned with `cond` representing the simulations with forced isotropy (when `COND_FORCE_ISO`=`YES` in
conduction coefficient *κ* in the isotropic heat flux `nirvanaUser.h`) only the `cond`-array has to be assigned with `cond`
**F**<sub>*C*</sub> = −*κ*∇*T*. representing now the conduction coefficient $\kappa$ in the isotropic
heat flux $\mathbf{F}_\mathrm{C}=-\kappa \nabla T$.
User-defined thermal conduction is enabled by appropriate choice in the Thermal conduction with user-defined conduction coefficients is enabled
parameter interface `nirvana.par` under the category by the appropriate choice of parameter `_C.conduction` in file
PHYSICS SPECIFICATIONS (code parameter: `_C.conduction`). `nirvana.par` below category PHYSICS SPECIFICATIONS.
#### Ohmic diffusion #### Ohmic diffusion
Ohmic diffusion enters the induction equation and energy equation as a Ohmic diffusion enters the induction equation and energy equation as a
field contribution given by field contribution given by
**E**<sub>*D*</sub> = *η*<sub>*D*</sub>∇×**B** $\mathbf{E}_\mathrm{D}=\eta_\mathrm{D} \nabla\times\mathbf{B}$
where *η*<sub>*D*</sub> \[m<sup>2</sup>⋅s<sup>−1</sup>\] is the where $\eta_\mathrm{D}$ in units \[m$^2\cdot$s$^{-1}$\] is the diffusion
diffusion coefficient. coefficient.
In the call of `diffusionCoeffUser()` function arguments are the A user-defined coefficients, $\eta_\mathrm{D}$, has to be assigned in
superblock pointer `g` and the array pointer `diff` of type `double***` module `diffusionCoeffUser.c` in the function
representing the diffusion coefficient *η*<sub>*D*</sub>:
diffusionCoeffUser(g,diff); diffusionCoeffUser(g,diff);
The user must assign the `diff`-array for the full index range of `g`, taking arguments `g`, a superblock pointer (not the mesh master
whence a loop like pointer), and the array pointer `diff` of type `double***` representing
the diffusion coefficient. Here is an example code fragment:
/* MY DIFFUSION COEFFICIENT */ /* MY DIFFUSION COEFFICIENT */
for(iz=0; iz<=g->nz; iz++) for(iz=0; iz<=g->nz; iz++)
...@@ -1425,39 +1455,39 @@ whence a loop like ...@@ -1425,39 +1455,39 @@ whence a loop like
diff[iz][iy][ix]=... diff[iz][iy][ix]=...
} }
must be added in `diffusionCoeffUser.c`. The coefficient should be Spatial indexing covers here the full domain of the superblock including
evaluated at cell-centroid coordinates ghost cells. The diffusion coefficient is assumed a cell-centroid
quantity, i.e., it should be defined at location
(`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`). (`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`).
User-defined Ohmic diffusion is enabled by appropriate choice in the Ohmic diffusion with user-defined diffusion coefficient is enabled by
parameter interface `nirvana.par` under the category appropriate choice of the parameter `_C.diffusion` in `nirvana.par`
PHYSICS SPECIFICATIONS (code parameter: `_C.diffusion`). below category PHYSICS SPECIFICATIONS.
## User-defined coefficient for ambipolar diffusion ### User-defined coefficient for ambipolar diffusion
The module `APdiffusionCoeffUser.c` serves as template for a The module `APdiffusionCoeffUser.c` serves as template for a
user-defined ambipolar diffusion coefficient. Ambipolar diffusion enters user-defined ambipolar diffusion coefficient.
the induction equation and energy equation as a field contribution given
by
**E**<sub>*AD*</sub> = *η*<sub>*AD*</sub>/*μ* **B**×[(∇×**B**)×**B**\] Ambipolar diffusion enters the induction equation and energy equation as
a field contribution given by
where *η*<sub>*AD*</sub> $\mathbf{E}_\mathrm{AD}=\eta_\mathrm{AD}/\mu\mathbf{B}\times\left[(\nabla\times\mathbf{B})
\[V⋅m⋅A<sup>−1</sup>⋅T<sup>−2</sup>\] denotes the \times\mathbf{B}\right]$
ambipolar diffusion coefficient.
*Note: The prefactor* *η*<sub>*AD*</sub>/*μ* *has units* where $\eta_\mathrm{AD}$ \[V$\cdot$m$\cdot$A$^{-1}\cdot$T$^{-2}$\]
m<sup>2</sup>⋅s<sup>−1</sup>⋅T<sup>−2</sup>. denotes the ambipolar diffusion coefficient. The prefactor
$\eta_\mathrm{AD}/\mu$ has units m$^2\cdot$s$^{-1}\cdot$T$^{-2}$.
In the call of `APdiffusionCoeffUser()` function arguments are the A user-defined coefficients, $\eta_\mathrm{AD}$, has to be assigned in
superblock pointer `g` and the array pointer `APdiff` of type module `APdiffusionCoeffUser.c` in the function
`double***` representing the ambipolar diffusion coefficient
*η*<sub>*AD*</sub>:
APdiffusionCoeffUser(g,APdiff); APdiffusionCoeffUser(g,APdiff);
The user must assign the `APdiff`-array for the full index range of `g`, taking arguments `g`, a superblock pointer (not the mesh master
whence a loop like pointer), and the array pointer `APdiff` of type `double***`
representing the ambipolar diffusion coefficient. Here is an example
code fragment:
/* MY AMBIPOLAR DIFFUSION COEFFICIENT */ /* MY AMBIPOLAR DIFFUSION COEFFICIENT */
for(iz=0; iz<=nz; iz++) for(iz=0; iz<=nz; iz++)
...@@ -1467,31 +1497,30 @@ whence a loop like ...@@ -1467,31 +1497,30 @@ whence a loop like
APdiff[iz][iy][ix]=... APdiff[iz][iy][ix]=...
} }
must be added in `APdiffusionCoeffUser.c`. The coefficient should be Spatial indexing covers here the full domain of the superblock including
evaluated at cell-centroid coordinates ghost cells. The ambipolar diffusion coefficient is assumed a
(`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`). An example definition for an cell-centroid quantity, i.e., it should be defined at location
ambipolar diffusion coefficient can be found in testproblem (`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`).
`/nirvana/testproblems/APDIFF/problem1`.
User-defined ambipolar diffusion is enabled by appropriate choice in the Ambipolar diffusion with user-defined ambipolar diffusion coefficient is
parameter interface `nirvana.par` under the category enabled by appropriate choice of the parameter `_C.APdiffusion` in
PHYSICS SPECIFICATIONS (code parameter: `_C.APdiffusion`). `nirvana.par` below category PHYSICS SPECIFICATIONS. An example can be
found in testproblem `/nirvana/testproblems/APDIFF/problem1`.
## User-defined body force ### User-defined body force
The module `forceUser.c` serves as template for coding a user-defined The module `forceUser.c` serves as template for a user-defined
*specific* body force **f** (force per mass in units *specific* body force $\mathbf{f}$ (force per mass in units
N⋅kg<sup>−1</sup>). The body force then enters the momentum N$\cdot$kg$^{-1}$). The body force enters the momentum equation and
equation and energy equation as source term. energy equation as source term.
In the call of `forceUser()` the function arguments are the superblock A user body force has to be defined in the function
pointer `g` and the array pointers `fx`,`fy`,`fz` of type `double***`
representing the force components:
forceUser(g,fx,fy,fz); forceUser(g,fx,fy,fz);
The user must assign the `fx`,`fy`,`fz`-arrays for active grid cells of taking arguments `g`, a superblock pointer (not the mesh master
`g`, whence a loop like pointer), and the array pointers `fx`,`fy`,`fz` of type `double***`
representing the force components. Here is an example code fragment:
/* MY BODY FORCE */ /* MY BODY FORCE */
for(iz=g->izs; iz<=g->ize; iz++) for(iz=g->izs; iz<=g->ize; iz++)
...@@ -1505,76 +1534,62 @@ The user must assign the `fx`,`fy`,`fz`-arrays for active grid cells of ...@@ -1505,76 +1534,62 @@ The user must assign the `fx`,`fy`,`fz`-arrays for active grid cells of
fz[iz][iy][ix]=... fz[iz][iy][ix]=...
} }
must be added in `forceUser.c`. Strictly speaking, the body force is to Spatial indexing covers the active domain of the superblock. The force
be understood as a cell-averaged quantity. Therefore, **f** should be components are assumed cell-centroid quantities, i.e., it should be
evaluated at cell-centroid coordinates (`g->xc[ix]`,`g->yc[iy]`,`g->zc[iz]`) defined at location (`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 A body force is enabled by specification of the parameter `_C.force` in
`nirvana.par` under the category PHYSICS SPECIFICATIONS (code file `nirvana.par` below category PHYSICS SPECIFICATIONS. An example can
parameter: `_C.force`). be found in testproblem `/nirvana/testproblems/MHD/problem21`.
## User-defined cooling and heating function ### User-defined cooling- and heating function
The modules `sourceCoolingUser.c` and `sourceHeatingUser.c` serve as The modules `sourceCoolingUser.c` and `sourceHeatingUser.c` serve as
templates for coding a user-defined cooling function, templates for coding a user-defined cooling function,
*L*<sub>*cool*</sub>≤0, and heating function, $L_\mathrm{cool}\le 0$, and heating function, $L_\mathrm{heat}\ge 0$,
*L*<sub>*heat*</sub>≥0, respectively. Both functions are allowed respectively. Both functions are allowed to depend on temperature $T$
to depend on temperature *T* and density 𝜚. The net heatloss, i.e. the and density $\varrho$. The net heatloss, i.e. the sum
sum $L_\mathrm{cool}(T,\varrho)+L_\mathrm{heat}(T,\varrho)$, enters as
*L*<sub>*cool*</sub>(*T*,𝜚)+*L*<sub>*heat*</sub>(*T*,𝜚) source term in the energy equation. $L_\mathrm{cool}$ and
of both functions, enters as a source term in the energy equation. $L_\mathrm{heat}$ are measured in units J$\cdot$s$^{-1}\cdot$m$^{-3}$.
*L*<sub>*cool*</sub> and *L*<sub>*heat*</sub> are measured
in units J⋅s<sup>−1</sup>⋅m<sup>−3</sup>. User-defined cooling/heating rates are to be defined in the functions
In the call of `sourceCoolingUser()` (`sourceHeatingUser()`) the f=sourceCoolingUser(T,rho,deriv,dfc)
function arguments are the temperature value `T`, density value `rho`,
the pointer `deriv` to the derivatives flag and the 2-element vector f=sourceHeatingUser(T,rho,deriv,dfh)
`dfc` (`dfh`):
taking arguments `T`, the temperature, `rho`, the gas density, the
f=sourceCoolingUser(T,rho,deriv,dfc); pointer `deriv`, a flag to tell the calling function whether the user
provides the derivatives $\partial L_\mathrm{cool}\partial\varrho$,
f=sourceHeatingUser(T,rho,deriv,dfh); $\partial L_\mathrm{cool}\partial T$ to be stored in the 2-element
vector `dfc[0], dfc[1]` and the derivatives
The user must define the return value $\partial L_\mathrm{heat}\partial\varrho$,
`f` = *L*<sub>*cool*</sub>(`T`,`rho`) $\partial L_\mathrm{heat}\partial T$ to be stored in the 2-element
(*L*<sub>*heat*</sub>(`T`,`rho`)) in `sourceCoolingUser.c` vector `dfh[0], dfh[1]`. The default value is preset to `deriv=YES` in
(`sourceHeatingUser.c`). The `deriv`-flag is thought to indicate the both functions. The user may avoid to explicitely calculating
calling function whether a user provides the derivatives of derivatives by setting `deriv=NO`. Then, derivatives are computed by the
*L*<sub>*cool*</sub>(*T*,𝜚) and code with a finite difference approximation. The functions return `f`,
*L*<sub>*heat*</sub>(*T*,𝜚) with respect to *T* and 𝜚 himself. the user-defined cooling/heating rates.
The `deriv`-flag is preset to `NO` when entering the code function. If
the user sets **Important:** Defining derivatives of the cooling/heating rates is a
prerequisite for the exponential Rosenbrock numerical method, the
*deriv=YES; Field-length-based mesh refinement criterion and for the cooling/heating
timestep computation.
the derivatives must be stored by the user in the 2-element vector `dfc`
(`dfh`) with `dfc[0]` (`dfh[0]`) storing the derivative User-defined cooling/heating is enabled by appropriately choosing the
∂*L*<sub>*cool*</sub>/∂*T* (∂*L*<sub>*heat*</sub>/∂*T*) and parameter `_C.heatloss` in file `nirvana.par` below category PHYSICS
`dfc[1]` (`dfh[1]`) storing the derivative SPECIFICATIONS. An example can be found in testproblem
∂*L*<sub>*cool*</sub>/∂𝜚 (∂*L*<sub>*heat*</sub>/∂𝜚). `/nirvana/testproblems/HEATLOSS/problem2`.
Otherwise derivatives are automatically computed by a finite difference
approximation. ### User-defined equation of state
*Note: Knowledge of the derivatives is a prerequisite in the numerical There are two types of an user-defined EOS: analytic and tabulated.
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 #### Analytic EOS
NIRVANA makes possible a simple user-defined implementation of an NIRVANA allows a simple user-defined implementation of an analytic EOS
analytic EOS. This requires the defintion of the macros by definition of the following macros:
| macro | thermodynamic quantity | | macro | thermodynamic quantity |
|:------------------|:-----------------------| |:------------------|:-----------------------|
...@@ -1583,15 +1598,16 @@ analytic EOS. This requires the defintion of the macros ...@@ -1583,15 +1598,16 @@ analytic EOS. This requires the defintion of the macros
| `TUSR(rho,eth)` | temperature | | `TUSR(rho,eth)` | temperature |
| `ETUSR(rho,eth)` | thermal energy density | | `ETUSR(rho,eth)` | thermal energy density |
as a function of mass density 𝜚 (`rho`) and thermal energy density *ε* The macros are in general functions of the mass density $\varrho$
(`eth`). The macros are declared in the header file `nirvanaUser.h`. The (`rho`) and thermal energy density $\varepsilon$ (`eth`). It are
following code fragment (taken from testproblem declared in the header file `nirvanaUser.h`. In order to give an
`/nirvana/testproblems/GRAVITY/problem4`) illustrates the definition of example, the following code fragment (taken from `nirvanaUser.h` in
a barotropic EOS of the form 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> $p(\varrho)=a^2\varrho \left[1+(\varrho/\varrho_0)^{4/3}\right]^{1/2}$
with *a*=200 and 𝜚<sub>0</sub>=10<sup>−10</sup>. with $a=200$ and $\varrho_0=10^{-10}$.
#define PUSR(rho,eth) (40000.*(rho)*sqrt(1.+pow((rho)/1.e-10,4./3.))+0*(eth)) #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.)) \ #define CS2USR(rho,eth) (40000.*sqrt(1.+pow((rho)/1.e-10,4./3.)) \
...@@ -1600,78 +1616,82 @@ with *a*=200 and 𝜚<sub>0</sub>=10<sup>−10</sup>. ...@@ -1600,78 +1616,82 @@ with *a*=200 and 𝜚<sub>0</sub>=10<sup>−10</sup>.
#define TUSR(rho,eth) (MOK*_C.mean_molecular_weight*PUSR(rho,eth)/rho) #define TUSR(rho,eth) (MOK*_C.mean_molecular_weight*PUSR(rho,eth)/rho)
#define ETUSR(rho,eth) (PUSR(rho,eth)/(_C.gamma-1.)) #define ETUSR(rho,eth) (PUSR(rho,eth)/(_C.gamma-1.))
Since *p* solely depends on 𝜚 the energy density parameter `eth` is Since $p$ solely depends on $\varrho$ the energy density parameter `eth`
irrelevant is this case. The `ETUSR` definition makes use of the prior is irrelevant is this case. The `ETUSR` definition makes use of the
definition of `PUSR` and assumes the usual caloric EOS with ratio of prior definition of `PUSR` and assumes the usual caloric EOS with ratio
specific heats *γ* (code variable: `_C.gamma`). No energy equation is of specific heats $\gamma$ (code variable `_C.gamma`). A barotropic EOS
solved in this testproblem. is consistent with solving no energy equation.
*Note 1: If the energy equation is solved together with an analytic EOS **Important:** If an energy equation is solved together with an analytic
the macro `ETUSR` must read* EOS the macro `ETUSR` must necessarily read
#define ETUSR(rho,eth) (eth) #define ETUSR(rho,eth) (eth)
*for obvious reasons of consistency.* for obvious reasons of consistency.
*Note 2: The macros are used in function `TDeos()` in module `utilTD.c` **Important:** In case of an analytic EOS the macros are used in the
and are needed in the numerical procedure.* function `TDeos()` in module `utilTD.c`, and it are also used at other
places in the numerical solver.
An analytic EOS is enabled by appropriate choice in the parameter An analytic EOS is enabled by appropriately choosing the parameter
interface `nirvana.par` under category PHYSICS SPECIFICATIONS (code `_C.eos` in file `nirvana.par` below category PHYSICS SPECIFICATIONS.
paramter: `_C.eos`).
#### Tabulated EOS #### Tabulated EOS
NIRVANA also permits the implementation of a tabulated EOS. It requires NIRVANA also permits the implementation of a tabulated EOS. However, the
the creation of look-up tables for the logarithm of pressure, log *p*, procedure is more involved and requires the creation of look-up tables
and logarithm of temperature, log *T*, both as functions of the for the logarithm of pressure, $\log p$, and logarithm of temperature,
logarithm of mass density 𝜚 and the logarithm of thermal energy density $\log T$, both as functions of the logarithm of mass density $\varrho$
*ε*. The look-up tables for a given (log 𝜚, log *ε*)-domain are and the logarithm of thermal energy density $\varepsilon$. The look-up
generated by cubic interpolation from user-specified sampling data. The tables for a prescribed $(\log\varrho, \log\varepsilon)$-range are
sampling data {(log *p*)<sub>*ij*</sub>} generated by cubic interpolation from user-specified data
({(log *T*)<sub>*ij*</sub>}) on the grid $\{(\log p)_{ij}\}$ and $\{(\log T)_{ij}\}$ on an rectilinear mesh
{(log 𝜚)<sub>*i*</sub>} × {(log *ε*)<sub>*j*</sub>} has to be defined by $\{(\log\varrho)_i\}\times\{(\log\varepsilon)_j\}$, $i=0,{\tt n1}-1$ and
the user in the interface `eosTabPressureUser.c` $j=0,{\tt n2}-1$ with `n1`$\times$`n2` the number of mesh points.
(`eosTabTemperatureUser.c`). The discretizations in log 𝜚-space and
log *ε*-space, {(log 𝜚)<sub>*i*</sub>; i=0,`n1`-1} and {(log *ε*)<sub>*j*</sub>; j=0,`n2`-1}, must be equidistant Data for the look-up tables has to be defined in the function
where `n1` and `n2` are the number of sampling points. `eosTabPressureUser()` in module `eosTabPressureUser.c` for $\log p$,
and in the function `eosTabTemperatureUser()` in module
Actually, in `eosTabPressureUser()` (`eosTabTemperatureUser()`) the user `eosTabTemperatureUser.c` for $\log T$.
must return a structure of type `UTAB` (declared in header `tabular.h`)
which collects all relevant sampling data. Denoting by `ut` an instance In both, `eosTabPressureUser()` and `eosTabTemperatureUser()`, the user
of `UTAB` the structure elements are accessed by `ut.V`. The following must assign and return a structure of type `UTAB` (declared in header
structure elements `V` must be defined by the user: `tabular.h`) which collects all relevant data to construct the tables.
Denoting by `ut` an instance of `UTAB` the structure is accessed by
`ut.V` with elements `V` summarized in the following table:
| `V` | meaning | | `V` | meaning |
|:----------------|:---------------------------------------------------------------------------------------| |:----------------|:---------------------------------------------------------------------------------------|
| `n1` | number of sampling points in log 𝜚 | | `n1` | number of mesh points in log 𝜚 |
| `n2` | number of sampling points in log *ε* | | `n2` | number of mesh points in log *ε* |
| `u1_lo`,`u1_up` | lower,upper value of log 𝜚-range | | `u1_lo`,`u1_up` | lower,upper value of log 𝜚-range |
| `u2_lo`,`u2_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>*ij*</sub>} ({(log *T*)<sub>*ij*</sub>}) | | `dt2[i][j]` | 2D data array for {(log *p*)<sub>*ij*</sub>} ({(log *T*)<sub>*ij*</sub>}) |
The subsequent example taken from testproblem Here is an example how to code `eosTabPressureUser()` taken from
`/nirvana/testproblem/MHD/problem1B` provides sampling data for testproblem `/nirvana/testproblem/MHD/problem1B` which provides data for
tabulation of the caloric EOS *p*=(*γ* − 1)*ε*. the standard EOS of caloric type, $p=(\gamma -1)\varepsilon$, in
tabulated form:
UTAB eosTabPressureUser(void) UTAB eosTabPressureUser(void)
/** /**
* Template for the definition of user-defined sampling data for ln(pressure) * Template for the definition of user-defined data for ln(pressure) as a
* on 2D [ln(density),ln(thermal_energy_density)]-space. The sample data is * function of (ln(density),ln(thermal_energy_density)). The returned structure
* used for creation of the 2D look-up table _TABLP in a tabulated EOS. * is used to create the look-up table _TABLP.
* *
* @retval ut UTAB structure for creating the 2D look-up table _TABLP * @retval ut UTAB structure for creating the look-up table _TABLP
* *
* @author Udo Ziegler, AIP * @author Udo Ziegler, AIP
* *
* @note The UTAB structure is declared in tabular.h and contains the * @note The UTAB structure is declared in tabular.h and contains the
* following elements which has to be assigned here by the user: * following elements which has to be assigned here by the user:
* ut.n1 = number of equidistant points for ln(density); * ut.n1 = number of equidistant points for ln(density)
* ut.n2 = number of equidistant points for * ut.n2 = number of equidistant points for
* ln(thermal_energy_density); * ln(thermal_energy_density)
* [ut.u1_lo,ut.u1_up] = ln(density) variable range; * [ut.u1_lo,ut.u1_up] = ln(density) range
* [ut.u2_lo,ut.u2_up] = ln(thermal_energy_density) variable range; * [ut.u2_lo,ut.u2_up] = ln(thermal_energy_density) range
* ut.dt2 = 2D array for ln(pressure) sampling data. * ut.dt2 = 2D data array for ln(pressure)
*
* Example implementation: tabulated ideal EOS for testproblem * Example implementation: tabulated ideal EOS for testproblem
* MHD/problem1B. * MHD/problem1B.
*/ */
...@@ -1702,7 +1722,7 @@ tabulation of the caloric EOS *p*=(*γ* − 1)*ε*. ...@@ -1702,7 +1722,7 @@ tabulation of the caloric EOS *p*=(*γ* − 1)*ε*.
for(i2=1; i2<n2; i2++) for(i2=1; i2<n2; i2++)
let[i2]=let[i2-1]+dlet; let[i2]=let[i2-1]+dlet;
/* LN(PRESSURE) SAMPLING DATA FOR ln(p)=ln(_C.gamma-1)+ln(et) */ /* LN(PRESSURE) DATA FOR ln(p)=ln(_C.gamma-1)+ln(et) */
ut.dt2=Array2(n1-1,n2-1); ut.dt2=Array2(n1-1,n2-1);
lc=log(_C.gamma-1.); lc=log(_C.gamma-1.);
...@@ -1715,57 +1735,56 @@ tabulation of the caloric EOS *p*=(*γ* − 1)*ε*. ...@@ -1715,57 +1735,56 @@ tabulation of the caloric EOS *p*=(*γ* − 1)*ε*.
return ut; return ut;
} }
First, the number of sampling points `ut.n1` in (log 𝜚)-direction and First, the number of mesh points, `ut.n1` and `ut.n2`, in
`ut.n2` in (log *ε*)-direction are set followed by the ranges $\log\varrho$-direction and $\log\varepsilon$-direction is set. It
\[`ut.u1_lo`,`ut.u1_up`\] for log 𝜚 and \[`ut.u2_lo`,`ut.u2_up`\] for follows the range \[`ut.u1_lo`,`ut.u1_up`\] for $\log\varrho$ and
log *ε*, respectively. Then, the 2D array `ut.dt2` is allocated and \[`ut.u2_lo`,`ut.u2_up`\] for $\log\varepsilon$, respectively. Then, the
sampling data {(log *p*)<sub>*ij*</sub>} is assigned. The returned 2D array `ut.dt2` is allocated and $\{(\log p)_{ij}\}$ data is assigned.
`ut` is used by the calling function to finally create the public `ut` is returned and used by the code to finally create the look-up
look-up table `_TABLP`. A corresponding user definition for log *T* in table `_TABLP`. Actually, `_TABLP` is a (global) pointer of struct type
`eosTabTemperatureUser.c` provides the input for the public look-up `TAB` (declared in the header file `tabular.h`) and represents the
table `_TABLT`. tabulated EOS $\log p(\log\varrho,\log\varepsilon)$. Correspondingly,
coding `eosTabTemperatureUser()` for $\log T$ would provide the input
*Note: `_TABLP` and `_TABLT` are actually pointers of struct type `TAB` for the look-up table `_TABLT`. `_TABLP` and `_TABLT` are used in the
(declared in the header `tabular.h`). It are used in the function code function `TDeos()` in module `utilTD.c`.
`TDeos()` in module `utilTD.c` and in the numerical procedure.*
A tabulated EOS is enabled by appropriate choice of the parameter
A tabulated EOS is enabled by appropriate choice in the parameter `_C.eos` in the file `nirvana.par` below category PHYSICS
interface `nirvana.par` under category PHYSICS SPECIFICATIONS (code SPECIFICATIONS.
paramter: `_C.eos`).
### User-defined initial mesh refinement and refinement control
## User-defined mesh refinement
For certain problems it may be an advantage to start a simulation with a
For certain problems it may be advantageous to start a simulation with a pre-refined mesh in parts of the computational domain, or to have some
pre-refined mesh in some parts of the computational domain or to retrict control over mesh refinement in certain domain parts during runtime.
mesh refinement to some parts of it during runtime. Such situations can This can be done with help of the functions `initDomainUser()` and
be realize with help of the user interfaces `initDomainUser.c` and `checkDomainUser()` in modules `initDomainUser.c` and
`checkDomainUser.c`. `checkDomainUser.c`.
#### Initial mesh refinement #### Initial mesh refinement
The module `initDomainUser.c` serves as template to initialize a In order to refine the mesh initially in a certain spatial domain the
user-defined refined mesh at simulation start. `initDomainUser()` is user must code the function
called by the mesh refinement algorithm. Function arguments are the
coordinates of a testpoint (`x`,`y`,`z`): level=initDomainUser(x,y,z)
level=initDomainUser(x,y,z); taking arguments (`x`,`y`,`z`) which are coordinates of a testpoint. The
user must check whether the testpoint lies within its specified domain
The user must return the mesh refinement level (`level`) with which the and return the requested refinement level `level` for that domain.
vicinity of the given testpoint should be resolved with the limitation
\|`level`\|≤`_C.imr` where code paramter `_C.imr` is `level` is limited by the constraint
the maximum allowed initial mesh refinement level specified by the user $|\mathtt{level}|\le\mathtt{\_C.imr}$ where `_C.imr` is the maximum
in the parameter interface `nirvana.par` under the category allowed initial mesh refinement level as specified in file `nirvana.par`
MESH REFINEMENT. In practice, the user defines spatial domains in below category MESH REFINEMENT. In practice, the user defines spatial
`initDomainUser.c` via mathematical relations, checks whether the domains in `initDomainUser()` via mathematical relations. Here is an
testpoint is contained and then returns the requested refinement level. example which illustrates the procedure. It refines the vicinity of a
The following code fragment illustrates the procedure for a spherical spherical interface with radius $0.1$ centered at $(x,y,z)=(0.3,0,0)$
interface at radius 0.1 relative to the center with refinement level 4:
(*x*, *y*, *z*)=(0.3, 0, 0) to be resolved with refinement level 4:
int initDomainUser(double x, double y, double z) int initDomainUser(double x, double y, double z)
/** /**
* Template for the specification of user-defined spatial domains for the * Template for the specification of user-defined spatial domains for initial
* purpose of initial mesh refinement. * mesh refinement.
* *
* @param x,y,z spatial coords of testpoint * @param x,y,z spatial coords of testpoint
* *
...@@ -1773,14 +1792,15 @@ interface at radius 0.1 relative to the center ...@@ -1773,14 +1792,15 @@ interface at radius 0.1 relative to the center
* *
* @author Udo Ziegler, AIP * @author Udo Ziegler, AIP
* *
* @note Domains must be defined via mathematical relations. The function * @note Domains within the computational grid are to be defined by the user
* then checks whether a testpoint (x,y,z) lies inside a specified * via mathematical relations. This function simply checks then whether
* domain. If this is the case the user must return the requested * a testpoint (x,y,z) lies inside the specified domain(s). If yes, the
* refinement level for this domain. Negative values (level<0) * user must return the requested refinement level for the domain.
* indicate that the domain should be statically refined (i.e. there * Especially, if negative values (level<0) are returned the domain is
* is no derefinement during runtime). For positive values the refined * assumed to be permanently refined with |level| (i.e., there
* domain is subject to subsequent derefinement according to the regular * is no derefinement during runtime). Otherwise, if level>0, the
* criteria. Multiple domains are allowed. * initially refined domain is subject to derefinement during runtime.
* Multiple domains are possible.
*/ */
{ {
int level=0; int level=0;
...@@ -1798,96 +1818,107 @@ interface at radius 0.1 relative to the center ...@@ -1798,96 +1818,107 @@ interface at radius 0.1 relative to the center
return level; return level;
} }
The coding here realizes a buffer zone approach by placing three hollow Note that a buffer zone is realized by placing three spherical shells of
globes of decreasing thickness around the interface. This ensures that decreasing thickness around the interface. This is needed because it
the interface is safely resolved with the requested degree of ensures that the interface can safely be resolved with the requested
refinement. degree of refinement level. Further examples can be found in the
testproblem `/nirvana/testproblems/MHD/problem16` and, for a more
*Note: The use of buffer zones is strongly recommended in the setup of complicated situation, in the testproblem
highly refined, localized regions since the refinement algorithm does `/nirvana/testproblems/GRAVITY/problem1`.
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
(`level`&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 `level`=0 has no effect on mesh refinement. A
return value `level`=−1 prohibits mesh refinement. A return
value `level`&gt;0 forces mesh refinement up to refinement
level `level` with the limitation `level`≤`_C.amr`
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 **Important:** Working with buffer zones is strongly recommended in the
`/nirvana/testproblems/GRAVITY/problem3`. setup of highly refined, localized regions since the refinement
algorithm does not allow jumps in refinement across two level. A
sufficient buffer zone thus ensures that the refinement hierarchy can be
established smoothly.
An initially refined domain is, by default, subject to derefinement
during runtime according to the standard refinement criteria. However,
it can also be marked as static by returning a negative value
$\mathtt{level}<0$ with $|\mathtt{level}|$ then the requested refinement
level. No derefinement occurs during runtime but, on the other hand,
further dynamic refinement on the initially refined domains is possible.
#### Mesh refinement control
## User-defined background magnetic field The user can control mesh refinement in spatial domains during runtime
by the function
The module `sourceB0User.c` serves as template level=checkDomainUser(x,y,z)
for coding a background magnetic field, **B**<sub>0</sub>, in the context
of the magnetic field splitting formalism.
**B**<sub>0</sub> must be a time-independent, divergence-free potential field, i.e.,
d**B**<sub>0</sub>/dt=0, taking arguments (`x`,`y`,`z`) which are the coordinates of a testpoint
as in `initDomainUser()`. As in `initDomainUser()` the user defines
spatial domains via mathematical relations, checks whether the testpoint
lies inside a certain domain or not, and returns the corresponding
control value `level` for that domain:
∇⋅**B**<sub>0</sub>=0 - $\mathtt{level}=0$: no restriction for mesh refinement
∇×**B**<sub>0</sub>=0. - $\mathtt{level}=-1$: prohibits mesh refinement on that domain
- $\mathtt{level}>0$: forces step-by-step mesh refinement on that
domain up to level `level` subject to the constraint
$\mathtt{level}\le \mathtt{\_C.amr}$ where parameter `_C.amr` is the
maximum allowed mesh refinement level as specified in file
`nirvana.par` below category MESH REFINEMENT.
An example can be found in the testproblem
`/nirvana/testproblems/GRAVITY/problem3`.
### User-defined background magnetic field
The module `sourceB0User.c` containing the function `sourceB0User()`
serves as template for defining a background magnetic field,
$\mathbf{B}_0$, in the context of the magnetic field splitting formalism
(cf. [physics
guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf)).
$\mathbf{B}_0$ must be a time-independent, divergence-free potential
field, i.e.,
$$\partial_t\mathbf{B}_0=0$$,
$$\nabla\!\cdot\!{\mathbf{B}_0}=0$$ and
$$\nabla\!\times\!{\mathbf{B}_0}=0$$.
In the call of `sourceB0User()` In the call of `sourceB0User()`
the function arguments are the pointer to the array `B0` of magnetic field the function arguments are the pointer to the array `B0` of magnetic
components and the coordinates `x,y,z`: field The function
sourceB0User(B0,x,y,z); sourceB0User(B0,x,y,z)
The user must define the components takes arguments `B0`, a pointer to a 3-element vector for the components
`B0[0]`, `B0[1]` and `B0[2]` in the x-,y- and z-direction, respectively, of the magnetic field, and the spatial coordinates `x,y,z` of a point.
as a function of (x,y,z). The user must define the components `B0[0]`=$B0_x$, `B0[1]`=$B0_y$ and
`B0[2]`=$B0_z$, respectively,
An example definition for a background field (magnetic dipole) can be found in An example definition for a background field (a magnetic dipole field)
testproblem `/nirvana/testproblems/MHD/problem30`. can be found in testproblem `/nirvana/testproblems/MHD/problem30`.
The B-splitting scheme is enabled The B-splitting scheme is enabled by setting the macro
by setting the macro `B_FIELD_SPLITTING=YES` in the user interface `B_FIELD_SPLITTING=YES` in the header file `nirvanaUser.h` (cf.
[nirvanaUser.h](3.2-User-interfaces#user-controllable-macros). [User-controllable macros](#user-controllable-macros)).
## Specification of NCCM parameters ### User-defined MHD driver
The parameter file `NCCM.par` serves as user interface to the The function `driverMHDUser()` in module `driverMHDUser.c` is thought to
multi-species framework/NCCM. `NCCM.par` is grouped into the category allow a modification of the MHD flow during the integration cycle.
sections Typically, this is a flow driver, or to realize special conditions
inside the active domain.
- **SPECIES SPECIFICATION** – selection of species. An example is given in testproblem
`/nirvana/testproblems/MHD/problem26`. Here, `driverMHDUser()`
implements a time-dependent velocity driver to initialize the
propagation of Alfven wings.
- **REACTIONS SPECIFICATION** – selection of chemical reactions. ### Specification of NCCM parameters
- **THERMAL PROCESSES SPECIFICATION** – selection of microphysical The parameter file `NCCM.par` serves as user interface to NIRVANA's
thermal processes. multi-species framework and chemistry and cooling model (NCCM).
`NCCM.par` is grouped into the sections
A section starts with the category name preceded by the ’`>`’-character. - **SPECIES SPECIFICATION** -- selection of species.
- **REACTIONS SPECIFICATION** -- selection of chemical reactions.
- **THERMAL PROCESSES SPECIFICATION** -- selection of microphysical
thermal processes.
**SPECIES SPECIFICATION**: **SPECIES SPECIFICATION**:
...@@ -1902,40 +1933,39 @@ A section starts with the category name preceded by the ’`>`’-character. ...@@ -1902,40 +1933,39 @@ A section starts with the category name preceded by the ’`>`’-character.
. . . .
. . . .
The following species are currently contained in the multi-species In the SPECIES SPECIFICATION section the user can select individual
framework: species. The following species are currently contained:
- chemical elements H, D, He, C, N, O, Ne, Mg, Si, Fe and their ions - chemical elements H, D, He, C, N, O, Ne, Mg, Si, Fe and their ions
- negatively charged elements H<sup>−</sup>, C<sup>−</sup>, - negatively charged elements H$^-$, C$^-$, O$^-$
O<sup>−</sup>
- molecules H$_2$, H$_2$+, H$_3^+$, HD, O$_2$, C$_2$, O$_2^+$, CH,
- molecules H<sub>2</sub>, H<sub>2</sub>+, H<sub>3</sub><sup>+</sup>, CH$^+$, CH$_2$, CH$_2^+$, CH$_3^+$, OH, OH$^+$, CO, CO$^+$, H$_2$O,
HD, O<sub>2</sub>, C<sub>2</sub>, O<sub>2</sub><sup>+</sup>, CH, H$_2$O$^+$, H$_3$O$^+$, HCO$^+$ and HOC$^+$
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>, Each species is defined by a unique label `LABEL` (eg., He+, Fe+16,
H<sub>2</sub>O, H<sub>2</sub>O<sup>+</sup>, HOC+), its chemical signature (H\_\_D\_\_He_C\_\_N\_\_O\_\_Ne_Mg_Si_Fe),
H<sub>3</sub>O<sup>+</sup>, HCO<sup>+</sup> and HOC<sup>+</sup> i.e., the chemical elements is made of, its charge `CHARGE` in units of
the elementary charge e, and its initial mass fraction `MASSFRAC`.
Each specie is defined by a unique label (`LABEL`), its chemical `LABEL`, the chemical signature and CHARGE should **NOT** be changed by
signature (H\_\_D\_\_He\_C\_\_N\_\_O\_\_Ne\_Mg\_Si\_Fe), i.e., the the user. Species are selected (unselected) by setting a `+` (`-`) sign
chemical elements a specie consists of, and its charge (`CHARGE`) in in front of the species label. Optionally, the user can assign each
units of the elementary charge e. These basic species properties should species its initial mass fraction which considerably simplifies the
*NOT* be changed here by the user. A specie is selected (unselected) by definition of IC for the species number densities in `configUser.c` (cf.
setting a `+` (`-`) sign in front of the specie label. Optionally, the [Defining initial conditions](#defining-initial-conditions)). The sum of
user can assign each specie its initial mass fraction (`MASSFRAC`) which mass fractions should be 1. If not, the mass fractions are automatically
considerably simplifies the definition of IC for the species number rescaled changing the given value. The species properties are collected
densities in the user interface `configUser.c` (cf. in a substructure `_C.X[is]` below the master control structure `_C`
\[Defining IC\](3.2-User-interfaces\#defining-initial-conditions)). The where `is` denotes the species index. Within the specified procedure, it
species properties are collected in a substructure `_C.X[is]` of the is possible to add further species.
public structure `_C` where `is` denotes the species index.
**Important:** A species 'electron' with label e- is automatically
*Note 1: If the sum of `MASSFRAC` over all species does not amount to generated by the code and not shown in `NCCM.par`.
unity the numbers are proportionally rescaled. This changes their
values.* If the multi-species framework is enabled (section PHYSICS
SPECIFICATIONS in file `nirvana.par`) an advection(-reaction) equation
*Note 2: A specie electron with label e- is automatically generated by is solved for each selected species.
the code.*
**REACTIONS SPECIFICATION**: **REACTIONS SPECIFICATION**:
...@@ -1951,12 +1981,13 @@ the code.* ...@@ -1951,12 +1981,13 @@ the code.*
. . . .
. . . .
The chemical reactions in section REACTIONS SPECIFICATION are grouped In the REACTIONS SPECIFICATIONS section chemical reactions between
into the classes: species are selected. The predefined chemical reactions so far are
grouped into the following categories:
- subnetworks for the ionization structure of each element - subnetworks for the ionization structure of each element
- subnetworks for H<sub>2</sub> formation and HD formation - subnetworks for H$_2$ formation and HD formation
- C-bearing and O-bearing chemical cycles - C-bearing and O-bearing chemical cycles
...@@ -1964,19 +1995,24 @@ into the classes: ...@@ -1964,19 +1995,24 @@ into the classes:
- dust-assisted reactions - dust-assisted reactions
Please consult the [chemistry
guide](https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/ChemistryGuide.pdf)
for more details on the current chemical network of NIRVANA.
Each reaction is characterized by an ID (a running number) followed by Each reaction is characterized by an ID (a running number) followed by
its definition according to the pattern its definition according to the pattern
STCOEFF REACTANT [& STCOEFF REACTANT] > STCOEFF PRODUCT [& STCOEFF PRODUCT] STCOEFF REACTANT [& STCOEFF REACTANT] > STCOEFF PRODUCT [& STCOEFF PRODUCT]
with stoichiometric coefficient STCOEFF, reactants REACTANT and products with stoichiometric coefficients STCOEFF, reactants REACTANT and
PRODUCT. A reaction is selected (unselected) by setting a `+ (-)` sign products PRODUCT. A reaction is selected (unselected) by setting a
in front of the reaction ID. `+ (-)` sign in front of the reaction ID. Within the specified
nomenclature it is possible to add further reactions. Note that the code
*Note: The code parser checks each reaction for charge conservation and parser checks each reaction for charge conservation and mass
mass conservation as well as the existence of the labels of involved conservation, as well as the existence of the involved species (by
species. This rection check can be disabled by setting the macro checking for `LABEL` in the SPECIES SPECIFICATIONS section). Reaction
`CHECK_REACTIONS` to `NO` in `nirvanaUser.h`.* checking can be disabled by setting the macro `CHECK_REACTIONS=NO` in
header file `nirvanaUser.h`.
**THERMAL PROCESSES SPECIFICATION**: **THERMAL PROCESSES SPECIFICATION**:
...@@ -2000,17 +2036,18 @@ species. This rection check can be disabled by setting the macro ...@@ -2000,17 +2036,18 @@ species. This rection check can be disabled by setting the macro
N >H2O: CE COOLING N >H2O: CE COOLING
N >OH: CE COOLING N >OH: CE COOLING
The NCCM currently contains the following microphysical thermal In the THERMAL PROCESSES SPECIFICATION section the user can select
processes: individual thermal process associated with the chemical network. The
NCCM currently includes the following processes:
- Bremsstrahlung for the H, He, O, Ne and Fe ions - Bremsstrahlung for H, He and He ions
- chemical- and atomic line cooling of H, D and He - 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 - fine-structure line cooling of the metals C, N, O, Ne, Mg, Si and Fe
- chemical-, rotovibrational line- and collision-induced cooling of - chemical-, rotovibrational line- and collision-induced cooling of
H<sub>2</sub> H$_2$
- rotovibrational line cooling of HD - rotovibrational line cooling of HD
...@@ -2018,36 +2055,49 @@ processes: ...@@ -2018,36 +2055,49 @@ processes:
- rotovibrational line cooling of CO - rotovibrational line cooling of CO
- rotovibrational line cooling of H<sub>2</sub>O - rotovibrational line cooling of H$_2$O
- rotational line cooling of OH - rotational line cooling of OH
A thermal process is activated (deactivated) by specifying the value Y A thermal process is selected (unselected) by specifying the value Y
(N). Every text in a line following the ’`>`’-character is ignored by (N).
the parser and serves for comments.
## User-controllable macros ### User-controllable macros
Besides the parameters collected in the user interfaces `nirvana.par` The files `nirvana.par` and `NCCM.par` serve for the specification of
and `NCCM.par` there are a number of further parameters which can be primary simulation parameters. controlled by the user. In practice,
controlled by the user. In practice, however, the user rarely has to however, the user rarely has to change them. There are further
change them. These parameters are defined in the header file controllable parameters defined as macros in the header file
`nirvanaUser.h` as macros and it are grouped into the categories: `nirvanaUser.h`. It are grouped into the categories/subcategories:
- NUMERICS-RELATED MACROS - NUMERICS-RELATED MACROS
- FV SCHEME
- MESH REFINEMENT
- GRAVITY/MULTIGRID
- RKL SCHEME
- HEATLOSS SOLVER
- CHEMICAL REACTIONS SOLVER
- PHYSICS-RELATED MACROS - PHYSICS-RELATED MACROS
- NCCM-RELATED MACROS - NCCM-RELATED MACROS
- MACROS RELATED TO ANALYTICAL EOS - ANALYTICAL EOS-RELATED MACROS
The following table compiles a complete list of macros including a short The following table compiles a complete list of macros including a short
description of their meaning and, in some cases, their permissable description of their meaning and, in some cases, their allowed range of
values (given in brackets). values (given in brackets).
**NUMERICS-RELATED MACROS**: **NUMERICS-RELATED MACROS**:
- FV SCHEME
- `SPACE_ORDER` (2): spatial order of the numerical scheme - `SPACE_ORDER` (2): spatial order of the numerical scheme
- `LIM(a)` ({MM(a),MC(a),VL(a)}): type of slope limiter in the - `LIM(a)` ({MM(a),MC(a),VL(a)}): type of slope limiter in the
...@@ -2058,65 +2108,111 @@ values (given in brackets). ...@@ -2058,65 +2108,111 @@ values (given in brackets).
- `PSI` (\[0.5,1\]): parameter in multi-d limiter - `PSI` (\[0.5,1\]): parameter in multi-d limiter
- `MAXLEVEL` (&lt;128): maximum refinement level - `HLLD_PRESSURE_CORR` ({YES,NO}): option to enable/disable the
low Mach pressure correction term in the HLLD Riemann solver
according to [@MM21]
- `BORIS_CORR` ({YES,NO,AUTO}): option to enable/disable the Boris
correction (cf. [Code features](#code-features))
- `BORIS_CRED`: reduced speed of light in Boris correction
- `BORIS_CRED_MIN`: lower limit of reduced speed of light in units
$2|{\bf v}|+c_s$
- `BORIS_AUTO_CA_MAX`: Alfven speed threshold in units
$2|{\bf v}|+c_s$ above which the Boris correction is avtivated
when `BORIS_CORR`=`AUTO`
- `MG_ITMAX`: maximum number of iterations allowed in the multigrid - MESH REFINEMENT
- `MAXLEVEL` ($<128$): maximum refinement level
- `NOBR` ({YES,NO}): option to prevent mesh refinement at domain
boundaries
- GRAVITY/MULTIGRID
- `GR_RES_NORM` (0,1,2): norm to be used in calculating the
residual in the multigrid solver:
- 0: 1-norm (normalized to initial 1-norm)
- 1: relative 1-norm
- 2: relative max-norm
- `MG_ITMAX`: maximum number of multigrid iterations before
termination
- `MG_TOL` (typical: $10^{-6}$): error tolerance in the multigrid
solver solver
- `MG_TYPE` (MULT): used type of multigrid solver - `MG_N_SUBLEVELS` (1,2): number of sublevels (below zero level)
to be used in the multigrid solver
- `MG_TOL` (typical: 10<sup>−6</sup>): error tolerance in the - `MG_SMOOTHER`
multigrid solver ({AUTO,POINT,XLINE,YLINE,ZLINE,XYPLANE,XZPLANE,YZPLANE}):
smoother type
- `RKL_COURANT_EXPL` (&lt;0.5): explicit Courant number used in the - `MG_ORDER` ({LEX,RB}): relaxation order
RKL solver
- `RKL_MAX_COURANT` (typical: &lt;1000): maximum Courant-like number - `MG_BASE_ITMAX`: base level solver maximum iterations
allowed in the RKL solver
- `RKL_DT_LIM`: RKL-related timestep limit as fraction of the - `MG_BASE_RES`: base level residual threshold
dynamical timestep
- `HEATLOSS_TOL` (typical: 10<sup>−5</sup>): relative error - `MG_BASE_ORDER` ({LEX,RB}): base level relaxation order
tolerance in heatloss solver
- `HEATLOSS_ATOL`: absolute error tolerance in heatloss solver - `MG_BASE_SOR` ({YES,NO}): option to enable overrelaxation in the
base level solver
- `REACTIONS_TOL` (typical: 10<sup>−5</sup>): relative error - `MG_BASE_CHEB_ACC` ({YES,NO}): option to enable Chebyshev
tolerance in NCCM solver acceleration in the
- `REACTIONS_ATOL_X`: absolute error tolerance for number densities in - `MG_MAX_ANISO_X` ($\ge 1$): allowed y,z-to-x anisotropy
NCCM solver
- `REACTIONS_ATOL_T`: absolute error tolerance for the temperature in - `MG_MAX_ANISO_Y` ($\ge 1$): allowed z-to-y anisotropy base level
NCCM solver solver
- `HLLD_PRESSURE_CORR` ({YES,NO}): option to enable/disable the low Mach - RKL SCHEME
pressure correction term in the HLLD Riemann solver according to
\[[MM21](3.2-User-interfaces#references)\]
- `BORIS_CORR` ({YES,NO,AUTO}): option to enable/disable the Boris correction - `RKL_COURANT_EXPL` ($<0.5$): explicit Courant number used in the
(cf. [Functionality overview](3.1-Code-basics#functionality-overview)) RKL solver
- `BORIS_CRED`: reduced speed of light in Boris correction - `RKL_MAX_COURANT` (typical: $<1000$): maximum Courant-like
number allowed in the RKL solver
- `BORIS_CRED_MIN`: lower limit of reduced speed of light in units - `RKL_DT_LIM`: RKL-related timestep limit as fraction of the
2|**v**|+c<sub>s</sub> in Boris correction dynamical timestep
- HEATLOSS SOLVER
- `HEATLOSS_TOL` (typical: $10^{-5}$): relative error tolerance in
heatloss solver
- `HEATLOSS_ATOL`: absolute error tolerance in heatloss solver
- CHEMICAL REACTIONS SOLVER
- `REACTIONS_TOL` (typical: $10^{-5}$): relative error tolerance
in NCCM solver
- `REACTIONS_ATOL_X`: absolute error tolerance for number
densities in NCCM solver
- `BORIS_AUTO_CA_MAX`: Alfven speed threshold in units 2|**v**|+c<sub>s</sub> - `REACTIONS_ATOL_T`: absolute error tolerance for the temperature
above which the Boris correction is avtivated when `BORIS_CORR=AUTO` in NCCM solver
**PHYSICS-RELATED MACROS** **PHYSICS-RELATED MACROS**
- `CENTRIFUGAL_FORCE` ({YES,NO}): option to enable/disable the - `CENTRIFUGAL_FORCE` ({YES,NO}): option to include the centrifugal
centrifugal force term force term
- `COND_FORCE_ISO` ({YES,NO}): option to force isotropy in thermal - `COND_FORCE_ISO` ({YES,NO}): option to force isotropy in thermal
conduction conduction
- `B_FIELD_SPLITTING` ({YES,NO}): option to enable/disable the magnetic - `B_FIELD_SPLITTING` ({YES,NO}): option to use the magnetic field
field splitting formalism splitting formalism (cf. [Code features](#code-features))
(cf. [Functionality overview](3.1-Code-basics#functionality-overview))
**NCCM-RELATED MACROS** **NCCM-RELATED MACROS**
...@@ -2126,28 +2222,28 @@ values (given in brackets). ...@@ -2126,28 +2222,28 @@ values (given in brackets).
- `ADUST`: dust abundance relative to Milky Way - `ADUST`: dust abundance relative to Milky Way
- `ISRF_G` (Milky Way: 1.13): Habing strength of interstellar - `ISRF_G`: Habing strength of interstellar radiation field (1.13)
radiation field
- `DBNCCM` ({YES,NO}): NCCM database availability flag - `DBNCCM` ({YES,NO}): NCCM database availability flag
- `CHECK_REACTIONS` ({YES,NO}): option to check chemical reactions - `CHECK_REACTIONS` ({YES,NO}): option to check chemical reactions
with regard to mass conservation, charge conservation and existence with regard to mass conservation, charge conservation and existence
of species labels of involved species names
**MACROS RELATED TO ANALYTICAL EOS** **ANALYTICAL EOS-RELATED MACROS**
- `PUSR(rho,eth)`: analytic expression for the pressure *p* as a - `PUSR(rho,eth)`: user-defined expression for the pressure $p$ as a
function of 𝜚 (`rho`) and *ε* (`eth`) function of $\varrho$ (`rho`) and $\varepsilon$ (`eth`)
- `CS2USR(rho,eth)`: analytic expression for the square of sound speed - `CS2USR(rho,eth)`: user-defined expression for the square of sound
*c*<sub>*s*</sub><sup>2</sup>(𝜚, *ε*) speed $c_s^2(\varrho,\varepsilon)$
- `TUSR(rho,eth)`: analytic expression for the temperature *T*(𝜚, *ε*) - `TUSR(rho,eth)`: user-defined expression for the temperature
$T(\varrho,\varepsilon)$
- `ETUSR(rho,eth)`: analytic expression of the thermal energy density - `ETUSR(rho,eth)`: user-defined expression of the thermal energy
*ε*(𝜚) as a function of 𝜚 in cases without energy equation. density $\varepsilon(\varrho)$ as a function of $\varrho$ in
Otherwise identity. simulations without energy equation. Otherwise identity.
### References ### References
......