Boundary conditions#

Boundary conditions for non-orthogonal grids#

If non-orthogonal grids are used (meaning that the x- and y-directions are not orthogonal, so g12 != 0.), then corner cells may be required. The boundary conditions are applied in corner cells[#disablecorners]_ by applying the y-boundary condition using x-boundary values. This requires that x-boundary conditions are applied before y-boundary conditions. The ordering is taken care of by the methods described in this section, but also needs to be respected by any custom boundary conditions in user code (e.g. sheath boundary conditions). Note that the iterators returned by the BoutMesh methods iterateBndryLowerY, iterateBndryLowerInnerY, iterateBndryLowerOuterY, iterateBndryUpperY, iterateBndryUpperInnerY, and iterateBndryUpperOuterY do include the corner cells at the domain boundary corners.

Parallel boundary conditions#

Unless using slab geometry (with ParallelTransformIdentity, see Parallel Transforms), some special handling is needed to apply parallel boundary conditions. The details depend on the parallel derivative scheme being used, see below. The default bndry_yup and bndry_ydown settings would apply boundary conditions in the poloidal, rather than parallel direction. As the poloidal direction generally has a coarse resolution, that is not sufficient to resolve perpendicular gradients, this would result in large numerical inaccuracies in the boundary conditions.

Shifted metric boundary conditions#

When using the ShiftedMetric implementation of ParallelTransform (by setting mesh:paralleltransform = shifted, see Shifted metric), the recommended method is to apply boundary conditions directly to the yup and ydown parallel slices. This can be done by setting bndry_par_yup and bndry_par_ydown, or bndry_par_all to set both at once. The possible values are parallel_dirichlet_o1, parallel_dirichlet_o2, parallel_dirichlet_o3, parallel_neumann_o1, parallel_neumann_o2 and parallel_neumann_o3. The stencils used are the same as for the standard boundary conditions without the parallel_ prefix, but are applied directly to parallel slices. The boundary condition can only be applied after the parallel slices are calculated, which is usually done during a call to Mesh::communicate(), so the applyBoundary() method must be called explicitly (when boundary conditions are applied automatically to evolving variables, they cannot set these parallel boundary conditions). For maximum efficiency, set bndry_yup and bndry_ydown to none to skip using any boundary condition to set the unused boundary cells of the base variable.

For example, for an evolving variable f, put a section in the BOUT.inp input file like

[f]
bndry_xin = dirichlet
bndry_xout = dirichlet
bndry_par_all = parallel_neumann_o2
bndry_ydown = none
bndry_yup = none

and in the PhysicsModel::rhs() function, before taking any derivatives of f, call

mesh->communicate(f);
f.applyBoundary();

The bndry_par_* options only provide a subset of boundary conditions. If others are required, they can be used with a different, slightly less optimised method. The modifier fromFieldAligned() applies a boundary condition by first transforming the variable to a globally field aligned grid, then applying the boundary condition, then transforming back to the toroidal grid. When this method is used, the boundary conditions must be applied before communicating, so that the parallel slices are calculated using the boundary cells of the base variable (for variables that have been added to the time solver, this will automatically be the case). For example, the settings in BOUT.inp for a Robin parallel boundary condition could be

[f]
bndry_xin = dirichlet
bndry_xout = dirichlet
bndry_yup = fromFieldAligned(robin(1, -1, 0))
bndry_ydown = fromFieldAligned(robin(1, 1, 0))

Aligned transform boundary conditions#

When using the ‘aligned transform’ method for parallel derivatives (see Aligned transform), the way to apply parallel boundary conditions depends on how the method was implemented.

For the ‘implicit transform’ version where the transformations to and from the field aligned grid are done within each parallel derivative or interpolation operator, the parallel boundary conditions must be applied to the base variable, so they must be applied using the fromFieldAligned() modifier, as described in the previous section (Shifted metric boundary conditions).

For the optimised method with separate objects for the field aligned versions of variables, it would be correct to apply boundary conditions using the fromFieldAligned() modifier before calculating the field aligned versions, but would add extra interpolations. Therefore the recommended way to apply parallel boundary conditions is to apply them directly to the field aligned versions of variables. Since the objects for the field aligned versions are not added to the time solver, it is necessary to load boundary conditions explicitly from the BOUT.inp input file during PhysicsModel::init(), for example by calling:

f_aligned.setBoundary("f_aligned")

where the argument to setBoundary() specifies the name of the section in BOUT.inp from which boundary conditions will be read. Then the boundary conditions must be applied explicitly after the field aligned object has been calculated in PhysicsModel::rhs(), for example:

f_aligned = fromFieldAligned(f);
f_aligned.applyBoundary();

The boundary condition should be applied directly to the array in f_aligned (not to parallel slices, which are not created for this scheme), so uses the ‘standard’ bndry_yup/bndry_ydown. Radial boundary points should never be used from the aligned object, so its x-boundaries should be set to none, and parallel boundary points should never be used from the base variable, so its y-boundaries should be set to none. The input sections for f and f_aligned might look like

[f]
bndry_xin = dirichlet
bndry_xout = dirichlet
bndry_yup = none
bndry_ydown = none

[f_aligned]
bndry_xin = none
bndry_xout = none
bndry_yup = free_o3
bndry_ydown = free_o3

FCI boundary conditions#

When using the FCI method (FCI method), parallel boundary conditions must be applied to the parallel slices using bndry_par_yup and bndry_par_ydown, or bndry_par_all to set both together. It is suggested, at least if there are boundaries in the y-direction of the grid, to set bndry_yup = none and bndry_down = none to skip unnecessary operations on y-boundary cells of the base variable. For example, for an evolving variable f, put a section in the BOUT.inp input file like

[f]
bndry_xin = dirichlet
bndry_xout = dirichlet
bndry_par_all = parallel_dirichlet_o2
bndry_ydown = none
bndry_yup = none

One should not that the parallel boundary conditions have to be applied after communication, while the perpendicular ones before:

f.applyBoundary();
mesh->communicate(f);
f.applyParallelBoundary("parallel_neumann_o2");

Note that during grid generation care has to be taken to ensure that there are no “short” connection lengths. Otherwise it can happen that for a point on a slice, both yup() and ydown() are boundary cells, and interpolation into the boundary can only use the single point on the given cell.

Relaxing boundaries#

All boundaries can be modified to be “relaxing” which are a combination of zero-gradient time-derivative, and whatever boundary condition they are applied to. The idea is that this prevents sharp discontinuities at boundaries during transients, whilst maintaining the desired boundary condition on longer time-scales. In some cases this can improve the numerical stability and timestep.

For example, relax(dirichlet) will make a field \(f\) at point \(i\) in the boundary follow a point \(i-1\) in the domain:

\[.{{\frac{\partial f}{\partial t}}}|_i = .{{\frac{\partial f}{\partial t}}}|_{i-1} - f_i / \tau\]

where \(\tau\) is a time-scale for the boundary (currently set to 0.1, but will be a global option). When the time-derivatives are slow close to the boundary, the boundary relaxes to the desired condition (Dirichlet in this case), but when the time-derivatives are large then the boundary approaches Neumann to reduce discontinuities.

By default, the relaxation rate is set to \(10\) (i.e. a time-scale of \(\tau=0.1\)). To change this, give the rate as the second argument e.g. relax(dirichlet, 2) would relax to a Dirichlet boundary condition at a rate of \(2\).

Changing the width of boundaries#

To change the width of a boundary region, the width modifier changes the width of a boundary region before applying the boundary condition, then changes the width back afterwards. To use, specify the boundary condition and the width, for example

bndry_core = width( neumann , 4 )

would apply a Neumann boundary condition on the innermost 4 cells in the core, rather than the usual 2. When combining with other boundary modifiers, this should be applied first e.g.

bndry_sol = width( relax( dirichlet ), 3)

would relax the last 3 cells towards zero, whereas

bndry_sol = relax( width( dirichlet, 3) )

would only apply to the usual 2, since relax didn’t use the updated width.

Limitations:

Because it modifies then restores a globally-used BoundaryRegion, this code is not thread safe.
Boundary conditions can’t be applied across processors, and no checks are done that the width asked for fits within a single processor.

Examples#

This example is taken from the UEDGE benchmark test (in examples/uedge-benchmark):

[All]
bndry_all = neumann # Default for all variables, boundaries

[Ni]
bndry_target = neumann
bndry_core = relax(dirichlet(1.))   # 1e13 cm^-3 on core boundary
bndry_all  = relax(dirichlet(0.1))  # 1e12 cm^-3 on other boundaries

[Vi]
bndry_ydown = relax(dirichlet(-1.41648))   # -3.095e4/Vi_x
bndry_yup   = relax(dirichlet( 1.41648))

The variable Ni (density) is set to a Neumann boundary condition on the targets (yup and ydown), relaxes towards \(1\) on the core boundary, and relaxes to \(0.1\) on all other boundaries. Note that the bndry_target = neumann needs to be in the Ni section: If we just had

[All]
bndry_all = neumann # Default for all variables, boundaries

[Ni]
bndry_core = relax(dirichlet(1.))   # 1e13 cm^-3 on core boundary
bndry_all  = relax(dirichlet(0.1))  # 1e12 cm^-3 on other boundaries

then the “target” boundary condition for Ni would first search in the [Ni] section for bndry_target, then for bndry_all in the [Ni] section. This is set to relax(dirichlet(0.1)), not the Neumann condition desired.

Boundary regions#

The boundary condition code needs ways to loop over the boundary regions, without needing to know the details of the mesh.

At the moment two mechanisms are provided: A RangeIterator over upper and lower Y boundaries, and a vector of BoundaryRegion objects.

// Boundary region iteration
virtual const RangeIterator iterateBndryLowerY() const = 0;
virtual const RangeIterator iterateBndryUpperY() const = 0;

bool hasBndryLowerY();
bool hasBndryUpperY();

bool BoundaryOnCell; // NB: DOESN'T REALLY BELONG HERE

The RangeIterator class is an iterator which allows looping over a set of indices. For example, in src/solver/solver.cxx to loop over the upper Y boundary of a 2D variable var:

for(RangeIterator xi = mesh->iterateBndryUpperY(); !xi.isDone(); xi++) {
  ...
}

The BoundaryRegion class is defined in include/boundary_region.hxx

Y-Boundaries#

The sheath boundaries are often implemented in the physics model. Previously of they where implemented using a RangeIterator:

class yboundary_example_legacy {
public:
  yboundary_example_legacy(Options* opt, const Field3D& N, const Field3D& V)
      : N(N), V(V) {
    Options& options = *opt;
    lower_y = options["lower_y"].doc("Boundary on lower y?").withDefault<bool>(lower_y);
    upper_y = options["upper_y"].doc("Boundary on upper y?").withDefault<bool>(upper_y);
  }

  void rhs() {
    BoutReal totalFlux = 0;
    if (lower_y) {
      for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) {
        for (int jz = 0; jz < mesh->LocalNz; jz++) {
          // Calculate flux through surface [normalised m^-2 s^-1],
          // should be positive since V < 0.0
          BoutReal flux =
              -0.5 * (N(r.ind, mesh->ystart, jz) + N(r.ind, mesh->ystart - 1, jz)) * 0.5
              * (V(r.ind, mesh->ystart, jz) + V(r.ind, mesh->ystart - 1, jz));
          totalFlux += flux;
        }
      }
    }
    if (upper_y) {
      for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) {
        for (int jz = 0; jz < mesh->LocalNz; jz++) {
          // Calculate flux through surface [normalised m^-2 s^-1],
          // should be positive since V < 0.0
          BoutReal flux = -0.5 * (N(r.ind, mesh->yend, jz) + N(r.ind, mesh->yend + 1, jz))
                          * 0.5
                          * (V(r.ind, mesh->yend, jz) + V(r.ind, mesh->yend + 1, jz));
          totalFlux += flux;
        }
      }
    }
  }

