File mesh.hxx

Interface for mesh classes. Contains standard variables and useful routines.

Changelog

2014-12 Ben Dudson bd512@york.ac.uk

  • Removing coordinate system into separate Coordinates class
  • Adding index derivative functions from derivs.cxx

2010-06 Ben Dudson, Sean Farley

  • Initial version, adapted from GridData class
  • Incorporates code from topology.cpp and Communicator

Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu

Contact: Ben Dudson, bd512@york.ac.uk

This file is part of BOUT++.

BOUT++ is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

BOUT++ is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with BOUT++. If not, see http://www.gnu.org/licenses/.

Typedefs

using comm_handle = void *

Type used to return pointers to handles.

class Mesh

Subclassed by BoutMesh

Public Functions

Mesh()

Constructor for a “bare”, uninitialised Mesh Only useful for testing

Mesh(GridDataSource *s, Options *options)

Constructor

Parameters
  • s: The source to be used for loading variables
  • options: The options section for settings

~Mesh()

Destructor.

virtual int load()

Loads the mesh values

Currently need to create and load mesh in separate calls because creating Fields uses the global “mesh” pointer which isn’t created until Mesh is constructed

virtual void outputVars(Datafile &file)

Add output variables to a data file These are used for post-processing

int get(std::string &sval, const std::string &name, const std::string &def = "")

Get a string from the input source

Return
zero if successful, non-zero on failure
Parameters
  • sval: The value will be put into this variable
  • name: The name of the variable to read
  • def: The default value if not found

int get(int &ival, const std::string &name, int def = 0)

Get an integer from the input source

Return
zero if successful, non-zero on failure
Parameters
  • ival: The value will be put into this variable
  • name: The name of the variable to read
  • def: The default value if not found

int get(BoutReal &rval, const std::string &name, BoutReal def = 0.0)

Get a BoutReal from the input source

Return
zero if successful, non-zero on failure
Parameters
  • rval: The value will be put into this variable
  • name: The name of the variable to read
  • def: The default value if not found

int get(bool &bval, const std::string &name, bool def = false)

Get a bool from the input source

Return
zero if successful, non-zero on failure
Parameters
  • bval: The value will be put into this variable
  • name: The name of the variable to read
  • def: The default value if not found

int get(Field2D &var, const std::string &name, BoutReal def = 0.0)

Get a Field2D from the input source including communicating guard cells

Return
zero if successful, non-zero on failure
Parameters
  • var: This will be set to the value. Will be allocated if needed
  • name: Name of the variable to read
  • def: The default value if not found

int get(Field3D &var, const std::string &name, BoutReal def = 0.0, bool communicate = true)

Get a Field3D from the input source

Return
zero if successful, non-zero on failure
Parameters
  • var: This will be set to the value. Will be allocated if needed
  • name: Name of the variable to read
  • def: The default value if not found
  • communicate: Should the field be communicated to fill guard cells?

int get(FieldPerp &var, const std::string &name, BoutReal def = 0.0, bool communicate = true)

Get a FieldPerp from the input source

Return
zero if successful, non-zero on failure
Parameters
  • var: This will be set to the value. Will be allocated if needed
  • name: Name of the variable to read
  • def: The default value if not found
  • communicate: Should the field be communicated to fill guard cells?

int get(Vector2D &var, const std::string &name, BoutReal def = 0.0)

Get a Vector2D from the input source. If var is covariant then this gets three Field2D variables with “_x”, “_y”, “_z” appended to name If var is contravariant, then “x”, “y”, “z” are appended to name

By default all fields revert to zero

Return
zero always.
Parameters
  • var: This will be set to the value read
  • name: The name of the vector. Individual fields are read based on this name by appending. See above
  • def: The default value if not found (used for all the components)

int get(Vector3D &var, const std::string &name, BoutReal def = 0.0)

Get a Vector3D from the input source. If var is covariant then this gets three Field3D variables with “_x”, “_y”, “_z” appended to name If var is contravariant, then “x”, “y”, “z” are appended to name

By default all fields revert to zero

Return
zero always.
Parameters
  • var: This will be set to the value read
  • name: The name of the vector. Individual fields are read based on this name by appending. See above
  • def: The default value if not found (used for all the components)

bool isDataSourceGridFile() const

Test if input source was a grid file.

bool sourceHasVar(const std::string &name)

Wrapper for GridDataSource::hasVar.

bool sourceHasXBoundaryGuards()

Wrapper for GridDataSource::hasXBoundaryGuards.

