File field3d.hxx

Defines

__FIELD3D_H__

Functions

FieldPerp operator+(const Field3D &lhs, const FieldPerp &rhs)
FieldPerp operator-(const Field3D &lhs, const FieldPerp &rhs)
FieldPerp operator*(const Field3D &lhs, const FieldPerp &rhs)
FieldPerp operator/(const Field3D &lhs, const FieldPerp &rhs)
Field3D operator+(const Field3D &lhs, const Field3D &rhs)
Field3D operator-(const Field3D &lhs, const Field3D &rhs)
Field3D operator*(const Field3D &lhs, const Field3D &rhs)
Field3D operator/(const Field3D &lhs, const Field3D &rhs)
Field3D operator+(const Field3D &lhs, const Field2D &rhs)
Field3D operator-(const Field3D &lhs, const Field2D &rhs)
Field3D operator*(const Field3D &lhs, const Field2D &rhs)
Field3D operator/(const Field3D &lhs, const Field2D &rhs)
Field3D operator+(const Field3D &lhs, BoutReal rhs)
Field3D operator-(const Field3D &lhs, BoutReal rhs)
Field3D operator*(const Field3D &lhs, BoutReal rhs)
Field3D operator/(const Field3D &lhs, BoutReal rhs)
Field3D operator+(BoutReal lhs, const Field3D &rhs)
Field3D operator-(BoutReal lhs, const Field3D &rhs)
Field3D operator*(BoutReal lhs, const Field3D &rhs)
Field3D operator/(BoutReal lhs, const Field3D &rhs)
Field3D operator-(const Field3D &f)

Unary minus. Returns the negative of given field, iterates over whole domain including guard/boundary cells.

Field3D pow(const Field3D &lhs, const Field2D &rhs, const std::string &rgn = "RGN_ALL")

Exponent: pow(lhs, lhs) is lhs raised to the power of rhs

Extra overloads not provided by the templates in field.hxx

This loops over the entire domain, including guard/boundary cells by default (can be changed using the rgn argument). If CHECK >= 3 then the result will be checked for non-finite numbers

Field3D pow(const Field3D &lhs, const Field2D &rhs, REGION rgn)
FieldPerp pow(const Field3D &lhs, const FieldPerp &rhs, const std::string &rgn = "RGN_ALL")
FieldPerp pow(const Field3D &lhs, const FieldPerp &rhs, REGION rgn)
void checkData(const Field3D &f, const std::string &region = "RGN_NOBNDRY")

