Use gridvar: namespace for input variables read from grid

Author: ZedThree

manual/sphinx/user_docs/input_grids.rst

@@ -28,7 +28,7 @@ tensor is the identity matrix), but this can be changed by specifying
 the metric tensor components.
 
 Integer quantities such as ``nx`` can be numbers (like “260”), or
-expressions (like “256 + 2\*MXG”). 
+expressions (like “256 + 2\*MXG”).
 A common use is to make ``x`` and ``z`` dimensions have the same
 number of points, when ``x`` has ``mxg`` boundary cells on each
 boundary but ``z`` does not (since it is usually periodic):
@@ -37,14 +37,14 @@ boundary but ``z`` does not (since it is usually periodic):
 
     [mesh]
     nx = nz + 2*mxg  # X grid size
-    nz = 256         # Z grid size            
-    mxg = 2            
+    nz = 256         # Z grid size
+    mxg = 2
 
 
 Note that the order of the defintion within a section isn't important,
 variables can be used before they are defined. All variables are first
 read, and only processed if they are used.
-    
+
 Expressions are always calculated in floating point; When expressions
 are used to set integer quantities (such as the number of grid
 points), the expressions are calculated in floating point and then
@@ -60,7 +60,7 @@ the ``round`` function:
     nx = 256.4   # Error!
     nx = round(256.4) # ok
 
-    
+
 Real (floating-point) values can also be expressions, allowing quite
 complicated analytic inputs. For example in the example ``test-griddata``:
 
@@ -158,7 +158,9 @@ Cartesian.
 You can read additional quantities from the grid and make them available in
 expressions in the input file by listing them in the ``input:grid_variables``
 section, with the key being the name in the grid file (``mesh:file``) and the
-value being the type (one of ``field3d``, ``field2d``, ``boutreal``):
+value being the type (one of ``field3d``, ``field2d``, ``boutreal``). Variables
+read in this way are made available under a ``gridvar:`` namespace to avoid
+collisions with other input variables:
 
 .. code-block:: cfg
 
@@ -168,7 +170,7 @@ value being the type (one of ``field3d``, ``field2d``, ``boutreal``):
    scale = boutreal
 
    [mesh]
-   B = (scale / rho) * cos(theta)
+   B = (gridvar:scale / gridvar:rho) * cos(gridvar:theta)
 
 This section describes how to generate inputs for tokamak equilibria. If
 you’re not interested in tokamaks then you can skip to the next section.
@@ -219,7 +221,7 @@ From EFIT files
 A separate tool (in python) called `Hypnotoad <https://github.com/boutproject/hypnotoad>`_
 has been developed to create BOUT++ input files from R-Z equilibria. This can read EFIT ’g’
 (geqdsk) files, find flux surfaces, and calculate metric
-coefficients. 
+coefficients.
 
 From GRIDUE files
 --------------
@@ -281,7 +283,7 @@ As in the above code, creating an output file consists of the following steps:
 
 1. Define a magnetic field
 2. Define the grid points. This can be broken down into:
-   
+
    a) Define 2D "poloidal" grids
    b) Form a 3D grid by putting 2D grids together along the Y direction
 
@@ -307,9 +309,9 @@ In this case with 10 points in y (second argument to ``rectangular_grid(nx,ny,nz
 the y locations are :math:`\left(0.5, 1.5, 2.5, \ldots, 9.5\right)`.
 
 At each of these y locations ``rectangular_grid`` defines a rectangular 2D poloidal grid in
-the X-Z coordinates, by default with a length of 1 in each direction and centred on :math:`x=0,z=0`. 
+the X-Z coordinates, by default with a length of 1 in each direction and centred on :math:`x=0,z=0`.
 These 2D poloidal grids are then put together into a 3D ``Grid``. This process can be customised
-by separating step 2 (the ``rectangular_grid`` call) into stages 2a) and 2b). 
+by separating step 2 (the ``rectangular_grid`` call) into stages 2a) and 2b).
 For example, to create a periodic rectangular grid we could use the following::
 
    import numpy as np
@@ -333,7 +335,7 @@ input file (this is in ``examples/zoidberg/tokamak.py``)::
 
    import numpy as np
    import zoidberg
