Mesh¶
The mesh is used in pretty much all parts of the code, and deals with
things like the geometry of the mesh (metric tensors etc.), and how
the mesh is divided between processors (communications). The Mesh
class defines an interface, and there is currently a single
implementation:
BoutMesh
(src/mesh/boutmesh.cxx
) which is backwards compatible with the BOUT and BOUT-06 codes. This is a logically rectangular mesh so the number of radial points (x) can’t change in the poloidal direction (y).
Grid data sources¶
All data sources inherit from GridDataSource
. They must supply a
method to test if a variable exists, GridDataSource::hasVar
:
bool hasVar(const string &name);
and then use the get
methods to get
integers or reals:
bool get(Mesh *m, <type> &variable, const string &name);
Loading a mesh¶
The Mesh constructor
takes GridDataSource
and
Options
objects. You can also call Mesh::create
with just one of
these objects, which will call out to the MeshFactory
singleton to
create a mesh “automatically”. This is the way that it is done in
bout++.cxx. Once you
have instantiated a Mesh
object, you can then call Mesh::load
to
read in all the appropriate variables from the GridDataSource
:
mesh = Mesh::create(); ///< Create the mesh
mesh->load(); ///< Load from sources. Required for Field initialisation
For post-processing of the results, it’s useful to have mesh
quantities in the dump files along with the results. To do this,
there’s the function Mesh::outputVars
(see also Datafile
and
Options
):
// Create an output file from an Options object
dump = Datafile(options->getSection("output"));
// Possibly add some other variables to the output file
...
// Save mesh configuration into output file
mesh->outputVars(dump);
which is called during BOUT++ initialisation.
Implementation: BoutMesh¶
BoutMesh
class uses the BOUT indices (which trace back to UEDGE):
int ixseps1, ixseps2, jyseps1_1, jyseps2_1, jyseps1_2, jyseps2_2;
ixseps1
and ixseps2
give
the X location of the separatrices, and are equal in the case of
single-null configurations. The indexing is such that all points 0
<= x < ixseps1
are inside the separatrix, whilst ixseps1 <= x <
LocalNx
are outside. See BOUT++ Topology for more details.
Index ranges¶
The Mesh
class includes several public members which describe the size
of the mesh, and are used all over BOUT++ to loop over variables:
/// Size of the mesh on this processor including guard/boundary cells
int LocalNx, LocalNy, LocalNz;
/// Local ranges of data (inclusive), excluding guard cells
int xstart, xend, ystart, yend;
Getting data¶
The Mesh::load
code above needs to read data for the mesh, and
physics codes usually need to read their initial profiles during
initialisation. To do this, Mesh provides an overloaded function
Mesh::get
:
int get(var, const char *name); // Request data from mesh file
where var
can be just about any BOUT++ datatype (Field2D
,
Vector3D
etc.).
Implementation: BoutMesh¶
For integers and BoutReals, the implementation is fairly trivial. Uses the Mesh protected functions to find a data source and read data from it:
GridDataSource* s = findSource(name); // Find a source of data
s->open(name); // Open the source
bool success = s->fetch(&ival, name); // Get the data
s->close(); // Close the source
To read 2D and 3D fields, the branch-cuts need to be taken into account.
Communications¶
The most common type of communication is to just exchange all guard cells with neighboring processors. Mesh provides the following commands for doing this:
template <typename... Ts>
int communicate(Ts&... ts); // Communicate one or more fields
int communicate(FieldGroup); // Communicate a group of fields
comm_handle send(FieldGroup); // Send data
int wait(comm_handle); // Receive data
Mesh::communicate
can be used to communicate any number of variables
together, and makes the code quite clear. For example in
examples/DriftInstability/2fluid.cxx
around line 360:
// Need to communicate jpar
mesh->communicate(jpar);
Since this uses the FieldData
interface like Datafile, this can be
used to communicate all BOUT++ field data types. You can also create a
FieldGroup
object to group fields together, then communicate them
all together:
FieldGroup comgrp; // Group of variables for communication
Field3D P;
Vector3D V;
comgrp.add(P); // Add the variables
comgrp.add(V); // Usually done in PhysicsModel::init
mesh->communicate(comgrp); // Communicate in PhysicsModel::rhs
Internally, this is how the templated Mesh::communicate
works.
If you want to overlap communications with calculations then use the
Mesh::send
and Mesh::wait
functions instead of
Mesh::communicate
:
comm_handle ch = mesh->send(comgrp); // Start the communications
// Calculations which don't need variables in comgrp
wait(ch); // Wait for all communications to finish
Implementation: BoutMesh¶
In BoutMesh
, the communication is controlled by the variables:
int UDATA_INDEST, UDATA_OUTDEST, UDATA_XSPLIT;
int DDATA_INDEST, DDATA_OUTDEST, DDATA_XSPLIT;
int IDATA_DEST, ODATA_DEST;
In the Y direction, each boundary region (Up and Down in Y)
can be split into two, with 0 <= x < UDATA_XSPLIT
going to the
processor index UDATA_INDEST
, and UDATA_INDEST <= x < LocalNx
going
to UDATA_OUTDEST
. Similarly for the Down boundary. Since there are
no branch-cuts in the X direction, there is just one destination for the
Inner and Outer boundaries. In all cases a negative
processor number means that there’s a domain boundary.
X communications¶
For parallel Laplacian inversions, communication is needed in the X direction only, and involves quantities which are not in Fields:
bool firstX(); // True if at the inner X boundary
bool lastX(); // True if at the outer X boundary
int NXPE, PE_XIND; // Number of processors in X, and X processor index
int sendXOut(BoutReal *buffer, int size, int tag);
sendXIn(BoutReal *buffer, int size, int tag);
comm_handle irecvXOut(BoutReal *buffer, int size, int tag);
comm_handle irecvXIn(BoutReal *buffer, int size, int tag);
The variables Mesh::NXPE
and Mesh::PE_XIND
shouldn’t really be
there, but are currently needed because the SPT algorithm in
LaplaceSPT
needs to know when it’s going to be next and so keep
track of which processor number is currently working. This logic to
pass a problem along a chain in X should really be moved into Mesh.
Y-Z surface communications¶
Some operations (like parallel inversions in
bout++/src/invert/invert_parderiv.cxx
) need to be performed on Y-Z
surfaces, i.e. slices at constant X. This needs to be able to handle
open and closed surfaces, and that closed surfaces may need a shift in
the Z direction to match one end onto the other (a twist-shift
condition).
The simplest operation is to average a quantity over Y with
averageY
.
To test if a particular surface is closed, there is the function
periodicY
.
The most general way to access data on surfaces is to use the
SurfaceIter
iterator, which can be created using
SurfaceIter::SurfaceIter
:
SurfaceIter* surface(mesh);
This then allows looping over the surfaces in the usual way:
for(surf->first(); !surf->isDone(); surf->next()) {
...
}
To test if the surface is closed, there’s the test SurfaceIter::closed
:
bool surf->closed(BoutReal &ts)
which returns true if the surface is closed, along with the twist-shift angle.
Initial profiles¶
The initial profiles code needs to construct a solution which is
smooth everywhere, with a form of perturbation specified in the input
file for each direction. In order to do this, it needs a continuous
function to use as an index. This is supplied by the functions
Mesh::GlobalX
and Mesh::GlobalY
:
BoutReal GlobalX(int jx); // Continuous X index between 0 and 1
BoutReal GlobalY(int jy); // Continuous Y index (0 -> 1)
which take a local x or y index and return a globally continuous x or y index.
Differencing¶
The mesh spacing is given by the public members Mesh::dx
, Mesh::dy
and Mesh::dx
:
// These used for differential operators
Field2D dx, dy;
Field2D d2x, d2y; // 2nd-order correction for non-uniform meshes
BoutReal zlength, dz; // Derived from options (in radians)
Metrics¶
While Mesh
handles the numerical details of the mesh, the “physical”
details are handled by Coordinates
. The contravariant and covariant
metric tensor components are public members of Coordinates
:
// Contravariant metric tensor (g^{ij})
Field2D g11, g22, g33, g12, g13, g23; // These are read in grid.cxx
// Covariant metric tensor
Field2D g_11, g_22, g_33, g_12, g_13, g_23;
int calcCovariant(); // Invert contravatiant metric to get covariant
int calcContravariant(); // Invert covariant metric to get contravariant
If only one of these sets is modified by an external code, then
Coordinates::calcCovariant
and Coordinates::calcContravariant
can
be used to calculate the other (uses Gauss-Jordan currently).
From the metric tensor components, Coordinates
calculates several
other useful quantities:
int jacobian(); // Calculate J and Bxy
Field2D J; // Jacobian
Field2D Bxy; // Magnitude of B = nabla z times nabla x
/// Calculate differential geometry quantities from the metric tensor
int geometry();
// Christoffel symbol of the second kind (connection coefficients)
Field2D G1_11, G1_22, G1_33, G1_12, G1_13;
Field2D G2_11, G2_22, G2_33, G2_12, G2_23;
Field2D G3_11, G3_22, G3_33, G3_13, G3_23;
Field2D G1, G2, G3;
These quantities are public and accessible everywhere, but this is because they are needed in a lot of the code. They shouldn’t change after initialisation, unless the physics model starts doing fancy things with deforming meshes.
Miscellaneous¶
There are some public members of Mesh
which are there for some
specific task and don’t really go anywhere else (yet).
To perform radial derivatives in tokamak geometry, interpolation is needed in the Z direction. This is done by shifting in Z by a phase factor, performing the derivatives, then shifting back. The following public variables are currently used for this:
bool ShiftXderivs; // Use shifted X derivatives
int ShiftOrder; // Order of shifted X derivative interpolation
Field2D zShift; // Z shift for each point (radians)
Field2D ShiftTorsion; // d <pitch angle> / dx. Needed for vector differentials (Curl)
Field2D IntShiftTorsion; // Integrated shear (I in BOUT notation)
bool IncIntShear; // Include integrated shear (if shifting X)
int TwistOrder; // Order of twist-shift interpolation
This determines what order method to use for the interpolation at the
twist-shift location, with 0
meaning FFT during communication. Since
this must be 0 at the moment it’s fairly redundant and should be
removed.
A (currently experimental) feature is:
bool StaggerGrids; ///< Enable staggered grids (Centre, Lower). Otherwise all vars are cell centred (default).