private:
  bool lower_y{true};
  bool upper_y{true};
  const Field3D& N;
  const Field3D& V;
}

This can be replaced using the YBoundary class, which not only simplifies the code, but also allows to have the same code working with non-field-aligned geometries, as flux coordinate independent (FCI) method:

#include <bout/yboundary_regions.hxx>

class yboundary_example {
public:
  yboundary_example(Options* opt, const Field3D& N, const Field3D& V) : N(N), V(V) {
    // Set what kind of yboundaries you want to include
    yboundary.init(opt);
  }

  void rhs() {
    BoutReal totalFlux = 0;
    yboundary.iter_pnts([&](auto& pnt) {
      BoutReal flux = pnt.interpolate_sheath_o2(N) * pnt.interpolate_sheath_o2(V);
    });
  }

private:
  YBoundary ybounday;
  const Field3D& N;
  const Field3D& V;
};

There are several member functions of pnt. pnt is of type BoundaryRegionParIterBase and BoundaryRegionIter, and both should provide the same interface. If they don’t that is a bug, as the above code is a template, that gets instantiated for both types, and thus requires both classes to provide the same interface, one for FCI-like boundaries and one for field aligned boundaries.

Here is a short summary of some members of pnt, where f is a :

Table 4 Members for boundary operation#
Function	Description
`pnt.ythis(f)`	Returns the value at the last point in the domain
`pnt.ynext(f)`	Returns the value at the first point in the domain
`pnt.yprev(f)`	Returns the value at the second to last point in the domain, if it is valid. NB: this point may not be valid.
`pnt.interpolate_sheath_o2(f)`	Returns the value at the boundary, assuming the bounday value has been set
`pnt.extrapolate_sheath_o1(f)`	Returns the value at the boundary, extrapolating from the bulk, first order
`pnt.extrapolate_sheath_o2(f)`	Returns the value at the boundary, extrapolating from the bulk, second order
`pnt.extrapolate_next_o{1,2}(f)`	Extrapolate into the boundary from the bulk, first or second order
`pnt.extrapolate_grad_o{1,2}(f)`	Extrapolate the gradient into the boundary, first or second order
`pnt.dirichlet_o{1,2,3}(f, v)`	Apply dirichlet boundary conditions with value `v` and given order
`pnt.neumann_o{1,2,3}(f, v)`	Applies a gradient of `v / dy` boundary condition.
`pnt.limitFree(f)`	Extrapolate into the boundary using only monotonic decreasing values. `f` needs to be positive.
`pnt.dir`	The direction of the boundary.

