Problem Setup setup.cpp
- The source code
setup.cpp
contains the code specific to the physical setup at hand. It serves several purposes: Get the setup parameters from the problem input file.
Init arrays and variables you will need
Set the initial conditions
Define and enroll setup-specific functions which will be called from the integration loop
Define and enroll setup-specific outputs
Most of these initialisations relies on the class Setup
which has to be implemented in your
setup.cpp
. At this stage, if you are not familliar with the idefix_loop
structures,
IdefixArray
data arrays and Idefix class structure, it is strongly recommended you have a
look at the Programming guide.
The Setup
class
The Setup
class is declared as follows:
class Setup {
public:
Setup(Input &, Grid &, DataBlock &, Output &);
void InitFlow(DataBlock &);
};
As it can be seen, this class consist of a constructor and one mandatory method: InitFlow
, which will handle
the initial condition.
The Setup
constructor
Let us start with the constructor Setup(Input &, Grid &, DataBlock &, Output &).
This constructor is called on the code startup and allows the user to load and set setup-specific parameters. The parameters are four objects
which have already been initialised when Setup
is called: Input
, Grid
, DataBlock
and Output
(see Useful classes).
A typical constructor first loads the setup parameters calling accessors from the Input
object (see The Input class). Then,
if there are some user-defined functions (for instance a user-defined potential, boundary condition or output),
the constructor also enrolls these functions before returning.
Function enrollment
The enrollment of user functions is required in Idefix whenever a parameter is set to “userdef” in
the input file, and for user-defined outputs. This can be seen as a way to link the user
setup to the main code at runtime, and avoid the need to pre-define tens of empty functions. Function enrollment
is achieved by calling one of the EnrollXXX
function of the class associated to it.
For instance, the Hydro
class provide the following list of enrollment functions (declared in hydro.hpp):
// Enroll user-defined boundary conditions
void EnrollUserDefBoundary(UserDefBoundaryFunc);
void EnrollInternalBoundary(InternalBoundaryFunc);
void EnrollEmfBoundary(EmfBoundaryFunc);
// Enroll user-defined gravitational potential
void EnrollGravPotential(GravPotentialFunc);
// Enroll user source terms
void EnrollUserSourceTerm(SrcTermFunc);
// Enroll user-defined ohmic, ambipolar and Hall diffusivities
void EnrollOhmicDiffusivity(DiffusivityFunc);
void EnrollAmbipolarDiffusivity(DiffusivityFunc);
void EnrollHallDiffusivity(DiffusivityFunc);
// Enroll user-defined isothermal sound speed
void EnrollIsoSoundSpeed(IsoSoundSpeedFunc);
When called, these function expects the address of the user-defined function. These user-defined function should have the following signatures (declared in hydro_defs.hpp):
using UserDefBoundaryFunc = void (*) (DataBlock &, int dir, BoundarySide side,
const real t);
using GravPotentialFunc = void (*) (DataBlock &, const real t, IdefixArray1D<real>&,
IdefixArray1D<real>&, IdefixArray1D<real>&,
IdefixArray3D<real> &);
using SrcTermFunc = void (*) (DataBlock &, const real t, const real dt);
using InternalBoundaryFunc = void (*) (DataBlock &, const real t);
using EmfBoundaryFunc = void (*) (DataBlock &, const real t);
using DiffusivityFunc = void (*) (DataBlock &, const real t, IdefixArray3D<real> &);
using IsoSoundSpeedFunc = void (*) (DataBlock &, const real t, IdefixArray3D<real> &);
Example
The following example have a user-defined gravitational potential, and defines a Setup
constructor which reads a parameter from the .ini file and enroll the user-defined potential.
// a global variable which stores the mass of some object
real Mass;
// user-defined potential
void Potential(DataBlock& data, const real t, IdefixArray1D<real>& x1, IdefixArray1D<real>& x2, IdefixArray1D<real>& x3, IdefixArray3D<real>& phi) {
idefix_for("Potential",0,data.np_tot[KDIR], 0, data.np_tot[JDIR], 0, data.np_tot[IDIR],
KOKKOS_LAMBDA (int k, int j, int i) {
phi(k,j,i) = -Mass/x1(i);
});
}
// Setup constructor
Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) {
// Read some parameter from the ini file
Mass = input.Get<real>("Setup","mass",0);
// Enroll the user-defined potential
data.hydro.EnrollGravPotential(&Potential);
}
User-defined boundaries
If one (or several) boundaries are set to userdef
in the input file, the user needs to
enroll a user-defined boundary function in the Setup
constructor as for the other user-def functions (see Function enrollment).
Note that even if several boundaries are userdef
in the input file, only one user-defined function
is required. When Idefix calls the user defined boundary function, it sets the direction of the boundary (dir=IDIR
, JDIR
,
or KDIR
) and the side of the bondary (side=left
or side=right
). If conveninent, one can use
the BoundaryFor
wrapper functions to automatically loop on the boundary specified by dir
and side
.
A typical user-defined boundary condition function looks like this:
void UserdefBoundary(DataBlock& data, int dir, BoundarySide side, real t) {
IdefixArray4D<real> Vc = data.hydro.Vc;
IdefixArray4D<real> Vs = data.hydro.Vs;
if(dir==IDIR) {
data.hydro.boundary.BoundaryFor("UserDefBoundary", dir, side,
KOKKOS_LAMBDA (int k, int j, int i) {
Vc(RHO,k,j,i) = 1.0;
Vc(VX1,k,j,i) = 0.0;
Vc(VX2,k,j,i) = 0.0;
Vc(VX3,k,j,i) = 0.0;
});
// For magnetic field (defined on cell sides), we need specific wrapper functions
// Note that we don't need to initialise the field component parallel to dir, as it is
// automatically reconstructed from the solenoidal condition and the tangential components
data.hydro.boundary.BoundaryForX2s("UserDefBoundaryBX2s", dir, side,
KOKKOS_LAMBDA (int k, int j, int i) {
Vs(BX2s,k,j,i) = 0.0;
});
data.hydro.boundary.BoundaryForX3s("UserDefBoundaryBX3s", dir, side,
KOKKOS_LAMBDA (int k, int j, int i) {
Vs(BX3s,k,j,i) = 0.0;
});
}
}
Setup::InitFlow
method
Basics of the Initflow method
InitFlow
is a method of the Setup
class and is called by Idefix after the Setup
constructor.
Its role is to define the initial conditions for the flow, initializing the Vc
(and Vs
in MHD)
arrays of the Hydro
class, for instance. Because this method does not have to be fast, since
it is called only once, it is customary to initialise the flow on the host, and then send it to the
device.
For this, it is useful to first define a mirror DataBlockHost
(see DataBlockHost class)
of the DataBlock
given in argument and initialse the flow in DataBlockHost
using a standard
C loop on the host, as in the example below.
void Setup::InitFlow(DataBlock &data) {
// Create a host copy of the DataBlock given in argument
DataBlockHost dataHost(data);
// Because we initialise the arrays in DataBlockHost,
// we can execute the loop on the host
for(int k = 0; k < dataHost.np_tot[KDIR] ; k++) {
for(int j = 0; j < dataHost.np_tot[JDIR] ; j++) {
for(int i = 0; i < dataHost.np_tot[IDIR] ; i++) {
real x = dataHost.x[IDIR](i);
real y = dataHost.x[JDIR](j);
real z = dataHost.x[KDIR](k);
dataHost.Vc(RHO,k,j,i) = 1.0;
dataHost.Vc(PRS,k,j,i) = 1.0;
dataHost.Vc(VX1,k,j,i) = -sin(y);
dataHost.Vc(VX2,k,j,i) = sin(x)+cos(z);
dataHost.Vc(VX3,k,j,i) = cos(x);
dataHost.Vs(BX1s,k,j,i) = -sin(y);
dataHost.Vs(BX2s,k,j,i) = sin(x);
dataHost.Vs(BX3s,k,j,i) = 0.0;
}
}
}
// Do not forget to send our initialisation to the parent dataBlock!
dataHost.SyncToDevice();
}
Warning
Do not forget to sync your DataBlockHost to its parent DataBlock using the
DataBlockHost::SyncToDevice()
method!
Initialising the magnetic field
When MHD is used, the face-centered magnetic field stored in Vs
should be initialised with a divergence-free
field at machine precision. This might not always be straightforward for some complex field geometry,
so dataBlockHost can also be initialised with a vector potential, from which the face-centered field
can be automatically derived using DataBlockHost::MakeVsFromAmag
as in the example below:
void Setup::InitFlow(DataBlock &data) {
// Create a host copy of the DataBlock given in argument
DataBlockHost dataHost(data);
// Allocate an array on host to store the vector potential (3 components are expected)
IdefixHostArray4D<real> A = IdefixHostArray4D<real>("Setup_VectorPotential", 3,
data.np_tot[KDIR],
data.np_tot[JDIR],
data.np_tot[IDIR]);
for(int k = 0; k < dataHost.np_tot[KDIR] ; k++) {
for(int j = 0; j < dataHost.np_tot[JDIR] ; j++) {
for(int i = 0; i < dataHost.np_tot[IDIR] ; i++) {
real x = dataHost.x[IDIR](i);
real y = dataHost.x[JDIR](j);
real z = dataHost.x[KDIR](k);
// Initialise Vc field (not shown)
// ...
// Initialise the 3 components of the vector potential
A(IDIR,k,j,i) = 0.0;
A(JDIR,k,j,i) = 0.0;
A(KDIR,k,j,i) = -y*B0;
}
}
}
// Compute the face centered Vs from the vector potential
dataHost.MakeVsFromAmag(A);
// Do not forget to send our initialisation to the parent dataBlock!
dataHost.SyncToDevice();
}
Initialising from a restart dump
In some cases, it can be useful to initialise the flow from a dump taken from a previous
simulation. While one can simply use the -restart
option on the commandline to resume
a simulation (see Command line options), there are some situation when one needs to create
a new initial condition by extrapolating or extanding a restart dump (such as in a resolution
test or a dimension change). In this case, one should use the DumpImage
class which provides
all the tools needed to read a restart dump (see also DumpImage class).
One typically first construct an instance of DumpImage
in the Setup
constructor, and then
use this instance to initialise the flow in Setup::InitFlow
. The procedure is examplified below,
assuming we want to create a dump from mydump.dmp
:
DumpImage *image; // Global pointer to our DumpImage
// Setup constructor
Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output) {
image = new DumpImage("mydump.dmp",output); // load the dump file and store it in a DumpImage
}
// Flow initialisation, read directly from the DumpImage
void Setup::InitFlow(DataBlock &data) {
// Create a host copy
DataBlockHost d(data);
for(int k = d.beg[KDIR]; k < d.end[KDIR] ; k++) {
for(int j = d.beg[JDIR]; j < d.end[JDIR] ; j++) {
for(int i = d.beg[IDIR]; i < d.end[IDIR] ; i++) {
// Note that the restart dump array only contains the full (global) active domain
// (i.e. it excludes the boundaries, but it is not decomposed accross MPI procs)
int iglob=i-2*d.beg[IDIR]+d.gbeg[IDIR];
int jglob=j-2*d.beg[JDIR]+d.gbeg[JDIR];
int kglob=k-2*d.beg[KDIR]+d.gbeg[KDIR];
d.Vc(RHO,k,j,i) = image->arrays["Vc-RHO"](kglob,jglob,iglob);
d.Vc(PRS,k,j,i) = image->arrays["Vc-PRS"](kglob,jglob,iglob);
d.Vc(VX1,k,j,i) = image->arrays["Vc-VX1"](kglob,jglob,iglob);
}}}
// For magnetic variable, we should fill the entire active domain, hence an additional
// point in the field direction
for(int k = d.beg[KDIR]; k < d.end[KDIR] ; k++) {
for(int j = d.beg[JDIR]; j < d.end[JDIR] ; j++) {
for(int i = d.beg[IDIR]; i < d.end[IDIR]+IOFFSET ; i++) {
int iglob=i-2*d.beg[IDIR]+d.gbeg[IDIR];
int jglob=j-2*d.beg[JDIR]+d.gbeg[JDIR];
int kglob=k-2*d.beg[KDIR]+d.gbeg[KDIR];
d.Vs(BX1s,k,j,i) = image->arrays["Vs-BX1s"](kglob,jglob,iglob);
}}}
// And so on for the other components
// ..
delete image; // don't forget to free the memory allocated for dumpImage!
// Send our datablock to the device
d.SyncToDevice();
}
Note
Note that the naming convention in DumpImage::arrays
combines the original array and variable names.
It is generically written XX-YYY
where XX
is the array name in the dataBlock
(e.g.
Vc
or Vs
) and YYY
is the variable name (e.g. VX2
or BX3s
).
User-defined analysis
User-defined analysis and outputs can be coded in the setup.cpp
file. Follow the
guidelines in Outputs.