Variable initialisation¶
Variables in BOUT++ are not initialised automatically, but must be
explicitly given a value. For example the following code declares a
Field3D
variable then attempts to access a particular element:
Field3D f; // Declare a variable
f(0,0,0) = 1.0; // Error!
This results in an error because the data array to store values in f
has not been allocated. Allocating data can be done in several ways:
Initialise with a value:
Field3D f = 0.0; // Allocates memory, fills with zeros f(0,0,0) = 1.0; // ok
This cannot be done at a global scope, since it requires the mesh to already exist and have a defined size.
Set to a scalar value:
Field3D f; f = 0.0; // Allocates memory, fills with zeros f(0,0,0) = 1.0; // ok
Note that setting a field equal to another field has the effect of making both fields share the same underlying data. This behaviour is similar to how NumPy arrays behave in Python.
Field3D g = 0.0; // Allocates memory, fills with zeros Field3D f = g; // f now shares memory with g f(0,0,0) = 1.0; // g also modified
To ensure that a field has a unique underlying memory array call the
Field3D::allocate()
method before writing to individual indices.Use
Field3D::allocate()
to allocate memory:Field3D f; f.allocate(); // Allocates memory, values undefined f(0,0,0) = 1.0; // ok
In a BOUT++ simulation some variables are typically evolved in time. The initialisation of these variables is handled by the time integration solver.
Initialisation of time evolved variables¶
Each variable being evolved has its own section, with the same name as
the output data. For example, the high\(\beta\) model has
variables “P”, “jpar”, and “U”, and so has sections [P]
, [jpar]
,
[U]
(not case sensitive).
Expressions¶
The recommended way to initialise a variable is to use the function
option for each variable:
[p]
function = 1 + gauss(x0.5)*gauss(y)*sin(z)
This evaluates an analytic expression to initialise the \(P\)
variable. Expressions can include the usual operators
(+
,
,*
,/
), including ^
for exponents. The
following values are also already defined:
Name 
Description 

x 
\(x\) position between \(0\) and \(1\) 
y 
\(y\) position between \(0\) and \(2\pi\) (excluding the last point) 
z 
\(z\) position between \(0\) and \(2\pi\) (excluding the last point) 
pi π 
\(3.1415\ldots\) 
By default, \(x\) is defined as i / (nx  2*MXG)
, where MXG
is the width of the boundary region, by default 2. Hence \(x\)
actually goes from 0 on the leftmost point to (nx1)/(nx4)
on the
rightmost point. This is not a particularly good definition, but for
most cases its sufficient to create some initial profiles. For some
problems like island reconnection simulations, it’s useful to define
\(x\) in a particular way which is more symmetric than the default.
To do this, set in BOUT.inp
[mesh]
symmetricGlobalX = true
This will change the definition of \(x\) to i / (nx  1)
, so
\(x\) is then between \(0\) and \(1\) everywhere.
By default the expressions are evaluated in a fieldaligned coordinate system,
i.e. if you are using the [mesh]
option paralleltransform = shifted
,
the input f
will have f = fromFieldAligned(f)
applied before being
returned. To switch off this behaviour and evaluate the input expressions in
coordinates with orthogonal xz (i.e. toroidal \(\{\psi,\theta,\phi\}\)
coordinates when using paralleltransform = shifted
), set in BOUT.inp
[input]
transform_from_field_aligned = false
The functions in Table 2 are also available in expressions.
Name 
Description 


Absolute value \(x\) 

Inverse trigonometric functions 



Ballooning transform, using \(n\) terms (default 3) 

Cosine 

Hyperbolic cosine 

Exponential 

Hyperbolic tangent 

Gaussian \(\exp(x^2/2) / \sqrt{2\pi}\) 

Gaussian \(\exp[x^2/(2w^2)] / (w\sqrt{2\pi})\) 

Heaviside function: \(1\) if \(x > 0\) otherwise \(0\) 

Natural logarithm 

Maximum (variable arguments) 

Minimum (variable arguments) 