-   
+
    field = zoidberg.field.GEQDSK("g014220.00200") # Read magnetic field
 
    grid = zoidberg.grid.rectangular_grid(100, 10, 100,
@@ -345,13 +347,13 @@ input file (this is in ``examples/zoidberg/tokamak.py``)::
 
    # Create the forward and backward maps
    maps = zoidberg.make_maps(grid, field)
-   
+
    # Save to file
    zoidberg.write_maps(grid, field, maps, gridfile="grid.fci.nc")
 
    # Plot grid points and the points they map to in the forward direction
    zoidberg.plot.plot_forward_map(grid, maps)
-   
+
 In the last example only one poloidal grid was created (a ``RectangularPoloidalGrid``)
 and then re-used for each y slice. We can instead define a different grid for each y
 position. For example, to define a grid which expands along y (for some reason) we could do::
@@ -387,10 +389,10 @@ One way to create this grid is to define the grid points manually e.g.::
    r,theta = np.meshgrid(np.linspace(1,2,10),
                          np.linspace(0,2*np.pi, 10),
                          indexing="ij")
-   
+
    R = r * np.sin(theta)
    Z = r * np.cos(theta)
-  
+
    poloidal_grid = zoidberg.poloidal_grid.StructuredPoloidalGrid(R,Z)
 
 For more complicated shapes than circles, Zoidberg comes with an
@@ -402,21 +404,21 @@ outer boundaries::
    inner = zoidberg.rzline.shaped_line(R0=3.0, a=0.5,
                             elong=1.0, triang=0.0, indent=1.0,
                             n=50)
-   
+
    outer = zoidberg.rzline.shaped_line(R0=2.8, a=1.5,
                             elong=1.0, triang=0.0, indent=0.2,
                             n=50)
-   
+
    poloidal_grid = zoidberg.poloidal_grid.grid_elliptic(inner, outer,
                                                  100, 100, show=True)
 
 which should produce the figure below:
 
 .. figure:: ../figs/zoidberg/elliptic_grid.png
    :name: elliptic
-   :alt: 
+   :alt:
    :scale: 50
-   
+
    A grid produced by ``grid_elliptic`` from shaped inner and outer lines
 
 
@@ -440,14 +442,14 @@ flux surface.
 At the moment this will not work correctly for slab geometries, but expects
 closed flux surfaces such as in a stellarator or tokamak. A simple test case
 is a straight stellarator::
-   
+
    import zoidberg
    field = zoidberg.field.StraightStellarator(I_coil=0.4, yperiod=10)
 
 By default ``StraightStellarator`` calculates the magnetic field due to four coils which spiral around
 the axis at a distance :math:`r=0.8` in a classical stellarator configuration. The ``yperiod``
 argument is the period in y after which the coils return to their starting locations.
-   
+
 To visualise the Poincare plot for this stellarator field, pass the ``MagneticField`` object
 to ``zoidberg.plot.plot_poincare``, together with start location(s) and periodicity information::
 
@@ -459,14 +461,14 @@ which should produce the following figure:
    :name: poincare
    :alt: Points on four oval shaped flux surfaces in x-z at three locations along the y direction
    :scale: 50
-   
+
    Poincare map of straight stellarator showing a single flux
    surface. Each colour corresponds to a different x-z plane
-   in the y direction. 
-           
+   in the y direction.
+
 The inputs here are the starting location :math:`\left(x,z\right) = \left(0.4, 0.0\right)`,
 and the periodicity in the y direction (10.0). By default this will
-integrate from this given starting location 40 times (``revs`` option) around the y domain (0 to 10). 
+integrate from this given starting location 40 times (``revs`` option) around the y domain (0 to 10).
 
 To create an ``RZline`` from these Poincare plots we need a
 list of points in order around the line. Since the points
@@ -476,14 +478,14 @@ this is a `known hard problem <https://en.wikipedia.org/wiki/Travelling_salesman
 but fortunately in this case the nearest neighbour algorithm seems to be quite robust provided there are enough points.
 
 An example of calculating a Poincare plot on a single y slice (y=0) and producing an ``RZline`` is::
-   
+
    from zoidberg.fieldtracer import trace_poincare
    rzcoord, ycoords = trace_poincare(field, 0.4, 0.0, 10.0,
                                      y_slices=[0])
-   
+
    R = rzcoord[:,0,0]
    Z = rzcoord[:,0,1]
-          
+
    line = zoidberg.rzline.line_from_points(R, Z)
 
    line.plot()
@@ -495,7 +497,7 @@ approximation to the flux surface, increase the number of points by setting the
 (y revolutions) in the ``trace_poincare`` call.
 
 In general the points along this line are not evenly
-distributed, but tend to cluster together in some regions and have large gaps in others. 
+distributed, but tend to cluster together in some regions and have large gaps in others.
 The elliptic grid generator places grid points on the boundaries
 which are uniform in the index of the ``RZline`` it is given.
 Passing a very uneven set of points will therefore result in
@@ -509,7 +511,7 @@ to create a grid file for a straight stellarator.
 
 Sections below now describe each part of Zoidberg in more detail. Further documentation
 of the API can be found in the docstrings and unit tests.
-   
+
 Magnetic fields
 ~~~~~~~~~~~~~~~
 
@@ -518,12 +520,12 @@ Magnetic fields can be defined in either cylindrical or Cartesian coordinates:
 
 * In Cartesian coordinates all (x,y,z) directions have the same units of length
 * In cylindrical coordinates the y coordinate is assumed to be an angle, so that
-  the distance in y is given by :math:`ds = R dy` where :math:`R` is the major radius.  
+  the distance in y is given by :math:`ds = R dy` where :math:`R` is the major radius.
 
 Which coordinate is used is controlled by the ``Rfunc`` method, which should return the
 major radius if using a cylindrical coordinate system.
-Should return ``None`` for a Cartesian coordinate system (the default). 
-  
+Should return ``None`` for a Cartesian coordinate system (the default).
+
 Several implementations inherit from ``MagneticField``, and provide:
 ``Bxfunc``, ``Byfunc``, ``Bzfunc`` which give the components of the magnetic field in
 the x,y and z directions respectively. These should be in the same units (e.g. Tesla) for
@@ -550,14 +552,14 @@ and has a given major radius (default 1)::
    field = zoidberg.field.CurvedSlab()
 
 Note that this uses a large aspect-ratio approximation, so the major radius
-is constant across the domain (independent of x). 
-    
+is constant across the domain (independent of x).
+
 Straight stellarator
 ++++++++++++++++++++
 
 This is generated by four coils with alternating currents arranged
 on the edge of a circle, which spiral around the axis::
-   
+
    import zoidberg
    field = zoidberg.field.StraightStellarator()
 
@@ -588,7 +590,7 @@ Plotting the magnetic field
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Routines to plot the magnetic field are in ``zoidberg.plot``. They include Poincare plots
-and 3D field line plots. 
+and 3D field line plots.
 
 For example, to make a Poincare plot from a MAST equilibrium::
 
@@ -634,7 +636,7 @@ To create a rectangular grid, pass the number of points and lengths in the x and
 to ``RectangularPoloidalGrid``::
 
    import zoidberg
-   
+
    rect = zoidberg.poloidal_grid.RectangularPoloidalGrid( nx, nz, Lx, Lz )
 
 By default the middle of the rectangle is at :math:`\left(R,Z\right) = \left(0,0\right)`
@@ -649,7 +651,7 @@ To create the structured curvilinear grids inner and outer lines are needed
 with the following formula:
 
 .. math::
-   
+
    R = R_0 - b + \left(a + b \cos\left(\theta\right)\cos\left(\theta + \delta\sin\left(\theta\right)\right)\right)
 
    Z = \left(1 + \epsilon\right)a\sin\left(\theta\right)

src/field/field_factory.cxx

@@ -36,6 +36,7 @@
 #include <bout/utils.hxx>
 
 #include "fieldgenerators.hxx"
+#include "fmt/format.h"
 
 #include <cmath>
 #include <memory>
@@ -94,7 +95,8 @@ class FieldIndirect : public FieldGenerator {
 // Read variables from the grid file and make them available in expressions
 template <class T>
 auto add_grid_variable(FieldFactory& factory, Mesh& mesh, const std::string& name) {
-  factory.addGenerator(name, std::make_shared<GridVariable<T>>(&mesh, name));
+  factory.addGenerator(fmt::format("gridvar:{}", name),
+                       std::make_shared<GridVariable<T>>(&mesh, name));
 }
 
 auto read_grid_variables(FieldFactory& factory, Mesh& mesh, Options& options) {

tests/unit/field/test_field_factory.cxx

@@ -1066,7 +1066,7 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField3D) {
     dynamic_cast<FakeMesh*>(mesh)->setGridDataSource(new GridFile{filename});
     auto factory = FieldFactory{mesh, &options};
 
-    const auto output = factory.create3D("rho * cos(theta)");
+    const auto output = factory.create3D("gridvar:rho * cos(gridvar:theta)");
     const auto x = factory.create3D("x");
     EXPECT_TRUE(IsFieldEqual(output, x, "RGN_NOBNDRY", 1e-14));
   }
@@ -1094,7 +1094,7 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField2D) {
     dynamic_cast<FakeMesh*>(mesh)->setGridDataSource(new GridFile{filename});
     auto factory = FieldFactory{mesh, &options};
 
-    const auto output = factory.create2D("rho * cos(theta)");
+    const auto output = factory.create2D("gridvar:rho * cos(gridvar:theta)");
     const auto x = factory.create2D("x");
     EXPECT_TRUE(IsFieldEqual(output, x, "RGN_ALL", 1e-14));
   }