Changes

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