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

template<class DerivedType>
using RegisterMesh = MeshFactory::RegisterInFactory<DerivedType>#
using comm_handle = void*#

Type used to return pointers to handles.

class MeshFactory : public Factory<Mesh, MeshFactory, GridDataSource*, Options*>#

Public Functions

ReturnType create(const std::string &type, Options *options = nullptr, GridDataSource *source = nullptr) const#
ReturnType create(Options *options = nullptr, GridDataSource *source = nullptr) const#

Public Static Attributes

static constexpr auto type_name = "Mesh"#
static constexpr auto section_name = "mesh"#
static constexpr auto option_name = "type"#
static constexpr auto default_type = "bout"#
class Mesh#

Subclassed by BoutMesh

Public Functions

inline Mesh()#

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

Mesh(GridDataSource *s, Options *options)#

Constructor

Parameters:
  • s[in] The source to be used for loading variables

  • options[in] The options section for settings

virtual ~Mesh()#

Destructor.

inline 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

inline virtual void outputVars([[maybe_unused]] Options &output_options)#

Add output variables to output_options 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

Parameters:
  • sval[out] The value will be put into this variable

  • name[in] The name of the variable to read

  • def[in] The default value if not found

Returns:

zero if successful, non-zero on failure

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

Get an integer from the input source

Parameters:
  • ival[out] The value will be put into this variable

  • name[in] The name of the variable to read

  • def[in] The default value if not found

Returns:

zero if successful, non-zero on failure

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

Get a BoutReal from the input source

Parameters:
  • rval[out] The value will be put into this variable

  • name[in] The name of the variable to read

  • def[in] The default value if not found

Returns:

zero if successful, non-zero on failure

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

Get a bool from the input source

Parameters:
  • bval[out] The value will be put into this variable

  • name[in] The name of the variable to read

  • def[in] The default value if not found

Returns:

zero if successful, non-zero on failure

int get(Field2D &var, const std::string &name, BoutReal def = 0.0, bool communicate = true, CELL_LOC location = CELL_DEFAULT)#

Get a Field2D from the input source including communicating guard cells

Parameters:
  • var[out] This will be set to the value. Will be allocated if needed

  • name[in] Name of the variable to read

  • def[in] The default value if not found

  • communicate[in] Should the field be communicated to fill guard cells?

Returns:

zero if successful, non-zero on failure

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

Get a Field3D from the input source

Parameters:
  • var[out] This will be set to the value. Will be allocated if needed

  • name[in] Name of the variable to read

  • def[in] The default value if not found

  • communicate[in] Should the field be communicated to fill guard cells?

Returns:

zero if successful, non-zero on failure

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

Get a FieldPerp from the input source

Parameters:
  • var[out] This will be set to the value. Will be allocated if needed

  • name[in] Name of the variable to read

  • def[in] The default value if not found

  • communicate[in] Should the field be communicated to fill guard cells?

Returns:

zero if successful, non-zero on failure

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

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

Parameters:
  • var[in] This will be set to the value read

  • name[in] The name of the vector. Individual fields are read based on this name by appending. See above

  • def[in] The default value if not found (used for all the components)

  • communicate[in] Should the field be communicated to fill guard cells?

Returns:

zero always.

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

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

Parameters:
  • var[in] This will be set to the value read

  • name[in] The name of the vector. Individual fields are read based on this name by appending. See above

  • def[in] The default value if not found (used for all the components)

  • communicate[in] Should the field be communicated to fill guard cells?

Returns:

