Add ability to use variables from grid file in input file expressions

examples/elm-pb/elm_pb.cxx

@@ -1356,7 +1356,8 @@ class ELMpb : public PhysicsModel {
           // Only update if simulation time has advanced
           // Uses an exponential decay of the weighting of the value in the boundary
           // so that the solution is well behaved for arbitrary steps
-          BoutReal const weight = exp(-(t - phi_boundary_last_update) / phi_boundary_timescale);
+          const BoutReal weight =
+              exp(-(t - phi_boundary_last_update) / phi_boundary_timescale);
           phi_boundary_last_update = t;
 
           if (mesh->firstX()) {
@@ -1385,11 +1386,11 @@ class ELMpb : public PhysicsModel {
               }
 
               // Old value of phi at boundary. Note: this is constant in Z
-              BoutReal const oldvalue =
+              const BoutReal oldvalue =
                   0.5 * (phi(mesh->xstart - 1, j, 0) + phi(mesh->xstart, j, 0));
 
               // New value of phi at boundary, relaxing towards phivalue
-              BoutReal const newvalue = weight * oldvalue + (1. - weight) * phivalue;
+              const BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue;
 
               // Set phi at the boundary to this value
               for (int k = mesh->zstart; k <= mesh->zend; k++) {
@@ -1412,7 +1413,7 @@ class ELMpb : public PhysicsModel {
                   0.5 * (phi(mesh->xend + 1, j, 0) + phi(mesh->xend, j, 0));
 
               // New value of phi at boundary, relaxing towards phivalue
-              BoutReal const newvalue = weight * oldvalue + (1. - weight) * phivalue;
+              const BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue;
 
               // Set phi at the boundary to this value
               for (int k = mesh->zstart; k <= mesh->zend; k++) {
@@ -1625,7 +1626,7 @@ class ELMpb : public PhysicsModel {
         for (int jz = 0; jz < mesh->LocalNz; jz++) {
 
           // Zero-gradient potential
-          BoutReal const phisheath = phi_fa(r.ind, mesh->ystart, jz);
+          const BoutReal phisheath = phi_fa(r.ind, mesh->ystart, jz);
 
           BoutReal jsheath = -(sqrt(mi_me) / (2. * sqrt(PI))) * phisheath;
 
@@ -1646,7 +1647,7 @@ class ELMpb : public PhysicsModel {
         for (int jz = 0; jz < mesh->LocalNz; jz++) {
 
           // Zero-gradient potential
-          BoutReal const phisheath = phi_fa(r.ind, mesh->yend, jz);
+          const BoutReal phisheath = phi_fa(r.ind, mesh->yend, jz);
 
           BoutReal jsheath = (sqrt(mi_me) / (2. * sqrt(PI))) * phisheath;
 
@@ -2068,7 +2069,8 @@ class ELMpb : public PhysicsModel {
     ddt(P).applyBoundary("neumann");
 
     Field3D U1 = ddt(U);
-    U1 += (gamma * B0 * B0) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(ddt(P));
+    U1 +=
+        (gamma * B0 * B0) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(ddt(P));
 
     // Second matrix, solving Alfven wave dynamics
     static std::unique_ptr<InvertPar> invU{nullptr};

include/bout/field2d.hxx

@@ -4,9 +4,9 @@
  * \brief Definition of 2D scalar field class
  *
  **************************************************************************
- * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu
+ * Copyright 2010 - 2026 BOUT++ developers
  *
- * Contact: Ben Dudson, bd512@york.ac.uk
+ * Contact: Ben Dudson, dudson2@llnl.gov
  *
  * This file is part of BOUT++.
  *
@@ -33,13 +33,16 @@ class Field2D;
 #define BOUT_FIELD2D_H
 
 #include "bout/array.hxx"
-#include "bout/build_config.hxx"
+#include "bout/bout_types.hxx"
+#include "bout/build_defines.hxx"
 #include "bout/field.hxx"
 #include "bout/field_data.hxx"
 #include "bout/fieldperp.hxx"
 #include "bout/region.hxx"
 
 #include <cstddef>
+#include <iostream>
+#include <string>
 
 #if BOUT_HAS_RAJA
 #include "RAJA/RAJA.hpp" // using RAJA lib
@@ -206,7 +209,7 @@ public:
     }
 #endif
 
-    return data[jx * ny + jy];
+    return data[(jx * ny) + jy];
   }
   inline const BoutReal& operator()(int jx, int jy) const {
 #if CHECK > 2 && !BOUT_HAS_CUDA
@@ -220,7 +223,7 @@ public:
     }
 #endif
 
-    return data[jx * ny + jy];
+    return data[(jx * ny) + jy];
   }
 
   /*!

include/bout/field3d.hxx

@@ -1,7 +1,7 @@
 /**************************************************************************
- * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu
+ * Copyright 2010 - 2026 BOUT++ contributors
  *
- * Contact: Ben Dudson, bd512@york.ac.uk
+ * Contact: Ben Dudson, dudson2@llnl.gov
  *
  * This file is part of BOUT++.
  *
@@ -33,6 +33,7 @@ class Field3D;
 #include "bout/bout_types.hxx"
 #include "bout/field.hxx"
 #include "bout/field2d.hxx"
+#include "bout/field_data.hxx"
 #include "bout/fieldperp.hxx"
 #include "bout/region.hxx"
 #include "bout/traits.hxx"

include/bout/invertable_operator.hxx

@@ -135,9 +135,9 @@ public:
       : operatorFunction(func), preconditionerFunction(func),
         opt(optIn == nullptr ? Options::getRoot()->getSection("invertableOperator")
                              : optIn),
-        localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn), lib(opt){
+        localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn), lib(opt) {
 
-                                                                               };
+        };
 
   /// Destructor just has to cleanup the PETSc owned objects.
   ~InvertableOperator() {

include/bout/mesh.hxx

@@ -597,7 +597,7 @@ public:
   virtual int getLocalZIndexNoBoundaries(int zglobal) const = 0;
 
   /// Size of the mesh on this processor including guard/boundary cells
-  int LocalNx, LocalNy, LocalNz;
+  int LocalNx{0}, LocalNy{0}, LocalNz{0};
 
   /// Local ranges of data (inclusive), excluding guard cells
   int xstart, xend, ystart, yend, zstart, zend;

include/bout/options.hxx

@@ -74,27 +74,27 @@ class Options;
  * which can be used as a map.
  *
  *     Options options;
- *     
+ *
  *     // Set values
  *     options["key"] = 1.0;
  *
  *     // Get values. Throws BoutException if not found
- *     int val = options["key"]; // Sets val to 1 
+ *     int val = options["key"]; // Sets val to 1
  *
  *     // Return as specified type. Throws BoutException if not found
  *     BoutReal var = options["key"].as<BoutReal>();
  *
  *     // A default value can be used if key is not found
  *     BoutReal value = options["pi"].withDefault(3.14);
- *    
+ *
  *     // Assign value with source label. Throws if already has a value from same source
  *     options["newkey"].assign(1.0, "some source");
  *
  *     // Force assign a new value
  *     options["newkey"].force(2.0, "some source");
  *
  * A legacy interface is also supported:
- * 
+ *
  *     options.set("key", 1.0, "code"); // Sets a key from source "code"
  *
  *     int val;
@@ -119,9 +119,9 @@ class Options;
  *
  * Each Options object can also contain any number of sections, which are
  * themselves Options objects.
- * 
+ *
  *     Options &section = options["section"];
- * 
+ *
  * which can be nested:
  *
  *     options["section"]["subsection"]["value"] = 3;
@@ -134,13 +134,13 @@ class Options;
  *
  * e.g.
  *     options->getSection("section")->getSection("subsection")->set("value", 3);
- * 
+ *
  * Options also know about their parents:
  *
  *     Options &parent = section.parent();
- *     
+ *
  * or
- * 
+ *
  *     Options *parent = section->getParent();
  *
  * Root options object
@@ -150,8 +150,8 @@ class Options;
  * there is a global singleton Options object which can be accessed with a static function
  *
  *    Options &root = Options::root();
- * 
- * or 
+ *
+ * or
  *
  *    Options *root = Options::getRoot();
  *
@@ -193,7 +193,7 @@ public:
   /// @param[in] parent        Parent object
   /// @param[in] sectionName   Name of the section, including path from the root
   Options(Options* parent_instance, std::string full_name)
-      : parent_instance(parent_instance), full_name(std::move(full_name)){};
+      : parent_instance(parent_instance), full_name(std::move(full_name)) {};
 
   /// Initialise with a value
   /// These enable Options to be constructed using initializer lists

include/bout/output.hxx

@@ -166,13 +166,13 @@ class ConditionalOutput : public Output {
 public:
   /// @param[in] base    The Output object which will be written to if enabled
   /// @param[in] enabled Should this be enabled by default?
-  ConditionalOutput(Output* base, bool enabled = true) : base(base), enabled(enabled){};
+  ConditionalOutput(Output* base, bool enabled = true) : base(base), enabled(enabled) {};
 
   /// Constuctor taking ConditionalOutput. This allows several layers of conditions
   ///
   /// @param[in] base    A ConditionalOutput which will be written to if enabled
   ///
-  ConditionalOutput(ConditionalOutput* base) : base(base), enabled(base->enabled){};
+  ConditionalOutput(ConditionalOutput* base) : base(base), enabled(base->enabled) {};
 
   /// If enabled, writes a string using fmt formatting
   /// by calling base->write
@@ -237,7 +237,7 @@ private:
 ///    output_debug << "debug message";
 /// compile but have no effect if BOUT_USE_OUTPUT_DEBUG is false
 template <typename T>
-DummyOutput& operator<<(DummyOutput& out, T const& UNUSED(t)) {
+DummyOutput& operator<<(DummyOutput& out, const T& UNUSED(t)) {
   return out;
 }
 
@@ -261,7 +261,7 @@ inline ConditionalOutput& operator<<(ConditionalOutput& out, stream_manipulator
 }
 
 template <typename T>
-ConditionalOutput& operator<<(ConditionalOutput& out, T const& t) {
+ConditionalOutput& operator<<(ConditionalOutput& out, const T& t) {
   if (out.isEnabled()) {
     *out.getBase() << t;
   }