bool sourceHasYBoundaryGuards()

Wrapper for GridDataSource::hasYBoundaryGuards.

template<typename ...Ts>
void communicate(Ts&... ts)

Communicate a list of FieldData objects Uses a variadic template (C++11) to pack all arguments into a FieldGroup

template<typename ...Ts>
void communicateXZ(Ts&... ts)
void communicate(FieldGroup &g)

Communicate a group of fields

void communicateXZ(FieldGroup &g)

Communcate guard cells in XZ only i.e. no Y communication

Parameters
  • g: The group of fields to communicate. Guard cells will be modified

void communicate(FieldPerp &f)

Communicate an X-Z field

This is a bit of a hack for now to get FieldPerp communications The FieldData class needs to be changed to accomodate FieldPerp objects

template<typename ...Ts>
comm_handle send(Ts&... ts)

Send a list of FieldData objects Packs arguments into a FieldGroup and passes to send(FieldGroup&).

virtual comm_handle send(FieldGroup &g) = 0

Perform communications without waiting for them to finish. Requires a call to wait() afterwards.

Return
handle to be used as input to wait()
Parameters
  • g: Group of fields to communicate

virtual int wait(comm_handle handle) = 0

Wait for the handle, return error code.

virtual MPI_Request sendToProc(int xproc, int yproc, BoutReal *buffer, int size, int tag) = 0

Low-level communication routine Send a buffer of data from this processor to another This must be matched by a corresponding call to receiveFromProc on the receiving processor

Parameters
  • xproc: X index of processor to send to
  • yproc: Y index of processor to send to
  • buffer: A buffer of data to send
  • size: The length of buffer
  • tag: A label, must be the same at receive

virtual comm_handle receiveFromProc(int xproc, int yproc, BoutReal *buffer, int size, int tag) = 0

Low-level communication routine Receive a buffer of data from another processor Must be matched by corresponding sendToProc call on the sending processor

Parameters
  • xproc: X index of sending processor
  • yproc: Y index of sending processor
  • buffer: The buffer to fill with data. Must already be allocated of length size
  • size: The length of buffer
  • tag: A label, must be the same as send

virtual int getNXPE() = 0

The number of processors in the X direction.

virtual int getNYPE() = 0

The number of processors in the Y direction.

virtual int getXProcIndex() = 0

This processor’s index in X direction.

virtual int getYProcIndex() = 0

This processor’s index in Y direction.

virtual bool firstX() = 0

Is this processor first in X? i.e. is there a boundary to the left in X?

virtual bool lastX() = 0

Is this processor last in X? i.e. is there a boundary to the right in X?

virtual int sendXOut(BoutReal *buffer, int size, int tag) = 0

Send a buffer of data to processor at X index +1

Parameters
  • buffer: The data to send. Must be at least length size
  • size: The number of BoutReals to send
  • tag: A label for the communication. Must be the same at receive

virtual int sendXIn(BoutReal *buffer, int size, int tag) = 0

Send a buffer of data to processor at X index -1

Parameters
  • buffer: The data to send. Must be at least length size
  • size: The number of BoutReals to send
  • tag: A label for the communication. Must be the same at receive

virtual comm_handle irecvXOut(BoutReal *buffer, int size, int tag) = 0

Receive a buffer of data from X index +1

Parameters
  • buffer: A buffer to put the data in. Must already be allocated of length size
  • size: The number of BoutReals to receive and put in buffer
  • tag: A label for the communication. Must be the same as sent

virtual comm_handle irecvXIn(BoutReal *buffer, int size, int tag) = 0

Receive a buffer of data from X index -1

Parameters
  • buffer: A buffer to put the data in. Must already be allocated of length size
  • size: The number of BoutReals to receive and put in buffer
  • tag: A label for the communication. Must be the same as sent

MPI_Comm getXcomm()

Return communicator containing all processors in X.

virtual MPI_Comm getXcomm(int jy) const = 0

Return X communicator.

virtual MPI_Comm getYcomm(int jx) const = 0

Return Y communicator.

bool periodicY(int jx) const

Is local X index jx periodic in Y?

Parameters
  • jx: The local (on this processor) index in X

virtual bool periodicY(int jx, BoutReal &ts) const = 0

Is local X index jx periodic in Y?

Parameters
  • jx: The local (on this processor) index in X
  • ts: The Twist-Shift angle if periodic

virtual std::pair<bool, BoutReal> hasBranchCutLower(int jx) const = 0

Is there a branch cut at this processor’s lower y-boundary?