Throw an exception if f is not allocated or if any elements are non-finite (for CHECK > 2). Loops over all points including the boundaries by default (can be changed using the rgn argument

void checkData(const Field3D &f, REGION region)
Field3D filter(const Field3D &var, int N0, const std::string &rgn = "RGN_ALL")

Fourier filtering, removes all except one mode

Parameters
  • var: Variable to apply filter to
  • N0: The component to keep
  • rgn: The region to calculate the result over

Field3D filter(const Field3D &var, int N0, REGION rgn)
Field3D lowPass(const Field3D &var, int zmax, bool keep_zonal, const std::string &rgn = "RGN_ALL")

Fourier low pass filtering. Removes modes higher than zmax and optionally the zonal component

Parameters
  • var: Variable to apply filter to
  • zmax: Maximum mode in Z
  • keep_zonal: Keep the zonal component if true
  • rgn: The region to calculate the result over

Field3D lowPass(const Field3D &var, int zmax, bool keep_zonal, REGION rgn)
Field3D lowPass(const Field3D &var, int zmax, int keep_zonal, REGION rgn = RGN_ALL)

The argument keep_zonal used to be integer “zmin” this was a misnomer. Please use the version above which uses a bool instead

Field3D lowPass(const Field3D &var, int zmax, const std::string rgn = "RGN_ALL")

Fourier low pass filtering. Removes modes higher than zmax

Parameters
  • var: Variable to apply filter to
  • zmax: Maximum mode in Z
  • rgn: The region to calculate the result over

Field3D lowPass(const Field3D &var, int zmax, REGION rgn)
void shiftZ(Field3D &var, int jx, int jy, double zangle)

Perform a shift by a given angle in Z

Parameters
  • var: The variable to be modified in-place
  • jx: X index
  • jy: Y index
  • zangle: The Z angle to apply

void shiftZ(Field3D &var, BoutReal zangle, const std::string &rgn = "RGN_ALL")

Apply a phase shift by a given angle zangle in Z to all points

Parameters
  • var: The variable to modify in-place
  • zangle: The angle to shift by in Z
  • rgn: The region to calculate the result over

void shiftZ(Field3D &var, BoutReal zangle, REGION rgn)
Field2D DC(const Field3D &f, const std::string &rgn = "RGN_ALL")

Average in the Z direction

Parameters
  • f: Variable to average
  • rgn: The region to calculate the result over

Field2D DC(const Field3D &f, REGION rgn)
void invalidateGuards(Field3D &var)

Force guard cells of passed field var to NaN.

Field3D &ddt(Field3D &f)

Returns a reference to the time-derivative of a field f

Wrapper around member function f.timeDeriv()

template<>
std::string toString(const Field3D &val)

toString template specialisation Defined in utils.hxx

bool operator==(const Field3D &a, const Field3D &b)

Test if two fields are the same, by calculating the minimum absolute difference between them

std::ostream &operator<<(std::ostream &out, const Field3D &value)

Output a string describing a Field3D to a stream.

class Field3D : public Field, public FieldData
#include <field3d.hxx>

Class for 3D X-Y-Z scalar fields.

This class represents a scalar field defined over the mesh. It handles memory management, and provides overloaded operators for operations on the data, iterators and access methods.

Initialisation

Fields can be declared in any scope (even global), but cannot be accessed by index or used until the data is allocated.

Field3D f;   // Declare variable, no data allocated
f(0,0,0) = 1.0; // Error !

f = 0.0;  // Allocates memory, fills with value (0.0)

Field3D g(1.0); // Declares, allocates memory, fills with value (1.0)

Field3D h;   // not allocated
h.allocate();  // Data array allocated, values undefined
f(0,0,0) = 1.0; // ok

Copy-on-Write

A field is a reference to the underlying data array, so setting one field equal to another has the effect of making both fields share the same underlying data

Field3D f(0.0);
Field3D g = f; // f and g now share data
f(0,0,0) = 1.0; // g is also modified

Setting the entire field equal to a new value changes the reference:

Field3D f(0.0);
Field3D g = f; // f and g now share data
g = 1.0;   // g and f are now separate

To ensure that a field is unique, call allocate() which will make a copy of the underlying data if it is shared.

Field3D f(0.0);
Field3D g = f; // f and g now share data
g.allocate();  // Data copied so g and f don't share data
f(0,0,0) = 1.0; // ok

Data access

Individual data indices can be accessed by index using round brackets:

Field3D f;
f(0,1,2) = 1.0;  // Set value
BoutReal val = f(2,1,3);  // Get value

If CHECK is greater than 2, this function will perform bounds checking. This will significantly slow calculations.

Some methods, such as FFT routines, need access to a pointer to memory. For the Z dimension this can be done by passing only the X and Y indices

BoutReal *data = f(0,1);

data now points to f(0,1,0) and can be incremented to move in Z.

Iteration

To loop over all points in a field, a for loop can be used to get the indices:

Field3D f(0.0); // Allocate, set to zero

for( const auto &i : f ) {  // Loop over all points, with index i
  f[i] = 1.0;
}

There is also more explicit looping over regions:

for( const auto &i : f.region(RGN_ALL) ) {  // Loop over all points, with index i
  f[i] = 1.0;
}

Parallel (y) derivatives

In several numerical schemes the mapping along magnetic fields (default y direction) is a relatively complex map. To accommodate this, the values of a field in the positive (up) and negative (down) directions can be stored in separate fields.

Field3D f(0.0); // f allocated, set to zero

f.yup() // error; f.yup not allocated

f.clearParallelSlices(); // f.yup_fields and f.ydown_fields are now empty
f.yup() // error; f.yup not allocated

To have separate fields for yup and ydown, first call

f.splitParallelSlices(); // f.yup() and f.ydown() separate

f.yup(); // ok
f.yup()(0,1,0) // error; f.yup not allocated

f.yup() = 1.0; // Set f.yup() field to 1.0

f.yup()(0,1,0) // ok

Unnamed Group

Field3D &operator=(const Field3D &rhs)

Assignment operators

Field3D &operator=(const Field2D &rhs)
void operator=(const FieldPerp &rhs)

return void, as only part initialised

Field3D &operator=(BoutReal val)

Unnamed Group

Field3D &operator+=(const Field3D &rhs)

Addition operators

Field3D &operator+=(const Field2D &rhs)
Field3D &operator+=(BoutReal rhs)

Unnamed Group

Field3D &operator-=(const Field3D &rhs)

Subtraction operators

Field3D &operator-=(const Field2D &rhs)
Field3D &operator-=(BoutReal rhs)

Unnamed Group

Field3D &operator*=(const Field3D &rhs)

Multiplication operators

Field3D &operator*=(const Field2D &rhs)
Field3D &operator*=(BoutReal rhs)

Unnamed Group

Field3D &operator/=(const Field3D &rhs)

Division operators

Field3D &operator/=(const Field2D &rhs)
Field3D &operator/=(BoutReal rhs)

Public Types

using ind_type = Ind3D
using value_type = BoutReal

Data type stored in this field.

Public Functions

Field3D::Field3D(Mesh * localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE, DirectionTypes directions_in = {YDirectionType::Standard, ZDirectionType::Standard})

Constructor.

Constructor

Note: the global “mesh” can’t be passed here because fields may be created before the mesh is.

Field3D(const Field3D &f)

Copy constructor

Doesn’t copy any data, just create a new reference to the same data (copy on change later)

Field3D(Field3D &&f)

Move constructor.

Field3D(const Field2D &f)

Constructor from 2D field.

Field3D(BoutReal val, Mesh *localmesh = nullptr)

Constructor from value.

Field3D::Field3D(Array < BoutReal > data, Mesh * localmesh, CELL_LOC location = CELL_CENTRE, DirectionTypes directions_in = {YDirectionType::Standard, ZDirectionType::Standard})

Constructor from Array and Mesh.

~Field3D()

Destructor.

Field3D &allocate()

Ensures that memory is allocated and unique

bool isAllocated() const

Test if data is allocated

Field3D *timeDeriv()

Return a pointer to the time-derivative field

The first time this is called, a new field will be allocated. Subsequent calls return the same field

int getNx() const

Return the number of nx points

int getNy() const

Return the number of ny points

int getNz() const

Return the number of nz points

Field3D &setLocation(CELL_LOC location)
Field3D &setDirectionY(YDirectionType d)
void splitParallelSlices()

Ensure that this field has separate fields for yup and ydown.

void splitYupYdown()
void clearParallelSlices()

Clear the parallel slices, yup and ydown

void mergeYupYdown()
bool hasParallelSlices() const

Check if this field has yup and ydown fields.

bool hasYupYdown() const
Field3D &yup(std::vector<Field3D>::size_type index = 0)

Check if this field has yup and ydown fields Return reference to yup field

const Field3D &yup(std::vector<Field3D>::size_type index = 0) const

Return const reference to yup field.

Field3D &ydown(std::vector<Field3D>::size_type index = 0)

Return reference to ydown field.

const Field3D &ydown(std::vector<Field3D>::size_type index = 0) const

Return const reference to ydown field.

Field3D &ynext(int offset)

Return the parallel slice at offset

offset of 0 returns the main field itself

const Field3D &ynext(int offset) const
bool requiresTwistShift(bool twist_shift_enabled)

If twist_shift_enabled is true, does this Field3D require a twist-shift at branch cuts on closed field lines?

const Region<Ind3D> &getRegion(REGION region) const

Return a Region<Ind3D> reference to use to iterate over this field

Example

This loops over the interior points, not the boundary and inside the loop the index is used to calculate the difference between the point one index up in x (i.xp()) and one index down in x (i.xm()), putting the result into a different field ‘g’

for(const auto &i : f.getRegion(RGN_NOBNDRY)) { g[i] = f[i.xp()] - f[i.xm()]; }

const Region<Ind3D> &getRegion(const std::string &region_name) const
const Region<Ind2D> &getRegion2D(REGION region) const

Return a Region<Ind2D> reference to use to iterate over the x- and y-indices of this field

const Region<Ind2D> &getRegion2D(const std::string &region_name) const
Region<Ind3D>::RegionIndices::const_iterator begin() const
Region<Ind3D>::RegionIndices::const_iterator end() const
BoutReal &operator[](const Ind3D &d)
const BoutReal &operator[](const Ind3D &d) const
BoutReal &operator()(const IndPerp &d, int jy)
const BoutReal &operator()(const IndPerp &d, int jy) const
BoutReal &operator()(const Ind2D &d, int jz)
const BoutReal &operator()(const Ind2D &d, int jz) const
BoutReal &operator()(int jx, int jy, int jz)

Direct access to the underlying data array

If CHECK > 2 then bounds checking is performed

If CHECK <= 2 then no checks are performed, to allow inlining and optimisation of inner loops

const BoutReal &operator()(int jx, int jy, int jz) const
const BoutReal *operator()(int jx, int jy) const

Direct access to the underlying data array

This version returns a pointer to a data array, and is intended for use with FFT routines. The data is guaranteed to be contiguous in Z index

BoutReal *operator()(int jx, int jy)
bool isReal() const

Returns true if field consists of BoutReal values.

bool is3D() const

True if variable is 3D.

int byteSize() const

Number of bytes for a single point.

int BoutRealSize() const

Number of BoutReals (not implemented if not BoutReal)

void accept(FieldVisitor &v)

Visitor pattern support.

void doneComms()
Field3D &calcParallelSlices()
void applyBoundary(bool init = false)
void applyBoundary(BoutReal t)
void applyBoundary(const std::string &condition)
void applyBoundary(const char *condition)
void applyBoundary(const std::string &region, const std::string &condition)
void applyTDerivBoundary()
void setBoundaryTo(const Field3D &f3d)

Copy the boundary values half-way between cells This uses 2nd order central differences to set the value on the boundary to the value on the boundary in field f3d. Note: does not just copy values in boundary region.

void applyParallelBoundary()
void applyParallelBoundary(BoutReal t)
void applyParallelBoundary(const std::string &condition)
void applyParallelBoundary(const char *condition)
void applyParallelBoundary(const std::string &region, const std::string &condition)
void applyParallelBoundary(const std::string &region, const std::string &condition, Field3D *f)

Private Members

const Field2D *background = {nullptr}

Boundary - add a 2D field.

int nx = {-1}

Array sizes (from fieldmesh). These are valid only if fieldmesh is not null.

int ny = {-1}
int nz = {-1}
Array<BoutReal> data

Internal data array. Handles allocation/freeing of memory.

Field3D *deriv = {nullptr}

Time derivative (may be nullptr)

std::vector<Field3D> yup_fields = {}

Fields containing values along Y.

std::vector<Field3D> ydown_fields = {}

Friends

friend Vector3D
void swap(Field3D &first, Field3D &second)