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).