Return
pair<bool, BoutReal> - bool is true if there is a branch cut, BoutReal gives the total zShift for a 2pi poloidal circuit if there is a branch cut
Parameters
  • jx: The local (on this processor) index in X

virtual std::pair<bool, BoutReal> hasBranchCutUpper(int jx) const = 0

Is there a branch cut at this processor’s upper y-boundary?

Return
pair<bool, BoutReal> - bool is true if there is a branch cut, BoutReal gives the total zShift for a 2pi poloidal circuit if there is a branch cut
Parameters
  • jx: The local (on this processor) index in X

int ySize(int jx) const

The number of points in Y at fixed X index jx.

virtual bool firstY() const = 0

Is this processor first in Y? i.e. is there a boundary at lower Y?

virtual bool lastY() const = 0

Is this processor last in Y? i.e. is there a boundary at upper Y?

virtual bool firstY(int xpos) const = 0

Is this processor first in Y? i.e. is there a boundary at lower Y?

virtual bool lastY(int xpos) const = 0

Is this processor last in Y? i.e. is there a boundary at upper Y?

virtual int UpXSplitIndex() = 0

If the upper Y guard cells are split in two, return the X index where the split occurs.

virtual int DownXSplitIndex() = 0

If the lower Y guard cells are split in two, return the X index where the split occurs.

virtual int sendYOutIndest(BoutReal *buffer, int size, int tag) = 0

Send data.

virtual int sendYOutOutdest(BoutReal *buffer, int size, int tag) = 0
virtual int sendYInIndest(BoutReal *buffer, int size, int tag) = 0
virtual int sendYInOutdest(BoutReal *buffer, int size, int tag) = 0
virtual comm_handle irecvYOutIndest(BoutReal *buffer, int size, int tag) = 0

Non-blocking receive. Must be followed by a call to wait()

Parameters
  • buffer: A buffer of length size which must already be allocated
  • size: The number of BoutReals expected
  • tag: The tag number of the expected message

virtual comm_handle irecvYOutOutdest(BoutReal *buffer, int size, int tag) = 0

Non-blocking receive. Must be followed by a call to wait()

Parameters
  • buffer: A buffer of length size which must already be allocated
  • size: The number of BoutReals expected
  • tag: The tag number of the expected message

virtual comm_handle irecvYInIndest(BoutReal *buffer, int size, int tag) = 0

Non-blocking receive. Must be followed by a call to wait()

Parameters
  • buffer: A buffer of length size which must already be allocated
  • size: The number of BoutReals expected
  • tag: The tag number of the expected message

virtual comm_handle irecvYInOutdest(BoutReal *buffer, int size, int tag) = 0

Non-blocking receive. Must be followed by a call to wait()

Parameters
  • buffer: A buffer of length size which must already be allocated
  • size: The number of BoutReals expected
  • tag: The tag number of the expected message

virtual const RangeIterator iterateBndryLowerY() const = 0

Iterate over the lower Y boundary.

virtual const RangeIterator iterateBndryUpperY() const = 0

Iterate over the upper Y boundary.

virtual const RangeIterator iterateBndryLowerOuterY() const = 0
virtual const RangeIterator iterateBndryLowerInnerY() const = 0
virtual const RangeIterator iterateBndryUpperOuterY() const = 0
virtual const RangeIterator iterateBndryUpperInnerY() const = 0
bool hasBndryLowerY()

Is there a boundary on the lower guard cells in Y?

bool hasBndryUpperY()

Is there a boundary on the upper guard cells in Y?

virtual std::vector<BoundaryRegion *> getBoundaries() = 0

Return a vector containing all the boundary regions on this processor.

virtual void addBoundary(BoundaryRegion *bndry)

Add a boundary region to this processor.

virtual std::vector<BoundaryRegionPar *> getBoundariesPar() = 0

Get all the parallel (Y) boundaries on this processor.

virtual void addBoundaryPar(BoundaryRegionPar *bndry)

Add a parallel(Y) boundary to this processor.

virtual const Field3D smoothSeparatrix(const Field3D &f)

Branch-cut special handling (experimental)

virtual BoutReal GlobalX(int jx) const = 0

Continuous X index between 0 and 1.

virtual BoutReal GlobalY(int jy) const = 0

Continuous Y index (0 -> 1)

virtual BoutReal GlobalX(BoutReal jx) const = 0

Continuous X index between 0 and 1.

virtual BoutReal GlobalY(BoutReal jy) const = 0