zero always.

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>
inline 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>
inline void communicateXZ(Ts&... ts)#
template<typename ...Ts>
inline void communicateYZ(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 communicateYZ(FieldGroup &g)#

Communcate guard cells in YZ only i.e. no X communication

Parameters:

g – The group of fields to communicate. Guard cells will be modified

virtual 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>
inline comm_handle send(Ts&... ts)#

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

template<typename ...Ts>
inline comm_handle sendX(Ts&... ts)#

Send guard cells from a list of FieldData objects in the x-direction Packs arguments into a FieldGroup and passes to send(FieldGroup&). Perform communications without waiting for them to finish. Requires a call to wait() afterwards.

template<typename ...Ts>
inline comm_handle sendY(Ts&... ts)#

Send guard cells from a list of FieldData objects in the y-direction 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.

Parameters:

g – Group of fields to communicate

Returns:

handle to be used as input to wait()

virtual comm_handle sendX(FieldGroup &g, comm_handle handle = nullptr, bool disable_corners = false) = 0#

Send only the x-guard cells.

virtual comm_handle sendY(FieldGroup &g, comm_handle handle = nullptr) = 0#

Send only the y-guard cells.

virtual int wait(comm_handle handle) = 0#

Wait for the handle, return error code.

Wait for the handle, return error code

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() const = 0#

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

virtual bool lastX() const = 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[in] The data to send. Must be at least length size

  • size[in] The number of BoutReals to send

  • tag[in] 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[in] The data to send. Must be at least length size

  • size[in] The number of BoutReals to send

  • tag[in] 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[in] A buffer to put the data in. Must already be allocated of length size

  • size[in] The number of BoutReals to receive and put in buffer

  • tag[in] 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[in] A buffer to put the data in. Must already be allocated of length size

  • size[in] The number of BoutReals to receive and put in buffer

  • tag[in] A label for the communication. Must be the same as sent

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

inline MpiWrapper &getMpi()#

Return pointer to the mesh’s MPI Wrapper object.

virtual bool periodicY(int jx) const#

Is local X index jx periodic in Y?

Parameters:

jx[in] 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[in] The local (on this processor) index in X

  • ts[out] The Twist-Shift angle if periodic

virtual int numberOfYBoundaries() const = 0#

Get number of boundaries in the y-direction, i.e. locations where there are boundary cells in the global grid

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

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

Parameters:

jx[in] The local (on this processor) index in X

Returns:

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

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

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

Parameters:

jx[in] The local (on this processor) index in X

Returns:

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

virtual 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? Note: First on the global grid, not necessarily at a boundary

virtual bool lastY() const = 0#

Is this processor last in Y? Note: Last on the global grid, not necessarily at a boundary

virtual bool firstY(int xpos) const = 0#

Is this processor first in Y? Note: Not necessarily at a boundary, but first in the Y communicator for the flux surface through local X index xpos

virtual bool lastY(int xpos) const = 0#

Is this processor last in Y? Note: Not necessarily at a boundary, but last in the Y communicator for the flux surface through local X index xpos

virtual RangeIterator iterateBndryLowerY() const = 0#

Iterate over the lower Y boundary.

virtual RangeIterator iterateBndryUpperY() const = 0#

Iterate over the upper Y boundary.

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

Is there a boundary on the lower guard cells in Y on any processor along the X direction?

bool hasBndryUpperY()#

Is there a boundary on the upper guard cells in Y on any processor along the X direction?

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

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

inline virtual std::set<std::string> getPossibleBoundaries() const#

Get the set of all possible boundaries in this configuration.

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

inline virtual void addBoundaryPar(BoundaryRegionPar *bndry)#

Add a parallel(Y) boundary to this processor.

inline virtual 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)

virtual int localSize3D()#

Returns the number of unique cells (i.e., ones not used for communication) on this processor for 3D fields. Boundaries are only included to a depth of 1.

virtual int localSize2D()#

Returns the number of unique cells (i.e., ones not used for communication) on this processor for 2D fields. Boundaries are only included to a depth of 1.

virtual int localSizePerp()#

Returns the number of unique cells (i.e., ones not used for communication) on this processor for perpendicular fields. Boundaries are only included to a depth of 1.

virtual int globalStartIndex3D()#

Get the value of the first global 3D index on this processor.

virtual int globalStartIndex2D()#

Get the value of the first global 2D index on this processor.

virtual int globalStartIndexPerp()#

Get the value of the first global perpendicular index on this processor.

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 getLocalXIndex(int xglobal) const = 0#

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