Boundary regions#

Different regions of the boundary such as “core”, “sol” etc. are labelled by the Mesh class (i.e. BoutMesh), which implements a member function defined in mesh.hxx:

// Boundary regions
virtual vector<BoundaryRegion*> getBoundaries() = 0;

This returns a vector of pointers to BoundaryRegion objects, each of which describes a boundary region with a label, a BndryLoc location (i.e. inner x, outer x, lower y, upper y or all), and iterator functions for looping over the points. This class is defined in boundary_region.hxx:

/// Describes a region of the boundary, and a means of iterating over it
class BoundaryRegion {
  public:
  BoundaryRegion();
  BoundaryRegion(const string &name, int xd, int yd);
  virtual ~BoundaryRegion();

  string label; // Label for this boundary region

  BndryLoc location; // Which side of the domain is it on?

  int x,y; // Indices of the point in the boundary
  int bx, by; // Direction of the boundary [x+dx][y+dy] is going outwards

  virtual void first() = 0;
  virtual void next() = 0; // Loop over every element from inside out (in X or
Y first)
  virtual void nextX() = 0; // Just loop over X
  virtual void nextY() = 0; // Just loop over Y
  virtual bool isDone() = 0; // Returns true if outside domain. Can use this
with nested nextX, nextY
};

Example: To loop over all points in BoundaryRegion *bndry , use