Continuous Y index (0 -> 1)

int XGLOBAL(int xloc) const

Returns the global X index given a local index If the local index includes the boundary cells, then so does the global.

virtual int YGLOBAL(int yloc) const

Returns the global Y index given a local index The local index must include the boundary, the global index does not.

virtual int XLOCAL(int xglo) const = 0

Returns the local X index given a global index If the global index includes the boundary cells, then so does the local.

virtual int YLOCAL(int yglo) const = 0

Returns the local Y index given a global index If the global index includes the boundary cells, then so does the local.

virtual int getGlobalXIndex(int xlocal) const = 0

Returns a global X index given a local index. Global index includes boundary cells, local index includes boundary or guard cells.

virtual int getGlobalXIndexNoBoundaries(int xlocal) const = 0

Returns a global X index given a local index. Global index excludes boundary cells, local index includes boundary or guard cells.

virtual int getGlobalYIndex(int ylocal) const = 0

Returns a global Y index given a local index. Global index includes boundary cells, local index includes boundary or guard cells.

virtual int getGlobalYIndexNoBoundaries(int ylocal) const = 0

Returns a global Y index given a local index. Global index excludes boundary cells, local index includes boundary or guard cells.

virtual int getGlobalZIndex(int zlocal) const = 0

Returns a global Z index given a local index. Global index includes boundary cells, local index includes boundary or guard cells.

virtual int getGlobalZIndexNoBoundaries(int zlocal) const = 0

Returns a global Z index given a local index. Global index excludes boundary cells, local index includes boundary or guard cells.

Coordinates *getCoordinates(const CELL_LOC location = CELL_CENTRE)

Coordinate system.

std::shared_ptr<Coordinates> getCoordinatesSmart(const CELL_LOC location = CELL_CENTRE)
CELL_LOC getAllowedStaggerLoc(DIRECTION direction) const

Returns the non-CELL_CENTRE location allowed as a staggered location

int getNpoints(DIRECTION direction) const

Returns the number of grid points in the particular direction

int getNguard(DIRECTION direction) const

Returns the number of guard points in the particular direction

void recalculateStaggeredCoordinates()

Re-calculate staggered Coordinates, useful if CELL_CENTRE Coordinates are changed.

STAGGER getStagger(const CELL_LOC inloc, const CELL_LOC outloc, const CELL_LOC allowedloc) const

Determines the resultant output stagger location in derivatives given the input and output location. Also checks that the combination of locations is allowed

STAGGER getStagger(const CELL_LOC vloc, const CELL_LOC inloc, const CELL_LOC outloc, const CELL_LOC allowedloc) const

Determines the resultant output stagger location in derivatives given the input and output location. Also checks that the combination of locations is allowed. This overload also checks the location of a second input field (velocity) is consistent.