virtual int getLocalXIndexNoBoundaries(int xglobal) const = 0#

Returns a local X index given a global 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 getLocalYIndex(int yglobal) const = 0#

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

virtual int getLocalYIndexNoBoundaries(int yglobal) const = 0#

Returns a local Y index given a global 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.

virtual int getLocalZIndex(int zglobal) const = 0#

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

virtual int getLocalZIndexNoBoundaries(int zglobal) const = 0#

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

inline Coordinates *getCoordinates(const CELL_LOC location = CELL_CENTRE)#

Coordinate system.

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

Returns the non-CELL_CENTRE location allowed as a staggered location

inline int getNpoints(DIRECTION direction) const#

Returns the number of grid points in the particular direction

inline 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<class T>
const Region<typename T::ind_type> &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

inline const Region &getRegion(const std::string &region_name) const#
const Region<Ind3D> &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#
inline 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

inline void addRegion(const std::string &region_name, const Region<Ind2D> &region)#
inline 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)#
inline Ind3D ind2Dto3D(const Ind2D &ind2D, int jz = 0)#

Converts an Ind2D to an Ind3D using calculation.

inline Ind2D ind3Dto2D(const Ind3D &ind3D)#

Converts an Ind3D to an Ind2D using calculation.

inline IndPerp ind3DtoPerp(const Ind3D &ind3D)#

Converts an Ind3D to an IndPerp using calculation.

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

Converts an IndPerp to an Ind3D using calculation.

inline BOUT_HOST_DEVICE Ind2D map3Dto2D (const Ind3D &ind3D)

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

void createDefaultRegions()#

Create the default regions for the data iterator

Creates RGN_{ALL,NOBNDRY,NOX,NOY}

std::optional<size_t> getCommonRegion(std::optional<size_t>, std::optional<size_t>)#
size_t getRegionID(const std::string &region) const#
inline const Region<Ind3D> &getRegion(size_t RegionID) const#
inline const Region<Ind3D> &getRegion(std::optional<size_t> RegionID) const#
template<>
inline const Region<Ind3D> &getRegion(const std::string &region_name) const#
template<>
inline const Region<Ind2D> &getRegion(const std::string &region_name) const#
template<>
inline const Region<IndPerp> &getRegion(const std::string &region_name) const#

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 GlobalNxNoBoundaries#

Size of the global arrays excluding boundary points.

int GlobalNyNoBoundaries#
int GlobalNzNoBoundaries#
int OffsetX#

Note: These offsets only correct if Y guards are not included in the global array and are corrected in gridfromfile.cxx

int OffsetY#
int OffsetZ#

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

int MapGlobalX#

Map between local and global indices (MapGlobalX, MapGlobalY, MapGlobalZ) in the global index space maps to (MapLocalX, MapLocalY, MapLocalZ) locally. Note that boundary cells are included in the global index space, but communication guard cells are not.

int MapGlobalY#
int MapGlobalZ#

Start global indices.

int MapLocalX#
int MapLocalY#
int MapLocalZ#

Start local indices.

int MapCountX#
int MapCountY#
int MapCountZ#

Size of the mapped region.

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 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 = {MAXREGIONBLOCKSIZE}#
bool StaggerGrids = {false}#

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

const bool include_corner_cells#

Public Static Functions

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

Create a Mesh object

Parameters:
  • source[in] The data source to use for loading variables

  • opt[in] The option section. By default this is “mesh”

static 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[in] 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)

MpiWrapper *mpi = nullptr#

Pointer to the global MPI wrapper, for convenience.

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, size_t> regionMap3D#
std::vector<Region<Ind3D>> region3D#
std::vector<std::optional<size_t>> region3Dintersect#
std::map<std::string, Region<Ind2D>> regionMap2D#
std::map<std::string, Region<IndPerp>> regionMapPerp#
Array<int> indexLookup3Dto2D#
int localNumCells3D = -1#
int localNumCells2D = -1#
int localNumCellsPerp = -1#