for(bndry->first(); !bndry->isDone(); bndry->next()) {
  ...
}

Inside the loop, bndry->x and bndry->y are the indices of the point, whilst bndry->bx and bndry->by are unit vectors out of the domain. The loop is over all the points from the domain outwards i.e. the point [bndry->x - bndry->bx][bndry->y - bndry->by] will always be defined.

Sometimes it’s useful to be able to loop over just one direction along the boundary. To do this, it is possible to use nextX() or nextY() rather than next(). It is also possible to loop over both dimensions using:

for(bndry->first(); !bndry->isDone(); bndry->nextX())
  for(; !bndry->isDone(); bndry->nextY()) {
    ...
  }

Boundary operations#

On each boundary, conditions must be specified for each variable. The different conditions are imposed by BoundaryOp objects. These set the values in the boundary region such that they obey e.g. Dirichlet or Neumann conditions. The BoundaryOp class is defined in boundary_op.hxx:

/// An operation on a boundary
class BoundaryOp {
 public:
  BoundaryOp() {bndry = NULL;}
  BoundaryOp(BoundaryRegion *region)

  // Note: All methods must implement clone, except for modifiers (see below)
  virtual BoundaryOp* clone(BoundaryRegion *region, const list<string> &args);

  /// Apply a boundary condition on field f
  virtual void apply(Field2D &f) = 0;
  virtual void apply(Field3D &f) = 0;

