diff --git a/CMakeLists.txt b/CMakeLists.txt index c45fca3b72..acbaef9551 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -105,7 +105,6 @@ set(BOUT_SOURCES ./include/bout/deriv_store.hxx ./include/bout/derivs.hxx ./include/bout/difops.hxx - ./include/bout/expr.hxx ./include/bout/fft.hxx ./include/bout/field.hxx ./include/bout/field2d.hxx diff --git a/include/bout/boundary_region.hxx b/include/bout/boundary_region.hxx index 58de12045e..77a0d2f37c 100644 --- a/include/bout/boundary_region.hxx +++ b/include/bout/boundary_region.hxx @@ -8,11 +8,6 @@ class BoundaryRegion; #include class Mesh; -namespace bout { -namespace globals { -extern Mesh* mesh; ///< Global mesh -} // namespace globals -} // namespace bout /// Location of boundary enum class BndryLoc { @@ -40,10 +35,9 @@ class BoundaryRegionBase { public: BoundaryRegionBase() = delete; BoundaryRegionBase(std::string name, Mesh* passmesh = nullptr) - : localmesh(passmesh ? passmesh : bout::globals::mesh), label(std::move(name)) {} + : localmesh(passmesh), label(std::move(name)) {} BoundaryRegionBase(std::string name, BndryLoc loc, Mesh* passmesh = nullptr) - : localmesh(passmesh ? passmesh : bout::globals::mesh), label(std::move(name)), - location(loc) {} + : localmesh(passmesh), label(std::move(name)), location(loc) {} virtual ~BoundaryRegionBase() = default; diff --git a/include/bout/expr.hxx b/include/bout/expr.hxx deleted file mode 100644 index 267af202ed..0000000000 --- a/include/bout/expr.hxx +++ /dev/null @@ -1,208 +0,0 @@ -/************************************************************************** - * - * Operators, and support for template expressions - * - * Originally based on article by Klaus Kreft & Angelika Langer - * http://www.angelikalanger.com/Articles/Cuj/ExpressionTemplates/ExpressionTemplates.htm - * - * Parts adapted from Blitz++ library - * - **************************************************************************/ - -#ifndef BOUT_EXPR_H -#define BOUT_EXPR_H - -#warning expr.hxx is deprecated. Do not use! - -#include -#include -#include - -/// Literal class to capture BoutReal values in expressions -class Literal { -public: - /// Type of this expression - using type = Literal; - - Literal(BoutReal v) : val(v) {} - ~Literal() {} - BoutReal operator()(int x, int y, int z) const { return val; } - -private: - const BoutReal val; -}; - -class Field3DExpr { -public: - using type = Field3D; - - Field3DExpr(const Field3D& f) : data(&f(0, 0, 0)) {} - const BoutReal& operator()(int x, int y, int z) const { - return data[(x * bout::globals::mesh->LocalNy + y) * bout::globals::mesh->LocalNz - + z]; - } - -private: - const BoutReal* data; -}; - -class Field2DExpr { -public: - using type = Field2D; - - Field2DExpr(const Field2D& f) : data(&f(0, 0)) {} - const BoutReal& operator()(int x, int y, int z) const { - return data[x * bout::globals::mesh->LocalNy + y]; - } - -private: - const BoutReal* data; -}; - -/// Expression traits, to convert doubles etc. to Literal - -template -struct exprTraits { - using expr_type = ExprT; -}; - -template <> -struct exprTraits { - using expr_type = Literal; -}; - -template <> -struct exprTraits { - using expr_type = Literal; -}; - -template <> -struct exprTraits { - using expr_type = Literal; -}; - -/////////////////////////////////////////////// -// asExpr: convert objects to expressions - -template -struct asExpr { - using type = T; - static const T& getExpr(const T& x) { return x; } -}; - -template <> -struct asExpr { - using type = Literal; - static const Literal getExpr(const int& x) { return Literal(x); } -}; - -template <> -struct asExpr { - using type = Literal; - static const Literal getExpr(const double& x) { return Literal(x); } -}; - -template <> -struct asExpr { - using type = Literal; - static const Literal getExpr(const float& x) { return Literal(x); } -}; - -template <> -struct asExpr { - using type = Field3DExpr; - static const Field3DExpr getExpr(const Field3D& x) { return Field3DExpr(x); } -}; - -///////////////////////////////////////////////////////////// -// Type promotion. Work out the type of a calculation, -// based on the type of the arguments - -template // If in doubt, convert to Field3D -struct PromoteType { - using type = Field3D; -}; - -///////////////////////////////////////////////////////////// -// Binary expressions - -template -class BinaryExpr { -public: - BinaryExpr(const ExprT1& e1, const ExprT2& e2) : _expr1(e1), _expr2(e2) {} - - // Work out the type of the inputs - using ltype = typename exprTraits::expr_type; - using rtype = typename exprTraits::expr_type; - - /// Type of the resulting expression - using type = typename PromoteType::type; - - BoutReal operator()(int x, int y, int z) const { - return BinOp::apply((_expr1)(x, y, z), (_expr2)(x, y, z)); - } - -private: - ltype const _expr1; - rtype const _expr2; -}; - -template -struct BinaryResult { - using arg1 = typename asExpr::type; - using arg2 = typename asExpr::type; - using type = BinaryExpr; -}; - -/// Binary operator classes - -#define DEFINE_BINARY_OP(name, op) \ - struct name { \ - template \ - static inline T apply(T a, T b) { \ - return a op b; \ - } \ - }; - -DEFINE_BINARY_OP(Add, +) -DEFINE_BINARY_OP(Subtract, -) -DEFINE_BINARY_OP(Multiply, *) -DEFINE_BINARY_OP(Divide, /) - -struct Power { - template - static inline T apply(T a, T b) { - return pow(a, b); - } -}; - -/// Define functions add, mul which use operator structs -#define DEFINE_OVERLOAD_FUNC(name, func) \ - template \ - typename BinaryResult::type func(const ExprT1& e1, \ - const ExprT2& e2) { \ - using type = typename BinaryResult::type; \ - return type(asExpr::getExpr(e1), asExpr::getExpr(e2)); \ - } - -/// Addition of two Expressions -DEFINE_OVERLOAD_FUNC(Add, add); -/// Multiplication of two Expressions -DEFINE_OVERLOAD_FUNC(Multiply, mul); - -/// A function to evaluate expressions -template -const Field3D eval3D(Expr e) { - Field3D result; - result.allocate(); - for (int i = 0; i < bout::globals::mesh->LocalNx; i++) { - for (int j = 0; j < bout::globals::mesh->LocalNy; j++) { - for (int k = 0; k < bout::globals::mesh->LocalNz; k++) { - result(i, j, k) = e(i, j, k); - } - } - } - return result; -} - -#endif // BOUT_EXPR_H diff --git a/include/bout/field.hxx b/include/bout/field.hxx index 61edff6723..10d13c3b20 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -3,10 +3,10 @@ * \brief field base class definition for differencing methods * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors + * + * Contact: Ben Dudson, dudson2@llnl.gov * - * Contact: Ben Dudson, bd512@york.ac.uk - * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -41,7 +41,6 @@ class Field; #include "bout/region.hxx" #include "bout/traits.hxx" #include "bout/utils.hxx" -#include #include class Mesh; diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx index 92658f1bbf..1ca4f6efb3 100644 --- a/include/bout/field2d.hxx +++ b/include/bout/field2d.hxx @@ -4,9 +4,9 @@ * \brief Definition of 2D scalar field class * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010-2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -55,16 +55,12 @@ class Field2D : public Field { public: using ind_type = Ind2D; /*! - * Constructor, taking an optional mesh pointer - * This mesh pointer is not used until the data is allocated, - * since Field2D objects can be globals, created before a mesh - * has been created. + * Constructor * * @param[in] localmesh The mesh which defines the field size. * - * By default the global Mesh pointer (mesh) is used. */ - Field2D(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE, + Field2D(Mesh* localmesh, CELL_LOC location_in = CELL_CENTRE, DirectionTypes directions_in = {YDirectionType::Standard, ZDirectionType::Average}); @@ -84,7 +80,7 @@ public: * allocates data, and assigns the value \p val to all points including * boundary cells. */ - Field2D(BoutReal val, Mesh* localmesh = nullptr); + Field2D(BoutReal val, Mesh* localmesh); /// Constructor from Array and Mesh Field2D(Array data, Mesh* localmesh, CELL_LOC location = CELL_CENTRE, diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index a75e38df36..ab6035ada5 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -1,7 +1,7 @@ /************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -161,10 +161,8 @@ public: /*! * Constructor * - * Note: the global "mesh" can't be passed here because - * fields may be created before the mesh is. */ - Field3D(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE, + Field3D(Mesh* localmesh, CELL_LOC location_in = CELL_CENTRE, DirectionTypes directions_in = {YDirectionType::Standard, ZDirectionType::Standard}); @@ -178,7 +176,7 @@ public: /// Constructor from 2D field Field3D(const Field2D& f); /// Constructor from value - Field3D(BoutReal val, Mesh* localmesh = nullptr); + Field3D(BoutReal val, Mesh* localmesh); /// Constructor from Array and Mesh Field3D(Array data, Mesh* localmesh, CELL_LOC location = CELL_CENTRE, DirectionTypes directions_in = {YDirectionType::Standard, diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 94007a57a2..38f809cca2 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -502,7 +502,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { n = toFieldAligned(n_in, "RGN_NOX"); Field3D vy = toFieldAligned(v.y, "RGN_NOX"); - Field3D yresult = 0.0; + Field3D yresult{0.0, mesh}; yresult.setDirectionY(YDirectionType::Aligned); BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { diff --git a/include/bout/globals.hxx b/include/bout/globals.hxx index 64b3a09ee3..363165bf4a 100644 --- a/include/bout/globals.hxx +++ b/include/bout/globals.hxx @@ -3,9 +3,9 @@ * * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -29,7 +29,6 @@ #include "bout/macro_for_each.hxx" -class Mesh; class MpiWrapper; namespace bout { @@ -42,54 +41,8 @@ namespace globals { #define SETTING(name, val) name = val #endif -SETTING(Mesh* mesh, nullptr); ///< The mesh object SETTING(MpiWrapper* mpi, nullptr); ///< The MPI wrapper object -/// Define for reading a variable from the grid -#define GRID_LOAD1(var) mesh->get(var, #var) -#define GRID_LOAD2(var1, var2) \ - { \ - mesh->get(var1, #var1); \ - mesh->get(var2, #var2); \ - } -#define GRID_LOAD3(var1, var2, var3) \ - { \ - mesh->get(var1, #var1); \ - mesh->get(var2, #var2); \ - mesh->get(var3, #var3); \ - } -#define GRID_LOAD4(var1, var2, var3, var4) \ - { \ - mesh->get(var1, #var1); \ - mesh->get(var2, #var2); \ - mesh->get(var3, #var3); \ - mesh->get(var4, #var4); \ - } -#define GRID_LOAD5(var1, var2, var3, var4, var5) \ - { \ - mesh->get(var1, #var1); \ - mesh->get(var2, #var2); \ - mesh->get(var3, #var3); \ - mesh->get(var4, #var4); \ - mesh->get(var5, #var5); \ - } -#define GRID_LOAD6(var1, var2, var3, var4, var5, var6) \ - { \ - mesh->get(var1, #var1); \ - mesh->get(var2, #var2); \ - mesh->get(var3, #var3); \ - mesh->get(var4, #var4); \ - mesh->get(var5, #var5); \ - mesh->get(var6, #var6); \ - } - -/// Read fields from the global mesh -/// The name of the variable will be used as the name -/// in the input. -/// This should accept up to 10 arguments -#define GRID_LOAD(...) \ - { MACRO_FOR_EACH_FN(GRID_LOAD1, __VA_ARGS__) } - /////////////////////////////////////////////////////////////// #undef GLOBAL diff --git a/include/bout/interpolation_xz.hxx b/include/bout/interpolation_xz.hxx index 52dc38f174..ef1117bfad 100644 --- a/include/bout/interpolation_xz.hxx +++ b/include/bout/interpolation_xz.hxx @@ -1,8 +1,7 @@ /************************************************************************** - * Copyright 2010-2020 B.D.Dudson, S.Farley, P. Hill, J.T. Omotani, J.T. Parker, - * M.V.Umansky, X.Q.Xu + * Copyright 2010-2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -48,8 +47,7 @@ protected: public: XZInterpolation(int y_offset = 0, Mesh* localmeshIn = nullptr) - : y_offset(y_offset), - localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn) {} + : y_offset(y_offset), localmesh(localmeshIn) {} XZInterpolation(const BoutMask& mask, int y_offset = 0, Mesh* mesh = nullptr) : XZInterpolation(y_offset, mesh) { region = regionFromMask(mask, localmesh); @@ -273,11 +271,12 @@ public: static constexpr auto option_name = "type"; static constexpr auto default_type = "hermitespline"; - ReturnType create(Options* options = nullptr, Mesh* mesh = nullptr) const { + ReturnType create(Mesh* mesh, Options* options = nullptr) const { return Factory::create(getType(options), mesh); } - ReturnType create(const std::string& type, [[maybe_unused]] Options* options) const { - return Factory::create(type, nullptr); + ReturnType create(const std::string& type, Mesh* mesh, + [[maybe_unused]] Options* options) const { + return Factory::create(type, mesh); } static void ensureRegistered(); diff --git a/include/bout/interpolation_z.hxx b/include/bout/interpolation_z.hxx index 68cf5b0b06..1309009071 100644 --- a/include/bout/interpolation_z.hxx +++ b/include/bout/interpolation_z.hxx @@ -1,7 +1,7 @@ /************************************************************************** - * Copyright 2020 P. Hill, J.T. Omotani, J.T. Parker + * Copyright 2020 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -74,16 +74,16 @@ public: static constexpr auto option_name = "type"; static constexpr auto default_type = "hermitespline"; - ReturnType create(Options* options, int y_offset = 0, Mesh* mesh = nullptr, + ReturnType create(Mesh* mesh, Options* options, int y_offset = 0, Region region_in = {}) const { return Factory::create(options, y_offset, mesh, region_in); } - ReturnType create(int y_offset = 0, Mesh* mesh = nullptr, - Region region_in = {}) const { + ReturnType create(Mesh* mesh, int y_offset = 0, Region region_in = {}) const { return Factory::create(getType(nullptr), y_offset, mesh, region_in); } - ReturnType create(const std::string& type, [[maybe_unused]] Options* options) const { - return Factory::create(type, 0, nullptr, Region{}); + ReturnType create(const std::string& type, Mesh* mesh, + [[maybe_unused]] Options* options) const { + return Factory::create(type, 0, mesh, Region{}); } static void ensureRegistered(); diff --git a/include/bout/invert/laplacexz.hxx b/include/bout/invert/laplacexz.hxx index 11f1c69330..cc249ef8fe 100644 --- a/include/bout/invert/laplacexz.hxx +++ b/include/bout/invert/laplacexz.hxx @@ -7,9 +7,9 @@ * * ************************************************************************** - * Copyright 2015 B.D.Dudson + * Copyright 2015 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -47,12 +47,12 @@ public: static constexpr auto option_name = "type"; static constexpr auto default_type = "cyclic"; - ReturnType create(Mesh* mesh = nullptr, Options* options = nullptr, + ReturnType create(Mesh* mesh, Options* options = nullptr, CELL_LOC loc = CELL_CENTRE) const { return Factory::create(getType(options), mesh, optionsOrDefaultSection(options), loc); } - ReturnType create(const std::string& type, Options* options) const { - return Factory::create(type, nullptr, options, CELL_CENTRE); + ReturnType create(const std::string& type, Mesh* mesh, Options* options) const { + return Factory::create(type, mesh, options, CELL_CENTRE); } static void ensureRegistered(); @@ -65,9 +65,8 @@ using RegisterUnavailableLaplaceXZ = LaplaceXZFactory::RegisterUnavailableInFact class LaplaceXZ { public: - LaplaceXZ(Mesh* m = nullptr, Options* UNUSED(options) = nullptr, - const CELL_LOC loc = CELL_CENTRE) - : localmesh(m == nullptr ? bout::globals::mesh : m), location(loc) {} + LaplaceXZ(Mesh* m, Options* UNUSED(options) = nullptr, const CELL_LOC loc = CELL_CENTRE) + : localmesh(m), location(loc) {} virtual ~LaplaceXZ() = default; virtual void setCoefs(const Field2D& A, const Field2D& B) = 0; diff --git a/include/bout/invert_laplace.hxx b/include/bout/invert_laplace.hxx index 187056d115..1cc30238c1 100644 --- a/include/bout/invert_laplace.hxx +++ b/include/bout/invert_laplace.hxx @@ -133,20 +133,20 @@ constexpr int INVERT_KX_ZERO = 16; */ class LaplaceFactory - : public Factory { + : public Factory { public: static constexpr auto type_name = "Laplacian"; static constexpr auto section_name = "laplace"; static constexpr auto option_name = "type"; static constexpr auto default_type = LAPLACE_CYCLIC; - ReturnType create(Options* options = nullptr, CELL_LOC loc = CELL_CENTRE, - Mesh* mesh = nullptr, Solver* solver = nullptr) { + ReturnType create(Mesh* mesh, Options* options = nullptr, CELL_LOC loc = CELL_CENTRE, + Solver* solver = nullptr) { options = optionsOrDefaultSection(options); - return Factory::create(getType(options), options, loc, mesh, solver); + return Factory::create(getType(options), mesh, options, loc, solver); } - ReturnType create(const std::string& type, Options* options) const { - return Factory::create(type, options, CELL_CENTRE, nullptr, nullptr); + ReturnType create(const std::string& type, Mesh* mesh, Options* options) const { + return Factory::create(type, mesh, options, CELL_CENTRE, nullptr); } }; @@ -169,8 +169,8 @@ class Solver; /// Base class for Laplacian inversion class Laplacian { public: - Laplacian(Options* options = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + Laplacian(Mesh* mesh_in, Options* options = nullptr, const CELL_LOC loc = CELL_CENTRE, + Solver* solver = nullptr); virtual ~Laplacian() = default; /// Set coefficients for inversion. Re-builds matrices if necessary @@ -256,24 +256,21 @@ public: virtual Field2D solve(const Field2D& b, const Field2D& x0); /// Coefficients in tridiagonal inversion - void tridagCoefs(int jx, int jy, int jz, dcomplex& a, dcomplex& b, dcomplex& c, - const Field2D* ccoef = nullptr, const Field2D* d = nullptr, - CELL_LOC loc = CELL_DEFAULT); + /// This function is static so it can be called by BoundaryZeroLaplace2, BoundaryConstLaplace + static void tridagCoefs(Mesh* mesh, int jx, int jy, int jz, dcomplex& a, dcomplex& b, + dcomplex& c, const Field2D* ccoef = nullptr, + const Field2D* d = nullptr, CELL_LOC loc = CELL_DEFAULT); /*! * Create a new Laplacian solver * * @param[in] opt The options section to use. By default "laplace" will be used */ - static std::unique_ptr create(Options* opts = nullptr, + static std::unique_ptr create(Mesh* mesh_in, Options* opts = nullptr, const CELL_LOC location = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr) { - return LaplaceFactory::getInstance().create(opts, location, mesh_in, solver); + return LaplaceFactory::getInstance().create(mesh_in, opts, location, solver); } - static Laplacian* defaultInstance(); ///< Return pointer to global singleton - - static void cleanup(); ///< Frees all memory /// Add any output variables to \p output_options, for example, /// performance information, with optional name for the time @@ -330,14 +327,15 @@ protected: /// and this the last proc in X direction bool isOuterBoundaryFlagSetOnLastX(int flag) const; - void tridagCoefs(int jx, int jy, BoutReal kwave, dcomplex& a, dcomplex& b, dcomplex& c, - const Field2D* ccoef = nullptr, const Field2D* d = nullptr, - CELL_LOC loc = CELL_DEFAULT) { - tridagCoefs(jx, jy, kwave, a, b, c, ccoef, ccoef, d, loc); + static void tridagCoefs(Mesh* mesh, int jx, int jy, BoutReal kwave, dcomplex& a, + dcomplex& b, dcomplex& c, const Field2D* ccoef = nullptr, + const Field2D* d = nullptr, CELL_LOC loc = CELL_DEFAULT) { + tridagCoefs(mesh, jx, jy, kwave, a, b, c, ccoef, ccoef, d, loc); } - void tridagCoefs(int jx, int jy, BoutReal kwave, dcomplex& a, dcomplex& b, dcomplex& c, - const Field2D* c1coef, const Field2D* c2coef, const Field2D* d, - CELL_LOC loc = CELL_DEFAULT); + static void tridagCoefs(Mesh* mesh, int jx, int jy, BoutReal kwave, dcomplex& a, + dcomplex& b, dcomplex& c, const Field2D* c1coef, + const Field2D* c2coef, const Field2D* d, + CELL_LOC loc = CELL_DEFAULT); void tridagMatrix(dcomplex* avec, dcomplex* bvec, dcomplex* cvec, dcomplex* bk, int jy, int kz, BoutReal kwave, const Field2D* a, const Field2D* ccoef, @@ -359,8 +357,6 @@ private: int inner_boundary_flags; ///< Flags to set inner boundary condition int outer_boundary_flags; ///< Flags to set outer boundary condition - /// Singleton instance - static std::unique_ptr instance; /// Name for writing performance infomation; default taken from /// constructing `Options` section std::string performance_name; @@ -386,12 +382,4 @@ protected: void setPerformanceName(std::string name) { performance_name = std::move(name); } }; -//////////////////////////////////////////// -// Legacy interface -// These will be removed at some point - -void laplace_tridag_coefs(int jx, int jy, int jz, dcomplex& a, dcomplex& b, dcomplex& c, - const Field2D* ccoef = nullptr, const Field2D* d = nullptr, - CELL_LOC loc = CELL_DEFAULT); - #endif // BOUT_LAPLACE_H diff --git a/include/bout/invert_parderiv.hxx b/include/bout/invert_parderiv.hxx index e9623e0f9f..e06d920c64 100644 --- a/include/bout/invert_parderiv.hxx +++ b/include/bout/invert_parderiv.hxx @@ -7,9 +7,9 @@ * A + B * Grad2_par2 + C*D2DYDZ + D*D2DZ2 + E*DDY * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -43,19 +43,19 @@ constexpr auto PARDERIVCYCLIC = "cyclic"; class InvertPar; class InvertParFactory - : public Factory { + : public Factory { public: static constexpr auto type_name = "InvertPar"; static constexpr auto section_name = "parderiv"; static constexpr auto option_name = "type"; static constexpr auto default_type = PARDERIVCYCLIC; - ReturnType create(Options* options = nullptr, CELL_LOC location = CELL_CENTRE, - Mesh* mesh = nullptr) const { - return Factory::create(getType(options), options, location, mesh); + ReturnType create(Mesh* mesh, Options* options = nullptr, + CELL_LOC location = CELL_CENTRE) const { + return Factory::create(getType(options), mesh, options, location); } - ReturnType create(const std::string& type, Options* options) const { - return Factory::create(type, options, CELL_CENTRE, nullptr); + ReturnType create(const std::string& type, Mesh* mesh, Options* options) const { + return Factory::create(type, mesh, options, CELL_CENTRE); } static void ensureRegistered(); }; @@ -88,18 +88,16 @@ public: * with pure virtual members, so can't be created directly. * To create an InvertPar object call the create() static function. */ - InvertPar(Options* UNUSED(opt), CELL_LOC location_in, Mesh* mesh_in = nullptr) - : location(location_in), - localmesh(mesh_in == nullptr ? bout::globals::mesh : mesh_in) {} + InvertPar(Mesh* mesh_in, Options* UNUSED(opt), CELL_LOC location_in) + : location(location_in), localmesh(mesh_in) {} virtual ~InvertPar() = default; /*! * Create an instance of InvertPar */ - static std::unique_ptr create(Options* opt_in = nullptr, - CELL_LOC location_in = CELL_CENTRE, - Mesh* mesh_in = nullptr) { - return InvertParFactory::getInstance().create(opt_in, location_in, mesh_in); + static std::unique_ptr create(Mesh* mesh_in, Options* opt_in = nullptr, + CELL_LOC location_in = CELL_CENTRE) { + return InvertParFactory::getInstance().create(mesh_in, opt_in, location_in); } /*! diff --git a/include/bout/invert_pardiv.hxx b/include/bout/invert_pardiv.hxx index 0153cc1987..bef439c98b 100644 --- a/include/bout/invert_pardiv.hxx +++ b/include/bout/invert_pardiv.hxx @@ -7,7 +7,7 @@ * A + Div_par( B * Grad_par ) * ************************************************************************** - * Copyright 2010-2022 BOUT++ contributors + * Copyright 2010-2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -43,19 +43,19 @@ constexpr auto PARDIVCYCLIC = "cyclic"; class InvertParDiv; class InvertParDivFactory - : public Factory { + : public Factory { public: static constexpr auto type_name = "InvertParDiv"; static constexpr auto section_name = "pardiv"; static constexpr auto option_name = "type"; static constexpr auto default_type = PARDIVCYCLIC; - ReturnType create(Options* options = nullptr, CELL_LOC location = CELL_CENTRE, - Mesh* mesh = nullptr) const { - return Factory::create(getType(options), options, location, mesh); + ReturnType create(Mesh* mesh, Options* options = nullptr, + CELL_LOC location = CELL_CENTRE) const { + return Factory::create(getType(options), mesh, options, location); } - ReturnType create(const std::string& type, Options* options) const { - return Factory::create(type, options, CELL_CENTRE, nullptr); + ReturnType create(const std::string& type, Mesh* mesh, Options* options) const { + return Factory::create(type, mesh, options, CELL_CENTRE); } static void ensureRegistered(); }; @@ -88,18 +88,16 @@ public: * with pure virtual members, so can't be created directly. * To create an InvertParDiv object call the create() static function. */ - InvertParDiv(Options* UNUSED(opt), CELL_LOC location_in, Mesh* mesh_in = nullptr) - : location(location_in), - localmesh(mesh_in == nullptr ? bout::globals::mesh : mesh_in) {} + InvertParDiv(Mesh* mesh_in, Options* UNUSED(opt), CELL_LOC location_in) + : location(location_in), localmesh(mesh_in) {} virtual ~InvertParDiv() = default; /*! * Create an instance of InvertParDiv */ - static std::unique_ptr create(Options* opt_in = nullptr, - CELL_LOC location_in = CELL_CENTRE, - Mesh* mesh_in = nullptr) { - return InvertParDivFactory::getInstance().create(opt_in, location_in, mesh_in); + static std::unique_ptr create(Mesh* mesh_in, Options* opt_in = nullptr, + CELL_LOC location_in = CELL_CENTRE) { + return InvertParDivFactory::getInstance().create(mesh_in, opt_in, location_in); } /*! diff --git a/include/bout/invertable_operator.hxx b/include/bout/invertable_operator.hxx index fe139986be..4f1f1f04f9 100644 --- a/include/bout/invertable_operator.hxx +++ b/include/bout/invertable_operator.hxx @@ -136,7 +136,7 @@ public: : operatorFunction(func), preconditionerFunction(func), opt(optIn == nullptr ? Options::getRoot()->getSection("invertableOperator") : optIn), - localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn), lib(opt) { + localmesh(localmeshIn), lib(opt) { AUTO_TRACE(); }; diff --git a/include/bout/mask.hxx b/include/bout/mask.hxx index fd90ae7345..b9398acd60 100644 --- a/include/bout/mask.hxx +++ b/include/bout/mask.hxx @@ -56,7 +56,7 @@ public: explicit BoutMask(const Mesh& mesh, bool value = false) : BoutMask(mesh.LocalNx, mesh.LocalNy, mesh.LocalNz, value) {} explicit BoutMask(const Mesh* mesh = nullptr, bool value = false) - : BoutMask(mesh == nullptr ? *bout::globals::mesh : *mesh, value) {} + : BoutMask(*mesh, value) {} // Assignment from bool BoutMask& operator=(bool value) { diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 02e2a23905..d7b2e5a857 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -80,6 +80,10 @@ public: ReturnType create(const std::string& type, Options* options = nullptr, GridDataSource* source = nullptr) const; ReturnType create(Options* options = nullptr, GridDataSource* source = nullptr) const; + + ReturnType create(const std::string& type, Mesh*, Options* options) const { + return create(type, options); + } }; BOUT_ENUM_CLASS(BoundaryParType, all, xin, xout, fwd, bwd, xin_fwd, xout_fwd, xin_bwd, @@ -879,4 +883,14 @@ Mesh::getRegion(const std::string& region_name) const { return getRegionPerp(region_name); } +/// Define for reading a variable from a mesh +#define GRID_LOAD1(var) mesh->get(var, #var) + +/// Read fields from a mesh +/// The name of the variable will be used as the name +/// in the input. +/// This should accept up to 10 arguments +#define GRID_LOAD(...) \ + { MACRO_FOR_EACH_FN(GRID_LOAD1, __VA_ARGS__) } + #endif // BOUT_MESH_H diff --git a/include/bout/physicsmodel.hxx b/include/bout/physicsmodel.hxx index 9fa25d8b0f..04a02b23d7 100644 --- a/include/bout/physicsmodel.hxx +++ b/include/bout/physicsmodel.hxx @@ -11,9 +11,9 @@ * * Initial version * ************************************************************************** - * Copyright 2013 B.D.Dudson + * Copyright 2013 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -134,13 +134,14 @@ public: typename = std::enable_if_t>> using ModelJacobianFunc = int (Model::*)(BoutReal t); - PhysicsModel(); + /// A PhysicsModel requires at least a Mesh + PhysicsModel(Mesh* mesh_); PhysicsModel(Mesh* mesh_, bool output_enabled_, bool restart_enabled_) : mesh(mesh_), output_enabled(output_enabled_), restart_enabled(restart_enabled_) {} virtual ~PhysicsModel() = default; - Mesh* mesh{nullptr}; + Mesh* mesh; bout::DataFileFacade dump{}; bout::DataFileFacade restart{}; @@ -426,6 +427,8 @@ private: return init_err; \ } \ try { \ + auto mesh = Mesh::create(); \ + mesh->load(); \ auto model = bout::utils::make_unique(); \ auto solver = Solver::create(); \ solver->setModel(model.get()); \ diff --git a/include/bout/rkscheme.hxx b/include/bout/rkscheme.hxx index ba818c04fe..3b5b5cd825 100644 --- a/include/bout/rkscheme.hxx +++ b/include/bout/rkscheme.hxx @@ -53,6 +53,10 @@ public: static constexpr auto section_name = "solver"; static constexpr auto option_name = "scheme"; static constexpr auto default_type = RKSCHEME_RKF45; + + ReturnType create(const std::string& type, Mesh*, Options* options) const { + return Factory::create(type, options); + } }; /// Simpler name for Factory registration helper class diff --git a/include/bout/snb.hxx b/include/bout/snb.hxx index 1644e0059c..59fce0e6e6 100644 --- a/include/bout/snb.hxx +++ b/include/bout/snb.hxx @@ -3,6 +3,8 @@ #include +class Mesh; + namespace bout { /// Calculate heat flux using the Shurtz-Nicolai-Busquet (SNB) model @@ -20,11 +22,11 @@ namespace bout { class HeatFluxSNB { public: /// Construct using the options in the "snb" section. - HeatFluxSNB() : HeatFluxSNB(Options::root()["snb"]) {} + HeatFluxSNB(Mesh* mesh) : HeatFluxSNB(mesh, Options::root()["snb"]) {} /// Construct using options in given section. - explicit HeatFluxSNB(Options& options) { - invertpardiv = std::unique_ptr{InvertParDiv::create()}; + explicit HeatFluxSNB(Mesh* mesh, Options& options) { + invertpardiv = std::unique_ptr{InvertParDiv::create(mesh)}; // Read options. Note that the defaults are initialised already r = options["r"] diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 446acefecb..c177f02fae 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -109,6 +109,10 @@ public: #else SOLVERPVODE; #endif + + ReturnType create(const std::string& type, Mesh*, Options* options) const { + return Factory::create(type, options); + } }; /// Simpler name for Factory registration helper class diff --git a/src/bout++.cxx b/src/bout++.cxx index 7f23cf5f91..75fbfcd21c 100644 --- a/src/bout++.cxx +++ b/src/bout++.cxx @@ -4,7 +4,7 @@ * Adapted from the BOUT code by B.Dudson, University of York, Oct 2007 * ************************************************************************** - * Copyright 2010-2023 BOUT++ contributors + * Copyright 2010-2025 BOUT++ contributors * * Contact Ben Dudson, dudson2@llnl.gov * @@ -42,6 +42,7 @@ const char DEFAULT_DIR[] = "data"; #include "bout/invert/laplacexz.hxx" #include "bout/invert_laplace.hxx" #include "bout/invert_parderiv.hxx" +#include "bout/mesh.hxx" #include "bout/mpi_wrapper.hxx" #include "bout/msg_stack.hxx" #include "bout/openmpwrap.hxx" @@ -203,11 +204,6 @@ int BoutInitialise(int& argc, char**& argv) { bout::globals::mpi = new MpiWrapper(); - // Create the mesh - bout::globals::mesh = Mesh::create(); - // Load from sources. Required for Field initialisation - bout::globals::mesh->load(); - // time_report options are used in BoutFinalise, i.e. after we // check for unused options Options::root()["time_report"].setConditionallyUsed(); @@ -316,10 +312,10 @@ template mesh_options["ny"] = 4; mesh_options["nz"] = 4; - // We might need a global mesh for some types, so best make one + // We might need a mesh for some types, so best make one bout::globals::mpi = new MpiWrapper(); - bout::globals::mesh = Mesh::create(); - bout::globals::mesh->load(); + auto mesh = Mesh::create(); + mesh->load(); // An empty Options that we'll later check for used values Options help_options; @@ -327,7 +323,7 @@ template // Most likely failure is typo in type name, so we definitely want // to print that try { - factory.create(type, &help_options[Factory::section_name]); + factory.create(type, mesh, &help_options[Factory::section_name]); } catch (const BoutException& error) { std::cout << error.what() << std::endl; std::exit(EXIT_FAILURE); @@ -758,15 +754,9 @@ int BoutFinalise(bool write_settings) { output.write("\n"); } - // Delete the mesh - delete bout::globals::mesh; - // Make sure all processes have finished writing before exit bout::globals::mpi->MPI_Barrier(BoutComm::get()); - // Laplacian inversion - Laplacian::cleanup(); - // Delete field memory Array::cleanup(); Array::cleanup(); diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx index c8b9ebb689..bff36d8c8a 100644 --- a/src/field/field2d.cxx +++ b/src/field/field2d.cxx @@ -30,8 +30,6 @@ #include #include -#include // for mesh - #include #include @@ -94,13 +92,8 @@ Field2D::~Field2D() { delete deriv; } Field2D& Field2D::allocate() { if (data.empty()) { - if (!fieldmesh) { - // 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; - } + // Must have a mesh + ASSERT0(fieldmesh); data.reallocate(nx * ny); #if CHECK > 2 invalidateGuards(*this); diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 0d2bc0694e..2442369ef1 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -28,7 +28,6 @@ #include "bout/build_config.hxx" #include -#include #include @@ -110,14 +109,8 @@ Field3D::~Field3D() { delete deriv; } Field3D& Field3D::allocate() { if (data.empty()) { - if (!fieldmesh) { - // 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; - } + // Must have a mesh + ASSERT0(fieldmesh); data.reallocate(nx * ny * nz); #if CHECK > 2 invalidateGuards(*this); diff --git a/src/field/field_data.cxx b/src/field/field_data.cxx index b95c0462d9..bfc5c31b96 100644 --- a/src/field/field_data.cxx +++ b/src/field/field_data.cxx @@ -5,7 +5,6 @@ #include #include #include -#include #include #include @@ -49,11 +48,11 @@ CELL_LOC normaliseLocation(CELL_LOC location, Mesh* mesh) { } // namespace bout FieldData::FieldData(Mesh* localmesh, CELL_LOC location_in) - : fieldmesh(localmesh == nullptr ? bout::globals::mesh : localmesh), + : fieldmesh(localmesh), location(bout::normaliseLocation( 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 + // nullptr if the mesh hasn't been initialized yet if (fieldmesh != nullptr) { // sets fieldCoordinates by getting Coordinates for our location from // fieldmesh @@ -204,15 +203,11 @@ FieldGeneratorPtr FieldData::getBndryGenerator(BndryLoc location) { } Mesh* FieldData::getMesh() const { - if (fieldmesh != nullptr) { - return fieldmesh; - } - // Don't set fieldmesh=mesh here, so that fieldmesh==nullptr until - // allocate() is called in one of the derived classes. fieldmesh==nullptr - // indicates that some initialization that would be done in the - // constructor if fieldmesh was a valid Mesh object still needs to be - // done. - return bout::globals::mesh; + // fieldmesh==nullptr indicates that some initialization that would + // be done in the constructor if fieldmesh was a valid Mesh object + // still needs to be done. + ASSERT0(fieldmesh); + return fieldmesh; } FieldData& FieldData::setLocation(CELL_LOC new_location) { diff --git a/src/field/field_factory.cxx b/src/field/field_factory.cxx index f65f2e7f55..bc6c3967b4 100644 --- a/src/field/field_factory.cxx +++ b/src/field/field_factory.cxx @@ -1,7 +1,7 @@ /************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -19,8 +19,6 @@ * along with BOUT++. If not, see . * **************************************************************************/ -#include - #include #include @@ -86,8 +84,7 @@ class FieldIndirect : public FieldGenerator { // FieldFactory public functions FieldFactory::FieldFactory(Mesh* localmesh, Options* opt) - : fieldmesh(localmesh == nullptr ? bout::globals::mesh : localmesh), - options(opt == nullptr ? Options::getRoot() : opt) { + : fieldmesh(localmesh), options(opt == nullptr ? Options::getRoot() : opt) { // Set options // Note: don't use 'options' here because 'options' is a 'const Options*' diff --git a/src/field/fieldperp.cxx b/src/field/fieldperp.cxx index ca9bdc0397..1e275598a4 100644 --- a/src/field/fieldperp.cxx +++ b/src/field/fieldperp.cxx @@ -24,7 +24,6 @@ **************************************************************************/ #include -#include #include @@ -58,13 +57,8 @@ FieldPerp::FieldPerp(Array data_in, Mesh* localmesh, CELL_LOC location FieldPerp& FieldPerp::allocate() { if (data.empty()) { - if (!fieldmesh) { - // 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; - nz = fieldmesh->LocalNz; - } + // Needs to have a mesh + ASSERT0(fieldmesh); data.reallocate(nx * nz); #if CHECK > 2 invalidateGuards(*this); diff --git a/src/field/globalfield.cxx b/src/field/globalfield.cxx index 6bd2fd4dbf..098ddd918d 100644 --- a/src/field/globalfield.cxx +++ b/src/field/globalfield.cxx @@ -2,6 +2,7 @@ #include #include #include +#include GlobalField::GlobalField(Mesh* m, int proc, int xsize, int ysize, int zsize) : mesh(m), data_on_proc(proc), nx(xsize), ny(ysize), nz(zsize) { diff --git a/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx b/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx index 5ce4e540b7..6ecd5a0c12 100644 --- a/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx +++ b/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx @@ -12,9 +12,9 @@ * * Added DST option * ************************************************************************** - * Copyright 2013 B.D.Dudson + * Copyright 2013 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, benjamin.dudson@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -51,9 +51,10 @@ #include -LaplaceCyclic::LaplaceCyclic(Options* opt, const CELL_LOC loc, Mesh* mesh_in, +LaplaceCyclic::LaplaceCyclic(Mesh* mesh_in, Options* opt, const CELL_LOC loc, Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), Acoef(0.0), C1coef(1.0), C2coef(1.0), Dcoef(1.0) { + : Laplacian(mesh_in, opt, loc), Acoef(0.0, localmesh), C1coef(1.0, localmesh), + C2coef(1.0, localmesh), Dcoef(1.0, localmesh) { Acoef.setLocation(location); C1coef.setLocation(location); diff --git a/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx b/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx index 54b300da12..81bd4cb9a1 100644 --- a/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx +++ b/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx @@ -5,9 +5,9 @@ * the cyclic reduction class to solve tridiagonal systems. * ************************************************************************** - * Copyright 2013 B.D.Dudson + * Copyright 2013 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -59,8 +59,8 @@ RegisterLaplace registerlaplacecycle(LAPLACE_CYCLIC); */ class LaplaceCyclic : public Laplacian { public: - LaplaceCyclic(Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplaceCyclic(Mesh* mesh_in, Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, + Solver* solver = nullptr); ~LaplaceCyclic(); using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx index c50be1db85..6584457120 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx @@ -42,8 +42,8 @@ #include -LaplaceHypre3d::LaplaceHypre3d(Options* opt, const CELL_LOC loc, Mesh* mesh_in, Solver*) - : Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0), Ex(0.0), Ez(0.0), +LaplaceHypre3d::LaplaceHypre3d(Mesh* mesh_in, Options* opt, const CELL_LOC loc, Solver*) + : Laplacian(mesh_in, opt, loc), A(0.0), C1(1.0), C2(1.0), D(1.0), Ex(0.0), Ez(0.0), opts(opt == nullptr ? Options::getRoot()->getSection("laplace") : opt), lowerY(localmesh->iterateBndryLowerY()), upperY(localmesh->iterateBndryUpperY()), indexer(std::make_shared>( diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx index d58ef6f688..2ee01d7877 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx @@ -52,8 +52,8 @@ RegisterLaplace registerlaplacehypre3d(LAPLACE_HYPRE3D); class LaplaceHypre3d : public Laplacian { public: - LaplaceHypre3d(Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplaceHypre3d(Mesh* mesh_in, Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, + Solver* solver = nullptr); ~LaplaceHypre3d() override; void setCoefA(const Field2D& val) override { diff --git a/src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.cxx b/src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.cxx index 5255ea939d..64c232490b 100644 --- a/src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.cxx +++ b/src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.cxx @@ -47,8 +47,8 @@ #include -LaplaceIPT::LaplaceIPT(Options* opt, CELL_LOC loc, Mesh* mesh_in, Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), +LaplaceIPT::LaplaceIPT(Mesh* mesh_in, Options* opt, CELL_LOC loc, Solver* UNUSED(solver)) + : Laplacian(mesh_in, opt, loc), rtol((*opt)["rtol"].doc("Relative tolerance").withDefault(1.e-7)), atol((*opt)["atol"].doc("Absolute tolerance").withDefault(1.e-20)), maxits((*opt)["maxits"].doc("Maximum number of iterations").withDefault(100)), diff --git a/src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.hxx b/src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.hxx index 7111d77e19..0cc7db1a0b 100644 --- a/src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.hxx +++ b/src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.hxx @@ -3,9 +3,9 @@ * and tridiagonal solver. * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -51,8 +51,8 @@ RegisterLaplace registerlaplaceipt(LAPLACE_IPT); class LaplaceIPT : public Laplacian { public: - LaplaceIPT(Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplaceIPT(Mesh* mesh_in, Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, + Solver* solver = nullptr); ~LaplaceIPT() = default; friend class Level; diff --git a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx index 04c8c7e0c9..093eaa2afc 100644 --- a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx +++ b/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx @@ -8,7 +8,7 @@ ************************************************************************** * Copyright 2016 K.S. Kang * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -42,9 +42,10 @@ BoutReal soltime = 0.0, settime = 0.0; -LaplaceMultigrid::LaplaceMultigrid(Options* opt, const CELL_LOC loc, Mesh* mesh_in, +LaplaceMultigrid::LaplaceMultigrid(Mesh* mesh_in, Options* opt, const CELL_LOC loc, Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0) { + : Laplacian(mesh_in, opt, loc), A(0.0, localmesh), C1(1.0, localmesh), + C2(1.0, localmesh), D(1.0, localmesh) { TRACE("LaplaceMultigrid::LaplaceMultigrid(Options *opt)"); @@ -391,7 +392,7 @@ 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); + snprintf(outfile, 256, "test_matF_%d.mat", kMG->rProcI); output << "Out file= " << outfile << endl; outf = fopen(outfile, "w"); int dim = (lxx + 2) * (lzz + 2); @@ -418,7 +419,7 @@ FieldPerp LaplaceMultigrid::solve(const FieldPerp& b_in, const FieldPerp& x0) { FILE* outf; char outfile[256]; - sprintf(outfile, "test_matC%1d_%d.mat", i, kMG->rProcI); + snprintf(outfile, 256, "test_matC%1d_%d.mat", i, kMG->rProcI); output << "Out file= " << outfile << endl; outf = fopen(outfile, "w"); int dim = (kMG->lnx[i - 1] + 2) * (kMG->lnz[i - 1] + 2); diff --git a/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx index c728902815..f2d724ebdd 100644 --- a/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx +++ b/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx @@ -137,8 +137,8 @@ private: class LaplaceMultigrid : public Laplacian { public: - LaplaceMultigrid(Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplaceMultigrid(Mesh* mesh_in, Options* opt = nullptr, + const CELL_LOC loc = CELL_CENTRE, Solver* solver = nullptr); ~LaplaceMultigrid(){}; using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/multigrid/multigrid_solver.cxx b/src/invert/laplace/impls/multigrid/multigrid_solver.cxx index 5ff616aa54..413f875b74 100644 --- a/src/invert/laplace/impls/multigrid/multigrid_solver.cxx +++ b/src/invert/laplace/impls/multigrid/multigrid_solver.cxx @@ -197,7 +197,7 @@ void Multigrid1DP::setMultigridC(int UNUSED(plag)) { for (int i = level; i >= 0; i--) { FILE* outf; char outfile[256]; - sprintf(outfile, "2DP_matC%1d_%d.mat", i, rMG->rProcI); + snprintf(outfile, 256, "2DP_matC%1d_%d.mat", i, rMG->rProcI); output << "Out file= " << outfile << endl; outf = fopen(outfile, "w"); int dim = (rMG->lnx[i] + 2) * (rMG->lnz[i] + 2); @@ -224,7 +224,7 @@ void Multigrid1DP::setMultigridC(int UNUSED(plag)) { for (int i = level; i >= 0; i--) { FILE* outf; char outfile[256]; - sprintf(outfile, "S1D_matC%1d_%d.mat", i, sMG->rProcI); + snprintf(outfile, 256, "S1D_matC%1d_%d.mat", i, sMG->rProcI); output << "Out file= " << outfile << endl; outf = fopen(outfile, "w"); int dim = (sMG->lnx[i] + 2) * (sMG->lnz[i] + 2); @@ -458,7 +458,7 @@ void Multigrid1DP::convertMatrixF2D(int level) { if (pcheck == 3) { FILE* outf; char outfile[256]; - sprintf(outfile, "2DP_CP_%d.mat", rProcI); + snprintf(outfile, 256, "2DP_CP_%d.mat", rProcI); output << "Out file= " << outfile << endl; outf = fopen(outfile, "w"); fprintf(outf, "dim = (%d, %d)\n", ggx, gnz[0]); @@ -478,7 +478,7 @@ void Multigrid1DP::convertMatrixF2D(int level) { if (pcheck == 3) { FILE* outf; char outfile[256]; - sprintf(outfile, "2DP_Conv_%d.mat", rProcI); + snprintf(outfile, 256, "2DP_Conv_%d.mat", rProcI); output << "Out file= " << outfile << endl; outf = fopen(outfile, "w"); fprintf(outf, "dim = (%d, %d)\n", ggx, gnz[0]); @@ -627,7 +627,7 @@ void Multigrid2DPf1D::setMultigridC(int UNUSED(plag)) { for (int i = level; i >= 0; i--) { FILE* outf; char outfile[256]; - sprintf(outfile, "S2D_matC%1d_%d.mat", i, sMG->rProcI); + snprintf(outfile, 256, "S2D_matC%1d_%d.mat", i, sMG->rProcI); outf = fopen(outfile, "w"); output << "Out file= " << outfile << endl; int dim = (sMG->lnx[i] + 2) * (sMG->lnz[i] + 2); diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index e6f68d850d..5dc9e83dc7 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -18,9 +18,9 @@ * ========= * ************************************************************************** - * Copyright 2018 B.D.Dudson, M. Loiten, J. Omotani + * Copyright 2018 - 2025 BOUT++ contributoris * - * Contact: Ben Dudson, benjamin.dudson@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -148,10 +148,11 @@ #include "naulin_laplace.hxx" -LaplaceNaulin::LaplaceNaulin(Options* opt, const CELL_LOC loc, Mesh* mesh_in, +LaplaceNaulin::LaplaceNaulin(Mesh* mesh_in, Options* opt, const CELL_LOC loc, Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), Acoef(0.0), C1coef(1.0), C2coef(0.0), Dcoef(1.0), - delp2solver(nullptr), naulinsolver_mean_its(0.), ncalls(0) { + : Laplacian(mesh_in, opt, loc), Acoef(0.0, localmesh), C1coef(1.0, localmesh), + C2coef(0.0, localmesh), Dcoef(1.0, localmesh), delp2solver(nullptr), + naulinsolver_mean_its(0.), ncalls(0) { ASSERT1(opt != nullptr); // An Options pointer should always be passed in by LaplaceFactory @@ -201,7 +202,7 @@ LaplaceNaulin::LaplaceNaulin(Options* opt, const CELL_LOC loc, Mesh* mesh_in, "if it has been decreased below initial_underrelax_factor.") .withDefault(1.1); ASSERT0(underrelax_recovery >= 1.); - delp2solver = create(opt->getSection("delp2solver"), location, localmesh); + delp2solver = create(localmesh, opt->getSection("delp2solver"), location); std::string delp2type; opt->getSection("delp2solver")->get("type", delp2type, "cyclic"); // Check delp2solver is using an FFT scheme, otherwise it will not exactly diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.hxx b/src/invert/laplace/impls/naulin/naulin_laplace.hxx index 1998a046d7..b3c9241b42 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.hxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.hxx @@ -2,9 +2,9 @@ * Iterative solver to handle non-constant-in-z coefficients * ************************************************************************** - * Copyright 2018 B.D.Dudson, M. Loiten, J. Omotani + * Copyright 2018 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -41,8 +41,8 @@ RegisterLaplace registerlaplacenaulin(LAPLACE_NAULIN); */ class LaplaceNaulin : public Laplacian { public: - LaplaceNaulin(Options* opt = NULL, const CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplaceNaulin(Mesh* mesh_in, Options* opt = NULL, const CELL_LOC loc = CELL_CENTRE, + Solver* solver = nullptr); ~LaplaceNaulin() = default; using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/pcr/pcr.cxx b/src/invert/laplace/impls/pcr/pcr.cxx index 48bbdbac4b..4f8fc956db 100644 --- a/src/invert/laplace/impls/pcr/pcr.cxx +++ b/src/invert/laplace/impls/pcr/pcr.cxx @@ -61,8 +61,8 @@ using namespace std; -LaplacePCR::LaplacePCR(Options* opt, CELL_LOC loc, Mesh* mesh_in, Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), Acoef(0.0, localmesh), C1coef(1.0, localmesh), +LaplacePCR::LaplacePCR(Mesh* mesh_in, Options* opt, CELL_LOC loc, Solver* UNUSED(solver)) + : Laplacian(mesh_in, opt, loc), Acoef(0.0, localmesh), C1coef(1.0, localmesh), C2coef(1.0, localmesh), Dcoef(1.0, localmesh), nmode(maxmode + 1), ncx(localmesh->LocalNx), ny(localmesh->LocalNy), avec(ny, nmode, ncx), bvec(ny, nmode, ncx), cvec(ny, nmode, ncx) { @@ -818,7 +818,6 @@ void LaplacePCR ::pcr_forward_single_row(Matrix& a, Matrix& const int nlevel = log2(nprocs); const int nhprocs = nprocs / 2; int dist_rank = 1; - int dist2_rank = 2; /// Parallel cyclic reduction continues until 2x2 matrix are made between a pair of /// rank, (myrank, myrank+nhprocs). @@ -936,7 +935,6 @@ void LaplacePCR ::pcr_forward_single_row(Matrix& a, Matrix& } dist_rank *= 2; - dist2_rank *= 2; } /// Solving 2x2 matrix. All pair of ranks, myrank and myrank+nhprocs, solves it diff --git a/src/invert/laplace/impls/pcr/pcr.hxx b/src/invert/laplace/impls/pcr/pcr.hxx index ec4637f56c..4d8111be31 100644 --- a/src/invert/laplace/impls/pcr/pcr.hxx +++ b/src/invert/laplace/impls/pcr/pcr.hxx @@ -3,9 +3,9 @@ * and tridiagonal solver. * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -40,8 +40,8 @@ RegisterLaplace registerlaplacepcr(LAPLACE_PCR); class LaplacePCR : public Laplacian { public: - LaplacePCR(Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplacePCR(Mesh* mesh_in, Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, + Solver* solver = nullptr); ~LaplacePCR() = default; using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/pcr_thomas/pcr_thomas.cxx b/src/invert/laplace/impls/pcr_thomas/pcr_thomas.cxx index 61c8f58694..8efec8d52b 100644 --- a/src/invert/laplace/impls/pcr_thomas/pcr_thomas.cxx +++ b/src/invert/laplace/impls/pcr_thomas/pcr_thomas.cxx @@ -58,9 +58,9 @@ using namespace std; -LaplacePCR_THOMAS::LaplacePCR_THOMAS(Options* opt, CELL_LOC loc, Mesh* mesh_in, +LaplacePCR_THOMAS::LaplacePCR_THOMAS(Mesh* mesh_in, Options* opt, CELL_LOC loc, Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), Acoef(0.0, localmesh), C1coef(1.0, localmesh), + : Laplacian(mesh_in, opt, loc), Acoef(0.0, localmesh), C1coef(1.0, localmesh), C2coef(1.0, localmesh), Dcoef(1.0, localmesh), nmode(maxmode + 1), ncx(localmesh->LocalNx), ny(localmesh->LocalNy), avec(ny, nmode, ncx), bvec(ny, nmode, ncx), cvec(ny, nmode, ncx) { diff --git a/src/invert/laplace/impls/pcr_thomas/pcr_thomas.hxx b/src/invert/laplace/impls/pcr_thomas/pcr_thomas.hxx index e12a647789..cebfb1e90c 100644 --- a/src/invert/laplace/impls/pcr_thomas/pcr_thomas.hxx +++ b/src/invert/laplace/impls/pcr_thomas/pcr_thomas.hxx @@ -3,9 +3,9 @@ * and tridiagonal solver. * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -40,8 +40,8 @@ RegisterLaplace registerlaplacepcrthomas(LAPLACE_PCR_THOMAS); class LaplacePCR_THOMAS : public Laplacian { public: - LaplacePCR_THOMAS(Options* opt = nullptr, CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplacePCR_THOMAS(Mesh* mesh_in, Options* opt = nullptr, CELL_LOC loc = CELL_CENTRE, + Solver* solver = nullptr); ~LaplacePCR_THOMAS() = default; using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.cxx b/src/invert/laplace/impls/petsc/petsc_laplace.cxx index 5f417d4faa..afd87f02f0 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.cxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.cxx @@ -61,9 +61,10 @@ static PetscErrorCode laplacePCapply(PC pc, Vec x, Vec y) { PetscFunctionReturn(laplace->precon(x, y)); // NOLINT } -LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, +LaplacePetsc::LaplacePetsc(Mesh* mesh_in, Options* opt, const CELL_LOC loc, Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0), Ex(0.0), Ez(0.0), + : Laplacian(mesh_in, opt, loc), A(0.0, localmesh), C1(1.0, localmesh), + C2(1.0, localmesh), D(1.0, localmesh), Ex(0.0, localmesh), Ez(0.0, localmesh), issetD(false), issetC(false), issetE(false), lib(opt == nullptr ? &(Options::root()["laplace"]) : opt) { A.setLocation(location); @@ -308,7 +309,7 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, rightprec = (*opts)["rightprec"].doc("Right preconditioning?").withDefault(true); // Options for preconditioner are in a subsection - pcsolve = Laplacian::create(opts->getSection("precon")); + pcsolve = Laplacian::create(localmesh, opts->getSection("precon")); } // Ensure that the matrix is constructed first time diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.hxx b/src/invert/laplace/impls/petsc/petsc_laplace.hxx index 1d56abd00b..6a4075022c 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.hxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.hxx @@ -57,8 +57,8 @@ RegisterLaplace registerlaplacepetsc(LAPLACE_PETSC); class LaplacePetsc : public Laplacian { public: - LaplacePetsc(Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplacePetsc(Mesh* mesh_in, Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, + Solver* solver = nullptr); ~LaplacePetsc() { KSPDestroy(&ksp); VecDestroy(&xs); diff --git a/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx b/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx index efca1d70be..466d425105 100644 --- a/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx +++ b/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx @@ -4,9 +4,9 @@ * Using PETSc Solvers * ************************************************************************** - * Copyright 2013 J. Buchanan, J.Omotani + * Copyright 2013 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -48,9 +48,10 @@ static constexpr auto DEFAULT_PC_TYPE = PCHYPRE; static constexpr auto DEFAULT_PC_TYPE = PCGAMG; #endif // PETSC_HAVE_HYPRE -LaplacePetsc3dAmg::LaplacePetsc3dAmg(Options* opt, const CELL_LOC loc, Mesh* mesh_in, +LaplacePetsc3dAmg::LaplacePetsc3dAmg(Mesh* mesh_in, Options* opt, const CELL_LOC loc, Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0), Ex(0.0), Ez(0.0), + : Laplacian(mesh_in, opt, loc), A(0.0, localmesh), C1(1.0, localmesh), + C2(1.0, localmesh), D(1.0, localmesh), Ex(0.0, localmesh), Ez(0.0, localmesh), opts(opt == nullptr ? &(Options::root()["laplace"]) : opt), lower_boundary_flags((*opts)["lower_boundary_flags"].withDefault(0)), upper_boundary_flags((*opts)["upper_boundary_flags"].withDefault(0)), @@ -120,8 +121,10 @@ LaplacePetsc3dAmg::LaplacePetsc3dAmg(Options* opt, const CELL_LOC loc, Mesh* mes // Set up boundary conditions in operator const bool inner_X_neumann = isInnerBoundaryFlagSet(INVERT_AC_GRAD); - const auto inner_X_BC = inner_X_neumann ? -1. / coords->dx / sqrt(coords->g_11) : 0.5; - const auto inner_X_BC_plus = inner_X_neumann ? -inner_X_BC : 0.5; + const auto inner_X_BC = inner_X_neumann ? -1. / coords->dx / sqrt(coords->g_11) + : Coordinates::FieldMetric(0.5, localmesh); + const auto inner_X_BC_plus = + inner_X_neumann ? -inner_X_BC : Coordinates::FieldMetric(0.5, localmesh); BOUT_FOR_SERIAL(i, indexer->getRegionInnerX()) { operator3D(i, i) = inner_X_BC[i]; @@ -129,8 +132,10 @@ LaplacePetsc3dAmg::LaplacePetsc3dAmg(Options* opt, const CELL_LOC loc, Mesh* mes } const bool outer_X_neumann = isOuterBoundaryFlagSet(INVERT_AC_GRAD); - const auto outer_X_BC = outer_X_neumann ? 1. / coords->dx / sqrt(coords->g_11) : 0.5; - const auto outer_X_BC_minus = outer_X_neumann ? -outer_X_BC : 0.5; + const auto outer_X_BC = outer_X_neumann ? 1. / coords->dx / sqrt(coords->g_11) + : Coordinates::FieldMetric(0.5, localmesh); + const auto outer_X_BC_minus = + outer_X_neumann ? -outer_X_BC : Coordinates::FieldMetric(0.5, localmesh); BOUT_FOR_SERIAL(i, indexer->getRegionOuterX()) { operator3D(i, i) = outer_X_BC[i]; @@ -138,8 +143,10 @@ LaplacePetsc3dAmg::LaplacePetsc3dAmg(Options* opt, const CELL_LOC loc, Mesh* mes } const bool lower_Y_neumann = flagSet(lower_boundary_flags, INVERT_AC_GRAD); - const auto lower_Y_BC = lower_Y_neumann ? -1. / coords->dy / sqrt(coords->g_22) : 0.5; - const auto lower_Y_BC_plus = lower_Y_neumann ? -lower_Y_BC : 0.5; + const auto lower_Y_BC = lower_Y_neumann ? -1. / coords->dy / sqrt(coords->g_22) + : Coordinates::FieldMetric(0.5, localmesh); + const auto lower_Y_BC_plus = + lower_Y_neumann ? -lower_Y_BC : Coordinates::FieldMetric(0.5, localmesh); BOUT_FOR_SERIAL(i, indexer->getRegionLowerY()) { operator3D(i, i) = lower_Y_BC[i]; @@ -147,8 +154,10 @@ LaplacePetsc3dAmg::LaplacePetsc3dAmg(Options* opt, const CELL_LOC loc, Mesh* mes } const bool upper_Y_neumann = flagSet(upper_boundary_flags, INVERT_AC_GRAD); - const auto upper_Y_BC = upper_Y_neumann ? 1. / coords->dy / sqrt(coords->g_22) : 0.5; - const auto upper_Y_BC_minus = upper_Y_neumann ? -upper_Y_BC : 0.5; + const auto upper_Y_BC = upper_Y_neumann ? 1. / coords->dy / sqrt(coords->g_22) + : Coordinates::FieldMetric(0.5, localmesh); + const auto upper_Y_BC_minus = + upper_Y_neumann ? -upper_Y_BC : Coordinates::FieldMetric(0.5, localmesh); BOUT_FOR_SERIAL(i, indexer->getRegionUpperY()) { operator3D(i, i) = upper_Y_BC[i]; @@ -168,7 +177,9 @@ void setBC(PetscVector& rhs, const Field3D& b_in, if (flagSet(boundary_flags, INVERT_RHS)) { BOUT_FOR(index, region) { ASSERT1(std::isfinite(b_in[index])); } } else { - const auto& outer_X_BC = (flagSet(boundary_flags, INVERT_SET)) ? x0 : 0.0; + const auto& outer_X_BC = (flagSet(boundary_flags, INVERT_SET)) + ? x0 + : Coordinates::FieldMetric(0.0, b_in.getMesh()); BOUT_FOR_SERIAL(index, region) { rhs(index) = outer_X_BC[index]; } } } @@ -272,9 +283,9 @@ PetscMatrix& LaplacePetsc3dAmg::getMatrix3D() { } void LaplacePetsc3dAmg::updateMatrix3D() { - const Field3D dc_dx = issetC ? DDX(C2) : Field3D(); - const Field3D dc_dy = issetC ? DDY(C2) : Field3D(); - const Field3D dc_dz = issetC ? DDZ(C2) : Field3D(); + const Field3D dc_dx = issetC ? DDX(C2) : Field3D(localmesh); + const Field3D dc_dy = issetC ? DDY(C2) : Field3D(localmesh); + const Field3D dc_dz = issetC ? DDZ(C2) : Field3D(localmesh); const auto dJ_dy = DDY(coords->J / coords->g_22); // Set up the matrix for the internal points on the grid. diff --git a/src/invert/laplace/impls/petsc3damg/petsc3damg.hxx b/src/invert/laplace/impls/petsc3damg/petsc3damg.hxx index 146a915ce3..bf8745d8df 100644 --- a/src/invert/laplace/impls/petsc3damg/petsc3damg.hxx +++ b/src/invert/laplace/impls/petsc3damg/petsc3damg.hxx @@ -59,8 +59,8 @@ RegisterLaplace registerlaplacepetsc3damg(LAPLACE_PETSC3DAMG) class LaplacePetsc3dAmg : public Laplacian { public: - LaplacePetsc3dAmg(Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplacePetsc3dAmg(Mesh* mesh_in, Options* opt = nullptr, + const CELL_LOC loc = CELL_CENTRE, Solver* solver = nullptr); ~LaplacePetsc3dAmg() override; void setCoefA(const Field2D& val) override { diff --git a/src/invert/laplace/impls/serial_band/serial_band.cxx b/src/invert/laplace/impls/serial_band/serial_band.cxx index d7b4ac5d3b..12aa403ba0 100644 --- a/src/invert/laplace/impls/serial_band/serial_band.cxx +++ b/src/invert/laplace/impls/serial_band/serial_band.cxx @@ -3,9 +3,9 @@ * and band-diagonal solver * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -42,9 +42,10 @@ //#define SECONDORDER // Define to use 2nd order differencing -LaplaceSerialBand::LaplaceSerialBand(Options* opt, const CELL_LOC loc, Mesh* mesh_in, +LaplaceSerialBand::LaplaceSerialBand(Mesh* mesh_in, Options* opt, const CELL_LOC loc, Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), Acoef(0.0), Ccoef(1.0), Dcoef(1.0) { + : Laplacian(mesh_in, opt, loc), Acoef(0.0, localmesh), Ccoef(1.0, localmesh), + Dcoef(1.0, localmesh) { Acoef.setLocation(location); Ccoef.setLocation(location); Dcoef.setLocation(location); diff --git a/src/invert/laplace/impls/serial_band/serial_band.hxx b/src/invert/laplace/impls/serial_band/serial_band.hxx index 82e73b82e4..0f885a4fd8 100644 --- a/src/invert/laplace/impls/serial_band/serial_band.hxx +++ b/src/invert/laplace/impls/serial_band/serial_band.hxx @@ -3,9 +3,9 @@ * and band-diagonal solver * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -51,8 +51,8 @@ RegisterLaplace registerlaplaceserialband(LAPLACE_BAND); class LaplaceSerialBand : public Laplacian { public: - LaplaceSerialBand(Options* opt = nullptr, const CELL_LOC = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplaceSerialBand(Mesh* mesh_in, Options* opt = nullptr, const CELL_LOC = CELL_CENTRE, + Solver* solver = nullptr); ~LaplaceSerialBand(){}; using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/serial_tri/serial_tri.cxx b/src/invert/laplace/impls/serial_tri/serial_tri.cxx index 0fb9294d76..e23e1997e5 100644 --- a/src/invert/laplace/impls/serial_tri/serial_tri.cxx +++ b/src/invert/laplace/impls/serial_tri/serial_tri.cxx @@ -38,9 +38,10 @@ #include -LaplaceSerialTri::LaplaceSerialTri(Options* opt, CELL_LOC loc, Mesh* mesh_in, +LaplaceSerialTri::LaplaceSerialTri(Mesh* mesh_in, Options* opt, CELL_LOC loc, Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), A(0.0), C(1.0), D(1.0) { + : Laplacian(mesh_in, opt, loc), A(0.0, localmesh), C(1.0, localmesh), + D(1.0, localmesh) { A.setLocation(location); C.setLocation(location); D.setLocation(location); diff --git a/src/invert/laplace/impls/serial_tri/serial_tri.hxx b/src/invert/laplace/impls/serial_tri/serial_tri.hxx index 5b0419fa27..a1117b1ddf 100644 --- a/src/invert/laplace/impls/serial_tri/serial_tri.hxx +++ b/src/invert/laplace/impls/serial_tri/serial_tri.hxx @@ -3,9 +3,9 @@ * and tridiagonal solver. * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -39,8 +39,8 @@ RegisterLaplace registerlaplaceserialtri(LAPLACE_TRI); class LaplaceSerialTri : public Laplacian { public: - LaplaceSerialTri(Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplaceSerialTri(Mesh* mesh_in, Options* opt = nullptr, + const CELL_LOC loc = CELL_CENTRE, Solver* solver = nullptr); ~LaplaceSerialTri(){}; using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/spt/spt.cxx b/src/invert/laplace/impls/spt/spt.cxx index 2e4c844c94..14f09c2d96 100644 --- a/src/invert/laplace/impls/spt/spt.cxx +++ b/src/invert/laplace/impls/spt/spt.cxx @@ -10,9 +10,9 @@ * * Removed static variables in functions, changing to class members. * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -42,9 +42,10 @@ #include "spt.hxx" -LaplaceSPT::LaplaceSPT(Options* opt, const CELL_LOC loc, Mesh* mesh_in, +LaplaceSPT::LaplaceSPT(Mesh* mesh_in, Options* opt, const CELL_LOC loc, Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), Acoef(0.0), Ccoef(1.0), Dcoef(1.0) { + : Laplacian(mesh_in, opt, loc), Acoef(0.0, localmesh), Ccoef(1.0, localmesh), + Dcoef(1.0, localmesh) { Acoef.setLocation(location); Ccoef.setLocation(location); Dcoef.setLocation(location); diff --git a/src/invert/laplace/impls/spt/spt.hxx b/src/invert/laplace/impls/spt/spt.hxx index a9d5b2583f..2e6d297dfd 100644 --- a/src/invert/laplace/impls/spt/spt.hxx +++ b/src/invert/laplace/impls/spt/spt.hxx @@ -15,9 +15,9 @@ * * Removed static variables in functions, changing to class members. * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -67,8 +67,8 @@ class LaplaceSPT; */ class LaplaceSPT : public Laplacian { public: - LaplaceSPT(Options* opt = nullptr, const CELL_LOC = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); + LaplaceSPT(Mesh* mesh_in, Options* opt = nullptr, const CELL_LOC = CELL_CENTRE, + Solver* solver = nullptr); using Laplacian::setCoefA; void setCoefA(const Field2D& val) override { diff --git a/src/invert/laplace/invert_laplace.cxx b/src/invert/laplace/invert_laplace.cxx index 4032499781..3dbe417317 100644 --- a/src/invert/laplace/invert_laplace.cxx +++ b/src/invert/laplace/invert_laplace.cxx @@ -10,9 +10,9 @@ * Flags control the boundary conditions (see header file) * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -34,7 +34,6 @@ #include #include #include -#include #include #include #include @@ -65,9 +64,9 @@ **********************************************************************************/ /// Laplacian inversion initialisation. Called once at the start to get settings -Laplacian::Laplacian(Options* options, const CELL_LOC loc, Mesh* mesh_in, +Laplacian::Laplacian(Mesh* mesh_in, Options* options, const CELL_LOC loc, Solver* UNUSED(solver)) - : location(loc), localmesh(mesh_in == nullptr ? bout::globals::mesh : mesh_in) { + : location(loc), localmesh(mesh_in) { if (options == nullptr) { // Use the default options @@ -139,17 +138,6 @@ Laplacian::Laplacian(Options* options, const CELL_LOC loc, Mesh* mesh_in, .withDefault(0); } -std::unique_ptr Laplacian::instance = nullptr; - -Laplacian* Laplacian::defaultInstance() { - if (instance == nullptr) { - instance = create(); - } - return instance.get(); -} - -void Laplacian::cleanup() { instance.reset(); } - /********************************************************************************** * Solve routines **********************************************************************************/ @@ -260,31 +248,29 @@ Field2D Laplacian::solve(const Field2D& b, const Field2D& x0) { * MATRIX ELEMENTS **********************************************************************************/ -void Laplacian::tridagCoefs(int jx, int jy, int jz, dcomplex& a, dcomplex& b, dcomplex& c, - const Field2D* ccoef, const Field2D* d, CELL_LOC loc) { - - if (loc == CELL_DEFAULT) { - loc = location; - } - +void Laplacian::tridagCoefs(Mesh* mesh, int jx, int jy, int jz, dcomplex& a, dcomplex& b, + dcomplex& c, const Field2D* ccoef, const Field2D* d, + CELL_LOC loc) { ASSERT1(ccoef == nullptr || ccoef->getLocation() == loc); ASSERT1(d == nullptr || d->getLocation() == loc); - BoutReal kwave = jz * 2.0 * PI / coords->zlength()(jx, jy); // wave number is 1/[rad] + BoutReal kwave = + jz * 2.0 * PI / mesh->getCoordinates()->zlength()(jx, jy); // wave number is 1/[rad] - tridagCoefs(jx, jy, kwave, a, b, c, ccoef, d, loc); + tridagCoefs(mesh, jx, jy, kwave, a, b, c, ccoef, d, loc); } #if BOUT_USE_METRIC_3D -void Laplacian::tridagCoefs(int /* jx */, int /* jy */, BoutReal /* kwave */, - dcomplex& /* a */, dcomplex& /* b */, dcomplex& /* c */, - const Field2D* /* c1coef */, const Field2D* /* c2coef */, - const Field2D* /* d */, CELL_LOC /* loc */) { +void Laplacian::tridagCoefs(Mesh* /* mesh */, int /* jx */, int /* jy */, + BoutReal /* kwave */, dcomplex& /* a */, dcomplex& /* b */, + dcomplex& /* c */, const Field2D* /* c1coef */, + const Field2D* /* c2coef */, const Field2D* /* d */, + CELL_LOC /* loc */) { throw BoutException("Laplacian::tridagCoefs() does not support 3d metrics."); } #else -void Laplacian::tridagCoefs(int jx, int jy, BoutReal kwave, dcomplex& a, dcomplex& b, - dcomplex& c, const Field2D* c1coef, const Field2D* c2coef, - const Field2D* d, CELL_LOC loc) { +void Laplacian::tridagCoefs(Mesh* mesh, int jx, int jy, BoutReal kwave, dcomplex& a, + dcomplex& b, dcomplex& c, const Field2D* c1coef, + const Field2D* c2coef, const Field2D* d, CELL_LOC loc) { /* Function: Laplacian::tridagCoef * Purpose: - Set the matrix components of A in Ax=b, solving * @@ -314,13 +300,7 @@ void Laplacian::tridagCoefs(int jx, int jy, BoutReal kwave, dcomplex& a, dcomple * here) */ - Coordinates* localcoords; - if (loc == CELL_DEFAULT) { - loc = location; - localcoords = coords; - } else { - localcoords = localmesh->getCoordinates(loc); - } + Coordinates* localcoords = mesh->getCoordinates(loc); BoutReal coef1, coef2, coef3, coef4, coef5; @@ -330,11 +310,8 @@ void Laplacian::tridagCoefs(int jx, int jy, BoutReal kwave, dcomplex& a, dcomple coef4 = 0.0; coef5 = 0.0; - // If global flag all_terms are set (true by default) - if (all_terms) { - coef4 = localcoords->G1(jx, jy); // X 1st derivative - coef5 = localcoords->G3(jx, jy); // Z 1st derivative - } + coef4 = localcoords->G1(jx, jy); // X 1st derivative + coef5 = localcoords->G3(jx, jy); // Z 1st derivative if (d != nullptr) { // Multiply Delp2 component by a factor @@ -345,19 +322,17 @@ void Laplacian::tridagCoefs(int jx, int jy, BoutReal kwave, dcomplex& a, dcomple coef5 *= (*d)(jx, jy); } - if (nonuniform) { - // non-uniform mesh correction - if ((jx != 0) && (jx != (localmesh->LocalNx - 1))) { - coef4 -= 0.5 - * ((localcoords->dx(jx + 1, jy) - localcoords->dx(jx - 1, jy)) - / SQ(localcoords->dx(jx, jy))) - * coef1; - } + // non-uniform mesh correction + if ((jx != 0) && (jx != (mesh->LocalNx - 1))) { + coef4 -= 0.5 + * ((localcoords->dx(jx + 1, jy) - localcoords->dx(jx - 1, jy)) + / SQ(localcoords->dx(jx, jy))) + * coef1; } if (c1coef != nullptr) { // First derivative terms - if ((jx > 0) && (jx < (localmesh->LocalNx - 1))) { + if ((jx > 0) && (jx < (mesh->LocalNx - 1))) { BoutReal dc2dx_over_c1 = ((*c2coef)(jx + 1, jy) - (*c2coef)(jx - 1, jy)) / (2. * localcoords->dx(jx, jy) * ((*c1coef)(jx, jy))); coef4 += localcoords->g11(jx, jy) * dc2dx_over_c1; @@ -365,7 +340,7 @@ void Laplacian::tridagCoefs(int jx, int jy, BoutReal kwave, dcomplex& a, dcomple } } - if (localmesh->IncIntShear) { + if (mesh->IncIntShear) { // d2dz2 term coef2 += localcoords->g11(jx, jy) * localcoords->IntShiftTorsion(jx, jy) * localcoords->IntShiftTorsion(jx, jy); @@ -479,7 +454,8 @@ void Laplacian::tridagMatrix(dcomplex* avec, dcomplex* bvec, dcomplex* cvec, dco // The boundaries will be set according to the if-statements below. for (int ix = 0; ix <= ncx; ix++) { // Actually set the metric coefficients - tridagCoefs(xs + ix, jy, kwave, avec[ix], bvec[ix], cvec[ix], c1coef, c2coef, d); + tridagCoefs(localmesh, xs + ix, jy, kwave, avec[ix], bvec[ix], cvec[ix], c1coef, + c2coef, d); if (a != nullptr) { // Add A to bvec (the main diagonal in the matrix) bvec[ix] += (*a)(xs + ix, jy); @@ -801,15 +777,3 @@ bool Laplacian::isInnerBoundaryFlagSetOnFirstX(int flag) const { bool Laplacian::isOuterBoundaryFlagSetOnLastX(int flag) const { return isOuterBoundaryFlagSet(flag) and localmesh->lastX(); } - -/********************************************************************************** - * LEGACY INTERFACE - * - * These functions are depreciated, and will be removed in future - **********************************************************************************/ - -/// Returns the coefficients for a tridiagonal matrix for laplace. Used by Delp2 too -void laplace_tridag_coefs(int jx, int jy, int jz, dcomplex& a, dcomplex& b, dcomplex& c, - const Field2D* ccoef, const Field2D* d, CELL_LOC loc) { - Laplacian::defaultInstance()->tridagCoefs(jx, jy, jz, a, b, c, ccoef, d, loc); -} diff --git a/src/invert/laplacexy/laplacexy.cxx b/src/invert/laplacexy/laplacexy.cxx index bdd22e29bd..d6846a7588 100644 --- a/src/invert/laplacexy/laplacexy.cxx +++ b/src/invert/laplacexy/laplacexy.cxx @@ -7,7 +7,6 @@ #include #include #include -#include #include #include #include @@ -28,8 +27,8 @@ static PetscErrorCode laplacePCapply(PC pc, Vec x, Vec y) { } LaplaceXY::LaplaceXY(Mesh* m, Options* opt, const CELL_LOC loc) - : lib(opt == nullptr ? &(Options::root()["laplacexy"]) : opt), - localmesh(m == nullptr ? bout::globals::mesh : m), location(loc), monitor(*this) { + : lib(opt == nullptr ? &(Options::root()["laplacexy"]) : opt), localmesh(m), + location(loc), indexXY(m), monitor(*this) { Timer timer("invert"); if (opt == nullptr) { @@ -1217,7 +1216,7 @@ const Field2D LaplaceXY::solve(const Field2D& rhs, const Field2D& x0) { ////////////////////////// // Copy data into result - Field2D result; + Field2D result(localmesh); result.allocate(); result.setLocation(location); diff --git a/src/invert/laplacexy2/laplacexy2.cxx b/src/invert/laplacexy2/laplacexy2.cxx index 713fd8e68f..67de77fc45 100644 --- a/src/invert/laplacexy2/laplacexy2.cxx +++ b/src/invert/laplacexy2/laplacexy2.cxx @@ -4,7 +4,6 @@ #include #include -#include #include #include #include @@ -22,9 +21,8 @@ Ind2D index2d(Mesh* mesh, int x, int y) { } // namespace LaplaceXY2::LaplaceXY2(Mesh* m, Options* opt, const CELL_LOC loc) - : localmesh(m == nullptr ? bout::globals::mesh : m), - indexConverter(std::make_shared>( - localmesh, squareStencil(localmesh))), + : localmesh(m), indexConverter(std::make_shared>( + localmesh, squareStencil(localmesh))), matrix(PetscMatrix(indexConverter)), location(loc) { Timer timer("invert"); diff --git a/src/invert/laplacexy2/laplacexy2_hypre.cxx b/src/invert/laplacexy2/laplacexy2_hypre.cxx index a5028169f8..91ff61a7ae 100644 --- a/src/invert/laplacexy2/laplacexy2_hypre.cxx +++ b/src/invert/laplacexy2/laplacexy2_hypre.cxx @@ -4,7 +4,6 @@ #include "bout/assert.hxx" #include "bout/boutcomm.hxx" -#include "bout/globals.hxx" #include "bout/mesh.hxx" #include "bout/output.hxx" #include "bout/sys/timer.hxx" @@ -25,9 +24,8 @@ inline void gpuAssert(cudaError_t code, const char* file, int line, bool abort = #endif LaplaceXY2Hypre::LaplaceXY2Hypre(Mesh* m, Options* opt, const CELL_LOC loc) - : localmesh(m == nullptr ? bout::globals::mesh : m), - indexConverter(std::make_shared>( - localmesh, squareStencil(localmesh))), + : localmesh(m), indexConverter(std::make_shared>( + localmesh, squareStencil(localmesh))), M(indexConverter), x(indexConverter), b(indexConverter), linearSystem(*localmesh, (opt == nullptr) ? Options::root()["laplacexy"] : *opt), location(loc) { diff --git a/src/invert/parderiv/impls/cyclic/cyclic.cxx b/src/invert/parderiv/impls/cyclic/cyclic.cxx index 5a798ca1d5..dfe2dc5998 100644 --- a/src/invert/parderiv/impls/cyclic/cyclic.cxx +++ b/src/invert/parderiv/impls/cyclic/cyclic.cxx @@ -7,16 +7,10 @@ * * Parallel algorithm, using Cyclic Reduction * - * Author: Ben Dudson, University of York, Oct 2011 - * - * Known issues: - * ------------ - * - * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -53,8 +47,9 @@ #include -InvertParCR::InvertParCR(Options* opt, CELL_LOC location, Mesh* mesh_in) - : InvertPar(opt, location, mesh_in), A(1.0), B(0.0), C(0.0), D(0.0), E(0.0) { +InvertParCR::InvertParCR(Mesh* mesh_in, Options* opt, CELL_LOC location) + : InvertPar(mesh_in, opt, location), A(1.0, localmesh), B(0.0, localmesh), + C(0.0, localmesh), D(0.0, localmesh), E(0.0, localmesh), sg(localmesh) { // Number of k equations to solve for each x location nsys = 1 + (localmesh->LocalNz) / 2; diff --git a/src/invert/parderiv/impls/cyclic/cyclic.hxx b/src/invert/parderiv/impls/cyclic/cyclic.hxx index 9e76d125b5..eac7e5a9c9 100644 --- a/src/invert/parderiv/impls/cyclic/cyclic.hxx +++ b/src/invert/parderiv/impls/cyclic/cyclic.hxx @@ -7,8 +7,6 @@ * * Parallel algorithm, using Cyclic Reduction * - * Author: Ben Dudson, University of York, Oct 2011 - * * Known issues: * ------------ * @@ -18,9 +16,9 @@ * staggered mesh requires one-sided derivative as boundary condition. * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -60,8 +58,7 @@ RegisterUnavailableInvertPar registerinvertparcyclic{ class InvertParCR : public InvertPar { public: - explicit InvertParCR(Options* opt, CELL_LOC location = CELL_CENTRE, - Mesh* mesh_in = bout::globals::mesh); + explicit InvertParCR(Mesh* mesh_in, Options* opt, CELL_LOC location = CELL_CENTRE); using InvertPar::solve; const Field3D solve(const Field3D& f) override; @@ -98,7 +95,7 @@ public: } private: - Field2D A{0.0}, B{0.0}, C{0.0}, D{0.0}, E{0.0}; + Field2D A, B, C, D, E; Field2D sg; // Coefficient of DDY contribution to Grad2_par2 int nsys; diff --git a/src/invert/pardiv/impls/cyclic/pardiv_cyclic.cxx b/src/invert/pardiv/impls/cyclic/pardiv_cyclic.cxx index be65d5ab0b..8fe30d52be 100644 --- a/src/invert/pardiv/impls/cyclic/pardiv_cyclic.cxx +++ b/src/invert/pardiv/impls/cyclic/pardiv_cyclic.cxx @@ -14,7 +14,7 @@ * * ************************************************************************** - * Copyright 2010 - 2022 BOUT++ contributors + * Copyright 2010 - 2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -52,8 +52,8 @@ #include -InvertParDivCR::InvertParDivCR(Options* opt, CELL_LOC location, Mesh* mesh_in) - : InvertParDiv(opt, location, mesh_in) { +InvertParDivCR::InvertParDivCR(Mesh* mesh_in, Options* opt, CELL_LOC location) + : InvertParDiv(mesh_in, opt, location), A(1.0, localmesh), B(0.0, localmesh) { // Number of k equations to solve for each x location nsys = 1 + (localmesh->LocalNz) / 2; } diff --git a/src/invert/pardiv/impls/cyclic/pardiv_cyclic.hxx b/src/invert/pardiv/impls/cyclic/pardiv_cyclic.hxx index d66a672cb2..3455382042 100644 --- a/src/invert/pardiv/impls/cyclic/pardiv_cyclic.hxx +++ b/src/invert/pardiv/impls/cyclic/pardiv_cyclic.hxx @@ -16,7 +16,7 @@ * staggered mesh requires one-sided derivative as boundary condition. * ************************************************************************** - * Copyright 2010 - 2022 BOUT++ contributors + * Copyright 2010 - 2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -58,8 +58,7 @@ RegisterUnavailableInvertParDiv registerinvertpardivcyclic{ class InvertParDivCR : public InvertParDiv { public: - explicit InvertParDivCR(Options* opt, CELL_LOC location = CELL_CENTRE, - Mesh* mesh_in = bout::globals::mesh); + explicit InvertParDivCR(Mesh* mesh_in, Options* opt, CELL_LOC location = CELL_CENTRE); using InvertParDiv::solve; Field3D solve(const Field3D& f) override; @@ -78,7 +77,7 @@ public: } private: - Field2D A{1.0}, B{0.0}; + Field2D A, B; int nsys; }; diff --git a/src/invert/pardiv/invert_pardiv.cxx b/src/invert/pardiv/invert_pardiv.cxx index 30e421edd8..ec72377015 100644 --- a/src/invert/pardiv/invert_pardiv.cxx +++ b/src/invert/pardiv/invert_pardiv.cxx @@ -6,7 +6,7 @@ * (A + Div_par( B * Grad_par ) x = r * ************************************************************************** - * Copyright 2010-2022 BOUT++ contributors + * Copyright 2010-2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * diff --git a/src/mesh/boundary_standard.cxx b/src/mesh/boundary_standard.cxx index c12901e53a..4fe05a8d21 100644 --- a/src/mesh/boundary_standard.cxx +++ b/src/mesh/boundary_standard.cxx @@ -2763,7 +2763,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { do { for (int jz = 0; jz <= ncz / 2; jz++) { dcomplex la, lb, lc; - laplace_tridag_coefs(x - bx, y, jz, la, lb, lc); + Laplacian::tridagCoefs(mesh, x - bx, y, jz, la, lb, lc); if (bx > 0) { // Outer boundary swap(la, lc); @@ -2808,17 +2808,19 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // Constant X second derivative int bx = bndry->bx; // Loop over the Y dimension + Mesh* mesh = f.getMesh(); for (bndry->first(); !bndry->isDone(); bndry->nextY()) { int x = bndry->x; int y = bndry->y; // Calculate the Laplacian on the last point dcomplex la, lb, lc; - laplace_tridag_coefs(x - 2 * bx, y, 0, la, lb, lc); + Laplacian::tridagCoefs(mesh, x - 2 * bx, y, 0, la, lb, lc); dcomplex val = la * f(x - bx - 1, y) + lb * f(x - 2 * bx, y) + lc * f(x - 2 * bx + 1, y); // Loop in X towards edge of domain + Mesh* mesh = f.getMesh(); do { - laplace_tridag_coefs(x - bx, y, 0, la, lb, lc); + Laplacian::tridagCoefs(mesh, x - bx, y, 0, la, lb, lc); if (bx < 0) { // Lower X f(x, y) = ((val - lb * f(x - bx, y) + lc * f(x - 2 * bx, y)) / la).real(); } else { // Upper X @@ -2864,7 +2866,7 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) { // Calculate Delp2 on point MXG+1 (and put into c1) for (int jz = 0; jz <= ncz / 2; jz++) { dcomplex la, lb, lc; - laplace_tridag_coefs(x - 2 * bx, y, jz, la, lb, lc); + Laplacian::tridagCoefs(mesh, x - 2 * bx, y, jz, la, lb, lc); if (bx < 0) { // Inner X c1[jz] = la * c0[jz] + lb * c1[jz] + lc * c2[jz]; } else { // Outer X diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 34c524d1e7..3925acdaec 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -375,11 +375,15 @@ Coordinates::Coordinates(Mesh* mesh, FieldMetric dx, FieldMetric dy, FieldMetric FieldMetric g_33, FieldMetric g_12, FieldMetric g_13, FieldMetric g_23, FieldMetric ShiftTorsion, FieldMetric IntShiftTorsion) - : dx(std::move(dx)), dy(std::move(dy)), dz(dz), J(std::move(J)), Bxy(std::move(Bxy)), - g11(std::move(g11)), g22(std::move(g22)), g33(std::move(g33)), g12(std::move(g12)), - g13(std::move(g13)), g23(std::move(g23)), g_11(std::move(g_11)), - g_22(std::move(g_22)), g_33(std::move(g_33)), g_12(std::move(g_12)), - g_13(std::move(g_13)), g_23(std::move(g_23)), ShiftTorsion(std::move(ShiftTorsion)), + : dx(std::move(dx)), dy(std::move(dy)), dz(dz), d1_dx(mesh), d1_dy(mesh), d1_dz(mesh), + J(std::move(J)), Bxy(std::move(Bxy)), g11(std::move(g11)), g22(std::move(g22)), + g33(std::move(g33)), g12(std::move(g12)), g13(std::move(g13)), g23(std::move(g23)), + g_11(std::move(g_11)), g_22(std::move(g_22)), g_33(std::move(g_33)), + g_12(std::move(g_12)), g_13(std::move(g_13)), g_23(std::move(g_23)), G1_11(mesh), + G1_22(mesh), G1_33(mesh), G1_12(mesh), G1_13(mesh), G1_23(mesh), G2_11(mesh), + G2_22(mesh), G2_33(mesh), G2_12(mesh), G2_13(mesh), G2_23(mesh), G3_11(mesh), + G3_22(mesh), G3_33(mesh), G3_12(mesh), G3_13(mesh), G3_23(mesh), G1(mesh), G2(mesh), + G3(mesh), ShiftTorsion(std::move(ShiftTorsion)), IntShiftTorsion(std::move(IntShiftTorsion)), nz(mesh->LocalNz), localmesh(mesh), location(CELL_CENTRE) {} @@ -1424,7 +1428,7 @@ void Coordinates::setParallelTransform(Options* options) { "using 3d metrics. You must provide zShift_ylow in the grid " "file."); } - Field2D zShift_centre; + Field2D zShift_centre(localmesh); if (localmesh->get(zShift_centre, "zShift", 0.0, false)) { // No zShift variable. Try qinty in BOUT grid files if (localmesh->get(zShift_centre, "qinty", 0.0, false)) { @@ -1699,7 +1703,8 @@ Field3D Coordinates::Delp2(const Field3D& f, CELL_LOC outloc, bool useFFT) { // Perform x derivative dcomplex a, b, c; - laplace_tridag_coefs(jx, jy, jz, a, b, c, nullptr, nullptr, outloc); + Laplacian::tridagCoefs(localmesh, jx, jy, jz, a, b, c, nullptr, nullptr, + outloc); delft(jx, jz) = a * ft(jx - 1, jz) + b * ft(jx, jz) + c * ft(jx + 1, jz); } @@ -1762,7 +1767,7 @@ FieldPerp Coordinates::Delp2(const FieldPerp& f, CELL_LOC outloc, bool useFFT) { // Perform x derivative dcomplex a, b, c; - laplace_tridag_coefs(jx, jy, jz, a, b, c); + Laplacian::tridagCoefs(localmesh, jx, jy, jz, a, b, c); delft(jx, jz) = a * ft(jx - 1, jz) + b * ft(jx, jz) + c * ft(jx + 1, jz); } @@ -1835,7 +1840,7 @@ Field2D Coordinates::Laplace_perpXY([[maybe_unused]] const Field2D& A, [[maybe_unused]] const Field2D& f) { TRACE("Coordinates::Laplace_perpXY( Field2D )"); #if not(BOUT_USE_METRIC_3D) - Field2D result; + Field2D result = emptyFrom(f); result.allocate(); for (auto i : result.getRegion(RGN_NOBNDRY)) { result[i] = 0.; @@ -1893,7 +1898,7 @@ Field2D Coordinates::Laplace_perpXY([[maybe_unused]] const Field2D& A, const Coordinates::FieldMetric& Coordinates::invSg() const { if (invSgCache == nullptr) { - auto ptr = std::make_unique(); + auto ptr = std::make_unique(localmesh); (*ptr) = 1.0 / sqrt(g_22); invSgCache = std::move(ptr); } @@ -1913,7 +1918,7 @@ Coordinates::Grad2_par2_DDY_invSg(CELL_LOC outloc, const std::string& method) co invSgCache->applyParallelBoundary("parallel_neumann_o2"); // cache - auto ptr = std::make_unique(); + auto ptr = std::make_unique(localmesh); *ptr = DDY(*invSgCache, outloc, method) * invSg(); Grad2_par2_DDY_invSgCache[method] = std::move(ptr); return *Grad2_par2_DDY_invSgCache[method]; diff --git a/src/mesh/data/gridfromfile.cxx b/src/mesh/data/gridfromfile.cxx index 5625522795..d09aac7951 100644 --- a/src/mesh/data/gridfromfile.cxx +++ b/src/mesh/data/gridfromfile.cxx @@ -208,7 +208,6 @@ bool GridFile::getField(Mesh* m, T& var, const std::string& name, BoutReal def, // Total number of y-boundary cells in grid file, used for check later. // Value depends on if we are double-null or not. - int total_grid_yguards = 2 * grid_yguards; if (m->numberOfXPoints > 1) { ASSERT1(m->numberOfXPoints == 2); // Need to check if we are before or after the target in the middle of the @@ -218,9 +217,6 @@ bool GridFile::getField(Mesh* m, T& var, const std::string& name, BoutReal def, // Note: neither ny_inner nor OffsetY include guard cells ys += 2 * grid_yguards; } - - // Add y-boundary guard cells at upper target - total_grid_yguards += 2 * grid_yguards; } // Index offsets into destination diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index fe5422b4d1..ad02480e0a 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -481,7 +481,7 @@ void communicateFluxes(Field3D& f) { Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc) { - Field3D result = 0.0; + Field3D result{zeroFrom(f)}; ////////////////////////////////////////// // X-Z diffusion diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 574902ea7b..90c0fa97cf 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -33,6 +33,7 @@ #include #include #include +#include #include #include #include diff --git a/src/mesh/interpolation/interpolation_z.cxx b/src/mesh/interpolation/interpolation_z.cxx index 8d39e6baa8..55776557bf 100644 --- a/src/mesh/interpolation/interpolation_z.cxx +++ b/src/mesh/interpolation/interpolation_z.cxx @@ -1,7 +1,7 @@ /************************************************************************** - * Copyright 2020 B.D.Dudson, P. Hill, J. Omotani, J. Parker + * Copyright 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -24,8 +24,7 @@ #include ZInterpolation::ZInterpolation(int y_offset, Mesh* mesh, Region region_in) - : localmesh(mesh == nullptr ? bout::globals::mesh : mesh), region(region_in), - y_offset(y_offset) { + : localmesh(mesh), region(region_in), y_offset(y_offset) { if (region.size() == 0) { // Construct region that skips calculating interpolation in y-boundary regions that // should be filled by boundary conditions diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 659df8600d..0b86d9eb6a 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -63,11 +63,11 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options& auto& interpolation_options = options["xzinterpolation"]; interp = - XZInterpolationFactory::getInstance().create(&interpolation_options, &map_mesh); + XZInterpolationFactory::getInstance().create(&map_mesh, &interpolation_options); interp->setYOffset(offset); interp_corner = - XZInterpolationFactory::getInstance().create(&interpolation_options, &map_mesh); + XZInterpolationFactory::getInstance().create(&map_mesh, &interpolation_options); interp_corner->setYOffset(offset); // Index-space coordinates of forward/backward points diff --git a/src/mesh/parallel/shiftedmetricinterp.cxx b/src/mesh/parallel/shiftedmetricinterp.cxx index 7f3637e79c..89e67e2b15 100644 --- a/src/mesh/parallel/shiftedmetricinterp.cxx +++ b/src/mesh/parallel/shiftedmetricinterp.cxx @@ -53,10 +53,10 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in, // offset, so we don't need to faff about after this for (int y_offset = 0; y_offset < mesh.ystart; ++y_offset) { parallel_slice_interpolators[yup_index + y_offset] = - ZInterpolationFactory::getInstance().create(&interp_options, y_offset + 1, &mesh); + ZInterpolationFactory::getInstance().create(&mesh, &interp_options, y_offset + 1); parallel_slice_interpolators[ydown_index + y_offset] = - ZInterpolationFactory::getInstance().create(&interp_options, -y_offset - 1, - &mesh); + ZInterpolationFactory::getInstance().create(&mesh, &interp_options, + -y_offset - 1); // Find the index positions where the magnetic field line intersects the x-z plane // y_offset points up @@ -87,9 +87,9 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in, // Set up interpolation to/from field-aligned coordinates interp_to_aligned = - ZInterpolationFactory::getInstance().create(&interp_options, 0, &mesh); + ZInterpolationFactory::getInstance().create(&mesh, &interp_options, 0); interp_from_aligned = - ZInterpolationFactory::getInstance().create(&interp_options, 0, &mesh); + ZInterpolationFactory::getInstance().create(&mesh, &interp_options, 0); Field3D zt_prime_to(&mesh), zt_prime_from(&mesh); zt_prime_to.allocate(); @@ -120,7 +120,7 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in, yvalid = 20; } // Create regions for parallel boundary conditions - Field2D dy; + Field2D dy(&mesh); mesh.get(dy, "dy", 1.); auto forward_boundary_xin = std::make_shared( "parallel_forward_xin", BNDRY_PAR_FWD_XIN, +1, &mesh); diff --git a/src/mesh/surfaceiter.cxx b/src/mesh/surfaceiter.cxx index d687ccd68e..5f0ca32789 100644 --- a/src/mesh/surfaceiter.cxx +++ b/src/mesh/surfaceiter.cxx @@ -2,6 +2,7 @@ #include #include +#include #include int SurfaceIter::ySize() { return m->ySize(xpos); } diff --git a/src/physics/gyro_average.cxx b/src/physics/gyro_average.cxx index 10f3e02ac6..63aa7e7e02 100644 --- a/src/physics/gyro_average.cxx +++ b/src/physics/gyro_average.cxx @@ -6,9 +6,9 @@ * * Initial version, simple averaging operator * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -40,14 +40,14 @@ Field3D gyroTaylor0(const Field3D& f, const Field3D& rho) { Field3D gyroPade0(const Field3D& f, BoutReal rho, int inner_boundary_flags, int outer_boundary_flags) { - const Field2D a = 1.0; - const Field2D d = -rho * rho; + const Field2D a{1.0, f.getMesh()}; + const Field2D d{-rho * rho, f.getMesh()}; // Invert, leaving boundaries unchanged Timer timer("invert"); - auto* lap = Laplacian::defaultInstance(); + auto lap = Laplacian::create(f.getMesh()); lap->setCoefA(a); lap->setCoefC(1.0); @@ -60,13 +60,13 @@ Field3D gyroPade0(const Field3D& f, BoutReal rho, int inner_boundary_flags, Field3D gyroPade0(const Field3D& f, const Field2D& rho, int inner_boundary_flags, int outer_boundary_flags) { - const Field2D a = 1.0; - const Field2D d = -rho * rho; + const Field2D a{1.0, f.getMesh()}; + const Field2D d{-rho * rho}; // Invert, leaving boundaries unchanged Timer timer("invert"); - auto* lap = Laplacian::defaultInstance(); + auto lap = Laplacian::create(f.getMesh()); lap->setCoefA(a); lap->setCoefC(1.0); @@ -85,13 +85,13 @@ Field3D gyroPade0(const Field3D& f, const Field3D& rho, int inner_boundary_flags Field3D gyroPade1(const Field3D& f, BoutReal rho, int inner_boundary_flags, int outer_boundary_flags) { - const Field2D a = 1.0; - const Field2D d = -0.5 * rho * rho; + const Field2D a{1.0, f.getMesh()}; + const Field2D d{-0.5 * rho * rho, f.getMesh()}; // Invert, leaving boundaries unchanged Timer timer("invert"); - auto* lap = Laplacian::defaultInstance(); + auto lap = Laplacian::create(f.getMesh()); lap->setCoefA(a); lap->setCoefC(1.0); @@ -104,13 +104,13 @@ Field3D gyroPade1(const Field3D& f, BoutReal rho, int inner_boundary_flags, Field3D gyroPade1(const Field3D& f, const Field2D& rho, int inner_boundary_flags, int outer_boundary_flags) { - const Field2D a = 1.0; - const Field2D d = -0.5 * rho * rho; + const Field2D a{1.0, f.getMesh()}; + const Field2D d{-0.5 * rho * rho}; // Invert, leaving boundaries unchanged Timer timer("invert"); - auto* lap = Laplacian::defaultInstance(); + auto lap = Laplacian::create(f.getMesh()); lap->setCoefA(a); lap->setCoefC(1.0); diff --git a/src/physics/physicsmodel.cxx b/src/physics/physicsmodel.cxx index f4be3bde2a..8aaaaa98e7 100644 --- a/src/physics/physicsmodel.cxx +++ b/src/physics/physicsmodel.cxx @@ -7,10 +7,10 @@ * * Initial version * ************************************************************************** - * Copyright 2013 B.D.Dudson + * Copyright 2013 - 2025 BOUT++ contributors + * + * Contact: Ben Dudson, dudson2@llnl.gov * - * Contact: Ben Dudson, bd512@york.ac.uk - * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -61,10 +61,10 @@ void DataFileFacade::add(Vector3D* value, const std::string& name, bool save_rep } } // namespace bout -PhysicsModel::PhysicsModel() - : mesh(bout::globals::mesh), output_enabled(Options::root()["output"]["enabled"] - .doc("Write output files") - .withDefault(true)), +PhysicsModel::PhysicsModel(Mesh* mesh_) + : mesh(mesh_), output_enabled(Options::root()["output"]["enabled"] + .doc("Write output files") + .withDefault(true)), restart_enabled(Options::root()["restart_files"]["enabled"] .doc("Write restart files") .withDefault(true)) @@ -173,7 +173,7 @@ int PhysicsModel::postInit(bool restarting) { if (restart_enabled) { restartVars(restart_options); solver->outputVars(restart_options, false); - bout::globals::mesh->outputVars(restart_options); + mesh->outputVars(restart_options); restart_options["BOUT_VERSION"].force(bout::version::as_double, "PhysicsModel"); diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx index 23883cc043..514ff8671e 100644 --- a/src/solver/impls/arkode/arkode.cxx +++ b/src/solver/impls/arkode/arkode.cxx @@ -229,8 +229,7 @@ int ArkodeSolver::init() { throw BoutException("ARKodeSetUserData failed\n"); } - if (ARKodeSetLinear(arkode_mem, static_cast(set_linear)) - != ARK_SUCCESS) { + if (ARKodeSetLinear(arkode_mem, static_cast(set_linear)) != ARK_SUCCESS) { throw BoutException("ARKodeSetLinear failed\n"); } diff --git a/src/solver/impls/cvode/cvode.cxx b/src/solver/impls/cvode/cvode.cxx index f0d42d39bc..a661b211db 100644 --- a/src/solver/impls/cvode/cvode.cxx +++ b/src/solver/impls/cvode/cvode.cxx @@ -3,7 +3,7 @@ * * ************************************************************************** - * Copyright 2010-2024 BOUT++ contributors + * Copyright 2010 - 2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 17050ad775..f695ba3718 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -121,9 +121,9 @@ private: Vec x1; ///< Previous solution BoutReal time1{-1.0}; ///< Time of previous solution - SNES snes; ///< SNES context - Mat Jmf; ///< Matrix Free Jacobian - Mat Jfd; ///< Finite Difference Jacobian + SNES snes; ///< SNES context + Mat Jmf; ///< Matrix Free Jacobian + Mat Jfd; ///< Finite Difference Jacobian MatFDColoring fdcoloring{nullptr}; ///< Matrix coloring context ///< Jacobian evaluation @@ -135,10 +135,10 @@ private: std::string pc_hypre_type; ///< Hypre preconditioner type std::string line_search_type; ///< Line search type - bool matrix_free; ///< Use matrix free Jacobian - bool matrix_free_operator; ///< Use matrix free Jacobian in the operator? - int lag_jacobian; ///< Re-use Jacobian - bool use_coloring; ///< Use matrix coloring + bool matrix_free; ///< Use matrix free Jacobian + bool matrix_free_operator; ///< Use matrix free Jacobian in the operator? + int lag_jacobian; ///< Re-use Jacobian + 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/solver/solver.cxx b/src/solver/solver.cxx index a3f6a874f6..0345771404 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -1026,10 +1026,9 @@ std::unique_ptr Solver::create(const SolverType& type, Options* opts) { /// Perform an operation at a given Ind2D (jx,jy) location, moving data between BOUT++ and CVODE void Solver::loop_vars_op(Ind2D i2d, BoutReal* udata, int& p, SOLVER_VAR_OP op, bool bndry) { - // Use global mesh: FIX THIS! - Mesh* mesh = bout::globals::mesh; - int nz = mesh->LocalNz; + // This assumes that all fields have the same nz + int nz = f3d.front().var->getMesh()->LocalNz; switch (op) { case SOLVER_VAR_OP::LOAD_VARS: { diff --git a/src/sys/options/options_adios.cxx b/src/sys/options/options_adios.cxx index 09797647ed..92e36a3a02 100644 --- a/src/sys/options/options_adios.cxx +++ b/src/sys/options/options_adios.cxx @@ -7,7 +7,6 @@ #include "bout/bout.hxx" #include "bout/boutexception.hxx" -#include "bout/globals.hxx" #include "bout/mesh.hxx" #include "bout/sys/timer.hxx" @@ -42,8 +41,6 @@ Options readVariable(adios2::Engine& reader, adios2::IO& io, const std::string& std::vector data; adios2::Variable variable = io.InquireVariable(name); - using bout::globals::mesh; - if (variable.ShapeID() == adios2::ShapeID::GlobalValue) { T value; reader.Get(variable, &value, adios2::Mode::Sync); @@ -245,7 +242,7 @@ bool readAttribute(adios2::IO& io, const std::string& name, const std::string& t return false; } -Options OptionsADIOS::read([[maybe_unused]] bool lazy) { +Options OptionsADIOS::read(Mesh* mesh, [[maybe_unused]] bool lazy) { Timer timer("io"); // Open file diff --git a/src/sys/options/options_io.cxx b/src/sys/options/options_io.cxx index 6717b6b07d..3ed7aa39fe 100644 --- a/src/sys/options/options_io.cxx +++ b/src/sys/options/options_io.cxx @@ -1,6 +1,5 @@ #include "bout/options_io.hxx" #include "bout/bout.hxx" -#include "bout/globals.hxx" #include "bout/mesh.hxx" #include "options_adios.hxx" @@ -46,11 +45,11 @@ OptionsIOFactory::ReturnType OptionsIOFactory::createFile(const std::string& fil return create(getDefaultType(), options); } -void writeDefaultOutputFile(Options& data) { +void writeDefaultOutputFile(Options& data, Mesh* mesh) { // Add BOUT++ version and flags bout::experimental::addBuildFlagsToOptions(data); // Add mesh information - bout::globals::mesh->outputVars(data); + mesh->outputVars(data); // Write to the default output file OptionsIOFactory::getInstance().createOutput()->write(data); }