template<typename T>
T indexDDX(const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexD2DX2(const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexD4DX4(const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexDDY(const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexD2DY2(const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexD4DY4(const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexDDZ(const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexD2DZ2(const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexD4DZ4(const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexVDDX(const T &vel, const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const

Advection operator in index space in [] direction

\[ v \frac{d}{di} f \]

Parameters
  • v: The velocity in the Y direction
  • f: The field being advected
  • outloc: The cell location where the result is desired. The default is the same as f
  • method: The differencing method to use
  • region: The region of the grid for which the result is calculated.

template<typename T>
T indexFDDX(const T &vel, const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexVDDY(const T &vel, const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexFDDY(const T &vel, const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexVDDZ(const T &vel, const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
template<typename T>
T indexFDDZ(const T &vel, const T &f, CELL_LOC outloc = CELL_DEFAULT, const std::string &method = "DEFAULT", REGION region = RGN_NOBNDRY) const
const Field3D toFieldAligned(const Field3D &f, const REGION region = RGN_ALL)
const Field3D fromFieldAligned(const Field3D &f, const REGION region = RGN_ALL)
const Field2D toFieldAligned(const Field2D &f, const REGION region = RGN_ALL)
const Field2D fromFieldAligned(const Field2D &f, const REGION region = RGN_ALL)
bool canToFromFieldAligned()
void setParallelTransform(std::unique_ptr<ParallelTransform> pt)
void setParallelTransform()
ParallelTransform &getParallelTransform()
const Region &getRegion(const std::string &region_name) const

Get the named region from the region_map for the data iterator

Throws if region_name not found

const Region &getRegion3D(const std::string &region_name) const
const Region<Ind2D> &getRegion2D(const std::string &region_name) const
const Region<IndPerp> &getRegionPerp(const std::string &region_name) const
bool hasRegion3D(const std::string &region_name) const

Indicate if named region has already been defined.

bool hasRegion2D(const std::string &region_name) const
bool hasRegionPerp(const std::string &region_name) const
void addRegion(const std::string &region_name, const Region<> &region)

Add a new region to the region_map for the data iterator

Outputs an error message if region_name already exists

void addRegion(const std::string &region_name, const Region<Ind2D> &region)
void addRegion(const std::string &region_name, const Region<IndPerp> &region)
void addRegion3D(const std::string &region_name, const Region<Ind3D> &region)
void addRegion2D(const std::string &region_name, const Region<Ind2D> &region)
void addRegionPerp(const std::string &region_name, const Region<IndPerp> &region)
Ind3D ind2Dto3D(const Ind2D &ind2D, int jz = 0)

Converts an Ind2D to an Ind3D using calculation.

Ind2D ind3Dto2D(const Ind3D &ind3D)

Converts an Ind3D to an Ind2D using calculation.

IndPerp ind3DtoPerp(const Ind3D &ind3D)

Converts an Ind3D to an IndPerp using calculation.

Ind3D indPerpto3D(const IndPerp &indPerp, int jy = 0)

Converts an IndPerp to an Ind3D using calculation.

Ind2D map3Dto2D(const Ind3D &ind3D)

Converts an Ind3D to an Ind2D representing a 2D index using a lookup to be used with care.

void createDefaultRegions()

Create the default regions for the data iterator

Creates RGN_{ALL,NOBNDRY,NOX,NOY}

Public Members

bool periodicX = {false}

Domain is periodic in X?

int NXPE
int PE_XIND

Number of processors in X, and X processor index.

int GlobalNx
int GlobalNy
int GlobalNz

Size of the global arrays. Note: can have holes.

int OffsetX
int OffsetY
int OffsetZ

Offset of this mesh within the global array so startx on this processor is OffsetX in global

int LocalNx

Size of the mesh on this processor including guard/boundary cells.

int LocalNy
int LocalNz
int xstart

Local ranges of data (inclusive), excluding guard cells.

int xend
int ystart
int yend
int zstart
int zend
bool StaggerGrids = {false}

Enable staggered grids (Centre, Lower). Otherwise all vars are cell centred (default).

bool IncIntShear = {false}

Include integrated shear (if shifting X)

int numberOfXPoints = {0}
BoutReal fft_derivs_filter = {0.0}

Fraction of modes to filter. This is set in derivs_init from option “ddz:fft_filter”.

int maxregionblocksize

Public Static Functions

Mesh *create(GridDataSource *source, Options *opt = nullptr)

Create a Mesh object

Parameters
  • source: The data source to use for loading variables
  • opt: The option section. By default this is “mesh”

Mesh *create(Options *opt = nullptr)

Create a Mesh object

The source is determined by 1) If “file” is set in the options, read that 2) If “grid” is set in global options, read that 3) Use options as data source

Parameters
  • opt: Input options. Default is “mesh” section

Protected Functions

const std::vector<int> readInts(const std::string &name, int n)

Read a 1D array of integers.

int msg_len(const std::vector<FieldData *> &var_list, int xge, int xlt, int yge, int ylt)

Calculates the size of a message for a given x and y range.

void derivs_init(Options *options)

Initialise derivatives.

Initialise the derivative methods. Must be called before any derivatives are used.

Protected Attributes

GridDataSource *source = {nullptr}

Source for grid data.

std::map<CELL_LOC, std::shared_ptr<Coordinates>> coords_map

Coordinate systems at different CELL_LOCs.

Options *options = {nullptr}

Mesh options section.

bool calcParallelSlices_on_communicate = {true}

Set whether to call calcParallelSlices on all communicated fields (true) or not (false)

Private Functions

std::shared_ptr<Coordinates> createDefaultCoordinates(const CELL_LOC location, bool force_interpolate_from_centre = false)

Allocates default Coordinates objects By default attempts to read staggered Coordinates from grid data source, interpolating from CELL_CENTRE if not present. Set force_interpolate_from_centre argument to true to always interpolate (useful if CELL_CENTRE Coordinates have been changed, so reading from file would not be correct).

Private Members

std::map<std::string, Region<Ind3D>> regionMap3D
std::map<std::string, Region<Ind2D>> regionMap2D
std::map<std::string, Region<IndPerp>> regionMapPerp
Array<int> indexLookup3Dto2D