  virtual void apply(Vector2D &f);

  virtual void apply(Vector3D &f);

  /// Apply a boundary condition on ddt(f)
  virtual void apply_ddt(Field2D &f);
  virtual void apply_ddt(Field3D &f);
  virtual void apply_ddt(Vector2D &f);
  virtual void apply_ddt(Vector3D &f);

  BoundaryRegion *bndry;
};

(where the implementations have been removed for clarity). Which has a pointer to a BoundaryRegion object specifying which region this boundary is operating on.

Boundary conditions need to be imposed on the initial conditions (after PhysicsModel::init()), and on the time-derivatives (after PhysicsModel::rhs()). The apply() functions are therefore called during initialisation and given the evolving variables, whilst the apply_ddt functions are passed the time-derivatives.

To implement a boundary operation, as a minimum the apply(Field2D), apply(Field2D) and clone() need to be implemented: By default the apply(Vector) will call the apply(Field) functions on each component individually, and the apply_ddt() functions just call the apply() functions.

Example: Neumann boundary conditions are defined in boundary_standard.hxx:

/// Neumann (zero-gradient) boundary condition
class BoundaryNeumann : public BoundaryOp {
 public:
  BoundaryNeumann() {}
 BoundaryNeumann(BoundaryRegion *region):BoundaryOp(region) { }
  BoundaryOp* clone(BoundaryRegion *region, const list<string> &args);
  void apply(Field2D &f);
  void apply(Field3D &f);
};

and implemented in boundary_standard.cxx

void BoundaryNeumann::apply(Field2D &f) {
  // Loop over all elements and set equal to the next point in
  for(bndry->first(); !bndry->isDone(); bndry->next())
    f[bndry->x][bndry->y] = f[bndry->x - bndry->bx][bndry->y - bndry->by];
}

void BoundaryNeumann::apply(Field3D &f) {
  for(bndry->first(); !bndry->isDone(); bndry->next())
    for(int z=0;z<mesh->LocalNz;z++)
      f[bndry->x][bndry->y][z] = f[bndry->x - bndry->bx][bndry->y -
bndry->by][z];
}

This is all that’s needed in this case since there’s no difference between applying Neumann conditions to a variable and to its time-derivative, and Neumann conditions for vectors are just Neumann conditions on each vector component.

To create a boundary condition, we need to give it a boundary region to operate over:

BoundaryRegion *bndry = ...
BoundaryOp op = new BoundaryOp(bndry);

The clone function is used to create boundary operations given a single object as a template in BoundaryFactory. This can take additional arguments as a vector of strings - see explanation in Boundary factory.

Boundary modifiers#

To create more complicated boundary conditions from simple ones (such as Neumann conditions above), boundary operations can be modified by wrapping them up in a BoundaryModifier object, defined in boundary_op.hxx:

class BoundaryModifier : public BoundaryOp {
 public:
  virtual BoundaryOp* clone(BoundaryOp *op, const list<string> &args) = 0;
 protected:
  BoundaryOp *op;
};

Since BoundaryModifier inherits from BoundaryOp, modified boundary operations are just a different boundary operation and can be treated the same (Decorator pattern). Boundary modifiers could also be nested inside each other to create even more complicated boundary operations. Note that the clone function is different to the BoundaryOp one: instead of a BoundaryRegion to operate on, modifiers are passed a BoundaryOp to modify.

Currently the only modifier is BoundaryRelax, defined in boundary_standard.hxx:

