Changes
Page history
version 4.2
authored
Oct 28, 2025
by
Udo Ziegler
Show whitespace changes
Inline
Side-by-side
3-NIRVANA-user-guide/3.1-Code-basics.md
View page @
750513d8
`4.
1`
`4.
2`
| Udo Ziegler, AIP
CONTENTS:
[
Functionality overview
](
3.1-Code-basics#functionality-overview
)
[
Dynamic arrays
](
3.1-Code-basics
#dynamic-arrays
)
[
Mesh data structure
](
3.1-Code-basics
#mesh-data-structure
)
[
Adaptive mesh refinement
](
3.1-Code-basics
#adaptive-mesh-refinement
)
[
Parallelism and load balancing
](
3.1-Code-basics#parallelism-and-
load-balancing
)
[
Code features
](
#code-features
)
[
Dynamic arrays
](
#dynamic-arrays
)
[
Mesh data structure
](
#mesh-data-structure
)
[
Adaptive mesh refinement
](
#adaptive-mesh-refinement
)
[
Workload balancing
](
#work
load-balancing
)
##
Functionality overview
##
# Code features
#### Physics
system
#### Physics
and chemistry
NIRVANA is a multi-physics, multi-species, multi-scale
simulation code.
It implements the following physics system:
NIRVANA is a multi-physics, multi-species, multi-scale
Finite-Volume
simulation code.
It implements the following physics system:
-
ideal hydrodynamics (HD) and magnetohydrodynamics (MHD) for
non-relativistic, compressible flows
...
...
@@ -24,42 +24,40 @@ It implements the following physics system:
-
selfgravity
-
body force
-
external
body force
-
external cooling/heating function
-
passive tracer advection
-
multi-species advection-reaction
-
advection of passive tracer
The physical processes and their underlying equations are described in
more detail in the
p
hysics
guide
(
[
PDF
](
https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf
)
).
The
multi-species advection framework
is complemented by
the NIRVANA
Chemistry&Cooling Model (NCCM) – a comprehensive network of chemical
reactions and microphysical processes designed to model the interstellar
medium (ISM). The NCCM is described in the chemistry guide
(
[
PDF
](
https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/ChemistryGuide.pdf
)
).
In full generality the multi-physics, multi-species system is governe
d
by hyperbolic-parabolic-elliptic Partial Differential Equations (PDE)
with source terms, and supplemented by a set of Ordinary Differential
Equations (ODE) representing the chemo-thermal rate equations in the
NCCM. This highly complex system demands for robust numerical schemes
and time integrators for the different mathematical parts of the
equations. The numerical methods applied
in the
NIRVANA code are
documented in the numerics guide
(
[
PDF
](
https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/NumericsGuide.pdf
)
)
.
more detail in the
(
[
P
hysics
guide
](
https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf
)
).
The
physics system
is complemented by
-
a multi-species advection-reaction framework
-
a chemistry & cooling model (NCCM)
The NCCM is a comprehensive network of chemical reactions an
d
microphysical processes orignally designed to model the chemistry of the
interstellar medium. The NCCM is described in detail in the (
[
Chemistry
guide
](
https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/ChemistryGuide.pdf
)
).
In full generality NIRVANA solves a set of hyperbolic-parabolic-elliptic
partial differential equations with source terms decrib
in
g
the
physics
system supplemented by a set of ordinary differential equations for
solving the chemo-thermal rate equations in the NCCM
.
#### Coordinate systems
NIRVANA works for Cartesian, cylindrical and spherical grids in two
space dimensions (2D) and in three space dimensions (3D). Algorithms are
formulated in a
covariant fashion with the vector of metric scale
factors, (
*h*
<sub>
*x*
</sub>
,
*h*
<sub>
*y*
</sub>
,
*h*
<sub>
*z*
</sub>
),
determining the coordinate system with canonical coordinates
(
*x*
,
*y*
,
*z*
) as summarized in the following
table
.
NIRVANA works for Cartesian, cylindrical and spherical grids in two
- or
three space dimensions. The majority of algorithms are formulated in a
covariant fashion with the vector of metric scale
factors,
$(h_x,h_y,h_z)$, determining the geometry. NIRVANA introduces canonical
coordinates $(x,y,z)$ with the meaning as summarized in the following
table
:
| | | | |
|:--------------------------|:------------|:----------------------|:--------------------------------------------------------------------------------------------------------------------|
...
...
@@ -67,92 +65,69 @@ determining the coordinate system with canonical coordinates
| | | (
*x*
,
*y*
,
*z*
) | (
*h*
<sub>
*x*
</sub>
,
*h*
<sub>
*y*
</sub>
(
*x*
),
*h*
<sub>
*z*
</sub>
=
*h*
<sub>
*y*
</sub>
(
*x*
) ⋅
*h*
<sub>
*zy*
</sub>
(
*y*
)) |
| !
[
CART
](
uploads/5f4052ddeef58a9db3b74c9904879e35/CART.png
)
| CARTESIAN | (
*x*
,
*y*
,
*z*
) | (1, 1, 1) |
| !
[
CYL
](
uploads/5a1aea6b07b2ec9aca7f36aac8da0200/CYL.png
)
| CYLINDRICAL | (
*z*
,
*R*
,
*ϕ*
) | (1, 1,
*R*
) |
| !
[
SPH
](
uploads/00d2f413938859f1cdb3cc8cd0c1ca71/SPH.png
)
| SPHERICAL | (
*r*
,
*θ*
,
*ϕ*
) | (1,
*r*
,
*r*
⋅ sin
*θ*
) |
The grid is assumed evenly spaced in canonical coordinates with cell
sizes (or resolutions)
*δx*
<sup>
(0)
</sup>
=
*L*
<sub>
*x*
</sub>
/
*N*
<sub>
*x*
</sub>
,
*δy*
<sup>
(0)
</sup>
=
*L*
<sub>
*y*
</sub>
/
*N*
<sub>
*y*
</sub>
and
*δz*
<sup>
(0)
</sup>
=
*L*
<sub>
*z*
</sub>
/
*N*
<sub>
*z*
</sub>
in the
*x*
-,
*y*
- and
*z*
-direction, respectively.
*L*
<sub>
*x*
</sub>
,
*L*
<sub>
*y*
</sub>
and
*L*
<sub>
*z*
</sub>
are the lengths of the
computation domain and
*N*
<sub>
*x*
</sub>
,
*N*
<sub>
*y*
</sub>
and
*N*
<sub>
*z*
</sub>
are the corresponding numbers of grid cells.
Non-equidistant grid spacings are currently not possible.
| !
[
SPH
](
uploads/00d2f413938859f1cdb3cc8cd0c1ca71/SPH.png
)
| SPHERICAL | (
*r*
,
*θ*
,
*ϕ*
) | (1,
*r*
,
*r*
⋅ sin
*θ*
) |
*
Note
1: In cylindrical geometry the first
canonical coordinate i
s the
coordinate
*z*
along the cylindrical axis.
*
Note
that the
canonical
$x$-
coordinate i
n cylindrical geometry, for
example, is the
coordinate
$z$
along the cylindrical axis.
*
Note 2: Generic 1D simulations are not possible with NIRVANA but can be
realized with a slab-symmetric 2D setup, in principle.
*
**Important:**
NIRVANA does not allow for generic 1D simulations. 1D
problems can be set up in 2D choosing minimal grid size in the invariant
direction.
#### Distributed
memory simulations
#### Distributed
-
memory simulations
NIRVANA allows for distributed-memory (parallel) simulations on compute
clusters based on the Message Passing Interface (MPI). The underlying
concepts of parallelism and load balancing are described in some more
detail
[
here
](
3.1-Code-basics#parallelism-and-load-balancing
)
.
Distributed-memory simulations in NIRVANA are based on the Message
Passing Interface (MPI) on CPUs. See
[
Workload
balancing
](
#workload-balancing
)
for more details.
#### Adaptive mesh
simulations
#### Adaptive mesh
refinement (AMR)
NIRVANA features the possibility of multi-scale simulations by applying
the technique of Adaptive Mesh Refinement (AMR). AMR works in all
offered geometries and in combination with MPI. A description of the
underlying ideas can be found
[
here
](
3.1-Code-basics#adaptive-mesh-refinement
)
.
NIRVANA is an AMR code based on a block-structured approach. AMR works
for all three coordinate systems and in distributed-memory simulations.
See
[
AMR
](
#adaptive-mesh-refinement
)
for more details.
#### The TESTFIELDS infrastructure
NIRVANA supports a special functionality called TESTFIELDS. TESTFIELDS
provides the AMR- and MPI infrastructure for the additional solution of
a number of induction-like equations as they appear in the so-called
testfield method. For that purpose TESTFIELDS allocates arrays for
necessary testfield variables. Hereby, a testfield variable consists of
a pair of vectors. The first vector describes the actual testfield
fluctuations which obey an induction equation subject to the
solenoidality condition. The components of the fluctuation vector are
assumed face-centroid analog to the magnetic field components. The
second vector is for the corresponding electromotive force (EMF) under
the curl operator in the induction equation. The EMF components are
assumed edge-centered analog to the electric field.
*
Note: The TESTFIELDS infrastructure does not implement any numerical
method to solve induction equations for testfield fluctuations.
*
NIRVANA supports a very special functionality called TESTFIELDS.
TESTFIELDS provides the AMR- and MPI infrastructure for the additional
solution of a number of induction-like equations as they appear, for
instance, in the so-called testfield method to compute mean-field dynamo
coefficients. TESTFIELDS allocates arrays for necessary testfield
variables. A pair of testfield variables consists of the vector of
testfield fluctuations and the associated vector of electromotive force.
TESTFIELDS does
*NOT*
implement numerical schemes to solve the
underlying testfield equations.
#### The Boris correction
NIRVANA supports the so-called Boris correction for reducing the Alfven speed,
avoiding a too small timestep in MHD simulations with strong magnetic field
and/or extremely low gas density.
See the physics guide
(
[
PDF
](
https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf
)
).
for more information on the underlying idea and equations.
The use of the Boris correction is controlled
by setting corresponding macros in the user interface
[
nirvanaUser.h
](
3.2-User-interfaces#user-controllable-macros
)
.
NIRVANA implements the so-called Boris correction in the MHD equations.
Boris MHD artificially reduces wave propagation speeds in regions with
strong magnetic field and/or very low gas density which allows larger
integration timesteps. More details:
[
Physics
guide
](
https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf
)
.
Use of the Boris correction is controlled in the header file
[
`nirvanaUser.h`
](
#user-controllable-macros
)
.
#### Magnetic field splitting
NIRVANA implements a magnetic field splitting strategy. In this strategy the
magnetic field is decomposed into an imposed background field and a residual field.
By solving for the (presumably weaker) residual field accuracy of the solution
can be increased since the background field is discarded
in the time evolution of the total energy density.
For details on the modified MHD equations which are solved consult the
physics guide
(
[
PDF
](
https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf
)
).
Use of the magnetic field splitting is controlled
by setting the corresponding macro in the user interface
[
nirvanaUser.h
](
3.2-User-interfaces#user-controllable-macros
)
and by the definition of the background magnetic field in the user interface
[
sourceB0User.c
](
3.2-User-interfaces#user-defined-background-magnetic-field
)
.
## Dynamic arrays
NIRVANA implements a magnetic field splitting strategy. The magnetic
field is decomposed into a (stronger) background field and a (weaker)
residual field. This allows a more accurate computation of the thermal
pressure in adiabatic simulations since the modified total energy
density only contains the residual magnetic field part. More details:
[
Physics
guide
](
https://gitlab.aip.de/ziegler/NIRVANA/-/tree/master/doc/pdf/PhysicsGuide.pdf
)
.
Use of the magnetic field splitting scheme is controlled in the header
file
[
`nirvanaUser.h`
](
#user-controllable-macros
)
and definition of the
background magnetic field in module
[
`sourceB0User.c`
](
#user-defined-background-magnetic-field
)
.
### Dynamic arrays
The NIRVANA infrastructure implements its own proprietary routines for
dynamic array allocation. The available array allocation- and
deallocation functions are collected in the module
`utilAF.c`
. Dynamic
arrays are essential in the code usage. In order to demonstrate how to
apply these functionality the code fragment
arrays are essential in code usage. The code fragment
int nx,ny,nz;
int *v;
...
...
@@ -175,58 +150,55 @@ apply these functionality the code fragment
free_Array2(A2);
free_Array3(A3);
for instance,
creates
demonstrates how to create,
for instance,
-
\(
i
\)
a vector
`v[ix]`
,
`ix`
=0
...
`nx`
, of
`int`
type with
`nx`
+1
elements
using the function
`Arrayi()`
,
-
\(
i
\)
a vector
`v[ix]`
,
\
.
..
`nx`
, of
`int`
type with
`nx`
+1
elements
using the function
`Arrayi()`
,
-
\(
ii
\)
a 2D array
`A2[iy][ix]`
,
`ix`
=0
...
`nx`
,
`iy`
=0
...
`ny`
, of
`double`
type
with (
`nx`
+1)
×
(
`ny`
+1) elements using the function
-
\(
ii
\)
a 2D array
`A2[iy][ix]`
,
\
.
..
`nx`
,
\
.
..
`ny`
, of
`double`
type
with (
`nx`
+1)
$
\t
imes$
(
`ny`
+1) elements using the function
`Array2()`
,
-
\(
iii
\)
a 3D array
`A3[iz][iy][ix]`
,
`
ix`
=0...
`nx`
,
`iy`
=0...
`ny`
,
`iz`
=0...
`nz`
, of
`double`
type with
(
`nx`
+1)×(
`ny`
+1)×(
`nz`
+1) elements
using the function
`Array3()`
-
\(
iii
\)
a 3D array
`A3[iz][iy][ix]`
,
\.
..
`nx`
,
\.
..
`ny`
,
\.
..
`nz`
, of
`
double`
type with (
`nx`
+1)$
\t
imes$(
`ny`
+1)$
\t
imes$(
`nz`
+1) elements
using the function
`Array3()`
and subsequently deallocates them.
and its subsequent deallocation. Note that in multi-d arrays the fastest
index representing the $x$-direction is rightmost.
*Note 1: In multi-d arrays the fastest index is rightmost.*
*
Note 2: Allocation of an array must always be followed by its
**Important.**
Allocation of an array must always be followed by its
deallocation when it is no longer used in order to free memory
resources.
*
*
Note 3: The module
`utilAF.c`
only implements array functions of such
dimensionality and type as really needed so far in NIRVANA.
*
## Mesh data structure
Knowledge of the data structure representing the numerical mesh is
important for coding intial conditions (IC) and other user interfaces
like non-standard boundary conditions (BC). In general, the mesh is a
set of individual, logically rectangular grid blocks (also called
superblocks in the following) of different size and resolution which are
organized in a system of linked lists. The system of linked lists,
representing the mesh, is led by a master mesh pointer,
`gm`
, which is
function argument to many user interfaces like
`configUser()`
(defining
IC). Technically,
`gm`
is a doubly pointer of struct type
`GRD`
, i.e.
`**GRD`
. The
`GRD`
struct type is declared in the header file
`nirvana.h`
. Dereferencing
`gm`
to
`gm[l]`
gives the first superblock in
a linked list which collects all superblocks belonging to mesh
refinement level
*l*
.
`gm[l]`
, like any grid block, is of
`*GRD`
type,
i.e., a pointer to struct
`GRD`
. The base level
*l*
=0 starting with
pointer
`gm[0]`
is special because it spans the computational domain.
The superblocks making up a refinement level
*l*
can be reached
consecutively going through the corresponding linked list. Starting with
`gm[l]`
the next grid block in the list is obtained with the ’next grid’
operator
`->next`
. The linked list ends if the next grid pointer gives
the NULL pointer. The following image illustrates the linked list
concept.
resources.
### Mesh data structure
Minimal knowledge about NIRVANA's mesh data structure is necessary for
defining intial conditions and user-defined boundary condition, and to
program other user interfaces. In general, the mesh on which the
equations are integrated is a set of (logically) rectangular grid blocks
(called superblocks in the following) of different size and, if AMR, of
different resolution. These superblocks are organized in a system of
linked lists. The system of linked lists representing the full mesh has
a master mesh pointer,
`gm`
. Technically,
`gm`
is a doubly pointer of
struct type
`GRD`
, i.e.
`*GRD`
. The
`GRD`
struct type is declared in the
header file
`nirvana.h`
and contains all grid information (attributes,
coordinates, variables arrays, etc.). Dereferencing
`gm`
to
`gm[l]`
,
gives the first superblock in a linked list collecting all superblocks
belonging to mesh refinement level $l$.
`gm[l]`
, like any superblock, is
of
`GRD`
type, i.e., a pointer to struct
`GRD`
. The base level, $l=0$,
spanning the computational domain starts with pointer
`gm[0]`
. In
uniform grid simulations that is the only list existing. The superblocks
of a refinement level $l$ are obtained by running through the
corresponding linked list. Starting with
`gm[l]`
as the first superblock
in that list the next superblock is given by the 'next grid' operator
`->next`
, i.e., the superblock
`g`
following
`g[l]`
is
`g=g[l]->next`
.
The list ends if the next grid pointer is a NULL pointer. The following
image illustrates the linked list concept.
Denoting the variable for a superblock pointer as
`g`
looping over mesh
refinement level
*l*
(index
`il`
) is as simple as
Assuming a serial run for the moment, looping over all superblocks
`g`
of mesh
refinement level
$l$
(index
`il`
) is as simple as
for(g=gm[il]; (g); g=g->next) /* LOOP OVER SUPERBLOCKS; STOPS IF g==NULL */
{
...
...
@@ -241,82 +213,69 @@ and looping over the full mesh is as simple as
...
}
where
`_C.level_all`
is a global parameter storing the highest
refinement level at a time.
*
Note 1: In uniform grid simulations there is just a base level
represented by the linked list
`gm[0]`
.
*
*
Note 2: In distributed-memory simulations each MPI thread manages its
own master mesh pointer
`gm`
representing the mesh portion of the
partition. In particular, this means that a thread’s
`gm[0]`
only covers
a subset of the domain-spanning base level.
*
*
Note 3: In addition to
`gm`
NIRVANA implements another grid hierarchy
imaging the numerical mesh. This dual hierarchy is represented by the
public master mesh pointer
`_G0`
(of
`**GRD`
type).
`_G0`
is of minor (or
none) importance for the user. It is important, however, for inherent
code tasks mainly in the mesh refinement-, mesh repartitioning and mesh
sychronization algorithms. In contrast to a superblock with variable
size in
`gm`
grid blocks in
`_G0`
are generic: it have fixed size of 4
grid cells per space dimension and do not provide extra ghost cells for
adding boundary values like superblocks. The
`_G0`
generic grid block
hierarchy resembles an (incomplete) oct tree structure. Both mesh
descriptions,
`gm`
and
`_G0`
, make use of the same struct type
`GRD`
.
*
*
Note 4: The
`gm`
-mesh is created before integration out of the
`_G0`
-mesh by refinement-level-wise clustering of generic grid blocks to
superblocks. Superblocks automatically get ghost cells to include
boundary values. After integration the solution is copied onto the
`_G0`
-mesh and the
`gm`
-mesh is destroyed.
*
The struct type
`GRD`
subsumes all information defining a grid block
(superblock or generic): the grid block size and cell spacings, metric
measures, various cell coordinates, array pointers for primary
where
`_C.level_all`
stores the highest refinement level at a time. In
distributed-memory simulations each MPI thread has its own master mesh
pointer
`gm`
and, thus, linked lists
`gm[l]`
represent only the mesh
portion located on a partition. In particular, in a uniform grid
simulation there is just one superblock,
`g=gm[0]`
, on each thread when
the workload balancer is regular block decomposition (cf.
[
Workload
balancing
](
#workload-balancing
)
).
**Note:**
In addition to the mesh master pointer
`gm`
, there exists a
dual master mesh pointer,
`_G0`
, of same type
`*GRD`
.
`_G0`
likewise
represents the full grid hierarchy but in terms of its basic
constituents: generic grid block units of fixed size ($4^3$ cubes in 3D,
$4^2$ squares in 2D). In consequence, the
`_G0`
mesh resembles an
(incomplete) oct tree structure. Actually,
`gm`
is constructed out of
`_G0`
by a block clustering procedure. Whereas
`_G0`
is of very minor
importance for the user, code internal processes like mesh refinement,
mesh repartitioning and mesh sychronization rely on it.
The struct type
`GRD`
subsumes all information defining a grid block --
superblock or generic block: grid block size and its cell spacings,
metric measures, various cell coordinates, array pointers for primary
variables, derived variables and others, as well as several pointers for
linkage to other grid blocks
or objects
. An element
`V`
in
`GRD`
for a
grid block instance
`g`
(of type
`GRD`
) is accessible through
dereferencing
`g->V`
, e.g
,
linkage to other grid blocks. An element
`V`
in
`GRD`
for a
superblock
`g`
(of type
`GRD`
) is accessible through
dereferencing
`g->V`
, for
example
,
| | |
|------------------------------------:|:--------------------------------------------------------------|
|
`g->next`
| is the
next grid block relative
to
`g`
in the
grid block
list |
|
`g->next`
| is the
superblock next
to
`g`
in the
linked
list |
|
`g->rho[iz][iy][ix]`
| is the mass density at
`g`
-index (
`ix`
,
`iy`
,
`iz`
) |
|
`g->x[ix]`
,
`g->y[iy]`
,
`g->z[iz]`
| are cell-nodal coordinates at
`g`
-index (
`ix`
,
`iy`
,
`iz`
) |
|
`g->xc[ix]`
,
`g->yc[iy]`
,
`g->zc[iz]`
| are cell-centroid coordinates at
`g`
-index (
`ix`
,
`iy`
,
`iz`
) |
The following is a complete listing of the elements V of struct
`GRD`
grouped into the classes: grid connectivity, physical variables (primary
and derived), auxiliary arrays, grid metrics and grid attributes. Vector
and array elements in
`GRD`
are dynamic structures which are allocated
when a grid block is created. The array dimensionality is seen by the
number of
`[]`
-brackets. The indicated vector- and array indices refer
to:
The following is a complete listing of elements V in struct type
`GRD`
.
It are grouped into the categories: grid connectivity, physical
variables (primary and derived), auxiliary arrays, grid metrics and grid
attributes. Vector and array elements in
`GRD`
are dynamic structures
which are allocated when a grid block is created. Referred indices in a
vector or array have the following meaning:
| index | meaning |
|---------------------:|:-----------------------------------------------------------------------|
|
`ix`
,
`iy`
,
`iz`
| grid index in
*x*
,
*y*
,
*z*
-direction
|
|
`ibx`
,
`iby`
,
`ibz`
| generic
grid
block index in
*x*
,
*y*
,
*z*
-direction on superblock |
|
`ibcx`
,
`ibcy`
,
`ibcz`
| child gri
d
block index in
*x*
,
*y*
,
*z*
-direction on parent
grid
block |
|---------------------:|:-----------------------------------------------------------------------
-
|
|
`ix`
,
`iy`
,
`iz`
| grid
block cell
index in
*x*
,
*y*
,
*z*
-direction |
|
`ibx`
,
`iby`
,
`ibz`
| generic block index in
*x*
,
*y*
,
*z*
-direction on superblock
|
|
`ibcx`
,
`ibcy`
,
`ibcz`
| child g
ene
ri
c
block index in
*x*
,
*y*
,
*z*
-direction on parent block
|
|
`iu`
| physical variable counter |
|
`is`
| species counter
*s*
= 0,
*N*
<sub>
*s*
</sub>
− 1 |
|
`ic`
| tracers counter
*c*
= 0,
*N*
<sub>
*c*
</sub>
− 1 |
|
`it`
| testfield fluctuation variables counter
*t*
= 0,
*N*
<sub>
*t*
</sub>
− 1 |
|
`it`
| testfield
s
fluctuation variables counter
*t*
= 0,
*N*
<sub>
*t*
</sub>
− 1 |
The important elements for the user are the primary physical variables
and the grid coordinates which are most relevant for implementing user
interfaces like
`configUser.c`
.
Not all elements in
`GRD`
are of relevance for the user. Obviously, the
more important ones are the physical variables and grid coordinates.
Elements V of struct
`GRD`
:
| V (grid connectivity) | description |
|:------------------------------------|:--------------------------------------------------------
-
|
|:------------------------------------|:--------------------------------------------------------|
|
`next`
| next grid block pointer in the linked list |
|
`gc[ibcz][ibcy][ibcx]`
| child gri
d
block
pointer
s of parent grid block
|
|
`gb[ibz][iby][ibx]`
| superblock-
associated
generic
grid
block
pointers
|
|
`gp`
| pointer to
the
parent
generic
grid block |
|
`gs`
| pointer to a generic
grid
block
’s
associated superblock |
|
`gxl`
,
`gxu`
,
`gyl`
,
`gyu`
,
`gzl`
,
`gzu`
| pointers to
adjacent generic grid blocks of a grid block
|
|
`gc[ibcz][ibcy][ibcx]`
|
pointers to
child g
ene
ri
c
blocks of
a
parent grid block |
|
`gb[ibz][iby][ibx]`
|
pointers to
superblock-
defining
generic block
s
|
|
`gp`
| pointer to parent grid block
|
|
`gs`
| pointer to a generic block associated superblock
|
|
`gxl`
,
`gxu`
,
`gyl`
,
`gyu`
,
`gzl`
,
`gzu`
| pointers to
spatial grid block neighbors
|
|
`ol`
| pointer to a grid block-connected object list |
| V (primary physical variables) | description |
...
...
@@ -399,105 +358,111 @@ Elements V of struct `GRD`:
|
`pos[3]`
| grid index of a generic block within parent block (AMR) |
|
`posc[3]`
| grid index of a generic block within its associated superblock |
A cell with index (
`ix`
,
`iy`
,
`iz`
) of
a grid block
`g`
spans the domain
\[
`g->x[ix]`
,
`g->x[ix+1]`
\]
×
\[
`g->y[iy]`
,
`g->y[iy+1]`
\]
×
\[
`g->z[iz]`
,
`g->z[iz+1]`
\]
where (
`g->x[ix]`
,
`g->y[iy]`
,
`g->z[iz]`
) are the coordinates of
the
lower cell
node. Besides cell-nodal coordinates cell-centroid
c
oordinates locating the center of cell volume and
face
-
cent
roid
coordinates locating the center of cell faces are of importance. Their
definition for cylindrical- and spherical geometry is given in the tabl
e
below where
*Δ*
<sub>
`ix`
</sub>
(
*Δ*
<sub>
`iy`
</sub>
) denotes the
difference of a quantity evaluated at
`ix`
+1 (
`iy`
+1) and
`ix`
(
`iy`
).
Formal half-index subscripts mean that a quantity is to be evaluated at
the geometric cell center. Corresponding coordinates are called
cell-centered
.
A cell with index (
,,) within
a grid block
`g`
spans the domain
\[
`g->x[ix]`
,
`g->x[ix+1]`
\]
$
\t
imes$
\[
`g->y[iy]`
,
`g->y[iy+1]`
\]
$
\t
imes$
\[
`g->z[iz]`
,
`g->z[iz+1]`
\]
where (
`g->x[ix]`
,
`g->y[iy]`
,
`g->z[iz]`
) are the
nodal
coordinates of
the
lower cell
corner. Cell-centroid coordinates locate the volume
c
enter of a cell and face-centroid coordinates locate the
face
cent
ers
of a cell. Their definition for cylindrical- and spherical geometry is
given in the table below where $
\D
_\mix$ ($\D_
\m
iy$) denotes th
e
difference operator, $
\D
_\mix U=U_
{{
\t
t ix}+1}-U_{
\t
t ix}$ and analog
for $
\D
_
\m
iy$. Half-index subscripts mean evaluation at the geometric
center, i.e., at cell-centered coordinates. At certain cell locations
centroid coordinates and centered coordinates coincide which, in
general, is true in Cartesian geometry
.
| location |
`g->`
| definition: cylindrical | definition: spherical |
|:--------------------------
----|:----
---------|:--------------------------------------------------------------------------------|:----------------------------------------------------------------
----------------
|
|:--------------------------
|:
---------|:--------------------------------------------------------------------------------|:----------------------------------------------------------------|
| cell-centroid x |
`xc[ix]`
|
*z*
<sub>
`ix`
+1/2
</sub>
| 3
*Δ*
<sub>
`ix`
</sub>
*r*
<sup>
4
</sup>
/(4
*Δ*
<sub>
`ix`
</sub>
*r*
<sup>
3
</sup>
) |
| cell-centroid y |
`yc[iy]`
| 2
*Δ*
<sub>
`iy`
</sub>
*R*
<sup>
3
</sup>
/(3
*Δ*
<sub>
`iy`
</sub>
*R*
<sup>
2
</sup>
) |
*Δ*
<sub>
`iy`
</sub>
(
*θ*
cos
*θ*
− sin
*θ*
)/
*Δ*
<sub>
`iy`
</sub>
cos
*θ*
|
| cell-centroid z |
`zc[iz]`
|
*ϕ*
<sub>
`iz`
+1/2
</sub>
|
*ϕ*
<sub>
`iz`
+1/2
</sub>
|
|
*xy*
/
*xz*
-face-centroid x |
`xf[ix]`
|
*z*
<sub>
`ix`
+1/2
</sub>
| 2
*Δ*
<sub>
`ix`
</sub>
*r*
<sup>
3
</sup>
/(3
*Δ*
<sub>
`ix`
</sub>
*r*
<sup>
2
</sup>
) |
|
*yx*
-face-centroid y | (≡)
`yc[iy]`
| 2
*Δ*
<sub>
`iy`
</sub>
*R*
<sup>
3
</sup>
/(3
*Δ*
<sub>
`iy`
</sub>
*R*
<sup>
2
</sup>
) |
*Δ*
<sub>
`iy`
</sub>
(
*θ*
cos
*θ*
− sin
*θ*
)/
*Δ*
<sub>
`iy`
</sub>
cos
*θ*
|
*
Note 5: Other cell/face-centroid coordinates not listed in the table
are identical to their cell/face-centered coordinates.
*
*
Note 6: In Cartesian geometry all cell/face-centroid coordinates are
identical to their cell/face-centered counterparts.
*
## Adaptive mesh refinement
The mesh refinement algorithm relies on the oct-tree data structure
`_G0`
(of type
`*GRD`
) which represents a hierarchy of nested generic
grid blocks with a generic grid block having a fixed number of 4 grid
cells per coordinate direction independent of the refinement level. A
generic grid block of refinement level
*l*
has grid spacings
(
*δx*
<sup>
(
*l*
)
</sup>
,
*δy*
<sup>
(
*l*
)
</sup>
,
*δz*
<sup>
(
*l*
)
</sup>
) = (
*δx*
<sup>
(0)
</sup>
,
*δy*
<sup>
(0)
</sup>
,
*δz*
<sup>
(0)
</sup>
)/2
<sup>
*l*
</sup>
where
(
*δx*
<sup>
(0)
</sup>
,
*δy*
<sup>
(0)
</sup>
,
*δz*
<sup>
(0)
</sup>
)
is the grid spacing of the base level, i.e., there is a doubling of
resolution from one refinement level to the next higher.
Starting on the base level (
*l*
=0) mesh refinements are realized by
overlaying child grid blocks at locations which need refinement
according to different refinement criteria discussed in more detail
below. A child grid block covers one octant (a logical cube of 2 grid
cell length) of its parent grid block. Grid blocks can be arbitrarily
nested up to some maximum refinement level given by the
user-controllable code parameter
`_C.amr`
.
|
*yz*
-face-centroid x |
`x[ix]`
|
*z*
<sub>
`ix`
+1/2
</sub>
|
*r*
<sub>
`ix`
+1/2
</sub>
) |
|
*xy*
/
*yz*
-face-centroid y |
`yc[iy]`
| 2
*Δ*
<sub>
`iy`
</sub>
*R*
<sup>
3
</sup>
/(3
*Δ*
<sub>
`iy`
</sub>
*R*
<sup>
2
</sup>
) |
*Δ*
<sub>
`iy`
</sub>
(
*θ*
cos
*θ*
− sin
*θ*
)/
*Δ*
<sub>
`iy`
</sub>
cos
*θ*
|
|
*xz*
-face-centroid y |
`y[iy]`
|
*R*
<sub>
`iy`
+1/2
</sub>
|
*θ*
<sub>
`iy`
+1/2
</sub>
) |
|
*xy*
-face-centroid z |
`z[iz]`
|
*ϕ*
<sub>
`iz`
</sub>
|
*ϕ*
<sub>
`iz`
</sub>
|
|
*xz*
/
*yz*
-face-centroid z |
`zc[iz]`
|
*ϕ*
<sub>
`iz`
+1/2
</sub>
|
*ϕ*
<sub>
`iz`
+1/2
</sub>
|
### Adaptive mesh refinement
The mesh refinement algorithm relies on the oct-tree data structure of
generic blocks (fixed size of 4 cells per direction) represented by the
master mesh pointer
`_G0`
(of type
`*GRD`
). A generic block of
refinement level $l$ has cell spacings
$(
\d
x^{(l)},
\d
y^{(l)},
\d
z^{(l)})=
(
\d
x^{(0)},
\d
y^{(0)},
\d
z^{(0)})/2^l$ where
$(
\d
x^{(0)},
\d
y^{(0)},
\d
z^{(0)})$ is the cell spacing of the base
level, i.e., resolution doubles each time when refined.
Starting on the base level ($l=0$) mesh refinements are realized by
introducing child blocks at locations which need refinement according to
different refinement criteria. A child block covers one octant of its
parent block. Blocks can be arbitrarily nested up to some maximum
refinement level given by the parameter
`_C.amr`
.
AMR is controlled by different types of refinement criteria which can be
individually selected by the user. The standard implementation covers
criteria which are
the following criteria:
-
derivatives-based (most important in practice):
The criterion is applied to primary (or related) physical variables
*U*
( = {𝜚,
**m**
,
*e*
,
**B**
,
*C*
<sub>
*c*
</sub>
}). Undivided first
(
*δU*
) and second (
*δ*
<sup>
2
</sup>
*U*
) differences are
computed along various spatial directions. The switch
*α*
∈
\[
0, 1
\]
(code parameter:
`_C.amr_d1`
) selects between a purely
gradient-based criterion (
*α*
=1) and a purely
second-derivatives-based criterion (
*α*
=0).
*U*
<sub>
*r**e**f*
</sub>
are reference values to avoid division by
zero for indefinite variables. The exponent
*ξ*
∈
\[
0, 2
\]
(code
parameter:
`_C.amr_exp`
) introduces a power law functional on the
grid spacing ratio
*δx*
<sup>
(
*l*
)
</sup>
/
*δx*
<sup>
(0)
</sup>
which allows to control the behavior with progressive mesh
refinement. The macro
`FILTER`
(preset to 10
<sup>
− 2
</sup>
) is a
filter to suppress refinement at small-scale solution wiggles. The
individual refinement thresholds are denoted by ℰ
<sub>
*U*
</sub>
(code parameter:
`_C.amr_eps[0-4]`
).
*Note 1: Mesh refinement takes place if the criterion for refinement
is fulfilled for at least ONE variable *U*. On the other hand, mesh
$$
\l
eft[
\a
lpha
\f
rac{|
\d
U|}{|U|+U_
\m
athrm{ref}}+(1-
\a
lpha)
\f
rac{|
\d
^2U|}{|
\d
U|+{
\r
m FILTER}
\c
dot(|U|+U_
\m
athrm{ref})}
\r
ight]
\l
eft(
\f
rac{
\d
x^{(l)}}{
\d
x^{(0)}}
\r
ight)^{
\x
i}
\l
eft
\{\b
egin{array}{ll}
>\mathcal{E}_{U} &\exists U \quad\hbox{refinement}\\
<0.8
\m
athcal{E}_{U} &
\f
orall U
\q
uad
\h
box{derefinement}
\\
\e
nd{array}
\r
ight.$$
The criterion is applied to physical variables $U$
($=\{\varrho,\mathbf{m},e,\mathbf{B},C_\mathrm{c}\}$). Undivided
first ($\d U$) and second ($\d^2U$) differences are computed along
various spatial directions. The switch $\alpha\in [0,1]$ (code
parameter: `_C.amr_d1`) selects between a purely gradient-based
criterion ($\alpha =1$) and a purely second-derivatives-based
criterion ($\alpha =0$). $U_\mathrm{ref}$ are reference values to
avoid division by zero for indefinite variables. The exponent
$\xi\in [0,2]$ (code parameter: `_C.amr_exp`) introduces a power law
functional on the grid spacing ratio $\d x^{(l)}/\d x^{(0)}$ which
allows to tune the behavior with progressive mesh refinement. The
macro `FILTER` (preset to $10^{-2}$) is a filter to suppress
refinement at small-scale solution wiggles. The individual
refinement thresholds are denoted by $\mathcal{E}_U$ (code
parameters: `_C.amr_eps[0-4]`).
Mesh refinement takes place if the criterion for refinement is
fulfilled for at least ONE variable $U$. On the other hand, mesh
derefinement takes place if the criterion for derefinement is
fulfilled for ALL variables
*U*.*
*Note 2: The criterion is evaluated in the parent grid block octant
to be refined and in a buffer zone around it in order to capture
flow features right in time.*
-
Jeans-length-based (important in simulations involving selfgravity)
where the first factor is the local Jeans length with
*c*
<sub>
*s*
</sub>
the sound speed and
*G*
the gravitational
constant, and where
*δs*
= min {
*δx*
,
*h*
<sub>
*y*
</sub>
*δy*
,
*h*
<sub>
*z*
</sub>
*δz*
}
.
ℰ
<sub>
*Jeans*
</sub>
is a user-specific threshold
with
the
reciprocal value givi
ng
th
e minimum number of grid cells the local
Jeans length should be resolved
.
fulfilled for ALL variables
$U$. The criterion is evaluated in the
vicinity of a block's octant (including a buffer of 2 grid cells per
side) to decide whether a child block is to be created.
-
Jeans-length-based (important in simulations involving selfgravity):
$$
\l
eft(
\f
rac{
\p
i}{G}
\f
rac{c_s^2}{
\v
arrho}
\r
ight)^{1/2}
\c
dot
\m
athcal{E}_
\m
athrm{Jeans}
\l
eft
\{\b
egin{array}{ll}
<
\d
s &
\m
box{refinement}
\\
>1.25\d s & \mbox{derefinement}\\
\e
nd{array}
\r
ight.$$ where the first factor is the local Jeans
length with $c_s$ the sound speed and $G$ the gravitational
constant, and where $
\d
s=
\m
in
\{\d
x,h_y
\d
y, h_z
\d
z
\}
$
.
$
\m
athcal{E}_
\m
athrm{Jeans}$
is a user-specific threshold
giving
the
fraction of the local Jeans le
ngth
to be resolved by at least 1 grid
cell
.
-
Field-length-based (potentially relevant for simulations involving a
heatloss source but rarely used in practice):
where the first factor is the Field length with
*L*
(𝜚,
*T*
) the
density- and temperature-dependent heatloss function and
*κ*
the
thermal conduction coefficient. ℰ
<sub>
*Field*
</sub>
is a
user-specific threshold with the reciprocal value giving the minimum
number of grid cells the Field length should be resolved.
## Parallelism and load balancing
$$2
\p
i
\l
eft(
\f
rac{
\k
appa T}{
\m
ax
\l
eft
\{
T
\l
eft(
\f
rac{
\p
artial L}{
\p
artial T}
\r
ight)_{
\v
arrho}
-
\v
arrho
\l
eft(
\f
rac{
\p
artial L}{
\p
artial
\v
arrho}
\r
ight)_{T},EPS
\r
ight
\}
}
\r
ight)^{1/2}
\c
dot
\m
athcal{E}_
\m
athrm{Field}
\l
eft
\{\b
egin{array}{ll}
<
\d
s &
\m
box{refinement}
\\
>1.25\d s & \mbox{derefinement}\\
\e
nd{array}
\r
ight.$$ where the first factor is the Field length with
$L(
\v
arrho, T)$ the density- and temperature-dependent heatloss
function and $
\k
appa$ the thermal conduction coefficient.
$
\m
athcal{E}_
\m
athrm{Field}$ is a user-specific threshold giving the
fraction of the local Field length to be resolved by at least 1 grid
cell.
### Workload balancing
The NIRVANA code is parallelized making use of the MPI library. The code
infrastructure provides routines for mesh partitioning in order to
...
...
@@ -508,59 +473,53 @@ process is first collected and then transmitted in two separate single
streams representing physical data and associated meta data. Whenever
possible one-to-one MPI communications are done in non-blocking mode.
Three different load balancing schemes are in use in NIRVANA. Two of the
schemes serve to distribute workload among MPI threads for the non-NCCM
part of the solution process (anything but chemo-thermal reactions).
Workload balancing is here approached by an equal distribution of
generic grid blocks assuming that each generic grid block produce the
same computational costs.
The first and simpler of the two non-NCCM schemes uses regular domain
decomposition applicable exclusively to uniform grid simulations. In
regular decomposition the mesh is divided into a number of same-sized
(logically) rectangular submeshes equal to the number of MPI threads.
The dimension of a submesh in any coordinate direction is subject to the
limitation that it must be a multiple of 4 – the size of a generic grid
block.
The second, non-NCCM load balancing scheme is applied to AMR
simulations. Here, mesh partitioning relies on space-filling-curve
mapping techniques. In the space-filling-curve-based scheme the mesh is
partitioned by an equal subdivision of the space-filling-curve interval.
The distributed mesh subdomains are polyhedral bounded but keep optimal
data compactness (surface-to-volume ratio minimized). At the same time,
NIRVANA features inter-level data locality by allowing smart grid block
overlapping (shared block technique) at partition interfaces. This
combination of space-filling-curve mesh partitioning and shared block
method means that the portion of grid hierarchy assigned to any thread
can be traversed without communication at the expense of additional code
complexity.
The above described load balancing strategies are not optimal for the
NCCM part of the solution process. This is because NIRVANA uses an
efficient ODE solver for evolving the thermo-chemical rate equations
based on error-controlled adaptive time-stepping. Hence, the
computational cost to solve the rate equations varies from grid cell to
grid cell which can result in strong workload imbalance between threads
if grid blocks were equally distributed and the NCCM solution step
dominates the total computational cost. To address this issue NIRVANA
switches to a third load balancing approach for the NCCM solution
process. The principle idea behind is to measure the computational cost
for the solution of the rate equations for a conglomerate of grid cells
representing an octant of a generic grid block (as smallest unit) with
help of the MPI timer. The a-posteriori measure thus obtained is then
used as an a-priori estimate for the load balancer in the next NCCM
solver call. The approach is only approximate and works the better the
less changes in density and temperature occur between two successive
NCCM updates. The measured workload residuals on threads (difference
between actual and desired workload per thread) are diminished by an
octant-wise redistribution of relevant data from threads with workload
excess to threads with workload deficiency before integration. After
integration the solution is copied back onto the original mesh.
*
Note: The NCCM load balancer comes to its limit when the computational
cost is highly concentrated on only a few grid cells.
*
The non-NCCM solution part of the solver (anything but chemo-thermal
rate equations) and the NCCM solution part of the solver (chemistry) use
different load balancing schemes. For the non-NCCM part there are 2
choices available:
-
**regular block decomposition.**
Workload balancing is achieved by
an equal distribution of generic blocks among threads where each
thread holds a same-sized (logically) rectangular submesh of the
mesh. The dimension of a thread's submesh in any coordinate
direction must be a multiple of 4 -- the size of a generic block.
Regular block decomposition is prefered for uniform grid simulations
and is not applicable to AMR simulations.
-
**space-filling curve (SFC) decomposition.**
Workload balancing is
achieved by an equal distribution of generic blocks among threads
where mapping of the mesh follows a space-filling-curve. A thread's
mesh partition thus corresponds to a subinterval of the
space-filling curve. SFC decomposition is applied in AMR simulations
but is also available for uniform grid simulations. SFC
decomposition achieves good data compactness with minimized
surface-to-volume ratio of the overall partition structure.
**Note:**
NIRVANA-AMR features inter-level data locality which means
that the AMR grid hierarchy on any thread can be traversed without
communication at the expense of additional code complexity due to smart
partition overlapping.
The load balancing strategy for the NCCM part of the solver is more
involved. Since the ODE solver for evolving the thermo-chemical rate
equations forward in time is based on error-controlled adaptive
time-stepping, the computational cost can vary from grid block to grid
block substantially. Strong workload imbalance may occur when using the
non-NCCM balancer above. Instead a load balancer based on a performance
measure is used in addition: the computational cost for the solution of
the rate equations on each cell in an octant of a generic block is
measured with help of the MPI timer and summed for individual octants.
This a-posteriori measure is then used as an a-priori estimate in the
load balancer before the next call of the NCCM solver. Workload
residuals, i.e. differences between actual and desired workload per
thread, are diminished by an octant-wise redistribution of chemical
reaction data from threads with workload excess to threads with workload
deficiency. This load balancer approach is only approximate and works
the better the less changes in density and temperature occur between two
successive NCCM updates. It definitely comes to its limit when the
computational cost is highly concentrated on only a few grid cells.
PREV:
[
2 Getting started
](
2-Getting-started
)
NEXT:
[
3.2 User interfaces
](
3.2-User-interfaces
)
PREV:
[
2 Getting started
](
https://gitlab.aip.de/ziegler/NIRVANA/-/wikis/2-Getting-started
)
NEXT:
[
3.2 User interfaces
](
3.2-User-interfaces
)