If value < low, return low; If value > high, return high; otherwise return value 

A mixture of Fourier modes 

seed determines random phase (default 0.5) 

Exponent \(x^y\) 

Sine 

Hyperbolic sine 

\(\sqrt{x}\) 

Tangent 

The error function 

The hat function \(\frac{1}{2}(\tanh[s (x[c\frac{w}{2}])]\) \( \tanh[s (x[c+\frac{w}{2}])] )\) 

The modulo operator, returns floating point remainder 
In addition there are some special functions which enable control flow
Name 
Description 


If the first 

Evaluate expression 

Define a new scope with variables whose value can be
accessed using braces 
For fieldaligned tokamak simulations, the Y direction is along the
field and in the core this will have a discontinuity at the twistshift
location where fieldlines are matched onto each other. To handle this,
the ballooning
function applies a truncated Ballooning
transformation to construct a smooth initial perturbation:
There is an example code testballooning
which compares methods of
setting initial conditions with the ballooning transform.
The mixmode(x)
function is a mixture of Fourier modes of the form:
where \(\phi\) is a random phase between \(\pi\) and \(+\pi\), which depends on the seed. The factor in front of each term is chosen so that the 4th harmonic (\(i=4\)) has the highest amplitude. This is useful mainly for initialising turbulence simulations, where a mixture of mode numbers is desired.
Context variables and scope¶
Expressions can use a form of local variables, by using []()
to
define new scopes:
var = [a = 2,
b = 3]( {a} + {b}^{a} )
Where here the braces {}
refer to context variables, to
distinguish them from variables in the options which have no
braces. One application of these is a (modest) performance
improvement: If {a}
is a large expression then in the above
example it would only be evaluated once, the value stored as {a}
and used twice in the expression.
Passing data into expressions¶
A second application of context variables is that they can be set by the calling C++ code, providing a way for data to be passed from BOUT++ into these expressions. The evaluation of expressions is currently not very efficient, but this provides a very flexible way for the input options to modify simulation behaviour.
This can be done by first parsing an expression and then passing values
to generate
in the Context
object.
Field3D shear = ...; // Value calculated in BOUT++
FieldFactory factory(mesh);
auto gen = factory>parse("model:viscosity");
Field3D viscosity;
viscosity.allocate();
BOUT_FOR(i, viscosity.region("RGN_ALL")) {
viscosity[i] = gen>generate(bout::generator::Context(i, CELL_CENTRE, mesh, 0.0)
.set("shear", shear[i]));
}
Note that the Context
constructor takes the index, the cell
location (e.g. staggered), a mesh, and then the time (set to 0.0
here). Additional variables can be set
, “shear” in this case. In
the input options file (or command line) the viscosity could now be a
function of {shear}
[model]
viscosity = 1 + {shear}
Defining functions in input options¶
Defining context variables in a new scope can be used to define and
call functions, as in the above example viscosity
is a function of
{shear}
. For example we could define a cosh function using
mycosh = 0.5 * (exp({arg}) + exp({arg}))
which uses {arg}
as the input value. We could then call this function:
result = [arg = x*2](mycosh)
Recursive functions¶
By default recursive expressions are not allowed in the input options,
and a ParseException
will be thrown if circular dependencies
occur. Recursive functions can however be enabled by setting
input:max_recursion_depth != 0
e.g.:
[input]
max_recursion_depth = 10 # 0 = none, 1 = unlimited
By putting a limit on the depth, expressions should (eventually)
terminate or fail with a BoutException
, rather than entering an
infinite loop. To remove this restriction max_recursion_depth
can
be set to 1 to allow arbitrary recursion (limited by stack, memory
sizes).
If recursion is allowed, then the where
special function and
Context
scopes can be (ab)used to define quite general
functions. For example the Fibonnacci sequence 1,1,2,3,5,8,...
can
be generated:
fib = where({n}  2.5,
[n={n}1](fib) + [n={n}2](fib),
1)
so if n
= 1 or 2 then fib
= 1, but if n = 3 or above then
recursion is used.
Note: Use of this facility in general is not encouraged, as it can easily lead to very inefficient and hard to understand code. It is here because occasionally it might be necessary, and because making the input language Turing complete was irresistible.
Initalising variables with the FieldFactory
class¶
This class provides a way to generate a field with a specified form. For
example to create a variable var
from options we could write
FieldFactory f(mesh);
Field2D var = f.create2D("var");
This will look for an option called “var”, and use that expression to
initialise the variable var
. This could then be set in the BOUT.inp
file or on the command line.
var = gauss(x0.5,0.2)*gauss(y)*sin(3*z)
To do this, FieldFactory
implements a recursive descent
parser to turn a string containing something like
"gauss(x0.5,0.2)*gauss(y)*sin(3*z)"
into values in a
Field3D
or Field2D
object. Examples are
given in the testfieldfactory
example:
FieldFactory f(mesh);
Field2D b = f.create2D("1  x");
Field3D d = f.create3D("gauss(x0.5,0.2)*gauss(y)*sin(z)");
This is done by creating a tree of FieldGenerator
objects
which then generate the field values:
class FieldGenerator {
public:
virtual ~FieldGenerator() { }
virtual FieldGeneratorPtr clone(const list<FieldGeneratorPtr> args) {return NULL;}
virtual BoutReal generate(const bout::generator::Context& ctx) = 0;
};
where FieldGeneratorPtr
is an alias for
std::shared_ptr
, a shared pointer to a
FieldGenerator
. The Context
input to generate
is an object
containing values which can be used in expressions, in particular x
,
y
, z
and t
coordinates. Additional values can be stored in the
Context
object, allowing data from BOUT++ to be used in expressions.
There are also ways to manipulate Context
objects for more complex
expressions and functions, see below for details.
All classes inheriting from FieldGenerator
must implement
a FieldGenerator::generate()
function, which returns the
value at the given (x,y,z,t)
position. Classes should also implement
a FieldGenerator::clone()
function, which takes a list of
arguments and creates a new instance of its class. This takes as input
a list of other FieldGenerator
objects, allowing a
variable number of arguments.
The simplest generator is a fixed numerical value, which is
represented by a FieldValue
object:
class FieldValue : public FieldGenerator {
public:
FieldValue(BoutReal val) : value(val) {}
BoutReal generate(const bout::generator::Context&) override { return value; }
private:
BoutReal value;
};
Adding a new function¶
To add a new function to the FieldFactory, a new
FieldGenerator
class must be defined. Here we will use
the example of the sinh
function, implemented using a class
FieldSinh
. This takes a single argument as input, but
FieldPI
takes no arguments, and
FieldGaussian
takes either one or two. Study these after
reading this to see how these are handled.
First, edit src/field/fieldgenerators.hxx
and add a class
definition:
class FieldSinh : public FieldGenerator {
public:
FieldSinh(FieldGeneratorPtr g) : gen(g) {}
FieldGeneratorPtr clone(const list<FieldGenerator*> args) override;
BoutReal generate(const bout::generator::Context& ctx) override;
private:
FieldGeneratorPtr gen;
};
The gen
member is used to store the input argument. The
constructor takes a single input, the FieldGenerator
argument to the
sinh
function, which is stored in the member gen
.
Next edit src/field/fieldgenerators.cxx
and add the implementation
of the clone
and generate
functions:
FieldGeneratorPtr FieldSinh::clone(const list<FieldGeneratorPtr> args) {
if (args.size() != 1) {
throw ParseException("Incorrect number of arguments to sinh function. Expecting 1, got %d", args.size());
}
return std::make_shared<FieldSinh>(args.front());
}
BoutReal FieldSinh::generate(const bout::generator::Context& ctx) {
return sinh(gen>generate(ctx));
}
The clone
function first checks the number of arguments using
args.size()
. This is used in FieldGaussian
to handle
different numbers of input, but in this case we throw a
ParseException
if the number of inputs isn’t
one. clone
then creates a new FieldSinh
object,
passing the first argument ( args.front()
) to the constructor
(which then gets stored in the gen
member variable).
Note that std::make_shared
is used to make a shared pointer.
The generate
function for sinh
just gets the value of the
input by calling gen>generate(ctx)
with the input Context
object ctx
, calculates sinh
of it and returns the result.
The clone
function means that the parsing code can make copies of
any FieldGenerator
class if it’s given a single instance to start
with. The final step is therefore to give the FieldFactory
class an
instance of this new generator. Edit the FieldFactory
constructor
FieldFactory::FieldFactory()
in src/field/field_factory.cxx
and
add the line:
addGenerator("sinh", std::make_shared<FieldSinh>(nullptr));
That’s it! This line associates the string "sinh"
with a
FieldGenerator
. Even though FieldFactory
doesn’t know what type of FieldGenerator
it is, it can
make more copies by calling the clone
member function. This is a
useful technique for polymorphic objects in C++ called the “Virtual
Constructor” idiom.
Parser internals¶
The basic expression parser is defined in
include/bout/sys/expressionparser.hxx
and the code in
src/sys/expressionparser.cxx
. The FieldFactory
adds the
function in table Table 2 on top of this basic
functionality, and also uses Options
to resolve unknown symbols to
Options
.
When a FieldGenerator
is added using the addGenerator
function, it is entered into a std::map
which maps strings to
FieldGenerator
objects (include/bout/sys/expressionparser.hxx
):
std::map<std::string, FieldGeneratorPtr> gen;
Parsing a string into a tree of FieldGenerator
objects is done by a
first splitting the string up into separate tokens like operators like
’*’, brackets ’(’, names like ’sinh’ and so on (Lexical analysis), then recognising
patterns in the stream of tokens (Parsing). Recognising tokens is done
in src/sys/expressionparser.cxx
:
char ExpressionParser::LexInfo::nextToken() {
...
This returns the next token, and setting the variable char curtok
to
the same value. This can be one of:
1 if the next token is a number. The variable
BoutReal curval
is set to the value of the token2 for a symbol (e.g. “sinh”, “x” or “pi”). This includes anything which starts with a letter, and contains only letters, numbers, and underscores. The string is stored in the variable
string curident
3 for a
Context
parameter which appeared surrounded by braces{}
.0 to mean end of input
The character if none of the above. Since letters and numbers are taken care of (see above), this includes brackets and operators like ’+’ and ’’.
The parsing stage turns these tokens into a tree of
FieldGenerator
objects, starting with the parse()
function:
FieldGenerator* FieldFactory::parse(const string &input) {
...
which puts the input string into a stream so that nextToken()
can
use it, then calls the parseExpression()
function to do the actual
parsing:
FieldGenerator* FieldFactory::parseExpression() {
...
This breaks down expressions in stages, starting with writing every expression as:
expression := primary [ op primary ]
i.e. a primary expression, and optionally an operator and another
primary expression. Primary expressions are handled by the
parsePrimary()
function, so first parsePrimary()
is called, and
then parseBinOpRHS
which checks if there is an operator, and if so
calls parsePrimary()
to parse it. This code also takes care of
operator precedence by keeping track of the precedence of the current
operator. Primary expressions are then further broken down and can
consist of either a number, a name (identifier), a minus sign and a
primary expression, or brackets around an expression:
primary := number
:= identifier
:= '' primary
:= '(' expression ')'
:= '[' expression ']'
The minus sign case is needed to handle the unary minus e.g. "x"
.
Identifiers are handled in parseIdentifierExpr()
which handles
either variable names, or functions
identifier := name
:= name '(' expression [ ',' expression [ ',' ... ] ] ')'
i.e. a name, optionally followed by brackets containing one or more
expressions separated by commas. names without brackets are treated the
same as those with empty brackets, so "x"
is the same as "x()"
.
A list of inputs (list<FieldGeneratorPtr> args;
) is created, the
gen
map is searched to find the FieldGenerator
object
corresponding to the name, and the list of inputs is passed to the
object’s clone
function.