From 58397e0a2af834afb2fbb200a23abb97ac40cd41 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 1 Dec 2025 15:00:20 +0000 Subject: [PATCH 01/19] Remove unused `Context` constructor --- include/bout/sys/generator_context.hxx | 3 --- src/sys/generator_context.cxx | 3 --- 2 files changed, 6 deletions(-) diff --git a/include/bout/sys/generator_context.hxx b/include/bout/sys/generator_context.hxx index bbe60fab65..528b96113d 100644 --- a/include/bout/sys/generator_context.hxx +++ b/include/bout/sys/generator_context.hxx @@ -27,9 +27,6 @@ public: /// Context(int ix, int iy, int iz, CELL_LOC loc, Mesh* msh, BoutReal t); - /// Specify the values directly - Context(BoutReal x, BoutReal y, BoutReal z, Mesh* msh, BoutReal t); - /// If constructed without parameters, contains no values (null). /// Requesting x,y,z or t throws an exception Context() = default; diff --git a/src/sys/generator_context.cxx b/src/sys/generator_context.cxx index 31a5662378..6088f99fbd 100644 --- a/src/sys/generator_context.cxx +++ b/src/sys/generator_context.cxx @@ -45,8 +45,5 @@ Context::Context(const BoundaryRegion* bndry, int iz, CELL_LOC loc, BoutReal t, parameters["t"] = t; } -Context::Context(BoutReal x, BoutReal y, BoutReal z, Mesh* msh, BoutReal t) - : localmesh(msh), parameters{{"x", x}, {"y", y}, {"z", z}, {"t", t}} {} - } // namespace generator } // namespace bout From ab96f260fde8fe6b958e8370349ad63687665037 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 2 Dec 2025 12:59:00 +0000 Subject: [PATCH 02/19] Add cell indices to `generator::Context` --- include/bout/sys/generator_context.hxx | 10 +++++++++- src/sys/generator_context.cxx | 19 ++++++++----------- 2 files changed, 17 insertions(+), 12 deletions(-) diff --git a/include/bout/sys/generator_context.hxx b/include/bout/sys/generator_context.hxx index 528b96113d..a5dd9d236a 100644 --- a/include/bout/sys/generator_context.hxx +++ b/include/bout/sys/generator_context.hxx @@ -24,7 +24,6 @@ public: : Context(i.x(), i.y(), 0, (loc == CELL_ZLOW) ? CELL_CENTRE : loc, msh, t) {} /// Specify a cell index, together with the cell location, mesh and time - /// Context(int ix, int iy, int iz, CELL_LOC loc, Mesh* msh, BoutReal t); /// If constructed without parameters, contains no values (null). @@ -41,6 +40,11 @@ public: BoutReal z() const { return get("z"); } BoutReal t() const { return get("t"); } + /// Cell indices + int ix() const { return ix_; } + int jy() const { return jy_; } + int kz() const { return kz_;} + /// Set the value of a parameter with given name Context& set(const std::string& name, BoutReal value) { parameters[name] = value; @@ -73,6 +77,10 @@ public: } private: + int ix_{0}; + int jy_{0}; + int kz_{0}; + Mesh* localmesh{nullptr}; ///< The mesh on which the position is defined /// Contains user-set values which can be set and retrieved diff --git a/src/sys/generator_context.cxx b/src/sys/generator_context.cxx index 6088f99fbd..01090daa51 100644 --- a/src/sys/generator_context.cxx +++ b/src/sys/generator_context.cxx @@ -9,7 +9,7 @@ namespace bout { namespace generator { Context::Context(int ix, int iy, int iz, CELL_LOC loc, Mesh* msh, BoutReal t) - : localmesh(msh) { + : ix_(ix), jy_(iy), kz_(iz), localmesh(msh) { parameters["x"] = (loc == CELL_XLOW) ? 0.5 * (msh->GlobalX(ix) + msh->GlobalX(ix - 1)) : msh->GlobalX(ix); @@ -24,20 +24,17 @@ Context::Context(int ix, int iy, int iz, CELL_LOC loc, Mesh* msh, BoutReal t) } Context::Context(const BoundaryRegion* bndry, int iz, CELL_LOC loc, BoutReal t, Mesh* msh) - : localmesh(msh) { - - // Add one to X index if boundary is in -x direction, so that XLOW is on the boundary - const int ix = (bndry->bx < 0) ? bndry->x + 1 : bndry->x; + : // Add one to X index if boundary is in -x direction, so that XLOW is on the boundary + ix_((bndry->bx < 0) ? bndry->x + 1 : bndry->x), + jy_((bndry->by < 0) ? bndry->y + 1 : bndry->y), kz_(iz), localmesh(msh) { parameters["x"] = ((loc == CELL_XLOW) || (bndry->bx != 0)) - ? 0.5 * (msh->GlobalX(ix) + msh->GlobalX(ix - 1)) - : msh->GlobalX(ix); - - const int iy = (bndry->by < 0) ? bndry->y + 1 : bndry->y; + ? 0.5 * (msh->GlobalX(ix_) + msh->GlobalX(ix_ - 1)) + : msh->GlobalX(ix_); parameters["y"] = ((loc == CELL_YLOW) || (bndry->by != 0)) - ? PI * (msh->GlobalY(iy) + msh->GlobalY(iy - 1)) - : TWOPI * msh->GlobalY(iy); + ? PI * (msh->GlobalY(jy_) + msh->GlobalY(jy_ - 1)) + : TWOPI * msh->GlobalY(jy_); parameters["z"] = (loc == CELL_ZLOW) ? PI * (msh->GlobalZ(iz) + msh->GlobalZ(iz - 1)) : TWOPI * msh->GlobalZ(iz); From c5ef3d1d40621a2278e294fba75b4594bd569d7a Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 2 Dec 2025 13:28:48 +0000 Subject: [PATCH 03/19] Add ability to use variables from grid files in input expressions --- src/field/field_factory.cxx | 48 ++++++++++++++++++++++--- src/field/fieldgenerators.hxx | 28 +++++++++++++++ tests/unit/field/test_field_factory.cxx | 37 +++++++++++++++++++ 3 files changed, 108 insertions(+), 5 deletions(-) diff --git a/src/field/field_factory.cxx b/src/field/field_factory.cxx index 188a64cf0c..3e671b5c93 100644 --- a/src/field/field_factory.cxx +++ b/src/field/field_factory.cxx @@ -19,20 +19,24 @@ * along with BOUT++. If not, see . * **************************************************************************/ -#include #include -#include - +#include #include +#include #include #include -#include "bout/constants.hxx" - #include "fieldgenerators.hxx" +#include +#include +#include +#include +#include +#include + using bout::generator::Context; /// Helper function to create a FieldValue generator from a BoutReal @@ -80,6 +84,16 @@ class FieldIndirect : public FieldGenerator { FieldGeneratorPtr target; }; + +// Split a string on commas and trim whitespace from the results +auto trimsplit(const std::string& str) -> std::vector { + auto split = strsplit(str, ','); + std::vector result{}; + result.reserve(split.size()); + std::transform(split.begin(), split.end(), std::back_inserter(result), + [](const std::string& element) { return trim(element); }); + return result; +} } // namespace ////////////////////////////////////////////////////////// @@ -122,6 +136,30 @@ FieldFactory::FieldFactory(Mesh* localmesh, Options* opt) nonconst_options["input"]["max_recursion_depth"].as()); } + const auto field_variables = + nonconst_options["input"]["field_variables"] + .doc("Variables to read from the grid file and make available in expressions. " + "Comma-separated list of variable names") + .withDefault(""); + + if (not field_variables.empty()) { + if (not localmesh->isDataSourceGridFile()) { + throw BoutException("A grid file ('mesh:file') is required for `field_variables`"); + } + + for (const auto& name : trimsplit(field_variables)) { + if (not localmesh->sourceHasVar(name)) { + const auto filename = Options::root()["mesh"]["file"].as(); + throw BoutException( + "Grid file '{}' missing `{}` specified in `input:field_variables`", filename, + name); + } + Field3D var; + localmesh->get(var, name); + addGenerator(name, std::make_shared(var, name)); + } + } + // Useful values addGenerator("pi", std::make_shared(PI)); addGenerator("π", std::make_shared(PI)); diff --git a/src/field/fieldgenerators.hxx b/src/field/fieldgenerators.hxx index 050e335448..9b46aea6d3 100644 --- a/src/field/fieldgenerators.hxx +++ b/src/field/fieldgenerators.hxx @@ -8,10 +8,14 @@ #define BOUT_FIELDGENERATORS_H #include +#include #include +#include #include #include +#include +#include ////////////////////////////////////////////////////////// // Generators from values @@ -377,4 +381,28 @@ public: } }; +/// A `Field3D` that can be used in expressions +class Field3DVariable : public FieldGenerator { +public: + Field3DVariable(Field3D var, std::string name) + : variable(std::move(var)), name(std::move(name)) {} + + double generate(const bout::generator::Context& ctx) override { + return variable(ctx.ix(), ctx.jy(), ctx.kz()); + } + + FieldGeneratorPtr clone(const std::list args) override { + if (args.size() != 0) { + throw ParseException("Variable '{}' takes no arguments but got {:d}", args.size()); + } + return std::make_shared(variable, name); + } + + std::string str() const override { return name; } + +private: + Field3D variable; + std::string name; +}; + #endif // BOUT_FIELDGENERATORS_H diff --git a/tests/unit/field/test_field_factory.cxx b/tests/unit/field/test_field_factory.cxx index b45206b979..8a02033f01 100644 --- a/tests/unit/field/test_field_factory.cxx +++ b/tests/unit/field/test_field_factory.cxx @@ -1,17 +1,21 @@ #include "gtest/gtest.h" +#include "fake_mesh.hxx" #include "test_extras.hxx" #include "bout/boutexception.hxx" #include "bout/constants.hxx" #include "bout/field2d.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include "bout/globals.hxx" #include "bout/mesh.hxx" +#include "bout/options_io.hxx" #include "bout/output.hxx" #include "bout/paralleltransform.hxx" #include "bout/traits.hxx" #include "fake_mesh_fixture.hxx" +#include "test_tmpfiles.hxx" // The unit tests use the global mesh using namespace bout::globals; @@ -1027,3 +1031,36 @@ TYPED_TEST(FieldFactoryCreationTest, CreatePeriodicYacrossSeparatrix) { EXPECT_TRUE(IsFieldEqual(output, expected)); } + +struct FieldFactoryFieldVariableTest : public FakeMeshFixture { + WithQuietOutput quiet{output_info}; +}; + +TEST_F(FieldFactoryFieldVariableTest, CreateField3D) { + bout::testing::TempFile filename; + + { + // Write some fields to a grid file + FieldFactory factory{mesh}; + const auto rho = factory.create3D("sqrt(x^2 + y^2)"); + const auto theta = factory.create3D("atan(y, x)"); + Options grid{{"rho", rho}, + {"theta", theta}, + {"nx", mesh->LocalNx}, + {"ny", mesh->LocalNy}, + {"nz", mesh->LocalNz}}; + bout::OptionsIO::create(filename)->write(grid); + } + + { + Options options{{"mesh", {{"file", filename.string()}}}, + {"input", {{"field_variables", "rho, theta"}}}}; + + dynamic_cast(mesh)->setGridDataSource(new GridFile{filename}); + auto factory = FieldFactory{mesh, &options}; + + const auto output = factory.create3D("rho * cos(theta)"); + const auto x = factory.create3D("x"); + EXPECT_TRUE(IsFieldEqual(output, x)); + } +} From a843556c656e0a4f2ebe982ab995ecdfc8e4425c Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 2 Dec 2025 13:29:18 +0000 Subject: [PATCH 04/19] tests: Only remove temp files if test passes --- tests/unit/test_tmpfiles.hxx | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/tests/unit/test_tmpfiles.hxx b/tests/unit/test_tmpfiles.hxx index 6c37e7c792..db7387bb45 100644 --- a/tests/unit/test_tmpfiles.hxx +++ b/tests/unit/test_tmpfiles.hxx @@ -1,5 +1,6 @@ #pragma once +#include #include #include #include @@ -30,7 +31,11 @@ public: TempFile& operator=(const TempFile&) = delete; TempFile& operator=(TempFile&&) = delete; - ~TempFile() { std::filesystem::remove_all(filename); } + ~TempFile() { + if (std::uncaught_exceptions() <= 0) { + std::filesystem::remove_all(filename); + } + } // Enable conversions to std::string / const char* operator std::string() const { return filename.string(); } From 7c44c66752b32d4b8f9c9f87b434566a3b4c0284 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 2 Dec 2025 16:06:24 +0000 Subject: [PATCH 05/19] Enable reading `Field2D`s as well as `Field3D`s for expressions --- manual/sphinx/user_docs/input_grids.rst | 13 ++---- src/field/field_factory.cxx | 61 +++++++++++++++---------- src/field/fieldgenerators.hxx | 15 +++--- tests/unit/field/test_field_factory.cxx | 60 +++++++++++++++++++++++- 4 files changed, 109 insertions(+), 40 deletions(-) diff --git a/manual/sphinx/user_docs/input_grids.rst b/manual/sphinx/user_docs/input_grids.rst index 402fb30d06..2c17411141 100644 --- a/manual/sphinx/user_docs/input_grids.rst +++ b/manual/sphinx/user_docs/input_grids.rst @@ -156,19 +156,16 @@ these are the only quantities specified, then the coordinates revert to 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``): +expressions in the input file using ``input:read_Field3Ds`` and +``input:read_Field2Ds``: .. code-block:: cfg - [input:grid_variables] - rho = field2d - theta = field2d - scale = boutreal + [input] + read_Field2Ds = rho, theta [mesh] - B = (scale / rho) * cos(theta) + B = (1 / rho) * cos(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. diff --git a/src/field/field_factory.cxx b/src/field/field_factory.cxx index 3e671b5c93..e9dbc0c804 100644 --- a/src/field/field_factory.cxx +++ b/src/field/field_factory.cxx @@ -24,8 +24,12 @@ #include #include +#include +#include +#include #include #include +#include #include #include "fieldgenerators.hxx" @@ -94,6 +98,35 @@ auto trimsplit(const std::string& str) -> std::vector { [](const std::string& element) { return trim(element); }); return result; } + +// Read variables from the grid file and make them available in expressions +template +auto read_grid_variables(FieldFactory& factory, Mesh& mesh, Options& options, + const std::string& option_name) { + const auto field_variables = + options["input"][option_name] + .doc("Variables to read from the grid file and make available in expressions. " + "Comma-separated list of variable names") + .withDefault(""); + + if (not field_variables.empty()) { + if (not mesh.isDataSourceGridFile()) { + throw BoutException("A grid file ('mesh:file') is required for `input:{}`", + option_name); + } + + for (const auto& name : trimsplit(field_variables)) { + if (not mesh.sourceHasVar(name)) { + const auto filename = Options::root()["mesh"]["file"].as(); + throw BoutException("Grid file '{}' missing `{}` specified in `input:{}`", + filename, name, option_name); + } + T var; + mesh.get(var, name); + factory.addGenerator(name, std::make_shared>(var, name)); + } + } +} } // namespace ////////////////////////////////////////////////////////// @@ -136,30 +169,6 @@ FieldFactory::FieldFactory(Mesh* localmesh, Options* opt) nonconst_options["input"]["max_recursion_depth"].as()); } - const auto field_variables = - nonconst_options["input"]["field_variables"] - .doc("Variables to read from the grid file and make available in expressions. " - "Comma-separated list of variable names") - .withDefault(""); - - if (not field_variables.empty()) { - if (not localmesh->isDataSourceGridFile()) { - throw BoutException("A grid file ('mesh:file') is required for `field_variables`"); - } - - for (const auto& name : trimsplit(field_variables)) { - if (not localmesh->sourceHasVar(name)) { - const auto filename = Options::root()["mesh"]["file"].as(); - throw BoutException( - "Grid file '{}' missing `{}` specified in `input:field_variables`", filename, - name); - } - Field3D var; - localmesh->get(var, name); - addGenerator(name, std::make_shared(var, name)); - } - } - // Useful values addGenerator("pi", std::make_shared(PI)); addGenerator("π", std::make_shared(PI)); @@ -217,6 +226,10 @@ FieldFactory::FieldFactory(Mesh* localmesh, Options* opt) // Periodic in the Y direction? addGenerator("is_periodic_y", std::make_shared()); + + // Variables from the grid file + read_grid_variables(*this, *localmesh, nonconst_options, "read_Field3Ds"); + read_grid_variables(*this, *localmesh, nonconst_options, "read_Field2Ds"); } Field2D FieldFactory::create2D(const std::string& value, const Options* opt, diff --git a/src/field/fieldgenerators.hxx b/src/field/fieldgenerators.hxx index 9b46aea6d3..49b2cb9690 100644 --- a/src/field/fieldgenerators.hxx +++ b/src/field/fieldgenerators.hxx @@ -8,9 +8,9 @@ #define BOUT_FIELDGENERATORS_H #include -#include #include #include +#include #include #include @@ -311,7 +311,7 @@ public: // Constructor FieldTanhHat(FieldGeneratorPtr xin, FieldGeneratorPtr widthin, FieldGeneratorPtr centerin, FieldGeneratorPtr steepnessin) - : X(xin), width(widthin), center(centerin), steepness(steepnessin){}; + : X(xin), width(widthin), center(centerin), steepness(steepnessin) {}; // Clone containing the list of arguments FieldGeneratorPtr clone(const std::list args) override; BoutReal generate(const bout::generator::Context& pos) override; @@ -326,7 +326,7 @@ private: class FieldWhere : public FieldGenerator { public: FieldWhere(FieldGeneratorPtr test, FieldGeneratorPtr gt0, FieldGeneratorPtr lt0) - : test(test), gt0(gt0), lt0(lt0){}; + : test(test), gt0(gt0), lt0(lt0) {}; FieldGeneratorPtr clone(const std::list args) override { if (args.size() != 3) { @@ -382,9 +382,10 @@ public: }; /// A `Field3D` that can be used in expressions -class Field3DVariable : public FieldGenerator { +template > +class GridVariable : public FieldGenerator { public: - Field3DVariable(Field3D var, std::string name) + GridVariable(T var, std::string name) : variable(std::move(var)), name(std::move(name)) {} double generate(const bout::generator::Context& ctx) override { @@ -395,13 +396,13 @@ public: if (args.size() != 0) { throw ParseException("Variable '{}' takes no arguments but got {:d}", args.size()); } - return std::make_shared(variable, name); + return std::make_shared>(variable, name); } std::string str() const override { return name; } private: - Field3D variable; + T variable; std::string name; }; diff --git a/tests/unit/field/test_field_factory.cxx b/tests/unit/field/test_field_factory.cxx index 8a02033f01..1f88d9ce7e 100644 --- a/tests/unit/field/test_field_factory.cxx +++ b/tests/unit/field/test_field_factory.cxx @@ -1054,7 +1054,7 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField3D) { { Options options{{"mesh", {{"file", filename.string()}}}, - {"input", {{"field_variables", "rho, theta"}}}}; + {"input", {{"read_Field3Ds", "rho, theta"}}}}; dynamic_cast(mesh)->setGridDataSource(new GridFile{filename}); auto factory = FieldFactory{mesh, &options}; @@ -1064,3 +1064,61 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField3D) { EXPECT_TRUE(IsFieldEqual(output, x)); } } + +TEST_F(FieldFactoryFieldVariableTest, CreateField2D) { + bout::testing::TempFile filename; + + { + // Write some fields to a grid file + FieldFactory factory{mesh}; + const auto rho = factory.create2D("sqrt(x^2 + y^2)"); + const auto theta = factory.create2D("atan(y, x)"); + Options grid{{"rho", rho}, + {"theta", theta}, + {"nx", mesh->LocalNx}, + {"ny", mesh->LocalNy}, + {"nz", mesh->LocalNz}}; + bout::OptionsIO::create(filename)->write(grid); + } + + { + Options options{{"mesh", {{"file", filename.string()}}}, + {"input", {{"read_Field2Ds", "rho, theta"}}}}; + + dynamic_cast(mesh)->setGridDataSource(new GridFile{filename}); + auto factory = FieldFactory{mesh, &options}; + + const auto output = factory.create2D("rho * cos(theta)"); + const auto x = factory.create2D("x"); + EXPECT_TRUE(IsFieldEqual(output, x)); + } +} + +TEST_F(FieldFactoryFieldVariableTest, NoMeshFile) { + Options options{{"input", {{"read_Field2Ds", "rho, theta"}}}}; + + EXPECT_THROW((FieldFactory(mesh, &options)), BoutException); +} + +TEST_F(FieldFactoryFieldVariableTest, MissingVariable) { + bout::testing::TempFile filename; + + { + // Write some fields to a grid file + FieldFactory factory{mesh}; + const auto rho = factory.create3D("sqrt(x^2 + y^2)"); + Options grid{{"rho", rho}, + {"nx", mesh->LocalNx}, + {"ny", mesh->LocalNy}, + {"nz", mesh->LocalNz}}; + bout::OptionsIO::create(filename)->write(grid); + } + + { + Options options{{"mesh", {{"file", filename.string()}}}, + {"input", {{"read_Field3Ds", "rho, theta"}}}}; + + dynamic_cast(mesh)->setGridDataSource(new GridFile{filename}); + EXPECT_THROW((FieldFactory{mesh, &options}), BoutException); + } +} From d3b1346c2d3fd573dbb0c172e9a3adf540fda5a8 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 2 Dec 2025 16:49:00 +0000 Subject: [PATCH 06/19] Refactor reading grid variables for input expressions --- manual/sphinx/user_docs/input_grids.rst | 13 +++-- src/field/field_factory.cxx | 72 +++++++++++++------------ tests/unit/field/test_field_factory.cxx | 42 ++++++++++++--- 3 files changed, 81 insertions(+), 46 deletions(-) diff --git a/manual/sphinx/user_docs/input_grids.rst b/manual/sphinx/user_docs/input_grids.rst index 2c17411141..402fb30d06 100644 --- a/manual/sphinx/user_docs/input_grids.rst +++ b/manual/sphinx/user_docs/input_grids.rst @@ -156,16 +156,19 @@ these are the only quantities specified, then the coordinates revert to Cartesian. You can read additional quantities from the grid and make them available in -expressions in the input file using ``input:read_Field3Ds`` and -``input:read_Field2Ds``: +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``): .. code-block:: cfg - [input] - read_Field2Ds = rho, theta + [input:grid_variables] + rho = field2d + theta = field2d + scale = boutreal [mesh] - B = (1 / rho) * cos(theta) + B = (scale / rho) * cos(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. diff --git a/src/field/field_factory.cxx b/src/field/field_factory.cxx index e9dbc0c804..19e6a2fd14 100644 --- a/src/field/field_factory.cxx +++ b/src/field/field_factory.cxx @@ -31,15 +31,15 @@ #include #include #include +#include +#include +#include #include "fieldgenerators.hxx" -#include #include -#include #include #include -#include using bout::generator::Context; @@ -53,6 +53,8 @@ FieldGeneratorPtr generator(BoutReal* ptr) { return std::make_shared(ptr); } +BOUT_ENUM_CLASS(GridVariableFunction, field3d, field2d, boutreal); + namespace { /// Provides a placeholder whose target can be changed after creation. /// This enables recursive FieldGenerator expressions to be generated @@ -89,41 +91,44 @@ class FieldIndirect : public FieldGenerator { FieldGeneratorPtr target; }; -// Split a string on commas and trim whitespace from the results -auto trimsplit(const std::string& str) -> std::vector { - auto split = strsplit(str, ','); - std::vector result{}; - result.reserve(split.size()); - std::transform(split.begin(), split.end(), std::back_inserter(result), - [](const std::string& element) { return trim(element); }); - return result; -} - // Read variables from the grid file and make them available in expressions template -auto read_grid_variables(FieldFactory& factory, Mesh& mesh, Options& options, - const std::string& option_name) { - const auto field_variables = - options["input"][option_name] - .doc("Variables to read from the grid file and make available in expressions. " - "Comma-separated list of variable names") - .withDefault(""); - - if (not field_variables.empty()) { +auto add_grid_variable(FieldFactory& factory, Mesh& mesh, const std::string& name) { + T var; + mesh.get(var, name); + factory.addGenerator(name, std::make_shared>(var, name)); +} + +auto read_grid_variables(FieldFactory& factory, Mesh& mesh, Options& options) { + auto& field_variables = options["input"]["grid_variables"].doc( + "Variables to read from the grid file and make available in expressions"); + + for (const auto& [name, value] : field_variables) { if (not mesh.isDataSourceGridFile()) { - throw BoutException("A grid file ('mesh:file') is required for `input:{}`", - option_name); + throw BoutException( + "A grid file ('mesh:file') is required for `input:grid_variables`"); } - for (const auto& name : trimsplit(field_variables)) { - if (not mesh.sourceHasVar(name)) { - const auto filename = Options::root()["mesh"]["file"].as(); - throw BoutException("Grid file '{}' missing `{}` specified in `input:{}`", - filename, name, option_name); - } - T var; + if (not mesh.sourceHasVar(name)) { + const auto filename = Options::root()["mesh"]["file"].as(); + throw BoutException( + "Grid file '{}' missing `{}` specified in `input:grid_variables`", filename, + name); + } + + const auto func = value.as(); + switch (func) { + case GridVariableFunction::field3d: + add_grid_variable(factory, mesh, name); + break; + case GridVariableFunction::field2d: + add_grid_variable(factory, mesh, name); + break; + case GridVariableFunction::boutreal: + BoutReal var{}; mesh.get(var, name); - factory.addGenerator(name, std::make_shared>(var, name)); + factory.addGenerator(name, std::make_shared(var)); + break; } } } @@ -228,8 +233,7 @@ FieldFactory::FieldFactory(Mesh* localmesh, Options* opt) addGenerator("is_periodic_y", std::make_shared()); // Variables from the grid file - read_grid_variables(*this, *localmesh, nonconst_options, "read_Field3Ds"); - read_grid_variables(*this, *localmesh, nonconst_options, "read_Field2Ds"); + read_grid_variables(*this, *localmesh, nonconst_options); } Field2D FieldFactory::create2D(const std::string& value, const Options* opt, diff --git a/tests/unit/field/test_field_factory.cxx b/tests/unit/field/test_field_factory.cxx index 1f88d9ce7e..b6e916dc97 100644 --- a/tests/unit/field/test_field_factory.cxx +++ b/tests/unit/field/test_field_factory.cxx @@ -1053,8 +1053,9 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField3D) { } { - Options options{{"mesh", {{"file", filename.string()}}}, - {"input", {{"read_Field3Ds", "rho, theta"}}}}; + Options options{ + {"mesh", {{"file", filename.string()}}}, + {"input", {{"grid_variables", {{"rho", "field3d"}, {"theta", "field3d"}}}}}}; dynamic_cast(mesh)->setGridDataSource(new GridFile{filename}); auto factory = FieldFactory{mesh, &options}; @@ -1082,8 +1083,9 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField2D) { } { - Options options{{"mesh", {{"file", filename.string()}}}, - {"input", {{"read_Field2Ds", "rho, theta"}}}}; + Options options{ + {"mesh", {{"file", filename.string()}}}, + {"input", {{"grid_variables", {{"rho", "field2d"}, {"theta", "field2d"}}}}}}; dynamic_cast(mesh)->setGridDataSource(new GridFile{filename}); auto factory = FieldFactory{mesh, &options}; @@ -1094,8 +1096,33 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField2D) { } } +TEST_F(FieldFactoryFieldVariableTest, ReadBoutReal) { + bout::testing::TempFile filename; + + { + Options grid{{"rho", 4}, + {"theta", 5}, + {"nx", mesh->LocalNx}, + {"ny", mesh->LocalNy}, + {"nz", mesh->LocalNz}}; + bout::OptionsIO::create(filename)->write(grid); + } + + { + Options options{ + {"mesh", {{"file", filename.string()}}}, + {"input", {{"grid_variables", {{"rho", "boutreal"}, {"theta", "boutreal"}}}}}}; + + dynamic_cast(mesh)->setGridDataSource(new GridFile{filename}); + auto factory = FieldFactory{mesh, &options}; + + const auto output = factory.create3D("rho * theta"); + EXPECT_TRUE(IsFieldEqual(output, 4 * 5)); + } +} + TEST_F(FieldFactoryFieldVariableTest, NoMeshFile) { - Options options{{"input", {{"read_Field2Ds", "rho, theta"}}}}; + Options options{{"input", {{"grid_variables", {{"rho", "field3d"}}}}}}; EXPECT_THROW((FieldFactory(mesh, &options)), BoutException); } @@ -1115,8 +1142,9 @@ TEST_F(FieldFactoryFieldVariableTest, MissingVariable) { } { - Options options{{"mesh", {{"file", filename.string()}}}, - {"input", {{"read_Field3Ds", "rho, theta"}}}}; + Options options{ + {"mesh", {{"file", filename.string()}}}, + {"input", {{"grid_variables", {{"rho", "field3d"}, {"theta", "field3d"}}}}}}; dynamic_cast(mesh)->setGridDataSource(new GridFile{filename}); EXPECT_THROW((FieldFactory{mesh, &options}), BoutException); From 42b88b447af9d77dd752796ed010b907f36d5876 Mon Sep 17 00:00:00 2001 From: ZedThree <1486942+ZedThree@users.noreply.github.com> Date: Tue, 2 Dec 2025 17:16:30 +0000 Subject: [PATCH 07/19] Apply clang-format changes --- include/bout/sys/generator_context.hxx | 2 +- src/field/field_factory.cxx | 6 +++--- src/field/fieldgenerators.hxx | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/include/bout/sys/generator_context.hxx b/include/bout/sys/generator_context.hxx index a5dd9d236a..e39eea1835 100644 --- a/include/bout/sys/generator_context.hxx +++ b/include/bout/sys/generator_context.hxx @@ -43,7 +43,7 @@ public: /// Cell indices int ix() const { return ix_; } int jy() const { return jy_; } - int kz() const { return kz_;} + int kz() const { return kz_; } /// Set the value of a parameter with given name Context& set(const std::string& name, BoutReal value) { diff --git a/src/field/field_factory.cxx b/src/field/field_factory.cxx index 19e6a2fd14..fd06d90bd7 100644 --- a/src/field/field_factory.cxx +++ b/src/field/field_factory.cxx @@ -22,6 +22,8 @@ #include +#include +#include #include #include #include @@ -29,11 +31,9 @@ #include #include #include +#include #include #include -#include -#include -#include #include "fieldgenerators.hxx" diff --git a/src/field/fieldgenerators.hxx b/src/field/fieldgenerators.hxx index 49b2cb9690..a777738d02 100644 --- a/src/field/fieldgenerators.hxx +++ b/src/field/fieldgenerators.hxx @@ -311,7 +311,7 @@ public: // Constructor FieldTanhHat(FieldGeneratorPtr xin, FieldGeneratorPtr widthin, FieldGeneratorPtr centerin, FieldGeneratorPtr steepnessin) - : X(xin), width(widthin), center(centerin), steepness(steepnessin) {}; + : X(xin), width(widthin), center(centerin), steepness(steepnessin){}; // Clone containing the list of arguments FieldGeneratorPtr clone(const std::list args) override; BoutReal generate(const bout::generator::Context& pos) override; @@ -326,7 +326,7 @@ private: class FieldWhere : public FieldGenerator { public: FieldWhere(FieldGeneratorPtr test, FieldGeneratorPtr gt0, FieldGeneratorPtr lt0) - : test(test), gt0(gt0), lt0(lt0) {}; + : test(test), gt0(gt0), lt0(lt0){}; FieldGeneratorPtr clone(const std::list args) override { if (args.size() != 3) { From 5491001cc94efe60b1c155aea276eb8006c101cc Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 26 Mar 2026 15:39:37 -0700 Subject: [PATCH 08/19] FieldFactory: Use fieldmesh not localmesh The constructor argument `localmesh` can be null but in that case `fieldmesh` is set to `bout::globals::mesh`. --- src/field/field_factory.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/field/field_factory.cxx b/src/field/field_factory.cxx index fd06d90bd7..3a396b29fd 100644 --- a/src/field/field_factory.cxx +++ b/src/field/field_factory.cxx @@ -233,7 +233,7 @@ FieldFactory::FieldFactory(Mesh* localmesh, Options* opt) addGenerator("is_periodic_y", std::make_shared()); // Variables from the grid file - read_grid_variables(*this, *localmesh, nonconst_options); + read_grid_variables(*this, *fieldmesh, nonconst_options); } Field2D FieldFactory::create2D(const std::string& value, const Options* opt, From 7901046f56ab4329ecd356b416a37abd8afa6996 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 26 Mar 2026 15:41:13 -0700 Subject: [PATCH 09/19] GridVariable: Load grid variables lazily `BoutMesh::load` creates a FieldFactory to read variables such as mesh sizes. This then tries to load grid variables, which in turn starts trying to create coordinates, use regions etc., none of which has been constructed at that point. This defers reading until it's needed, at which point it reads the variable once and re-uses in subsequent calls. --- src/field/field_factory.cxx | 6 ++---- src/field/fieldgenerators.hxx | 38 ++++++++++++++++++++++++++-------- src/mesh/data/gridfromfile.cxx | 5 ++++- 3 files changed, 35 insertions(+), 14 deletions(-) diff --git a/src/field/field_factory.cxx b/src/field/field_factory.cxx index 3a396b29fd..793ea9f278 100644 --- a/src/field/field_factory.cxx +++ b/src/field/field_factory.cxx @@ -1,5 +1,5 @@ /************************************************************************** - * Copyright 2010-2025 BOUT++ contributors + * Copyright 2010 - 2026 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -94,9 +94,7 @@ class FieldIndirect : public FieldGenerator { // Read variables from the grid file and make them available in expressions template auto add_grid_variable(FieldFactory& factory, Mesh& mesh, const std::string& name) { - T var; - mesh.get(var, name); - factory.addGenerator(name, std::make_shared>(var, name)); + factory.addGenerator(name, std::make_shared>(&mesh, name)); } auto read_grid_variables(FieldFactory& factory, Mesh& mesh, Options& options) { diff --git a/src/field/fieldgenerators.hxx b/src/field/fieldgenerators.hxx index a777738d02..9a2be94849 100644 --- a/src/field/fieldgenerators.hxx +++ b/src/field/fieldgenerators.hxx @@ -381,29 +381,49 @@ public: } }; -/// A `Field3D` that can be used in expressions +/// A `Field3D` or `Field2D` that can be used in expressions +/// +/// The variable is read from Mesh when first used and shared between +/// clones. This is to avoid circular dependencies in construction. template > class GridVariable : public FieldGenerator { +private: + struct LazyLoaded { + LazyLoaded(Mesh* mesh, std::string name) : mesh(mesh), name(std::move(name)) {} + + Field3D get() { + if (!var.isAllocated()) { + this->mesh->get(this->var, this->name); + } + return this->var; + } + + Mesh* mesh; + std::string name; + T var {}; + }; + + std::shared_ptr variable; public: - GridVariable(T var, std::string name) - : variable(std::move(var)), name(std::move(name)) {} + GridVariable(Mesh* mesh, std::string name) + : variable(std::make_shared(mesh, std::move(name))) {} + + GridVariable(std::shared_ptr variable) + : variable(std::move(variable)) {} double generate(const bout::generator::Context& ctx) override { - return variable(ctx.ix(), ctx.jy(), ctx.kz()); + return variable->get()(ctx.ix(), ctx.jy(), ctx.kz()); } FieldGeneratorPtr clone(const std::list args) override { if (args.size() != 0) { throw ParseException("Variable '{}' takes no arguments but got {:d}", args.size()); } - return std::make_shared>(variable, name); + return std::make_shared>(variable); } - std::string str() const override { return name; } + std::string str() const override { return variable->name; } -private: - T variable; - std::string name; }; #endif // BOUT_FIELDGENERATORS_H diff --git a/src/mesh/data/gridfromfile.cxx b/src/mesh/data/gridfromfile.cxx index 27ad9411c2..35007bd1ff 100644 --- a/src/mesh/data/gridfromfile.cxx +++ b/src/mesh/data/gridfromfile.cxx @@ -197,7 +197,10 @@ bool GridFile::getField(Mesh* m, T& var, const std::string& name, BoutReal def, ///Ghost region widths. const int mxg = (m->LocalNx - (m->xend - m->xstart + 1)) / 2; const int myg = (m->LocalNy - (m->yend - m->ystart + 1)) / 2; - ///Check that ghost region widths are in fact integers + // Check grid has cells + ASSERT1(m->LocalNx > 0); + ASSERT1(m->LocalNy > 0); + // Check that ghost region widths are in fact integers ASSERT1((m->LocalNx - (m->xend - m->xstart + 1)) % 2 == 0); ASSERT1((m->LocalNy - (m->yend - m->ystart + 1)) % 2 == 0); From e9475ff2038e14139186a29cb2b39e70ea91af23 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 26 Mar 2026 21:10:36 -0700 Subject: [PATCH 10/19] FieldData: Don't getCoordinates in constructor No need to get coordinates before they are needed. Calling this can result in problems if mesh is not fully constructed. --- src/field/field_data.cxx | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/field/field_data.cxx b/src/field/field_data.cxx index b925cb1426..8f6a808482 100644 --- a/src/field/field_data.cxx +++ b/src/field/field_data.cxx @@ -53,11 +53,6 @@ FieldData::FieldData(Mesh* localmesh, CELL_LOC location_in) location_in, fieldmesh)) { // Need to check for nullptr again, because the // fieldmesh might still be // nullptr if the global mesh hasn't been initialized yet - if (fieldmesh != nullptr) { - // sets fieldCoordinates by getting Coordinates for our location from - // fieldmesh - getCoordinates(); - } } FieldData::FieldData(const FieldData& other) { From ab229707045de4fe736d3817b5d2d91c77d7260d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 26 Mar 2026 21:11:59 -0700 Subject: [PATCH 11/19] GridVariable: Check mesh->get return Not strictly necessary because variables are checked before creating a GridVariable, but doesn't hurt. --- src/field/fieldgenerators.hxx | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/field/fieldgenerators.hxx b/src/field/fieldgenerators.hxx index 9a2be94849..e8d38b1573 100644 --- a/src/field/fieldgenerators.hxx +++ b/src/field/fieldgenerators.hxx @@ -311,7 +311,7 @@ public: // Constructor FieldTanhHat(FieldGeneratorPtr xin, FieldGeneratorPtr widthin, FieldGeneratorPtr centerin, FieldGeneratorPtr steepnessin) - : X(xin), width(widthin), center(centerin), steepness(steepnessin){}; + : X(xin), width(widthin), center(centerin), steepness(steepnessin) {}; // Clone containing the list of arguments FieldGeneratorPtr clone(const std::list args) override; BoutReal generate(const bout::generator::Context& pos) override; @@ -326,7 +326,7 @@ private: class FieldWhere : public FieldGenerator { public: FieldWhere(FieldGeneratorPtr test, FieldGeneratorPtr gt0, FieldGeneratorPtr lt0) - : test(test), gt0(gt0), lt0(lt0){}; + : test(test), gt0(gt0), lt0(lt0) {}; FieldGeneratorPtr clone(const std::list args) override { if (args.size() != 3) { @@ -393,23 +393,26 @@ private: Field3D get() { if (!var.isAllocated()) { - this->mesh->get(this->var, this->name); + // Read variable from mesh + if (this->mesh->get(this->var, this->name) != 0) { + throw BoutException("Couldn't read GridVariable '{}'", this->name); + } } return this->var; } Mesh* mesh; std::string name; - T var {}; + T var{}; }; std::shared_ptr variable; + public: GridVariable(Mesh* mesh, std::string name) - : variable(std::make_shared(mesh, std::move(name))) {} + : variable(std::make_shared(mesh, std::move(name))) {} - GridVariable(std::shared_ptr variable) - : variable(std::move(variable)) {} + GridVariable(std::shared_ptr variable) : variable(std::move(variable)) {} double generate(const bout::generator::Context& ctx) override { return variable->get()(ctx.ix(), ctx.jy(), ctx.kz()); @@ -423,7 +426,6 @@ public: } std::string str() const override { return variable->name; } - }; #endif // BOUT_FIELDGENERATORS_H From d0726d99dfeac75e00006940330f05ce18e377ac Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 27 Mar 2026 16:57:30 -0700 Subject: [PATCH 12/19] FieldGenerators: Fixes and tidying Missing argument to exception, misc. headers. --- src/field/fieldgenerators.hxx | 38 ++++++++++++++++++----------------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/src/field/fieldgenerators.hxx b/src/field/fieldgenerators.hxx index e8d38b1573..c84a35109f 100644 --- a/src/field/fieldgenerators.hxx +++ b/src/field/fieldgenerators.hxx @@ -7,15 +7,20 @@ #ifndef BOUT_FIELDGENERATORS_H #define BOUT_FIELDGENERATORS_H +#include #include #include #include +#include #include #include +#include #include +#include #include #include +#include ////////////////////////////////////////////////////////// // Generators from values @@ -43,7 +48,7 @@ template class FieldGenOneArg : public FieldGenerator { ///< Template for single-argument function public: FieldGenOneArg(FieldGeneratorPtr g, const std::string& name = "function") - : gen(g), name(name) {} + : gen(std::move(g)), name(name) {} FieldGeneratorPtr clone(const std::list args) override { if (args.size() != 1) { throw ParseException("Incorrect number of arguments to {:s}. Expecting 1, got {:d}", @@ -70,7 +75,7 @@ class FieldGenTwoArg : public FieldGenerator { ///< Template for two-argument fu public: FieldGenTwoArg(FieldGeneratorPtr a, FieldGeneratorPtr b, const std::string& name = "function") - : A(a), B(b), name(name) {} + : A(std::move(a)), B(std::move(b)), name(name) {} FieldGeneratorPtr clone(const std::list args) override { if (args.size() != 2) { throw ParseException("Incorrect number of arguments to {:s}. Expecting 2, got {:d}", @@ -93,11 +98,13 @@ private: /// Arc (Inverse) tangent. Either one or two argument versions class FieldATan : public FieldGenerator { public: - FieldATan(FieldGeneratorPtr a, FieldGeneratorPtr b = nullptr) : A(a), B(b) {} + FieldATan(FieldGeneratorPtr a, FieldGeneratorPtr b = nullptr) + : A(std::move(a)), B(std::move(b)) {} FieldGeneratorPtr clone(const std::list args) override { if (args.size() == 1) { return std::make_shared(args.front()); - } else if (args.size() == 2) { + } + if (args.size() == 2) { return std::make_shared(args.front(), args.back()); } throw ParseException( @@ -148,7 +155,7 @@ public: FieldMin() = default; FieldMin(const std::list args) : input(args) {} FieldGeneratorPtr clone(const std::list args) override { - if (args.size() == 0) { + if (args.empty()) { throw ParseException("min function must have some inputs"); } return std::make_shared(args); @@ -157,10 +164,7 @@ public: auto it = input.begin(); BoutReal result = (*it)->generate(pos); for (; it != input.end(); it++) { - BoutReal val = (*it)->generate(pos); - if (val < result) { - result = val; - } + result = std::min(result, (*it)->generate(pos)); } return result; } @@ -175,7 +179,7 @@ public: FieldMax() = default; FieldMax(const std::list args) : input(args) {} FieldGeneratorPtr clone(const std::list args) override { - if (args.size() == 0) { + if (args.empty()) { throw ParseException("max function must have some inputs"); } return std::make_shared(args); @@ -184,10 +188,7 @@ public: auto it = input.begin(); BoutReal result = (*it)->generate(pos); for (; it != input.end(); it++) { - BoutReal val = (*it)->generate(pos); - if (val > result) { - result = val; - } + result = std::max(result, (*it)->generate(pos)); } return result; } @@ -206,7 +207,7 @@ class FieldClamp : public FieldGenerator { public: FieldClamp() = default; FieldClamp(FieldGeneratorPtr value, FieldGeneratorPtr low, FieldGeneratorPtr high) - : value(value), low(low), high(high) {} + : value(std::move(value)), low(std::move(low)), high(std::move(high)) {} FieldGeneratorPtr clone(const std::list args) override { if (args.size() != 3) { throw ParseException( @@ -242,7 +243,7 @@ private: /// Generator to round to the nearest integer class FieldRound : public FieldGenerator { public: - FieldRound(FieldGeneratorPtr g) : gen(g) {} + FieldRound(FieldGeneratorPtr g) : gen(std::move(g)) {} FieldGeneratorPtr clone(const std::list args) override { if (args.size() != 1) { @@ -419,8 +420,9 @@ public: } FieldGeneratorPtr clone(const std::list args) override { - if (args.size() != 0) { - throw ParseException("Variable '{}' takes no arguments but got {:d}", args.size()); + if (!args.empty()) { + throw ParseException("Variable '{}' takes no arguments but got {:d}", + variable->name, args.size()); } return std::make_shared>(variable); } From 2e967789c084339eebd4121926c62f3dd0fdb9d1 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 27 Mar 2026 21:53:47 -0700 Subject: [PATCH 13/19] GridVariable unit tests: Fix unit test Unit tests now include Y guard cells. --- src/mesh/data/gridfromfile.cxx | 1 + src/mesh/impls/bout/boutmesh.cxx | 3 +- tests/unit/field/test_field_factory.cxx | 57 ++++++++++++------------- 3 files changed, 30 insertions(+), 31 deletions(-) diff --git a/src/mesh/data/gridfromfile.cxx b/src/mesh/data/gridfromfile.cxx index 35007bd1ff..813ee1e3d8 100644 --- a/src/mesh/data/gridfromfile.cxx +++ b/src/mesh/data/gridfromfile.cxx @@ -3,6 +3,7 @@ #include #include +#include #include #include #include diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index f56979785c..3c05425241 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -2573,9 +2573,8 @@ bool BoutMesh::periodicY(int jx, BoutReal& ts) const { int BoutMesh::numberOfYBoundaries() const { if (jyseps2_1 != jyseps1_2) { return 2; - } else { - return 1; } + return 1; } std::pair BoutMesh::hasBranchCutLower(int jx) const { diff --git a/tests/unit/field/test_field_factory.cxx b/tests/unit/field/test_field_factory.cxx index b6e916dc97..c7b5bcd9ae 100644 --- a/tests/unit/field/test_field_factory.cxx +++ b/tests/unit/field/test_field_factory.cxx @@ -12,11 +12,14 @@ #include "bout/options_io.hxx" #include "bout/output.hxx" #include "bout/paralleltransform.hxx" +#include "bout/sys/generator_context.hxx" #include "bout/traits.hxx" #include "fake_mesh_fixture.hxx" #include "test_tmpfiles.hxx" +#include + // The unit tests use the global mesh using namespace bout::globals; using bout::generator::Context; @@ -885,7 +888,7 @@ TEST_F(FieldFactoryCreateAndTransformTest, Create2D) { mesh->getCoordinates()->setParallelTransform( bout::utils::make_unique(*mesh, true)); - FieldFactory factory; + const FieldFactory factory; auto output = factory.create2D("x"); @@ -900,7 +903,7 @@ TEST_F(FieldFactoryCreateAndTransformTest, Create3D) { mesh->getCoordinates()->setParallelTransform( bout::utils::make_unique(*mesh, true)); - FieldFactory factory; + const FieldFactory factory; auto output = factory.create3D("x"); @@ -916,7 +919,7 @@ TEST_F(FieldFactoryCreateAndTransformTest, Create2DNoTransform) { Options options; options["input"]["transform_from_field_aligned"] = false; - FieldFactory factory{mesh, &options}; + const FieldFactory factory{mesh, &options}; auto output = factory.create2D("x"); @@ -933,7 +936,7 @@ TEST_F(FieldFactoryCreateAndTransformTest, Create3DNoTransform) { Options options; options["input"]["transform_from_field_aligned"] = false; - FieldFactory factory{mesh, &options}; + const FieldFactory factory{mesh, &options}; auto output = factory.create3D("x"); @@ -947,7 +950,7 @@ TEST_F(FieldFactoryCreateAndTransformTest, Create2DCantTransform) { mesh->getCoordinates()->setParallelTransform( bout::utils::make_unique(*mesh, false)); - FieldFactory factory{mesh}; + const FieldFactory factory{mesh}; auto output = factory.create2D("x"); @@ -1041,14 +1044,12 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField3D) { { // Write some fields to a grid file - FieldFactory factory{mesh}; + const FieldFactory factory{mesh}; const auto rho = factory.create3D("sqrt(x^2 + y^2)"); const auto theta = factory.create3D("atan(y, x)"); - Options grid{{"rho", rho}, - {"theta", theta}, - {"nx", mesh->LocalNx}, - {"ny", mesh->LocalNy}, - {"nz", mesh->LocalNz}}; + const Options grid{{"rho", rho}, {"theta", theta}, + {"nx", mesh->LocalNx}, {"ny", mesh->LocalNy - 2}, + {"nz", mesh->LocalNz}, {"y_boundary_guards", 1}}; bout::OptionsIO::create(filename)->write(grid); } @@ -1062,7 +1063,7 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField3D) { const auto output = factory.create3D("rho * cos(theta)"); const auto x = factory.create3D("x"); - EXPECT_TRUE(IsFieldEqual(output, x)); + EXPECT_TRUE(IsFieldEqual(output, x, "RGN_ALL", 1e-14)); } } @@ -1071,14 +1072,12 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField2D) { { // Write some fields to a grid file - FieldFactory factory{mesh}; + const FieldFactory factory{mesh}; const auto rho = factory.create2D("sqrt(x^2 + y^2)"); const auto theta = factory.create2D("atan(y, x)"); - Options grid{{"rho", rho}, - {"theta", theta}, - {"nx", mesh->LocalNx}, - {"ny", mesh->LocalNy}, - {"nz", mesh->LocalNz}}; + const Options grid{{"rho", rho}, {"theta", theta}, + {"nx", mesh->LocalNx}, {"ny", mesh->LocalNy - 2}, + {"nz", mesh->LocalNz}, {"y_boundary_guards", 1}}; bout::OptionsIO::create(filename)->write(grid); } @@ -1092,7 +1091,7 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField2D) { const auto output = factory.create2D("rho * cos(theta)"); const auto x = factory.create2D("x"); - EXPECT_TRUE(IsFieldEqual(output, x)); + EXPECT_TRUE(IsFieldEqual(output, x, "RGN_ALL", 1e-14)); } } @@ -1100,11 +1099,11 @@ TEST_F(FieldFactoryFieldVariableTest, ReadBoutReal) { bout::testing::TempFile filename; { - Options grid{{"rho", 4}, - {"theta", 5}, - {"nx", mesh->LocalNx}, - {"ny", mesh->LocalNy}, - {"nz", mesh->LocalNz}}; + const Options grid{{"rho", 4}, + {"theta", 5}, + {"nx", mesh->LocalNx}, + {"ny", mesh->LocalNy}, + {"nz", mesh->LocalNz}}; bout::OptionsIO::create(filename)->write(grid); } @@ -1132,12 +1131,12 @@ TEST_F(FieldFactoryFieldVariableTest, MissingVariable) { { // Write some fields to a grid file - FieldFactory factory{mesh}; + const FieldFactory factory{mesh}; const auto rho = factory.create3D("sqrt(x^2 + y^2)"); - Options grid{{"rho", rho}, - {"nx", mesh->LocalNx}, - {"ny", mesh->LocalNy}, - {"nz", mesh->LocalNz}}; + const Options grid{{"rho", rho}, + {"nx", mesh->LocalNx}, + {"ny", mesh->LocalNy}, + {"nz", mesh->LocalNz}}; bout::OptionsIO::create(filename)->write(grid); } From 5e18e164cb00d3a6f9a7f3f6e2c90a7fbc65c523 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 28 Mar 2026 09:12:18 -0700 Subject: [PATCH 14/19] Field2D/3D: Defer sizing until allocation FieldFactory can create Fields when bout::mesh::global is non-null but the size of the mesh (LocalNx etc.) is not set, because FieldFactory is created inside BoutMesh::load() via Options. --- include/bout/field2d.hxx | 2 +- include/bout/field3d.hxx | 2 +- include/bout/mesh.hxx | 2 +- src/field/field2d.cxx | 24 +++++++++-------------- src/field/field3d.cxx | 41 ++++++++++++---------------------------- 5 files changed, 24 insertions(+), 47 deletions(-) diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx index efa878ef0a..374945d47e 100644 --- a/include/bout/field2d.hxx +++ b/include/bout/field2d.hxx @@ -289,7 +289,7 @@ private: Array data; /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null - int nx{-1}, ny{-1}; + int nx{0}, ny{0}; /// Time-derivative, can be nullptr Field2D* deriv{nullptr}; diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index fabe2b646d..422523ed2e 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -551,7 +551,7 @@ public: protected: /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null - int nx{-1}, ny{-1}, nz{-1}; + int nx{0}, ny{0}, nz{0}; /// Internal data array. Handles allocation/freeing of memory Array data; diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 1453f51f5f..ba2a1c640f 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -605,7 +605,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; diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx index b363eeef07..e246ad4c63 100644 --- a/src/field/field2d.cxx +++ b/src/field/field2d.cxx @@ -46,12 +46,8 @@ Field2D::Field2D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in, [[maybe_unused]] std::optional regionID) : Field(localmesh, location_in, directions_in) { - - if (fieldmesh) { - nx = fieldmesh->LocalNx; - ny = fieldmesh->LocalNy; - } - + // Note: Even if fieldmesh is not null, LocalNx and LocalNy may not + // be initialised. #if BOUT_USE_TRACK name = ""; #endif @@ -62,11 +58,6 @@ Field2D::Field2D(const Field2D& f) : Field(f), data(f.data) { #if BOUT_USE_TRACK name = f.name; #endif - - if (fieldmesh) { - nx = fieldmesh->LocalNx; - ny = fieldmesh->LocalNy; - } } Field2D::Field2D(BoutReal val, Mesh* localmesh) : Field2D(localmesh) { *this = val; } @@ -89,13 +80,16 @@ Field2D::~Field2D() { delete deriv; } Field2D& Field2D::allocate() { if (data.empty()) { - if (!fieldmesh) { + if (fieldmesh == nullptr) { // fieldmesh was not initialized when this field was initialized, so use - // the global mesh and set some members to default values + // the global mesh fieldmesh = bout::globals::mesh; - nx = fieldmesh->LocalNx; - ny = fieldmesh->LocalNy; } + // Get size from the mesh. + nx = fieldmesh->LocalNx; + ny = fieldmesh->LocalNy; + ASSERT1(nx > 0); + ASSERT1(ny > 0); data.reallocate(nx * ny); #if CHECK > 2 invalidateGuards(*this); diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 633eb42e9b..a449335f29 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -63,38 +63,21 @@ Field3D::Field3D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes direction #if BOUT_USE_TRACK name = ""; #endif - - if (fieldmesh) { - nx = fieldmesh->LocalNx; - ny = fieldmesh->LocalNy; - nz = fieldmesh->LocalNz; - } } /// Doesn't copy any data, just create a new reference to the same data (copy on change /// later) Field3D::Field3D(const Field3D& f) : Field(f), data(f.data), yup_fields(f.yup_fields), ydown_fields(f.ydown_fields), - regionID(f.regionID) { + regionID(f.regionID) {} - if (fieldmesh) { - nx = fieldmesh->LocalNx; - ny = fieldmesh->LocalNy; - nz = fieldmesh->LocalNz; - } -} - -Field3D::Field3D(const Field2D& f) : Field(f) { - - nx = fieldmesh->LocalNx; - ny = fieldmesh->LocalNy; - nz = fieldmesh->LocalNz; +Field3D::Field3D(const Field2D& f) + : Field(f), nx(fieldmesh->LocalNx), ny(fieldmesh->LocalNy), nz(fieldmesh->LocalNz) { *this = f; } Field3D::Field3D(const BoutReal val, Mesh* localmesh) : Field3D(localmesh) { - *this = val; } @@ -113,11 +96,8 @@ Field3DParallel::Field3DParallel(const BoutReal val, Mesh* localmesh) Field3D::Field3D(Array data_in, Mesh* localmesh, CELL_LOC datalocation, DirectionTypes directions_in) - : Field(localmesh, datalocation, directions_in), data(std::move(data_in)) { - - nx = fieldmesh->LocalNx; - ny = fieldmesh->LocalNy; - nz = fieldmesh->LocalNz; + : Field(localmesh, datalocation, directions_in), nx(fieldmesh->LocalNx), + ny(fieldmesh->LocalNy), nz(fieldmesh->LocalNz), data(std::move(data_in)) { ASSERT1(data.size() == nx * ny * nz); } @@ -126,14 +106,17 @@ Field3D::~Field3D() { delete deriv; } Field3D& Field3D::allocate() { if (data.empty()) { - if (!fieldmesh) { + if (fieldmesh == nullptr) { // fieldmesh was not initialized when this field was initialized, so use // the global mesh and set some members to default values fieldmesh = bout::globals::mesh; - nx = fieldmesh->LocalNx; - ny = fieldmesh->LocalNy; - nz = fieldmesh->LocalNz; } + nx = fieldmesh->LocalNx; + ny = fieldmesh->LocalNy; + nz = fieldmesh->LocalNz; + ASSERT1(nx > 0); + ASSERT1(ny > 0); + ASSERT1(nz > 0); data.reallocate(nx * ny * nz); #if CHECK > 2 invalidateGuards(*this); From 95b01d8bfca8ccbb9e140f18f44b68a55ae86a11 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 28 Mar 2026 09:15:05 -0700 Subject: [PATCH 15/19] Miscellaneous tidying Addressing clang-tidy comments and Clang compiler warnings (sprintf use). --- include/bout/field2d.hxx | 1 + src/field/field2d.cxx | 1 + src/field/field3d.cxx | 2 + .../impls/multigrid/multigrid_laplace.cxx | 29 +++++++------- .../impls/multigrid/multigrid_solver.cxx | 38 ++++++++----------- src/invert/laplace/impls/pcr/pcr.cxx | 23 ++++++----- 6 files changed, 46 insertions(+), 48 deletions(-) diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx index 374945d47e..5d07ea9aba 100644 --- a/include/bout/field2d.hxx +++ b/include/bout/field2d.hxx @@ -43,6 +43,7 @@ class Field2D; #include #include +#include #include #if BOUT_HAS_RAJA diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx index e246ad4c63..a0b2b440e5 100644 --- a/src/field/field2d.cxx +++ b/src/field/field2d.cxx @@ -29,6 +29,7 @@ #include "bout/build_defines.hxx" #include "bout/unused.hxx" +#include #include #include #include diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index a449335f29..1241dee465 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -25,8 +25,10 @@ * **************************************************************************/ +#include "bout/array.hxx" #include "bout/bout_types.hxx" #include "bout/build_defines.hxx" +#include "bout/field2d.hxx" #include "bout/index_derivs_interface.hxx" #include diff --git a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx index 4be5bfbad3..c06a83585b 100644 --- a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx +++ b/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx @@ -1,5 +1,5 @@ /************************************************************************** - * Perpendicular Laplacian inversion. + * Perpendicular Laplacian inversion. * Using Geometrical Multigrid Solver * * Equation solved is: @@ -9,7 +9,7 @@ * Copyright 2016 K.S. Kang * * Contact: Ben Dudson, bd512@york.ac.uk - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -32,8 +32,11 @@ #if not BOUT_USE_METRIC_3D +#include +#include #include #include +#include #if BOUT_USE_OPENMP #include @@ -107,11 +110,7 @@ LaplaceMultigrid::LaplaceMultigrid(Options* opt, const CELL_LOC loc, Mesh* mesh_ } Nz_global = localmesh->GlobalNz; Nz_local = Nz_global; // No parallelization in z-direction (for now) - // - //else { - // Nz_local = localmesh->zend - localmesh->zstart + 1; // excluding guard cells - // Nz_global = localmesh->GlobalNz - 2*localmesh->zstart; // excluding guard cells - // } + if (mgcount == 0) { output << "Nz=" << Nz_global << "(" << Nz_local << ")" << endl; } @@ -214,7 +213,9 @@ LaplaceMultigrid::LaplaceMultigrid(Options* opt, const CELL_LOC loc, Mesh* mesh_ } BOUT_OMP_SAFE(parallel) BOUT_OMP_SAFE(master) - { output << "Num threads = " << omp_get_num_threads() << endl; } + { + output << "Num threads = " << omp_get_num_threads() << endl; + } } } @@ -384,11 +385,9 @@ FieldPerp LaplaceMultigrid::solve(const FieldPerp& b_in, const FieldPerp& x0) { } if ((pcheck == 3) && (mgcount == 0)) { - FILE* outf; - char outfile[256]; - sprintf(outfile, "test_matF_%d.mat", kMG->rProcI); + std::string outfile = fmt::format("test_matF_{:d}.mat", kMG->rProcI); output << "Out file= " << outfile << endl; - outf = fopen(outfile, "w"); + FILE* outf = fopen(outfile.c_str(), "w"); int dim = (lxx + 2) * (lzz + 2); fprintf(outf, "dim = %d (%d, %d)\n", dim, lxx, lzz); @@ -411,11 +410,9 @@ FieldPerp LaplaceMultigrid::solve(const FieldPerp& b_in, const FieldPerp& x0) { output << i << "dimension= " << kMG->lnx[i - 1] << "(" << kMG->gnx[i - 1] << ")," << kMG->lnz[i - 1] << endl; - FILE* outf; - char outfile[256]; - sprintf(outfile, "test_matC%1d_%d.mat", i, kMG->rProcI); + std::string outfile = fmt::format("test_matC{:1d}_{:d}.mat", i, kMG->rProcI); output << "Out file= " << outfile << endl; - outf = fopen(outfile, "w"); + FILE* outf = fopen(outfile.c_str(), "w"); int dim = (kMG->lnx[i - 1] + 2) * (kMG->lnz[i - 1] + 2); fprintf(outf, "dim = %d (%d,%d)\n", dim, kMG->lnx[i - 1], kMG->lnz[i - 1]); diff --git a/src/invert/laplace/impls/multigrid/multigrid_solver.cxx b/src/invert/laplace/impls/multigrid/multigrid_solver.cxx index 5ff616aa54..074e482320 100644 --- a/src/invert/laplace/impls/multigrid/multigrid_solver.cxx +++ b/src/invert/laplace/impls/multigrid/multigrid_solver.cxx @@ -1,5 +1,5 @@ /************************************************************************** - * Perpendicular Laplacian inversion. + * Perpendicular Laplacian inversion. * Using Geometrical Multigrid Solver * * Equation solved is: @@ -9,7 +9,7 @@ * Copyright 2015 K.S. Kang * * Contact: Ben Dudson, bd512@york.ac.uk - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -35,6 +35,9 @@ #include "bout/unused.hxx" #include +#include +#include + Multigrid1DP::Multigrid1DP(int level, int lx, int lz, int gx, int dl, int merge, MPI_Comm comm, int check) : MultigridAlg(level, lx, lz, gx, lz, comm, check) { @@ -195,11 +198,9 @@ void Multigrid1DP::setMultigridC(int UNUSED(plag)) { if (pcheck == 2) { for (int i = level; i >= 0; i--) { - FILE* outf; - char outfile[256]; - sprintf(outfile, "2DP_matC%1d_%d.mat", i, rMG->rProcI); + std::string outfile = fmt::format("2DP_matC{:1d}_{:d}.mat", i, rMG->rProcI); output << "Out file= " << outfile << endl; - outf = fopen(outfile, "w"); + FILE* outf = fopen(outfile.c_str(), "w"); int dim = (rMG->lnx[i] + 2) * (rMG->lnz[i] + 2); fprintf(outf, "dim = %d (%d, %d)\n", dim, rMG->lnx[i], rMG->lnz[i]); @@ -222,11 +223,9 @@ void Multigrid1DP::setMultigridC(int UNUSED(plag)) { } if (pcheck == 3) { for (int i = level; i >= 0; i--) { - FILE* outf; - char outfile[256]; - sprintf(outfile, "S1D_matC%1d_%d.mat", i, sMG->rProcI); + std::string outfile = fmt::format("S1D_matC{:1d}_{:d}.mat", i, sMG->rProcI); output << "Out file= " << outfile << endl; - outf = fopen(outfile, "w"); + FILE* outf = fopen(outfile.c_str(), "w"); int dim = (sMG->lnx[i] + 2) * (sMG->lnz[i] + 2); fprintf(outf, "dim = %d\n", dim); @@ -456,11 +455,9 @@ void Multigrid1DP::convertMatrixF2D(int level) { } } if (pcheck == 3) { - FILE* outf; - char outfile[256]; - sprintf(outfile, "2DP_CP_%d.mat", rProcI); + std::string outfile = fmt::format("2DP_CP_{}.mat", rProcI); output << "Out file= " << outfile << endl; - outf = fopen(outfile, "w"); + FILE* outf = fopen(outfile.c_str(), "w"); fprintf(outf, "dim = (%d, %d)\n", ggx, gnz[0]); for (int ii = 0; ii < dim; ii++) { @@ -476,11 +473,10 @@ void Multigrid1DP::convertMatrixF2D(int level) { MPI_SUM, comm2D); if (pcheck == 3) { - FILE* outf; - char outfile[256]; - sprintf(outfile, "2DP_Conv_%d.mat", rProcI); + + std::string outfile = fmt::format("2DP_Conv_{}.mat", rProcI); output << "Out file= " << outfile << endl; - outf = fopen(outfile, "w"); + FILE* outf = fopen(outfile.c_str(), "w"); fprintf(outf, "dim = (%d, %d)\n", ggx, gnz[0]); for (int ii = 0; ii < dim; ii++) { @@ -625,10 +621,8 @@ void Multigrid2DPf1D::setMultigridC(int UNUSED(plag)) { } if (pcheck == 2) { for (int i = level; i >= 0; i--) { - FILE* outf; - char outfile[256]; - sprintf(outfile, "S2D_matC%1d_%d.mat", i, sMG->rProcI); - outf = fopen(outfile, "w"); + std::string outfile = fmt::format("S2D_matC{:1d}_{:d}.mat", i, sMG->rProcI); + FILE* outf = fopen(outfile.c_str(), "w"); output << "Out file= " << outfile << endl; int dim = (sMG->lnx[i] + 2) * (sMG->lnz[i] + 2); fprintf(outf, "dim = %d\n", dim); diff --git a/src/invert/laplace/impls/pcr/pcr.cxx b/src/invert/laplace/impls/pcr/pcr.cxx index 93fbf09cdb..325242a3b2 100644 --- a/src/invert/laplace/impls/pcr/pcr.cxx +++ b/src/invert/laplace/impls/pcr/pcr.cxx @@ -45,11 +45,14 @@ #include "../../common_transform.hxx" #include +#include +#include #include #include #include #include #include +#include #include #include #include @@ -124,8 +127,8 @@ LaplacePCR::LaplacePCR(Options* opt, CELL_LOC loc, Mesh* mesh_in, Solver* UNUSED // (unless periodic in x) xe = localmesh->LocalNx - 1; } - int n = xe - xs + 1; // Number of X points on this processor, - // including boundaries but not guard cells + const int n = xe - xs + 1; // Number of X points on this processor, + // including boundaries but not guard cells a.reallocate(nmode, n); b.reallocate(nmode, n); @@ -148,7 +151,7 @@ FieldPerp LaplacePCR::solve(const FieldPerp& rhs, const FieldPerp& x0) { FieldPerp x{emptyFrom(rhs)}; // Result - int jy = rhs.getIndex(); // Get the Y index + const int jy = rhs.getIndex(); // Get the Y index x.setIndex(jy); // Get the width of the boundary @@ -215,7 +218,7 @@ Field3D LaplacePCR::solve(const Field3D& rhs, const Field3D& x0) { outbndry = 1; } - int nx = xe - xs + 1; // Number of X points on this processor + const int nx = xe - xs + 1; // Number of X points on this processor // Get range of Y indices int ys = localmesh->ystart; @@ -313,7 +316,7 @@ void LaplacePCR ::cr_pcr_solver(Matrix& a_mpi, Matrix& b_mpi // don't want to copy them. // xs = xstart if a proc has no boundary points // xs = 0 if a proc has boundary points - int offset = localmesh->xstart - xs; + const int offset = localmesh->xstart - xs; aa(kz, ix + 1) = a_mpi(kz, ix + offset); bb(kz, ix + 1) = b_mpi(kz, ix + offset); cc(kz, ix + 1) = c_mpi(kz, ix + offset); @@ -400,7 +403,7 @@ void LaplacePCR ::apply_boundary_conditions(const Matrix& a, } } if (localmesh->lastX()) { - int n = xe - xs + 1; // actual length of array + const int n = xe - xs + 1; // actual length of array for (int kz = 0; kz < nsys; kz++) { for (int ix = n - localmesh->xstart; ix < n; ix++) { x(kz, ix) = (r(kz, ix) - a(kz, ix) * x(kz, ix - 1)) / b(kz, ix); @@ -419,7 +422,7 @@ void LaplacePCR ::cr_forward_multiple_row(Matrix& a, Matrix& Matrix& r) const { const int nsys = std::get<0>(a.shape()); - MPI_Comm comm = BoutComm::get(); + const MPI_Comm comm = BoutComm::get(); Array alpha(nsys); Array gamma(nsys); Array sbuf(4 * nsys); @@ -494,7 +497,7 @@ void LaplacePCR ::cr_backward_multiple_row(Matrix& a, Matrix Matrix& x) const { const int nsys = std::get<0>(a.shape()); - MPI_Comm comm = BoutComm::get(); + const MPI_Comm comm = BoutComm::get(); MPI_Status status; MPI_Request request[2]; @@ -533,7 +536,7 @@ void LaplacePCR ::cr_backward_multiple_row(Matrix& a, Matrix dist_row = dist_row / 2; } if (xproc < nprocs - 1) { - MPI_Wait(request + 1, &status); + MPI_Wait(&request[1], &status); } } @@ -555,7 +558,7 @@ void LaplacePCR ::pcr_forward_single_row(Matrix& a, Matrix& MPI_Status status; Array request(4); - MPI_Comm comm = BoutComm::get(); + const MPI_Comm comm = BoutComm::get(); const int nlevel = log2(nprocs); const int nhprocs = nprocs / 2; From fe00722b5bf016bb919b6cb3b929ccef93a4c245 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 28 Mar 2026 10:38:56 -0700 Subject: [PATCH 16/19] Field2D/3D: Partly revert 8ebc0a0 Set field sizes in the same way, but get and check sizes on allocation. --- include/bout/field2d.hxx | 2 +- include/bout/field3d.hxx | 2 +- src/field/field2d.cxx | 15 +++++++++++---- src/field/field3d.cxx | 14 +++++++++++++- 4 files changed, 26 insertions(+), 7 deletions(-) diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx index 5d07ea9aba..0b437e1bd8 100644 --- a/include/bout/field2d.hxx +++ b/include/bout/field2d.hxx @@ -290,7 +290,7 @@ private: Array data; /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null - int nx{0}, ny{0}; + int nx{-1}, ny{-1}; /// Time-derivative, can be nullptr Field2D* deriv{nullptr}; diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 422523ed2e..fabe2b646d 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -551,7 +551,7 @@ public: protected: /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null - int nx{0}, ny{0}, nz{0}; + int nx{-1}, ny{-1}, nz{-1}; /// Internal data array. Handles allocation/freeing of memory Array data; diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx index a0b2b440e5..9cad2652a0 100644 --- a/src/field/field2d.cxx +++ b/src/field/field2d.cxx @@ -4,7 +4,7 @@ * Class for 2D X-Y profiles * ************************************************************************** - * Copyright 2010 - 2025 BOUT++ developers + * Copyright 2010 - 2026 BOUT++ developers * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -47,15 +47,22 @@ Field2D::Field2D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in, [[maybe_unused]] std::optional regionID) : Field(localmesh, location_in, directions_in) { - // Note: Even if fieldmesh is not null, LocalNx and LocalNy may not - // be initialised. + if (fieldmesh != nullptr) { + // Note: Even if fieldmesh is not null, LocalNx and LocalNy may + // not be initialised. + nx = fieldmesh->LocalNx; + ny = fieldmesh->LocalNy; + } #if BOUT_USE_TRACK name = ""; #endif } Field2D::Field2D(const Field2D& f) : Field(f), data(f.data) { - + if (fieldmesh != nullptr) { + nx = fieldmesh->LocalNx; + ny = fieldmesh->LocalNy; + } #if BOUT_USE_TRACK name = f.name; #endif diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 1241dee465..61bc9ce86a 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -65,13 +65,25 @@ Field3D::Field3D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes direction #if BOUT_USE_TRACK name = ""; #endif + + if (fieldmesh != nullptr) { + nx = fieldmesh->LocalNx; + ny = fieldmesh->LocalNy; + nz = fieldmesh->LocalNz; + } } /// Doesn't copy any data, just create a new reference to the same data (copy on change /// later) Field3D::Field3D(const Field3D& f) : Field(f), data(f.data), yup_fields(f.yup_fields), ydown_fields(f.ydown_fields), - regionID(f.regionID) {} + regionID(f.regionID) { + if (fieldmesh != nullptr) { + nx = fieldmesh->LocalNx; + ny = fieldmesh->LocalNy; + nz = fieldmesh->LocalNz; + } +} Field3D::Field3D(const Field2D& f) : Field(f), nx(fieldmesh->LocalNx), ny(fieldmesh->LocalNy), nz(fieldmesh->LocalNz) { From cd04c41650c3d5527cf1fbbb61adbb8756de4484 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 28 Mar 2026 11:22:00 -0700 Subject: [PATCH 17/19] FieldFactory: Don't test values in Y boundary --- tests/unit/field/test_field_factory.cxx | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/tests/unit/field/test_field_factory.cxx b/tests/unit/field/test_field_factory.cxx index c7b5bcd9ae..7f31c83a3b 100644 --- a/tests/unit/field/test_field_factory.cxx +++ b/tests/unit/field/test_field_factory.cxx @@ -2,8 +2,10 @@ #include "fake_mesh.hxx" #include "test_extras.hxx" +#include "bout/bout_types.hxx" #include "bout/boutexception.hxx" #include "bout/constants.hxx" +#include "bout/coordinates.hxx" #include "bout/field2d.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" @@ -12,13 +14,16 @@ #include "bout/options_io.hxx" #include "bout/output.hxx" #include "bout/paralleltransform.hxx" +#include "bout/sys/expressionparser.hxx" #include "bout/sys/generator_context.hxx" #include "bout/traits.hxx" +#include "bout/utils.hxx" #include "fake_mesh_fixture.hxx" #include "test_tmpfiles.hxx" #include +#include // The unit tests use the global mesh using namespace bout::globals; @@ -777,7 +782,7 @@ TEST_F(FieldFactoryTest, Recursion) { opt["input"]["max_recursion_depth"] = 4; // Should be sufficient for n=6 // Create a factory with a max_recursion_depth != 0 - FieldFactory factory_rec(nullptr, &opt); + const FieldFactory factory_rec(nullptr, &opt); // Fibonacci sequence: 1 1 2 3 5 8 opt["fib"] = "where({n} - 2.5, [n={n}-1](fib) + [n={n}-2](fib), 1)"; @@ -809,7 +814,7 @@ TEST_F(FieldFactoryTest, ResolveLocalOptions) { options["f"] = "2 + 2"; options["g"] = "f * f"; - FieldFactoryExposer factory_local(mesh, &options); + const FieldFactoryExposer factory_local(mesh, &options); auto g = factory_local.resolve("g"); EXPECT_EQ(g->generate({}), 16); @@ -1040,7 +1045,7 @@ struct FieldFactoryFieldVariableTest : public FakeMeshFixture { }; TEST_F(FieldFactoryFieldVariableTest, CreateField3D) { - bout::testing::TempFile filename; + const bout::testing::TempFile filename; { // Write some fields to a grid file @@ -1063,12 +1068,12 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField3D) { const auto output = factory.create3D("rho * cos(theta)"); const auto x = factory.create3D("x"); - EXPECT_TRUE(IsFieldEqual(output, x, "RGN_ALL", 1e-14)); + EXPECT_TRUE(IsFieldEqual(output, x, "RGN_NOBNDRY", 1e-14)); } } TEST_F(FieldFactoryFieldVariableTest, CreateField2D) { - bout::testing::TempFile filename; + const bout::testing::TempFile filename; { // Write some fields to a grid file @@ -1096,7 +1101,7 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField2D) { } TEST_F(FieldFactoryFieldVariableTest, ReadBoutReal) { - bout::testing::TempFile filename; + const bout::testing::TempFile filename; { const Options grid{{"rho", 4}, @@ -1127,7 +1132,7 @@ TEST_F(FieldFactoryFieldVariableTest, NoMeshFile) { } TEST_F(FieldFactoryFieldVariableTest, MissingVariable) { - bout::testing::TempFile filename; + const bout::testing::TempFile filename; { // Write some fields to a grid file From 787538420396167dafd700b302a50899dc4f1274 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 28 Mar 2026 11:55:11 -0700 Subject: [PATCH 18/19] Clang formatting --- include/bout/field2d.hxx | 2 +- include/bout/invertable_operator.hxx | 4 +- include/bout/options.hxx | 24 +++++----- include/bout/output.hxx | 8 ++-- include/bout/physicsmodel.hxx | 47 +++++++++---------- include/bout/sys/generator_context.hxx | 2 +- include/bout/utils.hxx | 18 +++---- src/invert/laplace/common_transform.cxx | 24 +++++----- .../laplace/impls/serial_tri/serial_tri.hxx | 2 +- src/invert/laplace/impls/spt/spt.hxx | 14 +++--- .../laplacexy/impls/hypre/laplacexy-hypre.cxx | 6 ++- src/mesh/coordinates_accessor.cxx | 6 ++- src/mesh/parallel/fci_comm.hxx | 4 +- src/mesh/parallel/shiftedmetric.cxx | 4 +- src/solver/impls/petsc/petsc.cxx | 20 ++++---- src/solver/impls/snes/snes.hxx | 6 +-- src/sys/msg_stack.cxx | 8 +++- src/sys/options/options_netcdf.cxx | 40 ++++++++-------- 18 files changed, 123 insertions(+), 116 deletions(-) diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx index 0b437e1bd8..5143889d4c 100644 --- a/include/bout/field2d.hxx +++ b/include/bout/field2d.hxx @@ -42,8 +42,8 @@ class Field2D; #include "bout/region.hxx" #include -#include #include +#include #include #if BOUT_HAS_RAJA diff --git a/include/bout/invertable_operator.hxx b/include/bout/invertable_operator.hxx index d61c258654..6dcd92897a 100644 --- a/include/bout/invertable_operator.hxx +++ b/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() { diff --git a/include/bout/options.hxx b/include/bout/options.hxx index 03e95488a8..48840a7b6f 100644 --- a/include/bout/options.hxx +++ b/include/bout/options.hxx @@ -74,19 +74,19 @@ 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(); * * // 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"); * @@ -94,7 +94,7 @@ class Options; * 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 §ion = 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 diff --git a/include/bout/output.hxx b/include/bout/output.hxx index 9416ab411b..43fc6c59bf 100644 --- a/include/bout/output.hxx +++ b/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 -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 -ConditionalOutput& operator<<(ConditionalOutput& out, T const& t) { +ConditionalOutput& operator<<(ConditionalOutput& out, const T& t) { if (out.isEnabled()) { *out.getBase() << t; } diff --git a/include/bout/physicsmodel.hxx b/include/bout/physicsmodel.hxx index d653d96d68..8e8e3814cd 100644 --- a/include/bout/physicsmodel.hxx +++ b/include/bout/physicsmodel.hxx @@ -1,20 +1,20 @@ /*!************************************************************************ * \file physicsmodel.hxx - * + * * @brief Base class for Physics Models - * - * + * + * * * Changelog: - * + * * 2013-08 Ben Dudson * * Initial version - * + * ************************************************************************** * Copyright 2013 B.D.Dudson * * Contact: Ben Dudson, bd512@york.ac.uk - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -167,9 +167,9 @@ public: * * Output * ------ - * + * * The time derivatives will be put in the ddt() variables - * + * * Returns a flag: 0 indicates success, non-zero an error flag */ int runRHS(BoutReal time, bool linear = false); @@ -203,7 +203,7 @@ public: bool hasPrecon(); /*! - * Run the preconditioner. The system state should be in the + * Run the preconditioner. The system state should be in the * evolving variables, and the vector to be solved in the ddt() variables. * The result will be put in the ddt() variables. * @@ -219,7 +219,7 @@ public: /*! * Run the Jacobian-vector multiplication function - * + * * Note: this is usually only called by the Solver */ int runJacobian(BoutReal t); @@ -241,10 +241,10 @@ protected: // The init and rhs functions are implemented by user code to specify problem /*! * @brief This function is called once by the solver at the start of a simulation. - * + * * A valid PhysicsModel must implement this function - * - * Variables should be read from the inputs, and the variables to + * + * Variables should be read from the inputs, and the variables to * be evolved should be specified. */ virtual int init(bool restarting) = 0; @@ -258,7 +258,7 @@ protected: /*! * @brief This function is called by the time integration solver * at least once per time step - * + * * Variables being evolved will be set by the solver * before the call, and this function must calculate * and set the time-derivatives. @@ -278,10 +278,10 @@ protected: /// Add additional variables other than the evolving variables to the restart files virtual void restartVars(Options& options); - /* + /* If split operator is set to true, then convective() and diffusive() are called instead of rhs() - + For implicit-explicit schemes, convective() will typically be treated explicitly, whilst diffusive() will be treated implicitly. For unsplit methods, both convective and diffusive will be called @@ -334,7 +334,7 @@ protected: * * @param[in] var The variable to evolve * @param[in] name The name to use for variable initialisation and output - * + * * Note that the variable must not be destroyed (e.g. go out of scope) * after this call, since a pointer to \p var is stored in the solver. * @@ -358,11 +358,11 @@ protected: * Specify a constrained variable \p var, which will be * adjusted to make \p F_var equal to zero. * If the solver does not support constraints then this will throw an exception - * + * * @param[in] var The variable the solver should modify * @param[in] F_var The control variable, which the user will set * @param[in] name The name to use for initialisation and output - * + * */ bool bout_constrain(Field3D& var, Field3D& F_var, const char* name); @@ -491,8 +491,7 @@ private: /// Add fields to the solver. /// This should accept up to ten arguments -#define SOLVE_FOR(...) \ - { MACRO_FOR_EACH(SOLVE_FOR1, __VA_ARGS__) } +#define SOLVE_FOR(...) {MACRO_FOR_EACH(SOLVE_FOR1, __VA_ARGS__)} /// Write this variable once to the grid file #define SAVE_ONCE1(var) dump.addOnce(var, #var); @@ -532,8 +531,7 @@ private: dump.addOnce(var6, #var6); \ } -#define SAVE_ONCE(...) \ - { MACRO_FOR_EACH(SAVE_ONCE1, __VA_ARGS__) } +#define SAVE_ONCE(...) {MACRO_FOR_EACH(SAVE_ONCE1, __VA_ARGS__)} /// Write this variable every timestep #define SAVE_REPEAT1(var) dump.addRepeat(var, #var); @@ -573,7 +571,6 @@ private: dump.addRepeat(var6, #var6); \ } -#define SAVE_REPEAT(...) \ - { MACRO_FOR_EACH(SAVE_REPEAT1, __VA_ARGS__) } +#define SAVE_REPEAT(...) {MACRO_FOR_EACH(SAVE_REPEAT1, __VA_ARGS__)} #endif // BOUT_PHYSICS_MODEL_H diff --git a/include/bout/sys/generator_context.hxx b/include/bout/sys/generator_context.hxx index e39eea1835..080b4d8ed7 100644 --- a/include/bout/sys/generator_context.hxx +++ b/include/bout/sys/generator_context.hxx @@ -33,7 +33,7 @@ public: /// The location on the boundary Context(const BoundaryRegion* bndry, int iz, CELL_LOC loc, BoutReal t, Mesh* msh); Context(const BoundaryRegion* bndry, CELL_LOC loc, BoutReal t, Mesh* msh) - : Context(bndry, 0, loc, t, msh){}; + : Context(bndry, 0, loc, t, msh) {}; BoutReal x() const { return get("x"); } BoutReal y() const { return get("y"); } diff --git a/include/bout/utils.hxx b/include/bout/utils.hxx index 5088a025c1..eae3b27730 100644 --- a/include/bout/utils.hxx +++ b/include/bout/utils.hxx @@ -8,7 +8,7 @@ * Copyright 2010 - 2026 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -484,7 +484,7 @@ inline bool is_pow2(int x) { return x && !((x - 1) & x); } /*! * Return the sign of a number \p a - * by testing if a > 0 + * by testing if a > 0 */ template T SIGN(T a) { // Return +1 or -1 (0 -> +1) @@ -511,7 +511,7 @@ inline void checkData(BoutReal f) { } #else /// Ignored with disabled CHECK; Throw an exception if \p f is not finite -inline void checkData(BoutReal UNUSED(f)){}; +inline void checkData(BoutReal UNUSED(f)) {}; #endif /*! @@ -587,7 +587,7 @@ BoutReal stringToReal(const std::string& s); /*! * Convert a string to an int - * + * * Throws BoutException if can't be done */ int stringToInt(const std::string& s); @@ -604,7 +604,7 @@ std::list& strsplit(const std::string& s, char delim, /*! * Split a string on a given delimiter - * + * * @param[in] s The string to split (not modified by call) * @param[in] delim The delimiter to split on (single char) */ @@ -612,7 +612,7 @@ std::list strsplit(const std::string& s, char delim); /*! * Strips leading and trailing spaces from a string - * + * * @param[in] s The string to trim (not modified) * @param[in] c Collection of characters to remove */ @@ -620,7 +620,7 @@ std::string trim(const std::string& s, const std::string& c = " \t\r"); /*! * Strips leading spaces from a string - * + * * @param[in] s The string to trim (not modified) * @param[in] c Collection of characters to remove */ @@ -628,7 +628,7 @@ std::string trimLeft(const std::string& s, const std::string& c = " \t"); /*! * Strips leading spaces from a string - * + * * @param[in] s The string to trim (not modified) * @param[in] c Collection of characters to remove */ @@ -636,7 +636,7 @@ std::string trimRight(const std::string& s, const std::string& c = " \t\r"); /*! * Strips the comments from a string - * + * * @param[in] s The string to trim (not modified) * @param[in] c Collection of characters to remove */ diff --git a/src/invert/laplace/common_transform.cxx b/src/invert/laplace/common_transform.cxx index 847add4fea..98571624a1 100644 --- a/src/invert/laplace/common_transform.cxx +++ b/src/invert/laplace/common_transform.cxx @@ -28,8 +28,8 @@ FFTTransform::FFTTransform(const Mesh& mesh, int nmode, int xs, int xe, int ys, auto FFTTransform::forward(const Laplacian& laplacian, const Field3D& rhs, const Field3D& x0, const Field2D& Acoef, const Field2D& C1coef, - const Field2D& C2coef, - const Field2D& Dcoef) const -> Matrices { + const Field2D& C2coef, const Field2D& Dcoef) const + -> Matrices { Matrices result(nsys, nx); @@ -79,8 +79,8 @@ auto FFTTransform::forward(const Laplacian& laplacian, const Field3D& rhs, return result; } -auto FFTTransform::backward(const Field3D& rhs, - const Matrix& xcmplx3D) const -> Field3D { +auto FFTTransform::backward(const Field3D& rhs, const Matrix& xcmplx3D) const + -> Field3D { Field3D x{emptyFrom(rhs)}; // FFT back to real space @@ -162,8 +162,8 @@ auto FFTTransform::forward(const Laplacian& laplacian, const FieldPerp& rhs, return result; } -auto FFTTransform::backward(const FieldPerp& rhs, - const Matrix& xcmplx) const -> FieldPerp { +auto FFTTransform::backward(const FieldPerp& rhs, const Matrix& xcmplx) const + -> FieldPerp { FieldPerp x{emptyFrom(rhs)}; // FFT back to real space @@ -209,8 +209,8 @@ DSTTransform::DSTTransform(const Mesh& mesh, int nmode, int xs, int xe, int ys, auto DSTTransform::forward(const Laplacian& laplacian, const Field3D& rhs, const Field3D& x0, const Field2D& Acoef, const Field2D& C1coef, - const Field2D& C2coef, - const Field2D& Dcoef) const -> Matrices { + const Field2D& C2coef, const Field2D& Dcoef) const + -> Matrices { Matrices result(nsys, nx); @@ -260,8 +260,8 @@ auto DSTTransform::forward(const Laplacian& laplacian, const Field3D& rhs, return result; } -auto DSTTransform::backward(const Field3D& rhs, - const Matrix& xcmplx3D) const -> Field3D { +auto DSTTransform::backward(const Field3D& rhs, const Matrix& xcmplx3D) const + -> Field3D { Field3D x{emptyFrom(rhs)}; // DST back to real space @@ -346,8 +346,8 @@ auto DSTTransform::forward(const Laplacian& laplacian, const FieldPerp& rhs, return result; } -auto DSTTransform::backward(const FieldPerp& rhs, - const Matrix& xcmplx) const -> FieldPerp { +auto DSTTransform::backward(const FieldPerp& rhs, const Matrix& xcmplx) const + -> FieldPerp { FieldPerp x{emptyFrom(rhs)}; // DST back to real space diff --git a/src/invert/laplace/impls/serial_tri/serial_tri.hxx b/src/invert/laplace/impls/serial_tri/serial_tri.hxx index 4aed777b7c..332c272f3f 100644 --- a/src/invert/laplace/impls/serial_tri/serial_tri.hxx +++ b/src/invert/laplace/impls/serial_tri/serial_tri.hxx @@ -52,7 +52,7 @@ class LaplaceSerialTri : public Laplacian { public: LaplaceSerialTri(Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh* mesh_in = nullptr, Solver* solver = nullptr); - ~LaplaceSerialTri(){}; + ~LaplaceSerialTri() {}; using Laplacian::setCoefA; void setCoefA(const Field2D& val) override { diff --git a/src/invert/laplace/impls/spt/spt.hxx b/src/invert/laplace/impls/spt/spt.hxx index a3b1acb7c3..e951f431e7 100644 --- a/src/invert/laplace/impls/spt/spt.hxx +++ b/src/invert/laplace/impls/spt/spt.hxx @@ -1,16 +1,16 @@ /************************************************************************** - * Perpendicular Laplacian inversion. + * Perpendicular Laplacian inversion. * PARALLEL CODE - SIMPLE ALGORITHM - * + * * I'm just calling this Simple Parallel Tridag. Naive parallelisation of * the serial code. For use as a reference case. - * + * * Overlap calculation / communication of poloidal slices to achieve some * parallelism. * * Changelog * --------- - * + * * 2014-06 Ben Dudson * * Removed static variables in functions, changing to class members. * @@ -18,7 +18,7 @@ * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu * * Contact: Ben Dudson, bd512@york.ac.uk - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -142,8 +142,8 @@ private: Array buffer; }; - int ys, ye; // Range of Y indices - SPT_data slicedata; // Used to solve for a single FieldPerp + int ys, ye; // Range of Y indices + SPT_data slicedata; // Used to solve for a single FieldPerp Array alldata; // Used to solve a Field3D Array dc1d; ///< 1D in Z for taking FFTs diff --git a/src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx b/src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx index 439d07a4e5..61632a332d 100644 --- a/src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx +++ b/src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx @@ -23,8 +23,10 @@ #include #if BOUT_HAS_CUDA && defined(__CUDACC__) -#define gpuErrchk(ans) \ - { gpuAssert((ans), __FILE__, __LINE__); } +#define gpuErrchk(ans) \ + { \ + gpuAssert((ans), __FILE__, __LINE__); \ + } inline void gpuAssert(cudaError_t code, const char* file, int line, bool abort = true) { if (code != cudaSuccess) { fprintf(stderr, "GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line); diff --git a/src/mesh/coordinates_accessor.cxx b/src/mesh/coordinates_accessor.cxx index 0ce4b664b5..c874ea3ac2 100644 --- a/src/mesh/coordinates_accessor.cxx +++ b/src/mesh/coordinates_accessor.cxx @@ -45,8 +45,10 @@ CoordinatesAccessor::CoordinatesAccessor(const Coordinates* coords) { data[stripe_size * ind.ind + static_cast(Offset::symbol)] = coords->symbol[ind]; // Implement copy for each argument -#define COPY_STRIPE(...) \ - { MACRO_FOR_EACH(COPY_STRIPE1, __VA_ARGS__) } +#define COPY_STRIPE(...) \ + { \ + MACRO_FOR_EACH(COPY_STRIPE1, __VA_ARGS__) \ + } // Iterate over all points in the field // Note this could be 2D or 3D, depending on FieldMetric type diff --git a/src/mesh/parallel/fci_comm.hxx b/src/mesh/parallel/fci_comm.hxx index 324cae8a22..34512b18c4 100644 --- a/src/mesh/parallel/fci_comm.hxx +++ b/src/mesh/parallel/fci_comm.hxx @@ -67,7 +67,7 @@ struct ProcLocal { struct GlobalToLocal1D { GlobalToLocal1D(int mg, int npe, int localwith, bool periodic) : mg(mg), npe(npe), localwith(localwith), local(localwith - (2 * mg)), - global(local * npe), globalwith(global + (2 * mg)), periodic(periodic){}; + global(local * npe), globalwith(global + (2 * mg)), periodic(periodic) {}; ProcLocal convert(int id) const; int getLocalWith() const { return localwith; } int getGlobalWith() const { return globalwith; } @@ -104,7 +104,7 @@ public: const BoutReal& operator[](IndG3D ind) const; GlobalField3DAccessInstance(const GlobalField3DAccess* gfa, std::vector&& data) - : gfa(gfa), data(std::move(data)){}; + : gfa(gfa), data(std::move(data)) {}; private: const GlobalField3DAccess* gfa; diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx index 64c6d9a2ce..f167824541 100644 --- a/src/mesh/parallel/shiftedmetric.cxx +++ b/src/mesh/parallel/shiftedmetric.cxx @@ -39,8 +39,8 @@ void ShiftedMetric::checkInputGrid() { "Should be 'shiftedmetric'."); } } // else: parallel_transform variable not found in grid input, indicates older input - // file or grid from options so must rely on the user having ensured the type is - // correct + // file or grid from options so must rely on the user having ensured the type is + // correct } void ShiftedMetric::outputVars(Options& output_options) { diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 1e1e05596b..cc884c240c 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -64,10 +64,10 @@ class ColoringStencil { private: - bool static isInSquare(int const i, int const j, int const n_square) { + bool static isInSquare(const int i, const int j, const int n_square) { return std::abs(i) <= n_square && std::abs(j) <= n_square; } - bool static isInCross(int const i, int const j, int const n_cross) { + bool static isInCross(const int i, const int j, const int n_cross) { if (i == 0) { return std::abs(j) <= n_cross; } @@ -76,7 +76,7 @@ class ColoringStencil { } return false; } - bool static isInTaxi(int const i, int const j, int const n_taxi) { + bool static isInTaxi(const int i, const int j, const int n_taxi) { return std::abs(i) + std::abs(j) <= n_taxi; } @@ -614,7 +614,7 @@ int PetscSolver::init() { n_taxi = 2; } - auto const xy_offsets = ColoringStencil::getOffsets(n_square, n_taxi, n_cross); + const auto xy_offsets = ColoringStencil::getOffsets(n_square, n_taxi, n_cross); { // This is ugly but can't think of a better and robust way to // count the non-zeros for some arbitrary stencil @@ -734,7 +734,7 @@ int PetscSolver::init() { // Mark non-zero entries output_progress.write("Marking non-zero Jacobian entries\n"); - PetscScalar const val = 1.0; + const PetscScalar val = 1.0; for (int x = mesh->xstart; x <= mesh->xend; x++) { for (int y = mesh->ystart; y <= mesh->yend; y++) { @@ -753,21 +753,21 @@ int PetscSolver::init() { continue; } - int const ind2 = ROUND(index(xi, yi, 0)); + const int ind2 = ROUND(index(xi, yi, 0)); if (ind2 < 0) { continue; // A boundary point } // Depends on all variables on this cell for (int j = 0; j < n2d; j++) { - PetscInt const col = ind2 + j; + const PetscInt col = ind2 + j; PetscCall(MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES)); } } } // 3D fields for (int z = mesh->zstart; z <= mesh->zend; z++) { - int const ind = ROUND(index(x, y, z)); + const int ind = ROUND(index(x, y, z)); for (int i = 0; i < n3d; i++) { PetscInt row = ind + i; @@ -777,7 +777,7 @@ int PetscSolver::init() { // Depends on 2D fields for (int j = 0; j < n2d; j++) { - PetscInt const col = ind0 + j; + const PetscInt col = ind0 + j; PetscCall(MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES)); } @@ -802,7 +802,7 @@ int PetscSolver::init() { // 3D fields on this cell for (int j = 0; j < n3d; j++) { - PetscInt const col = ind2 + j; + const PetscInt col = ind2 + j; int ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); if (ierr != 0) { diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index fbcf3baddf..28f9881cab 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -202,8 +202,8 @@ private: BoutReal kI; ///< (0.2 - 0.4) Integral parameter (smooths history of changes) BoutReal kD; ///< (0.1 - 0.3) Derivative (dampens oscillation - optional) bool pid_consider_failures; ///< Reduce timestep increases if recent solves have failed - BoutReal recent_failure_rate; ///< Rolling average of recent failure rate - BoutReal last_failure_weight; ///< 1 / number of recent solves + BoutReal recent_failure_rate; ///< Rolling average of recent failure rate + BoutReal last_failure_weight; ///< 1 / number of recent solves BoutReal nl_its_prev; BoutReal nl_its_prev2; @@ -245,7 +245,7 @@ private: bool matrix_free_operator; ///< Use matrix free Jacobian in the operator? int lag_jacobian; ///< Re-use Jacobian bool jacobian_persists; ///< Re-use Jacobian and preconditioner across nonlinear solves - bool use_coloring; ///< Use matrix coloring + bool use_coloring; ///< Use matrix coloring bool jacobian_recalculated; ///< Flag set when Jacobian is recalculated bool prune_jacobian; ///< Remove small elements in the Jacobian? diff --git a/src/sys/msg_stack.cxx b/src/sys/msg_stack.cxx index 3dbd7c2797..bd14a766ed 100644 --- a/src/sys/msg_stack.cxx +++ b/src/sys/msg_stack.cxx @@ -59,7 +59,9 @@ void MsgStack::pop() { return; } BOUT_OMP_SAFE(single) - { --position; } + { + --position; + } } void MsgStack::pop(int id) { @@ -87,7 +89,9 @@ void MsgStack::clear() { void MsgStack::dump() { BOUT_OMP_SAFE(single) - { output << this->getDump(); } + { + output << this->getDump(); + } } std::string MsgStack::getDump() { diff --git a/src/sys/options/options_netcdf.cxx b/src/sys/options/options_netcdf.cxx index f21f461c09..9170f7c284 100644 --- a/src/sys/options/options_netcdf.cxx +++ b/src/sys/options/options_netcdf.cxx @@ -132,25 +132,27 @@ void readGroup(const std::string& filename, const NcGroup& group, Options& resul {s2i(dims[0].getSize()), s2i(dims[1].getSize()), s2i(dims[2].getSize())}); // We need to explicitly copy file, so that there is a pointer to the file, and // the file does not get closed, which would prevent us from reading. - result[var_name].setLazyLoad(std::make_unique( - int, int, int, int, int, int)>>( - [file, var](int xstart, int xend, int ystart, int yend, int zstart, - int zend) { - const auto i2s = [](int i) { - if (i < 0) { - throw BoutException("BadCast {} < 0", i); - } - return static_cast(i); - }; - Tensor value(xend - xstart + 1, yend - ystart + 1, - zend - zstart + 1); - const std::vector index{i2s(xstart), i2s(ystart), i2s(zstart)}; - const std::vector count{i2s(xend - xstart + 1), - i2s(yend - ystart + 1), - i2s(zend - zstart + 1)}; - var.getVar(index, count, value.begin()); - return value; - })); + result[var_name].setLazyLoad( + std::make_unique< + std::function(int, int, int, int, int, int)>>( + [file, var](int xstart, int xend, int ystart, int yend, int zstart, + int zend) { + const auto i2s = [](int i) { + if (i < 0) { + throw BoutException("BadCast {} < 0", i); + } + return static_cast(i); + }; + Tensor value(xend - xstart + 1, yend - ystart + 1, + zend - zstart + 1); + const std::vector index{i2s(xstart), i2s(ystart), + i2s(zstart)}; + const std::vector count{i2s(xend - xstart + 1), + i2s(yend - ystart + 1), + i2s(zend - zstart + 1)}; + var.getVar(index, count, value.begin()); + return value; + })); } else { Tensor value(static_cast(dims[0].getSize()), static_cast(dims[1].getSize()), From f94a1e12b5ef7daf3995828c07b1e15f5af0face Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 19 Jun 2026 14:22:11 +0100 Subject: [PATCH 19/19] Use `gridvar:` namespace for input variables read from grid --- manual/sphinx/user_docs/input_grids.rst | 82 +++++++++++++------------ src/field/field_factory.cxx | 4 +- tests/unit/field/test_field_factory.cxx | 4 +- 3 files changed, 47 insertions(+), 43 deletions(-) diff --git a/manual/sphinx/user_docs/input_grids.rst b/manual/sphinx/user_docs/input_grids.rst index 402fb30d06..337fb16938 100644 --- a/manual/sphinx/user_docs/input_grids.rst +++ b/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 `_ 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,11 +404,11 @@ 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) @@ -414,9 +416,9 @@ 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 #include "fieldgenerators.hxx" +#include "fmt/format.h" #include #include @@ -94,7 +95,8 @@ class FieldIndirect : public FieldGenerator { // Read variables from the grid file and make them available in expressions template auto add_grid_variable(FieldFactory& factory, Mesh& mesh, const std::string& name) { - factory.addGenerator(name, std::make_shared>(&mesh, name)); + factory.addGenerator(fmt::format("gridvar:{}", name), + std::make_shared>(&mesh, name)); } auto read_grid_variables(FieldFactory& factory, Mesh& mesh, Options& options) { diff --git a/tests/unit/field/test_field_factory.cxx b/tests/unit/field/test_field_factory.cxx index 7f31c83a3b..25177c3169 100644 --- a/tests/unit/field/test_field_factory.cxx +++ b/tests/unit/field/test_field_factory.cxx @@ -1066,7 +1066,7 @@ TEST_F(FieldFactoryFieldVariableTest, CreateField3D) { dynamic_cast(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(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)); }