/// Convert a boundary condition to a relaxing one
class BoundaryRelax : public BoundaryModifier {
 public:
  BoundaryRelax(BoutReal rate) {r = fabs(rate);}
  BoundaryOp* clone(BoundaryOp *op, const list<string> &args);

  void apply(Field2D &f);
  void apply(Field3D &f);

  void apply_ddt(Field2D &f);
  void apply_ddt(Field3D &f);
 private:
  BoundaryRelax() {} // Must be initialised with a rate
  BoutReal r;
};

Boundary factory#

The boundary factory creates new boundary operations from input strings, for example turning “relax(dirichlet)” into a relaxing Dirichlet boundary operation on a given region. It is defined in boundary_factory.hxx as a Singleton, so to get a pointer to the boundary factory use

BoundaryFactory *bfact = BoundaryFactory::getInstance();

and to delete this singleton, free memory and clean-up at the end use:

BoundaryFactory::cleanup();

Because users should be able to add new boundary conditions during PhysicsModel::init(), boundary conditions are not hard-wired into BoundaryFactory. Instead, boundary conditions must be registered with the factory, passing an instance which can later be cloned. This is done in bout++.cxx for the standard boundary conditions:

BoundaryFactory* bndry = BoundaryFactory::getInstance();
bndry->add(new BoundaryDirichlet(), "dirichlet");
...
bndry->addMod(new BoundaryRelax(10.), "relax");

where the add function adds BoundaryOp objects, whereas addMod adds BoundaryModifier objects. Note: The objects passed to BoundaryFactory will be deleted when cleanup() is called.

When a boundary operation is added, it is given a name such as “dirichlet”, and similarly for the modifiers (“relax” above). These labels and object pointers are stored internally in BoundaryFactory in maps defined in boundary_factory.hxx:

// Database of available boundary conditions and modifiers
map<string, BoundaryOp*> opmap;
map<string, BoundaryModifier*> modmap;

These are then used by BoundaryFactory::create():

/// Create a boundary operation object
BoundaryOp* create(const string &name, BoundaryRegion *region);
BoundaryOp* create(const char* name, BoundaryRegion *region);

to turn a string such as “relax(dirichlet)” and a BoundaryRegion pointer into a BoundaryOp object. These functions are implemented in boundary_factory.cxx, starting around line 42. The parsing is done recursively by matching the input string to one of:

modifier(<expression>, arg1, ...)
modifier(<expression>)
operation(arg1, ...)
operation

the <expression> variable is then resolved into a BoundaryOp object by calling create(<expression>, region).

When an operator or modifier is found, it is created from the pointer stored in the opmap or modmap maps using the clone method, passing a list<string> reference containing any arguments. It’s up to the operation implementation to ensure that the correct number of arguments are passed, and to parse them into floats or other types.

Example: The Dirichlet boundary condition can take an optional argument to change the value the boundary’s set to. In boundary_standard.cxx:

BoundaryOp* BoundaryDirichlet::clone(BoundaryRegion *region, const list<string>
&args) {
  if(!args.empty()) {
    // First argument should be a value
    stringstream ss;
    ss << args.front();

    BoutReal val;
    ss >> val;
    return new BoundaryDirichlet(region, val);
  }
  return new BoundaryDirichlet(region);
}

If no arguments are passed i.e. the string was “dirichlet” or “dirichlet()” then the args list is empty, and the default value (0.0) is used. If one or more arguments is used then the first argument is parsed into a BoutReal type and used to create a new BoundaryDirichlet object. If more arguments are passed then these are just ignored; probably a warning should be printed.

To set boundary conditions on a field, FieldData methods are defined in field_data.hxx:

// Boundary conditions
  void setBoundary(const string &name); ///< Set the boundary conditions
  void setBoundary(const string &region, BoundaryOp *op); ///< Manually set
  virtual void applyBoundary() {}
  virtual void applyTDerivBoundary() {};
 protected:
  vector<BoundaryOp*> bndry_op; // Boundary conditions

The FieldData::setBoundary() method is implemented in field_data.cxx. It first gets a vector of pointers to BoundaryRegions from the mesh, then loops over these calling BoundaryFactory::createFromOptions() for each one and adding the resulting boundary operator to the FieldData::bndry_op vector.

Boundary conditions

Contents

Boundary conditions#