diff --git a/include/boundary_region.hxx b/include/boundary_region.hxx index cc53f1b54e..09bdc3aa6c 100644 --- a/include/boundary_region.hxx +++ b/include/boundary_region.hxx @@ -8,7 +8,11 @@ class BoundaryRegion; #include class Mesh; -extern Mesh* mesh; +namespace bout { +namespace globals { + extern Mesh* mesh; ///< Global mesh +} // namespace bout +} // namespace globals /// Location of boundary enum BndryLoc {BNDRY_XIN=1, @@ -24,9 +28,9 @@ public: BoundaryRegionBase() = delete; BoundaryRegionBase(std::string name, Mesh *passmesh = nullptr) - : localmesh(passmesh ? passmesh : mesh), label(std::move(name)) {} + : localmesh(passmesh ? passmesh : bout::globals::mesh), label(std::move(name)) {} BoundaryRegionBase(std::string name, BndryLoc loc, Mesh *passmesh = nullptr) - : localmesh(passmesh ? passmesh : mesh), label(std::move(name)), location(loc) {} + : localmesh(passmesh ? passmesh : bout::globals::mesh), label(std::move(name)), location(loc) {} virtual ~BoundaryRegionBase() {} diff --git a/include/bout.hxx b/include/bout.hxx index 921d81bc7f..278a2a7fa2 100644 --- a/include/bout.hxx +++ b/include/bout.hxx @@ -38,6 +38,7 @@ #include "boutcomm.hxx" +#include #include "globals.hxx" #include "field2d.hxx" @@ -65,6 +66,13 @@ const BoutReal BOUT_VERSION = BOUT_VERSION_DOUBLE; ///< Version number +#ifndef BOUT_NO_USING_NAMESPACE_BOUTGLOBALS +// Include using statement by default in user code. +// Macro allows us to include bout.hxx or physicsmodel.hxx without the using +// statement in library code. +using namespace bout::globals; +#endif // BOUT_NO_USING_NAMESPACE_BOUTGLOBALS + // BOUT++ main functions /*! diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index f1bbf8bd6e..8015d97f0f 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -10,6 +10,7 @@ #include "../vector2d.hxx" #include "../utils.hxx" +#include namespace FV { /*! @@ -182,6 +183,10 @@ namespace FV { ASSERT2(f_in.getLocation() == v_in.getLocation()); + Mesh* mesh = f_in.getMesh(); + ASSERT1(mesh == v_in.getMesh()); + ASSERT1(mesh == wave_speed.getMesh()); + CellEdges cellboundary; Field3D f = mesh->toFieldAligned(f_in); @@ -341,6 +346,9 @@ namespace FV { const Field3D Div_f_v(const Field3D &n_in, const Vector3D &v, bool bndry_flux) { ASSERT2(n_in.getLocation() == v.getLocation()); + Mesh* mesh = n_in.getMesh(); + ASSERT1(mesh == v.x.getMesh()); + CellEdges cellboundary; Coordinates *coord = n_in.getCoordinates(); diff --git a/include/bout/invert/laplacexz.hxx b/include/bout/invert/laplacexz.hxx index 5bdb3902be..34200a648e 100644 --- a/include/bout/invert/laplacexz.hxx +++ b/include/bout/invert/laplacexz.hxx @@ -39,8 +39,8 @@ class LaplaceXZ { public: LaplaceXZ(Mesh* m = nullptr, Options* UNUSED(options) = nullptr, - const CELL_LOC loc = CELL_CENTRE) - : localmesh(m==nullptr ? mesh : m), location(loc) {} + const CELL_LOC loc = CELL_CENTRE) + : localmesh(m == nullptr ? bout::globals::mesh : m), location(loc) {} virtual ~LaplaceXZ() {} virtual void setCoefs(const Field2D &A, const Field2D &B) = 0; diff --git a/include/bout/invertable_operator.hxx b/include/bout/invertable_operator.hxx index 6886cb4a13..e5650f70c3 100644 --- a/include/bout/invertable_operator.hxx +++ b/include/bout/invertable_operator.hxx @@ -35,6 +35,7 @@ class InvertableOperator; #ifdef BOUT_HAS_PETSC +#include #include #include #include @@ -133,7 +134,7 @@ public: : operatorFunction(func), preconditionerFunction(func), opt(optIn == nullptr ? optIn : Options::getRoot()->getSection("invertableOperator")), - localmesh(localmeshIn == nullptr ? mesh : localmeshIn), doneSetup(false) { + localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn), doneSetup(false) { AUTO_TRACE(); }; diff --git a/include/bout/physicsmodel.hxx b/include/bout/physicsmodel.hxx index 6e98ecebbf..8f4f3266d1 100644 --- a/include/bout/physicsmodel.hxx +++ b/include/bout/physicsmodel.hxx @@ -43,6 +43,7 @@ class PhysicsModel; #include "solver.hxx" #include "unused.hxx" #include "bout/macro_for_each.hxx" + /*! Base class for physics models */ @@ -310,7 +311,7 @@ private: solver->setModel(model); \ Monitor * bout_monitor = new BoutMonitor(); \ solver->addMonitor(bout_monitor, Solver::BACK); \ - solver->outputVars(dump); \ + solver->outputVars(bout::globals::dump); \ solver->solve(); \ delete model; \ delete solver; \ diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 411011f82e..5d0108f588 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -71,7 +71,9 @@ typedef int (*TimestepMonitorFunc)(Solver *solver, BoutReal simtime, BoutReal la #include "vector2d.hxx" #include "vector3d.hxx" +#define BOUT_NO_USING_NAMESPACE_BOUTGLOBALS #include "physicsmodel.hxx" +#undef BOUT_NO_USING_NAMESPACE_BOUTGLOBALS #include #include @@ -290,7 +292,7 @@ class Solver { /// /// @param[inout] outputfile The file to add variable to /// @param[in] save_repeat If true, add variables with time dimension - void outputVars(Datafile &outputfile, bool save_repeat=true); + virtual void outputVars(Datafile &outputfile, bool save_repeat=true); /*! * Create a Solver object. This uses the "type" option diff --git a/include/datafile.hxx b/include/datafile.hxx index ecfd369312..22cbc2c3be 100644 --- a/include/datafile.hxx +++ b/include/datafile.hxx @@ -15,10 +15,6 @@ class Datafile; #define __DATAFILE_H__ #include "bout_types.hxx" -#include "field2d.hxx" -#include "field3d.hxx" -#include "vector2d.hxx" -#include "vector3d.hxx" #include "options.hxx" #include "bout/macro_for_each.hxx" @@ -26,6 +22,11 @@ class Datafile; #include #include +class Mesh; +class Field2D; +class Field3D; +class Vector2D; +class Vector3D; #include #include @@ -37,7 +38,7 @@ class Datafile; */ class Datafile { public: - Datafile(Options *opt = nullptr); + Datafile(Options *opt = nullptr, Mesh* mesh_in = nullptr); Datafile(Datafile &&other) noexcept; ~Datafile(); // need to delete filename @@ -78,6 +79,7 @@ class Datafile { void setAttribute(const std::string &varname, const std::string &attrname, BoutReal value); private: + Mesh* mesh; bool parallel; // Use parallel formats? bool flush; // Flush after every write? bool guards; // Write guard cells? @@ -137,63 +139,63 @@ class Datafile { }; /// Write this variable once to the grid file -#define SAVE_ONCE1(var) dump.add(var, #var, 0); +#define SAVE_ONCE1(var) bout::globals::dump.add(var, #var, 0); #define SAVE_ONCE2(var1, var2) { \ - dump.add(var1, #var1, 0); \ - dump.add(var2, #var2, 0);} + bout::globals::dump.add(var1, #var1, 0); \ + bout::globals::dump.add(var2, #var2, 0);} #define SAVE_ONCE3(var1, var2, var3) {\ - dump.add(var1, #var1, 0); \ - dump.add(var2, #var2, 0); \ - dump.add(var3, #var3, 0);} + bout::globals::dump.add(var1, #var1, 0); \ + bout::globals::dump.add(var2, #var2, 0); \ + bout::globals::dump.add(var3, #var3, 0);} #define SAVE_ONCE4(var1, var2, var3, var4) { \ - dump.add(var1, #var1, 0); \ - dump.add(var2, #var2, 0); \ - dump.add(var3, #var3, 0); \ - dump.add(var4, #var4, 0);} + bout::globals::dump.add(var1, #var1, 0); \ + bout::globals::dump.add(var2, #var2, 0); \ + bout::globals::dump.add(var3, #var3, 0); \ + bout::globals::dump.add(var4, #var4, 0);} #define SAVE_ONCE5(var1, var2, var3, var4, var5) {\ - dump.add(var1, #var1, 0); \ - dump.add(var2, #var2, 0); \ - dump.add(var3, #var3, 0); \ - dump.add(var4, #var4, 0); \ - dump.add(var5, #var5, 0);} + bout::globals::dump.add(var1, #var1, 0); \ + bout::globals::dump.add(var2, #var2, 0); \ + bout::globals::dump.add(var3, #var3, 0); \ + bout::globals::dump.add(var4, #var4, 0); \ + bout::globals::dump.add(var5, #var5, 0);} #define SAVE_ONCE6(var1, var2, var3, var4, var5, var6) {\ - dump.add(var1, #var1, 0); \ - dump.add(var2, #var2, 0); \ - dump.add(var3, #var3, 0); \ - dump.add(var4, #var4, 0); \ - dump.add(var5, #var5, 0); \ - dump.add(var6, #var6, 0);} + bout::globals::dump.add(var1, #var1, 0); \ + bout::globals::dump.add(var2, #var2, 0); \ + bout::globals::dump.add(var3, #var3, 0); \ + bout::globals::dump.add(var4, #var4, 0); \ + bout::globals::dump.add(var5, #var5, 0); \ + bout::globals::dump.add(var6, #var6, 0);} #define SAVE_ONCE(...) \ { MACRO_FOR_EACH(SAVE_ONCE1, __VA_ARGS__) } /// Write this variable every timestep -#define SAVE_REPEAT1(var) dump.add(var, #var, 1); +#define SAVE_REPEAT1(var) bout::globals::dump.add(var, #var, 1); #define SAVE_REPEAT2(var1, var2) { \ - dump.add(var1, #var1, 1); \ - dump.add(var2, #var2, 1);} + bout::globals::dump.add(var1, #var1, 1); \ + bout::globals::dump.add(var2, #var2, 1);} #define SAVE_REPEAT3(var1, var2, var3) {\ - dump.add(var1, #var1, 1); \ - dump.add(var2, #var2, 1); \ - dump.add(var3, #var3, 1);} + bout::globals::dump.add(var1, #var1, 1); \ + bout::globals::dump.add(var2, #var2, 1); \ + bout::globals::dump.add(var3, #var3, 1);} #define SAVE_REPEAT4(var1, var2, var3, var4) { \ - dump.add(var1, #var1, 1); \ - dump.add(var2, #var2, 1); \ - dump.add(var3, #var3, 1); \ - dump.add(var4, #var4, 1);} + bout::globals::dump.add(var1, #var1, 1); \ + bout::globals::dump.add(var2, #var2, 1); \ + bout::globals::dump.add(var3, #var3, 1); \ + bout::globals::dump.add(var4, #var4, 1);} #define SAVE_REPEAT5(var1, var2, var3, var4, var5) {\ - dump.add(var1, #var1, 1); \ - dump.add(var2, #var2, 1); \ - dump.add(var3, #var3, 1); \ - dump.add(var4, #var4, 1); \ - dump.add(var5, #var5, 1);} + bout::globals::dump.add(var1, #var1, 1); \ + bout::globals::dump.add(var2, #var2, 1); \ + bout::globals::dump.add(var3, #var3, 1); \ + bout::globals::dump.add(var4, #var4, 1); \ + bout::globals::dump.add(var5, #var5, 1);} #define SAVE_REPEAT6(var1, var2, var3, var4, var5, var6) {\ - dump.add(var1, #var1, 1); \ - dump.add(var2, #var2, 1); \ - dump.add(var3, #var3, 1); \ - dump.add(var4, #var4, 1); \ - dump.add(var5, #var5, 1); \ - dump.add(var6, #var6, 1);} + bout::globals::dump.add(var1, #var1, 1); \ + bout::globals::dump.add(var2, #var2, 1); \ + bout::globals::dump.add(var3, #var3, 1); \ + bout::globals::dump.add(var4, #var4, 1); \ + bout::globals::dump.add(var5, #var5, 1); \ + bout::globals::dump.add(var6, #var6, 1);} #define SAVE_REPEAT(...) \ { MACRO_FOR_EACH(SAVE_REPEAT1, __VA_ARGS__) } diff --git a/include/dataformat.hxx b/include/dataformat.hxx index b2ebd5d468..462bf712c0 100644 --- a/include/dataformat.hxx +++ b/include/dataformat.hxx @@ -40,9 +40,12 @@ class DataFormat; #include #include +class Mesh; + // Can't copy, to control access to file class DataFormat { public: + DataFormat(Mesh* mesh_in = nullptr); virtual ~DataFormat() { } // File opening routines virtual bool openr(const char *name) = 0; @@ -195,6 +198,9 @@ class DataFormat { /// ------- /// value A BoutReal attribute of the variable virtual bool getAttribute(const std::string &varname, const std::string &attrname, BoutReal &value) = 0; + + protected: + Mesh* mesh; }; // For backwards compatability. In formatfactory.cxx diff --git a/include/field.hxx b/include/field.hxx index 216c12a2a8..fa0710e0d4 100644 --- a/include/field.hxx +++ b/include/field.hxx @@ -33,6 +33,7 @@ class Field; #include "bout_types.hxx" #include "boutexception.hxx" +#include #include "msg_stack.hxx" #include "stencils.hxx" #include @@ -41,7 +42,6 @@ class Field; class Mesh; class Coordinates; -extern Mesh * mesh; ///< Global mesh #ifdef TRACK #include @@ -94,7 +94,7 @@ class Field { if (fieldmesh){ return fieldmesh; } else { - return mesh; + return bout::globals::mesh; } } diff --git a/include/globals.hxx b/include/globals.hxx index e98fe36345..88bab32731 100644 --- a/include/globals.hxx +++ b/include/globals.hxx @@ -27,10 +27,13 @@ #ifndef __GLOBALS_H__ #define __GLOBALS_H__ -#include "bout/mesh.hxx" #include "datafile.hxx" #include "bout/macro_for_each.hxx" +class Mesh; + +namespace bout { +namespace globals { #ifndef GLOBALORIGIN #define GLOBAL extern #define SETTING(name, val) extern name @@ -83,5 +86,7 @@ GLOBAL Datafile dump; #undef GLOBAL #undef SETTING +} // namespace globals +} // namespace bout #endif // __GLOBALS_H__ diff --git a/include/interpolation.hxx b/include/interpolation.hxx index ad72870878..6c688afb69 100644 --- a/include/interpolation.hxx +++ b/include/interpolation.hxx @@ -243,7 +243,7 @@ protected: public: Interpolation(int y_offset = 0, Mesh* localmeshIn = nullptr) - : localmesh(localmeshIn == nullptr ? mesh : localmeshIn), + : localmesh(localmeshIn == nullptr ? bout::globals::mesh : localmeshIn), skip_mask(*localmesh, false), y_offset(y_offset) {} Interpolation(const BoutMask &mask, int y_offset = 0, Mesh *mesh = nullptr) : Interpolation(y_offset, mesh) { diff --git a/include/invert_laplace.hxx b/include/invert_laplace.hxx index bab11d7ec6..e6a4417208 100644 --- a/include/invert_laplace.hxx +++ b/include/invert_laplace.hxx @@ -105,7 +105,7 @@ const int INVERT_OUT_RHS = 32768; ///< Use input value in RHS at outer boundary /// Base class for Laplacian inversion class Laplacian { public: - Laplacian(Options *options = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh* mesh_in = mesh); + Laplacian(Options *options = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh* mesh_in = nullptr); virtual ~Laplacian() {} /// Set coefficients for inversion. Re-builds matrices if necessary @@ -192,7 +192,7 @@ public: * * @param[in] opt The options section to use. By default "laplace" will be used */ - static Laplacian *create(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); + static Laplacian *create(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = nullptr); static Laplacian* defaultInstance(); ///< Return pointer to global singleton static void cleanup(); ///< Frees all memory diff --git a/include/invert_parderiv.hxx b/include/invert_parderiv.hxx index c46dfe817f..6eff81f6b1 100644 --- a/include/invert_parderiv.hxx +++ b/include/invert_parderiv.hxx @@ -63,8 +63,8 @@ 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), Mesh *mesh_in = mesh) - : localmesh(mesh_in) {} + InvertPar(Options *UNUSED(opt), Mesh *mesh_in = nullptr) + : localmesh(mesh_in==nullptr ? bout::globals::mesh : mesh_in) {} virtual ~InvertPar() {} /*! @@ -72,7 +72,7 @@ public: * * Note: For consistency this should be renamed "create" and take an Options* argument */ - static InvertPar* Create(Mesh *mesh_in = mesh); + static InvertPar* Create(Mesh *mesh_in = nullptr); /*! * Solve the system of equations diff --git a/include/mask.hxx b/include/mask.hxx index cc1539e550..60b82c0fa9 100644 --- a/include/mask.hxx +++ b/include/mask.hxx @@ -61,7 +61,7 @@ public: BoutMask(Mesh& mesh, bool value=false) : BoutMask(mesh.LocalNx, mesh.LocalNy, mesh.LocalNz, value) {} // Default constructor uses global mesh - BoutMask() : BoutMask(*mesh) {} + BoutMask() : BoutMask(*bout::globals::mesh) {} // Assignment from bool BoutMask& operator=(bool value) { diff --git a/src/bout++.cxx b/src/bout++.cxx index 6e4ef0b5cd..53a5e45c39 100644 --- a/src/bout++.cxx +++ b/src/bout++.cxx @@ -44,13 +44,16 @@ const char DEFAULT_LOG[] = "BOUT.log"; #include "mpi.h" #include -#include #include #include #include #include #include +#define BOUT_NO_USING_NAMESPACE_BOUTGLOBALS +#include +#undef BOUT_NO_USING_NAMESPACE_BOUTGLOBALS + #include #include @@ -475,9 +478,9 @@ int BoutInitialise(int &argc, char **&argv) { try { ///////////////////////////////////////////// - mesh = Mesh::create(); ///< Create the mesh - mesh->load(); ///< Load from sources. Required for Field initialisation - mesh->setParallelTransform(); ///< Set the parallel transform from options + bout::globals::mesh = Mesh::create(); ///< Create the mesh + bout::globals::mesh->load(); ///< Load from sources. Required for Field initialisation + bout::globals::mesh->setParallelTransform(); ///< Set the parallel transform from options ///////////////////////////////////////////// /// Get some settings @@ -493,23 +496,23 @@ int BoutInitialise(int &argc, char **&argv) { // Set up the "dump" data output file output << "Setting up output (dump) file\n"; - dump = Datafile(options->getSection("output")); + bout::globals::dump = Datafile(options->getSection("output"), bout::globals::mesh); /// Open a file for the output if(append) { - dump.opena("%s/BOUT.dmp.%s", data_dir.c_str(), dump_ext.c_str()); + bout::globals::dump.opena("%s/BOUT.dmp.%s", data_dir.c_str(), dump_ext.c_str()); }else { - dump.openw("%s/BOUT.dmp.%s", data_dir.c_str(), dump_ext.c_str()); + bout::globals::dump.openw("%s/BOUT.dmp.%s", data_dir.c_str(), dump_ext.c_str()); } /// Add book-keeping variables to the output files - dump.add(const_cast(BOUT_VERSION), "BOUT_VERSION", false); - dump.add(simtime, "t_array", true); // Appends the time of dumps into an array - dump.add(iteration, "iteration", false); + bout::globals::dump.add(const_cast(BOUT_VERSION), "BOUT_VERSION", false); + bout::globals::dump.add(simtime, "t_array", true); // Appends the time of dumps into an array + bout::globals::dump.add(iteration, "iteration", false); //////////////////////////////////////////// - mesh->outputVars(dump); ///< Save mesh configuration into output file + bout::globals::mesh->outputVars(bout::globals::dump); ///< Save mesh configuration into output file }catch(BoutException &e) { output_error.write(_("Error encountered during initialisation: %s\n"), e.what()); @@ -556,10 +559,10 @@ int BoutFinalise() { } // Delete the mesh - delete mesh; + delete bout::globals::mesh; // Close the output file - dump.close(); + bout::globals::dump.close(); // Make sure all processes have finished writing before exit MPI_Barrier(BoutComm::get()); @@ -621,7 +624,7 @@ int BoutMonitor::call(Solver *solver, BoutReal t, int iter, int NOUT) { iteration = iter; /// Write dump file - dump.write(); + bout::globals::dump.write(); /// Collect timing information BoutReal wtime = Timer::resetTime("run"); diff --git a/src/field/field.cxx b/src/field/field.cxx index 3e4fd46987..4d133f722d 100644 --- a/src/field/field.cxx +++ b/src/field/field.cxx @@ -32,10 +32,8 @@ #include #include -Field::Field(Mesh *localmesh) : fieldmesh(localmesh) { - if (fieldmesh == nullptr) { - fieldmesh = mesh; - } +Field::Field(Mesh *localmesh) + : fieldmesh(localmesh==nullptr ? bout::globals::mesh : localmesh) { // Note we would like to do `fieldCoordinates = getCoordinates();` here but can't // currently as this would lead to circular/recursive behaviour (getCoordinates would diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx index 7da30dacf4..5fcdef90df 100644 --- a/src/field/field2d.cxx +++ b/src/field/field2d.cxx @@ -39,6 +39,7 @@ #include #include +#include #include #include @@ -89,7 +90,7 @@ void Field2D::allocate() { if(data.empty()) { if(!fieldmesh) { /// If no mesh, use the global - fieldmesh = mesh; + fieldmesh = bout::globals::mesh; nx = fieldmesh->LocalNx; ny = fieldmesh->LocalNy; } @@ -549,8 +550,6 @@ void checkData(const Field2D &f, REGION region) { #if CHECK > 2 void invalidateGuards(Field2D &var) { - Mesh *localmesh = var.getMesh(); - BOUT_FOR(i, var.getRegion("RGN_GUARDS")) { var[i] = BoutNaN; } } #endif diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 302419bbec..ead57190c2 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -123,7 +123,7 @@ void Field3D::allocate() { if(data.empty()) { if(!fieldmesh) { /// If no mesh, use the global - fieldmesh = mesh; + fieldmesh = bout::globals::mesh; nx = fieldmesh->LocalNx; ny = fieldmesh->LocalNy; nz = fieldmesh->LocalNz; @@ -1082,8 +1082,6 @@ Field2D DC(const Field3D &f, REGION rgn) { #if CHECK > 2 void invalidateGuards(Field3D &var) { - Mesh *localmesh = var.getMesh(); - BOUT_FOR(i, var.getRegion("RGN_GUARDS")) { var[i] = BoutNaN; } } #endif diff --git a/src/field/field_data.cxx b/src/field/field_data.cxx index 4251ace1b4..8282f9cda4 100644 --- a/src/field/field_data.cxx +++ b/src/field/field_data.cxx @@ -1,4 +1,5 @@ +#include #include #include #include @@ -20,7 +21,7 @@ void FieldData::setBoundary(const std::string &name) { output_info << "Setting boundary for variable " << name << endl; /// Loop over the mesh boundary regions - for(const auto& reg : mesh->getBoundaries()) { + for(const auto& reg : bout::globals::mesh->getBoundaries()) { BoundaryOp* op = static_cast(bfact->createFromOptions(name, reg)); if (op != nullptr) bndry_op.push_back(op); @@ -28,9 +29,9 @@ void FieldData::setBoundary(const std::string &name) { } /// Get the mesh boundary regions - std::vector par_reg = mesh->getBoundariesPar(); + std::vector par_reg = bout::globals::mesh->getBoundariesPar(); /// Loop over the mesh parallel boundary regions - for(const auto& reg : mesh->getBoundariesPar()) { + for(const auto& reg : bout::globals::mesh->getBoundariesPar()) { BoundaryOpPar* op = static_cast(bfact->createFromOptions(name, reg)); if (op != nullptr) bndry_op_par.push_back(op); @@ -43,7 +44,7 @@ void FieldData::setBoundary(const std::string &name) { void FieldData::setBoundary(const std::string &UNUSED(region), BoundaryOp *op) { /// Get the mesh boundary regions - std::vector reg = mesh->getBoundaries(); + std::vector reg = bout::globals::mesh->getBoundaries(); /// Find the region @@ -72,7 +73,7 @@ void FieldData::addBndryFunction(FuncPtr userfunc, BndryLoc location){ void FieldData::addBndryGenerator(FieldGeneratorPtr gen, BndryLoc location) { if(location == BNDRY_ALL){ - for(const auto& reg : mesh->getBoundaries()) { + for(const auto& reg : bout::globals::mesh->getBoundaries()) { bndry_generator[reg->location] = gen; } } else { diff --git a/src/field/fieldperp.cxx b/src/field/fieldperp.cxx index 07686b90f3..8ac0a2906d 100644 --- a/src/field/fieldperp.cxx +++ b/src/field/fieldperp.cxx @@ -28,6 +28,7 @@ #include +#include #include #include #include @@ -50,7 +51,7 @@ void FieldPerp::allocate() { if (data.empty()) { if (!fieldmesh) { /// If no mesh, use the global - fieldmesh = mesh; + fieldmesh = bout::globals::mesh; nx = fieldmesh->LocalNx; nz = fieldmesh->LocalNz; } @@ -538,8 +539,6 @@ void checkData(const FieldPerp &f, REGION region) { #if CHECK > 2 void invalidateGuards(FieldPerp &var) { - Mesh *localmesh = var.getMesh(); - BOUT_FOR(i, var.getRegion("RGN_GUARDS")) { var[i] = BoutNaN; } } #endif diff --git a/src/field/initialprofiles.cxx b/src/field/initialprofiles.cxx index 16b0fbc156..5fa35f98bb 100644 --- a/src/field/initialprofiles.cxx +++ b/src/field/initialprofiles.cxx @@ -36,6 +36,9 @@ * **************************************************************************/ +#include +#include +#include #include #include #include diff --git a/src/field/vecops.cxx b/src/field/vecops.cxx index 306317d2a9..e61a711e64 100644 --- a/src/field/vecops.cxx +++ b/src/field/vecops.cxx @@ -24,6 +24,7 @@ * **************************************************************************/ +#include #include #include @@ -31,6 +32,7 @@ #include #include #include +#include /************************************************************************** * Gradient operators diff --git a/src/field/vector3d.cxx b/src/field/vector3d.cxx index 1959ce3bbe..80081d1931 100644 --- a/src/field/vector3d.cxx +++ b/src/field/vector3d.cxx @@ -472,7 +472,9 @@ const Vector3D Vector3D::operator/(const Field3D &rhs) const { ////////////////// DOT PRODUCT /////////////////// const Field3D Vector3D::operator*(const Vector3D &rhs) const { - Field3D result(x.getMesh()); + Mesh* mesh = x.getMesh(); + + Field3D result(mesh); ASSERT2(location == rhs.getLocation()) if(rhs.covariant ^ covariant) { diff --git a/src/fileio/datafile.cxx b/src/fileio/datafile.cxx index ba6b259395..80310d20d3 100644 --- a/src/fileio/datafile.cxx +++ b/src/fileio/datafile.cxx @@ -36,9 +36,11 @@ //#include "mpi.h" // For MPI_Wtime() #include +#include #include #include #include +#include #include #include #include @@ -46,10 +48,12 @@ #include #include "formatfactory.hxx" -Datafile::Datafile(Options *opt) : parallel(false), flush(true), guards(true), - floats(false), openclose(true), enabled(true), shiftOutput(false), - shiftInput(false), flushFrequencyCounter(0), flushFrequency(1), - file(nullptr), writable(false), appending(false), first_time(true) +Datafile::Datafile(Options *opt, Mesh* mesh_in) + : mesh(mesh_in==nullptr ? bout::globals::mesh : mesh_in), + parallel(false), flush(true), guards(true), floats(false), openclose(true), + enabled(true), shiftOutput(false), shiftInput(false), + flushFrequencyCounter(0), flushFrequency(1), file(nullptr), writable(false), + appending(false), first_time(true) { filenamelen=FILENAMELEN; filename=new char[filenamelen]; @@ -91,7 +95,7 @@ Datafile::Datafile(Datafile &&other) noexcept } Datafile::Datafile(const Datafile &other) : - parallel(other.parallel), flush(other.flush), guards(other.guards), + mesh(other.mesh), parallel(other.parallel), flush(other.flush), guards(other.guards), floats(other.floats), openclose(other.openclose), Lx(other.Lx), Ly(other.Ly), Lz(other.Lz), enabled(other.enabled), shiftOutput(other.shiftOutput), shiftInput(other.shiftInput), flushFrequencyCounter(other.flushFrequencyCounter), flushFrequency(other.flushFrequency), file(nullptr), writable(other.writable), appending(other.appending), first_time(other.first_time), @@ -106,6 +110,7 @@ Datafile::Datafile(const Datafile &other) : } Datafile& Datafile::operator=(Datafile &&rhs) noexcept { + mesh = rhs.mesh; parallel = rhs.parallel; flush = rhs.flush; guards = rhs.guards; @@ -189,7 +194,7 @@ bool Datafile::openw(const char *format, ...) { bout_vsnprintf(filename, filenamelen, format); // Get the data format - file = FormatFactory::getInstance()->createDataFormat(filename, parallel); + file = FormatFactory::getInstance()->createDataFormat(filename, parallel, mesh); if(!file) throw BoutException("Datafile::open: Factory failed to create a DataFormat!"); diff --git a/src/fileio/dataformat.cxx b/src/fileio/dataformat.cxx index 2fa1d9dd0d..cf09ea170b 100644 --- a/src/fileio/dataformat.cxx +++ b/src/fileio/dataformat.cxx @@ -1,8 +1,12 @@ +#include #include #include #include +DataFormat::DataFormat(Mesh* mesh_in) + : mesh(mesh_in==nullptr ? bout::globals::mesh : mesh_in) {} + bool DataFormat::openr(const std::string &name, int mype) { // Split into base name and extension size_t pos = name.find_last_of('.'); diff --git a/src/fileio/formatfactory.cxx b/src/fileio/formatfactory.cxx index 244b4da6e8..d2ee7fd8f8 100644 --- a/src/fileio/formatfactory.cxx +++ b/src/fileio/formatfactory.cxx @@ -12,6 +12,7 @@ #include #include +#include #include FormatFactory *FormatFactory::instance = nullptr; @@ -25,27 +26,29 @@ FormatFactory* FormatFactory::getInstance() { } // Work out which data format to use for given filename -std::unique_ptr FormatFactory::createDataFormat(const char *filename, bool parallel) { +std::unique_ptr FormatFactory::createDataFormat(const char *filename, + bool parallel, + Mesh* mesh_in) { if ((filename == nullptr) || (strcasecmp(filename, "default") == 0)) { // Return default file format if (parallel) { #ifdef PNCDF - return bout::utils::make_unique(); + return bout::utils::make_unique(mesh_in); #else } #ifdef NCDF4 - return bout::utils::make_unique(); + return bout::utils::make_unique(mesh_in); #else #ifdef NCDF - return bout::utils::make_unique(); + return bout::utils::make_unique(mesh_in); #else #ifdef HDF5 - return bout::utils::make_unique(); + return bout::utils::make_unique(mesh_in); #else #error No file format available; aborting. diff --git a/src/fileio/formatfactory.hxx b/src/fileio/formatfactory.hxx index c6830bb771..361773871c 100644 --- a/src/fileio/formatfactory.hxx +++ b/src/fileio/formatfactory.hxx @@ -14,7 +14,8 @@ public: static FormatFactory* getInstance(); std::unique_ptr createDataFormat(const char *filename = nullptr, - bool parallel = true); + bool parallel = true, + Mesh* mesh_in = nullptr); private: static FormatFactory* instance; ///< The only instance of this class (Singleton) diff --git a/src/fileio/impls/hdf5/h5_format.cxx b/src/fileio/impls/hdf5/h5_format.cxx index f6260982de..b938a52202 100644 --- a/src/fileio/impls/hdf5/h5_format.cxx +++ b/src/fileio/impls/hdf5/h5_format.cxx @@ -30,11 +30,12 @@ #include #include +#include #include #include #include -H5Format::H5Format(bool parallel_in) { +H5Format::H5Format(bool parallel_in, Mesh* mesh_in) : DataFormat(mesh_in) { parallel = parallel_in; x0 = y0 = z0 = t0 = 0; lowPrecision = false; @@ -69,7 +70,8 @@ H5Format::H5Format(bool parallel_in) { throw BoutException("Failed to set error stack to not print errors"); } -H5Format::H5Format(const char *name, bool parallel_in) { +H5Format::H5Format(const char *name, bool parallel_in, Mesh* mesh_in) + : DataFormat(mesh_in) { parallel = parallel_in; x0 = y0 = z0 = t0 = 0; lowPrecision = false; diff --git a/src/fileio/impls/hdf5/h5_format.hxx b/src/fileio/impls/hdf5/h5_format.hxx index 0853f695b6..b319aac461 100644 --- a/src/fileio/impls/hdf5/h5_format.hxx +++ b/src/fileio/impls/hdf5/h5_format.hxx @@ -54,9 +54,10 @@ class H5Format; class H5Format : public DataFormat { public: - H5Format(bool parallel_in = false); - H5Format(const char *name, bool parallel_in = false); - H5Format(const std::string &name, bool parallel_in = false) : H5Format(name.c_str(), parallel_in) {} + H5Format(bool parallel_in = false, Mesh* mesh_in = nullptr); + H5Format(const char *name, bool parallel_in = false, Mesh* mesh_in = nullptr); + H5Format(const std::string &name, bool parallel_in = false, Mesh* mesh_in = nullptr) + : H5Format(name.c_str(), parallel_in, mesh_in) {} ~H5Format(); using DataFormat::openr; diff --git a/src/fileio/impls/netcdf/nc_format.cxx b/src/fileio/impls/netcdf/nc_format.cxx index b648c638eb..7d9e2aa2a4 100644 --- a/src/fileio/impls/netcdf/nc_format.cxx +++ b/src/fileio/impls/netcdf/nc_format.cxx @@ -28,6 +28,7 @@ #include #include +#include #include #include @@ -37,7 +38,7 @@ using std::vector; // Define this to see loads of info messages //#define NCDF_VERBOSE -NcFormat::NcFormat() { +NcFormat::NcFormat(Mesh* mesh_in) : DataFormat(mesh_in) { dataFile = nullptr; x0 = y0 = z0 = t0 = 0; recDimList = new const NcDim*[4]; @@ -50,7 +51,7 @@ NcFormat::NcFormat() { fname = nullptr; } -NcFormat::NcFormat(const char *name) { +NcFormat::NcFormat(const char *name, Mesh* mesh_in) : DataFormat(mesh_in) { dataFile = nullptr; x0 = y0 = z0 = t0 = 0; recDimList = new const NcDim*[4]; diff --git a/src/fileio/impls/netcdf/nc_format.hxx b/src/fileio/impls/netcdf/nc_format.hxx index 4b65329e21..e1ebfaa964 100644 --- a/src/fileio/impls/netcdf/nc_format.hxx +++ b/src/fileio/impls/netcdf/nc_format.hxx @@ -55,9 +55,10 @@ class NcFormat; class NcFormat : public DataFormat { public: - NcFormat(); - NcFormat(const char *name); - NcFormat(const std::string &name) : NcFormat(name.c_str()) {} + NcFormat(Mesh* mesh_in = nullptr); + NcFormat(const char *name, Mesh* mesh_in = nullptr); + NcFormat(const std::string &name, Mesh* mesh_in = nullptr) + : NcFormat(name.c_str(), mesh_in) {} ~NcFormat(); using DataFormat::openr; diff --git a/src/fileio/impls/netcdf4/ncxx4.cxx b/src/fileio/impls/netcdf4/ncxx4.cxx index c6c084c7ef..a1c5c0e6e3 100644 --- a/src/fileio/impls/netcdf4/ncxx4.cxx +++ b/src/fileio/impls/netcdf4/ncxx4.cxx @@ -24,6 +24,7 @@ #ifdef NCDF4 +#include #include #include #include @@ -38,7 +39,7 @@ using namespace netCDF; // Define this to see loads of info messages //#define NCDF_VERBOSE -Ncxx4::Ncxx4() { +Ncxx4::Ncxx4(Mesh* mesh_in) : DataFormat(mesh_in) { dataFile = nullptr; x0 = y0 = z0 = t0 = 0; recDimList = new const NcDim*[4]; @@ -51,7 +52,7 @@ Ncxx4::Ncxx4() { fname = nullptr; } -Ncxx4::Ncxx4(const char *name) { +Ncxx4::Ncxx4(const char *name, Mesh* mesh_in) : DataFormat(mesh_in) { dataFile = nullptr; x0 = y0 = z0 = t0 = 0; recDimList = new const NcDim*[4]; diff --git a/src/fileio/impls/netcdf4/ncxx4.hxx b/src/fileio/impls/netcdf4/ncxx4.hxx index c0a9bdc6df..497b32df0b 100644 --- a/src/fileio/impls/netcdf4/ncxx4.hxx +++ b/src/fileio/impls/netcdf4/ncxx4.hxx @@ -56,9 +56,10 @@ class Ncxx4; class Ncxx4 : public DataFormat { public: - Ncxx4(); - Ncxx4(const char *name); - Ncxx4(const std::string &name) : Ncxx4(name.c_str()) {} + Ncxx4(Mesh* mesh_in = nullptr); + Ncxx4(const char *name, Mesh* mesh_in = nullptr); + Ncxx4(const std::string &name, Mesh* mesh_in = nullptr) + : Ncxx4(name.c_str(), mesh_in) {} ~Ncxx4(); using DataFormat::openr; diff --git a/src/fileio/impls/pnetcdf/pnetcdf.cxx b/src/fileio/impls/pnetcdf/pnetcdf.cxx index 2ba919bd66..4234de61b4 100644 --- a/src/fileio/impls/pnetcdf/pnetcdf.cxx +++ b/src/fileio/impls/pnetcdf/pnetcdf.cxx @@ -41,7 +41,7 @@ // Define this to see loads of info messages //#define NCDF_VERBOSE -PncFormat::PncFormat() { +PncFormat::PncFormat(Mesh* mesh_in) : DataFormat(mesh_in) { x0 = y0 = z0 = t0 = 0; lowPrecision = false; dimList = recDimList+1; @@ -51,7 +51,7 @@ PncFormat::PncFormat() { fname = nullptr; } -PncFormat::PncFormat(const char *name) { +PncFormat::PncFormat(const char *name, Mesh* mesh_in) : DataFormat(mesh_in) { x0 = y0 = z0 = t0 = 0; lowPrecision = false; dimList = recDimList+1; diff --git a/src/fileio/impls/pnetcdf/pnetcdf.hxx b/src/fileio/impls/pnetcdf/pnetcdf.hxx index 19f23d2917..dd355f22df 100644 --- a/src/fileio/impls/pnetcdf/pnetcdf.hxx +++ b/src/fileio/impls/pnetcdf/pnetcdf.hxx @@ -52,9 +52,10 @@ class PncFormat; class PncFormat : public DataFormat { public: - PncFormat(); - PncFormat(const char *name); - PncFormat(const std::string &name) : PncFormat(name.c_str()) {} + PncFormat(Mesh* mesh_in = nullptr); + PncFormat(const char *name, Mesh* mesh_in = nullptr); + PncFormat(const std::string &name, Mesh* mesh_in = nullptr) + : PncFormat(name.c_str(), mesh_in) {} ~PncFormat(); bool openr(const char *name) override; diff --git a/src/invert/lapack_routines.cxx b/src/invert/lapack_routines.cxx index 3efbcf1bd8..ce000e7d0b 100644 --- a/src/invert/lapack_routines.cxx +++ b/src/invert/lapack_routines.cxx @@ -38,6 +38,7 @@ #include #include #include +#include #ifdef LAPACK diff --git a/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx b/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx index af950b7cbe..a0d2d523ae 100644 --- a/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx +++ b/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx @@ -35,6 +35,7 @@ #include #include +#include #include #include #include diff --git a/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx b/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx index b007542569..d4419dd5f2 100644 --- a/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx +++ b/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx @@ -44,7 +44,7 @@ class LaplaceCyclic; */ class LaplaceCyclic : public Laplacian { public: - LaplaceCyclic(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); + LaplaceCyclic(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = nullptr); ~LaplaceCyclic(); using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx index 1ff40e9b0b..86b912154f 100644 --- a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx +++ b/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx @@ -28,6 +28,7 @@ **************************************************************************/ #include "multigrid_laplace.hxx" +#include #include #include diff --git a/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx index 0067e1bb51..e0e05328eb 100644 --- a/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx +++ b/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx @@ -132,7 +132,8 @@ private: class LaplaceMultigrid : public Laplacian { public: - LaplaceMultigrid(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); + LaplaceMultigrid(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, + Mesh *mesh_in = nullptr); ~LaplaceMultigrid() {}; void setCoefA(const Field2D &val) override { diff --git a/src/invert/laplace/impls/mumps/mumps_laplace.hxx b/src/invert/laplace/impls/mumps/mumps_laplace.hxx index 6c37e63fa9..82fc082b7e 100644 --- a/src/invert/laplace/impls/mumps/mumps_laplace.hxx +++ b/src/invert/laplace/impls/mumps/mumps_laplace.hxx @@ -36,7 +36,7 @@ class LaplaceMumps; class LaplaceMumps : public Laplacian { public: - LaplaceMumps(Options *UNUSED(opt) = nullptr, const CELL_LOC UNUSED(loc) = CELL_CENTRE, Mesh *UNUSED(mesh_in) = mesh) { + LaplaceMumps(Options *UNUSED(opt) = nullptr, const CELL_LOC UNUSED(loc) = CELL_CENTRE, Mesh *UNUSED(mesh_in) = nullptr) { throw BoutException("Mumps library not available"); } @@ -76,7 +76,7 @@ public: class LaplaceMumps : public Laplacian { public: - LaplaceMumps(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); + LaplaceMumps(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = nullptr); ~LaplaceMumps() { mumps_struc.job = -2; dmumps_c(&mumps_struc); diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.hxx b/src/invert/laplace/impls/naulin/naulin_laplace.hxx index b58d63f678..d1088722ab 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.hxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.hxx @@ -37,7 +37,7 @@ class LaplaceNaulin; */ class LaplaceNaulin : public Laplacian { public: - LaplaceNaulin(Options *opt = NULL, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); + LaplaceNaulin(Options *opt = NULL, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = nullptr); ~LaplaceNaulin(); // ACoef is not implemented because the delp2solver that we use can probably diff --git a/src/invert/laplace/impls/pdd/pdd.hxx b/src/invert/laplace/impls/pdd/pdd.hxx index 0dff1dd6e0..62c2031afb 100644 --- a/src/invert/laplace/impls/pdd/pdd.hxx +++ b/src/invert/laplace/impls/pdd/pdd.hxx @@ -34,13 +34,14 @@ class LaplacePDD; #ifndef __LAPLACE_PDD_H__ #define __LAPLACE_PDD_H__ +#include #include #include #include class LaplacePDD : public Laplacian { public: - LaplacePDD(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh) + LaplacePDD(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = nullptr) : Laplacian(opt, loc, mesh_in), Acoef(0.0), Ccoef(1.0), Dcoef(1.0), PDD_COMM_XV(123), PDD_COMM_Y(456) { Acoef.setLocation(location); diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.cxx b/src/invert/laplace/impls/petsc/petsc_laplace.cxx index c986b25910..fe9a4f2c62 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.cxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.cxx @@ -27,6 +27,7 @@ #include "petsc_laplace.hxx" +#include #include #include #include diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.hxx b/src/invert/laplace/impls/petsc/petsc_laplace.hxx index 61b094ace2..b71b504f59 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.hxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.hxx @@ -37,7 +37,7 @@ class LaplacePetsc; class LaplacePetsc : public Laplacian { public: - LaplacePetsc(Options *UNUSED(opt) = nullptr, const CELL_LOC UNUSED(loc) = CELL_CENTRE, Mesh *UNUSED(mesh_in) = mesh) { + LaplacePetsc(Options *UNUSED(opt) = nullptr, const CELL_LOC UNUSED(loc) = CELL_CENTRE, Mesh *UNUSED(mesh_in) = nullptr) { throw BoutException("No PETSc solver available"); } @@ -68,7 +68,7 @@ public: class LaplacePetsc : public Laplacian { public: - LaplacePetsc(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); + LaplacePetsc(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = nullptr); ~LaplacePetsc() { KSPDestroy( &ksp ); VecDestroy( &xs ); diff --git a/src/invert/laplace/impls/serial_band/serial_band.cxx b/src/invert/laplace/impls/serial_band/serial_band.cxx index f773f90996..38c0e3ffa3 100644 --- a/src/invert/laplace/impls/serial_band/serial_band.cxx +++ b/src/invert/laplace/impls/serial_band/serial_band.cxx @@ -27,6 +27,7 @@ #include #include "serial_band.hxx" +#include #include #include #include diff --git a/src/invert/laplace/impls/serial_band/serial_band.hxx b/src/invert/laplace/impls/serial_band/serial_band.hxx index 05a7649e32..41c29b8875 100644 --- a/src/invert/laplace/impls/serial_band/serial_band.hxx +++ b/src/invert/laplace/impls/serial_band/serial_band.hxx @@ -36,7 +36,7 @@ class LaplaceSerialBand; class LaplaceSerialBand : public Laplacian { public: - LaplaceSerialBand(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE, Mesh *mesh_in = mesh); + LaplaceSerialBand(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE, Mesh *mesh_in = 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 f2d98e3f22..efc40e7126 100644 --- a/src/invert/laplace/impls/serial_tri/serial_tri.cxx +++ b/src/invert/laplace/impls/serial_tri/serial_tri.cxx @@ -27,6 +27,7 @@ #include "globals.hxx" #include "serial_tri.hxx" +#include #include #include #include diff --git a/src/invert/laplace/impls/serial_tri/serial_tri.hxx b/src/invert/laplace/impls/serial_tri/serial_tri.hxx index 519dc9a361..770b71050d 100644 --- a/src/invert/laplace/impls/serial_tri/serial_tri.hxx +++ b/src/invert/laplace/impls/serial_tri/serial_tri.hxx @@ -35,7 +35,7 @@ class LaplaceSerialTri; class LaplaceSerialTri : public Laplacian { public: - LaplaceSerialTri(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); + LaplaceSerialTri(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = nullptr); ~LaplaceSerialTri(){}; using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/shoot/shoot_laplace.cxx b/src/invert/laplace/impls/shoot/shoot_laplace.cxx index 405fef96be..6a813409ef 100644 --- a/src/invert/laplace/impls/shoot/shoot_laplace.cxx +++ b/src/invert/laplace/impls/shoot/shoot_laplace.cxx @@ -32,6 +32,7 @@ */ #include "shoot_laplace.hxx" +#include #include #include #include diff --git a/src/invert/laplace/impls/shoot/shoot_laplace.hxx b/src/invert/laplace/impls/shoot/shoot_laplace.hxx index db6e1bd617..a31780a032 100644 --- a/src/invert/laplace/impls/shoot/shoot_laplace.hxx +++ b/src/invert/laplace/impls/shoot/shoot_laplace.hxx @@ -37,7 +37,7 @@ class LaplaceShoot; class LaplaceShoot : public Laplacian { public: - LaplaceShoot(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE, Mesh *mesh_in = mesh); + LaplaceShoot(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE, Mesh *mesh_in = nullptr); ~LaplaceShoot(){}; using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/spt/spt.cxx b/src/invert/laplace/impls/spt/spt.cxx index 705f198c29..9b4d9e7452 100644 --- a/src/invert/laplace/impls/spt/spt.cxx +++ b/src/invert/laplace/impls/spt/spt.cxx @@ -32,6 +32,7 @@ */ #include +#include #include #include #include diff --git a/src/invert/laplace/impls/spt/spt.hxx b/src/invert/laplace/impls/spt/spt.hxx index 44b7d98eb8..1bdc064b98 100644 --- a/src/invert/laplace/impls/spt/spt.hxx +++ b/src/invert/laplace/impls/spt/spt.hxx @@ -66,7 +66,7 @@ class LaplaceSPT; */ class LaplaceSPT : public Laplacian { public: - LaplaceSPT(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE, Mesh *mesh_in = mesh); + LaplaceSPT(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE, Mesh *mesh_in = nullptr); ~LaplaceSPT(); using Laplacian::setCoefA; diff --git a/src/invert/laplace/invert_laplace.cxx b/src/invert/laplace/invert_laplace.cxx index 038d7adc6c..5d5fde2518 100644 --- a/src/invert/laplace/invert_laplace.cxx +++ b/src/invert/laplace/invert_laplace.cxx @@ -31,6 +31,7 @@ * */ +#include #include #include #include @@ -52,7 +53,7 @@ /// Laplacian inversion initialisation. Called once at the start to get settings Laplacian::Laplacian(Options *options, const CELL_LOC loc, Mesh *mesh_in) - : location(loc), localmesh(mesh_in) { + : location(loc), localmesh(mesh_in==nullptr ? bout::globals::mesh : mesh_in) { if (options == nullptr) { // Use the default options diff --git a/src/invert/laplace/laplacefactory.cxx b/src/invert/laplace/laplacefactory.cxx index 777cecb2dc..e112249d57 100644 --- a/src/invert/laplace/laplacefactory.cxx +++ b/src/invert/laplace/laplacefactory.cxx @@ -41,6 +41,10 @@ Laplacian* LaplaceFactory::createLaplacian(Options *options, const CELL_LOC loc, if (options == nullptr) options = Options::getRoot()->getSection("laplace"); + if (mesh_in == nullptr) { + mesh_in = bout::globals::mesh; + } + std::string type; if(mesh_in->firstX() && mesh_in->lastX()) { diff --git a/src/invert/laplace/laplacefactory.hxx b/src/invert/laplace/laplacefactory.hxx index 215ee8790d..780e8c2410 100644 --- a/src/invert/laplace/laplacefactory.hxx +++ b/src/invert/laplace/laplacefactory.hxx @@ -11,7 +11,8 @@ class LaplaceFactory { /// Return a pointer to the only instance static LaplaceFactory* getInstance(); - Laplacian *createLaplacian(Options *options = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); + Laplacian *createLaplacian(Options *options = nullptr, const CELL_LOC loc = CELL_CENTRE, + Mesh *mesh_in = nullptr); private: LaplaceFactory() {} // Prevent instantiation of this class diff --git a/src/invert/laplacexy/laplacexy.cxx b/src/invert/laplacexy/laplacexy.cxx index 6fb32d5f14..dc7ec8b97e 100644 --- a/src/invert/laplacexy/laplacexy.cxx +++ b/src/invert/laplacexy/laplacexy.cxx @@ -8,6 +8,7 @@ #include #include +#include #include #include @@ -26,7 +27,7 @@ static PetscErrorCode laplacePCapply(PC pc,Vec x,Vec y) { } LaplaceXY::LaplaceXY(Mesh *m, Options *opt, const CELL_LOC loc) - : localmesh(m==nullptr ? mesh : m), location(loc) { + : localmesh(m==nullptr ? bout::globals::mesh : m), location(loc) { Timer timer("invert"); if (opt == nullptr) { diff --git a/src/invert/laplacexz/laplacexz.cxx b/src/invert/laplacexz/laplacexz.cxx index 1703c0b076..05bba728ca 100644 --- a/src/invert/laplacexz/laplacexz.cxx +++ b/src/invert/laplacexz/laplacexz.cxx @@ -9,7 +9,7 @@ LaplaceXZ* LaplaceXZ::create(Mesh *m, Options *options, const CELL_LOC loc) { if (m == nullptr) { // use global mesh - m = mesh; + m = bout::globals::mesh; } if (options == nullptr) { diff --git a/src/invert/parderiv/impls/cyclic/cyclic.hxx b/src/invert/parderiv/impls/cyclic/cyclic.hxx index 328c213a64..390228b57f 100644 --- a/src/invert/parderiv/impls/cyclic/cyclic.hxx +++ b/src/invert/parderiv/impls/cyclic/cyclic.hxx @@ -41,11 +41,12 @@ #include "invert_parderiv.hxx" #include "dcomplex.hxx" +#include #include "utils.hxx" class InvertParCR : public InvertPar { public: - InvertParCR(Options *opt, Mesh *mesh_in = mesh); + InvertParCR(Options *opt, Mesh *mesh_in = bout::globals::mesh); ~InvertParCR(); using InvertPar::solve; diff --git a/src/invert/parderiv/parderiv_factory.hxx b/src/invert/parderiv/parderiv_factory.hxx index d5a3a5cb17..5701f5ef48 100644 --- a/src/invert/parderiv/parderiv_factory.hxx +++ b/src/invert/parderiv/parderiv_factory.hxx @@ -11,9 +11,9 @@ class ParDerivFactory { /// Return a pointer to the only instance static ParDerivFactory* getInstance(); - InvertPar* createInvertPar(Mesh* mesh_in = mesh); - InvertPar *createInvertPar(const char *type, Options *opt = nullptr, Mesh* mesh_in = mesh); - InvertPar* createInvertPar(Options *opts, Mesh* mesh_in = mesh); + InvertPar* createInvertPar(Mesh* mesh_in = bout::globals::mesh); + InvertPar *createInvertPar(const char *type, Options *opt = nullptr, Mesh* mesh_in = bout::globals::mesh); + InvertPar* createInvertPar(Options *opts, Mesh* mesh_in = bout::globals::mesh); private: ParDerivFactory() {} // Prevent instantiation of this class static ParDerivFactory* instance; ///< The only instance of this class (Singleton) diff --git a/src/mesh/boundary_factory.cxx b/src/mesh/boundary_factory.cxx index 29f7df07a0..e2454d7371 100644 --- a/src/mesh/boundary_factory.cxx +++ b/src/mesh/boundary_factory.cxx @@ -1,6 +1,7 @@ #include #include #include +#include #include #include diff --git a/src/mesh/boundary_region.cxx b/src/mesh/boundary_region.cxx index bee26d439e..e1bdd88201 100644 --- a/src/mesh/boundary_region.cxx +++ b/src/mesh/boundary_region.cxx @@ -1,4 +1,5 @@ +#include #include #include #include diff --git a/src/mesh/boundary_standard.cxx b/src/mesh/boundary_standard.cxx index ba06f23feb..f80876d7fc 100644 --- a/src/mesh/boundary_standard.cxx +++ b/src/mesh/boundary_standard.cxx @@ -1,3 +1,4 @@ +#include #include #include #include @@ -29,6 +30,7 @@ void verifyNumPoints(BoundaryRegion *region, int ptsRequired) { int ptsAvailGlobal, ptsAvailLocal, ptsAvail; std::string side, gridType; + Mesh* mesh = region->localmesh; //Initialise var in case of no match and CHECK<=2 ptsAvail = ptsRequired; //Ensures test passes without exception @@ -121,6 +123,8 @@ void BoundaryDirichlet::apply(Field2D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); // Decide which generator to use @@ -308,6 +312,8 @@ void BoundaryDirichlet::apply(Field3D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); // Decide which generator to use @@ -545,6 +551,8 @@ void BoundaryDirichlet::apply_ddt(Field2D &f) { } void BoundaryDirichlet::apply_ddt(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) for(int z=0;zLocalNz;z++) @@ -573,6 +581,8 @@ void BoundaryDirichlet_O3::apply(Field2D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); // Decide which generator to use @@ -758,6 +768,8 @@ void BoundaryDirichlet_O3::apply(Field3D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); // Decide which generator to use @@ -957,6 +969,8 @@ void BoundaryDirichlet_O3::apply_ddt(Field2D &f) { } void BoundaryDirichlet_O3::apply_ddt(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Field3D *dt = f.timeDeriv(); bndry->first() ; @@ -988,6 +1002,8 @@ void BoundaryDirichlet_O4::apply(Field2D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); // Decide which generator to use @@ -1185,6 +1201,8 @@ void BoundaryDirichlet_O4::apply(Field3D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); // Decide which generator to use @@ -1389,6 +1407,8 @@ void BoundaryDirichlet_O4::apply_ddt(Field2D &f) { } void BoundaryDirichlet_O4::apply_ddt(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) for(int z=0;zLocalNz;z++) @@ -1417,6 +1437,8 @@ void BoundaryDirichlet_4thOrder::apply(Field2D &f) { } void BoundaryDirichlet_4thOrder::apply(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); // Set (at 4th order) the value at the mid-point between the guard cell and the grid cell to be val for(bndry->first(); !bndry->isDone(); bndry->next1d()) for(int z=0;zLocalNz;z++) { @@ -1432,6 +1454,8 @@ void BoundaryDirichlet_4thOrder::apply_ddt(Field2D &f) { } void BoundaryDirichlet_4thOrder::apply_ddt(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) for(int z=0;zLocalNz;z++) @@ -1452,6 +1476,8 @@ BoundaryOp* BoundaryNeumann_NonOrthogonal::clone(BoundaryRegion *region, const s } void BoundaryNeumann_NonOrthogonal::apply(Field2D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Coordinates *metric = f.getCoordinates(); // Calculate derivatives for metric use mesh->communicate(f); @@ -1492,6 +1518,8 @@ void BoundaryNeumann_NonOrthogonal::apply(Field2D &f) { } void BoundaryNeumann_NonOrthogonal::apply(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Coordinates *metric = f.getCoordinates(); // Calculate derivatives for metric use mesh->communicate(f); @@ -1557,6 +1585,8 @@ void BoundaryNeumann::apply(Field2D &f,BoutReal t) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Coordinates *metric = f.getCoordinates(); bndry->first(); @@ -1754,6 +1784,8 @@ void BoundaryNeumann::apply(Field3D &f) { void BoundaryNeumann::apply(Field3D &f,BoutReal t) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Coordinates *metric = f.getCoordinates(); bndry->first(); @@ -1957,6 +1989,8 @@ void BoundaryNeumann::apply_ddt(Field2D &f) { } void BoundaryNeumann::apply_ddt(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) for(int z=0;zLocalNz;z++) @@ -1979,6 +2013,8 @@ void BoundaryNeumann_O4::apply(Field2D &f) { } void BoundaryNeumann_O4::apply(Field2D &f,BoutReal t) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); // Set (at 4th order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used @@ -2036,6 +2072,8 @@ void BoundaryNeumann_O4::apply(Field3D &f) { } void BoundaryNeumann_O4::apply(Field3D &f,BoutReal t) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); // Decide which generator to use @@ -2091,6 +2129,8 @@ void BoundaryNeumann_O4::apply_ddt(Field2D &f) { } void BoundaryNeumann_O4::apply_ddt(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) for(int z=0;zLocalNz;z++) @@ -2121,6 +2161,8 @@ void BoundaryNeumann_4thOrder::apply(Field2D &f) { } void BoundaryNeumann_4thOrder::apply(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Coordinates *metric = f.getCoordinates(); // Set (at 4th order) the gradient at the mid-point between the guard cell and the grid cell to be val // This sets the value of the co-ordinate derivative, i.e. DDX/DDY not Grad_par/Grad_perp.x @@ -2139,6 +2181,8 @@ void BoundaryNeumann_4thOrder::apply_ddt(Field2D &f) { } void BoundaryNeumann_4thOrder::apply_ddt(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) for(int z=0;zLocalNz;z++) @@ -2164,6 +2208,8 @@ void BoundaryNeumannPar::apply(Field2D &f) { } void BoundaryNeumannPar::apply(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Coordinates *metric = f.getCoordinates(); for(bndry->first(); !bndry->isDone(); bndry->next()) for(int z=0;zLocalNz;z++) @@ -2219,6 +2265,8 @@ void BoundaryRobin::apply(Field2D &f) { } void BoundaryRobin::apply(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); if(fabs(bval) < 1.e-12) { for(bndry->first(); !bndry->isDone(); bndry->next()) for(int z=0;zLocalNz;z++) @@ -2242,6 +2290,8 @@ void BoundaryConstGradient::apply(Field2D &f){ } void BoundaryConstGradient::apply(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); for(bndry->first(); !bndry->isDone(); bndry->next()) for(int z=0;zLocalNz;z++) f(bndry->x, bndry->y, z) = 2.*f(bndry->x - bndry->bx, bndry->y - bndry->by, z) - f(bndry->x - 2*bndry->bx,bndry->y - 2*bndry->by,z); @@ -2290,6 +2340,8 @@ void BoundaryZeroLaplace::apply(Field2D &f) { } void BoundaryZeroLaplace::apply(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); int ncz = mesh->LocalNz; Coordinates *metric = f.getCoordinates(); @@ -2380,6 +2432,8 @@ void BoundaryZeroLaplace2::apply(Field2D &f) { } void BoundaryZeroLaplace2::apply(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); int ncz = mesh->LocalNz; ASSERT0(ncz % 2 == 0); // Allocation assumes even number @@ -2476,6 +2530,8 @@ void BoundaryConstLaplace::apply(Field3D &f) { throw BoutException("ERROR: Can't apply Zero Laplace condition to non-X boundaries\n"); } + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Coordinates *metric = f.getCoordinates(); int ncz = mesh->LocalNz; @@ -2544,6 +2600,9 @@ void BoundaryDivCurl::apply(Vector2D &UNUSED(f)) { } void BoundaryDivCurl::apply(Vector3D &var) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = var.x.getMesh()); + int jx, jy, jz, jzp, jzm; BoutReal tmp; @@ -2652,6 +2711,8 @@ void BoundaryFree_O2::apply(Field2D &f) { // Set (at 2nd order) the value at the mid-point between the guard cell and the grid cell to be val // N.B. Only first guard cells (closest to the grid) should ever be used + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); // Check for staggered grids @@ -2749,9 +2810,10 @@ void BoundaryFree_O2::apply(Field2D &f) { void BoundaryFree_O2::apply(Field3D &f) { // Extrapolate from the last evolved simulation cells into the guard cells at 3rd order. + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); - // Check for staggered grids CELL_LOC loc = f.getLocation(); @@ -2868,6 +2930,8 @@ void BoundaryFree_O2::apply_ddt(Field2D &f) { } void BoundaryFree_O2::apply_ddt(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) for(int z=0;zLocalNz;z++) @@ -2889,6 +2953,8 @@ BoundaryOp* BoundaryFree_O3::clone(BoundaryRegion *region, const std::listlocalmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); // Check for staggered grids @@ -2987,6 +3053,8 @@ void BoundaryFree_O3::apply(Field2D &f) { void BoundaryFree_O3::apply(Field3D &f) { // Extrapolate from the last evolved simulation cells into the guard cells at 3rd order. + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); @@ -3114,6 +3182,8 @@ void BoundaryFree_O3::apply_ddt(Field2D &f) { } void BoundaryFree_O3::apply_ddt(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); Field3D *dt = f.timeDeriv(); for(bndry->first(); !bndry->isDone(); bndry->next()) for(int z=0;zLocalNz;z++) @@ -3165,6 +3235,9 @@ void BoundaryRelax::apply_ddt(Field2D &f) { void BoundaryRelax::apply_ddt(Field3D &f) { TRACE("BoundaryRelax::apply_ddt(Field3D)"); + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); + // Make a copy of f Field3D g = f; // NOTE: This is not very efficient... copying entire field // Apply the boundary to g @@ -3238,6 +3311,9 @@ void BoundaryToFieldAligned::apply(Field2D &f, BoutReal t) { } void BoundaryToFieldAligned::apply(Field3D &f, BoutReal t) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); + //NOTE: This is not very efficient... updating entire field f = mesh->fromFieldAligned(f); @@ -3257,6 +3333,8 @@ void BoundaryToFieldAligned::apply_ddt(Field2D &f) { } void BoundaryToFieldAligned::apply_ddt(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); f = mesh->fromFieldAligned(f); ddt(f) = mesh->fromFieldAligned(ddt(f)); op->apply_ddt(f); @@ -3281,6 +3359,8 @@ void BoundaryFromFieldAligned::apply(Field2D &f, BoutReal t) { } void BoundaryFromFieldAligned::apply(Field3D &f, BoutReal t) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); //NOTE: This is not very efficient... shifting entire field f = mesh->toFieldAligned(f); @@ -3300,6 +3380,8 @@ void BoundaryFromFieldAligned::apply_ddt(Field2D &f) { } void BoundaryFromFieldAligned::apply_ddt(Field3D &f) { + Mesh* mesh = bndry->localmesh; + ASSERT1(mesh = f.getMesh()); f = mesh->toFieldAligned(f); ddt(f) = mesh->toFieldAligned(ddt(f)); op->apply_ddt(f); diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 974737da8a..7a4121a2d2 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -799,7 +799,7 @@ const Field3D Coordinates::Grad2_par2(const Field3D &f, CELL_LOC outloc, const s #include // Delp2 uses same coefficients as inversion code -const Field2D Coordinates::Delp2(const Field2D& f, CELL_LOC outloc, bool useFFT) { +const Field2D Coordinates::Delp2(const Field2D& f, CELL_LOC outloc, bool UNUSED(useFFT)) { TRACE("Coordinates::Delp2( Field2D )"); ASSERT1(location == outloc || outloc == CELL_DEFAULT); diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index c924ae11e0..b3ba1cc95f 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -24,7 +24,7 @@ **************************************************************************/ #include -#include +#include #include #include #include @@ -449,7 +449,9 @@ const Field2D Div_par_LtoC(const Field2D &var) { } const Field3D Div_par_LtoC(const Field3D &var) { - Field3D result; + Mesh* mesh = var.getMesh(); + + Field3D result(mesh); result.allocate(); Coordinates *metric = var.getCoordinates(CELL_CENTRE); @@ -478,7 +480,9 @@ const Field2D Div_par_CtoL(const Field2D &var) { } const Field3D Div_par_CtoL(const Field3D &var) { - Field3D result; + Mesh* mesh = var.getMesh(); + + Field3D result(mesh); result.allocate(); Coordinates *metric = var.getCoordinates(CELL_CENTRE); diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index d67fe767a6..f65b93911b 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -180,6 +180,9 @@ namespace FV { const Field3D D4DY4(const Field3D &d_in, const Field3D &f_in) { ASSERT2(d_in.getLocation() == f_in.getLocation()); + Mesh* mesh = d_in.getMesh(); + ASSERT1(mesh = f_in.getMesh()); + Field3D result = 0.0; result.setLocation(f_in.getLocation()); @@ -233,6 +236,8 @@ namespace FV { Field3D result = 0.0; result.setLocation(f_in.getLocation()); + Mesh* mesh = f_in.getMesh(); + // Convert to field aligned coordinates Field3D f = mesh->toFieldAligned(f_in); @@ -339,6 +344,7 @@ namespace FV { } void communicateFluxes(Field3D &f) { + Mesh* mesh = f.getMesh(); // Use X=0 as temporary buffer if (mesh->xstart != 2) diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 8ed15a1cc6..81ea690561 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -2223,7 +2223,7 @@ void BoutMesh::addBoundaryRegions() { all_boundaries.emplace_back("RGN_UPPER_Y"); // Inner X - if(mesh->firstX() && !mesh->periodicX) { + if(firstX() && !periodicX) { addRegion3D("RGN_INNER_X", Region(0, xstart-1, ystart, yend, 0, LocalNz-1, LocalNy, LocalNz, maxregionblocksize)); addRegion2D("RGN_INNER_X", Region(0, xstart-1, ystart, yend, 0, 0, @@ -2240,7 +2240,7 @@ void BoutMesh::addBoundaryRegions() { } // Outer X - if(mesh->firstX() && !mesh->periodicX) { + if(firstX() && !periodicX) { addRegion3D("RGN_OUTER_X", Region(xend+1, LocalNx-1, ystart, yend, 0, LocalNz-1, LocalNy, LocalNz, maxregionblocksize)); addRegion2D("RGN_OUTER_X", Region(xend+1, LocalNx-1, ystart, yend, 0, 0, diff --git a/src/mesh/interpolation/interpolation_factory.cxx b/src/mesh/interpolation/interpolation_factory.cxx index fb57cb1e75..f77d2f2658 100644 --- a/src/mesh/interpolation/interpolation_factory.cxx +++ b/src/mesh/interpolation/interpolation_factory.cxx @@ -47,16 +47,19 @@ Interpolation* InterpolationFactory::create(Options *options, Mesh *mesh) { Interpolation* InterpolationFactory::create(const std::string &name, Options *options, Mesh *localmesh) { // If no options section passed (e.g. for a variable), then use the // "interpolation" section - if (options == nullptr) + if (options == nullptr) { options = Options::getRoot()->getSection("interpolation"); + } // Use the global mesh if none passed - if (localmesh == nullptr) - localmesh = mesh; + if (localmesh == nullptr) { + localmesh = bout::globals::mesh; + } auto interp = findInterpolation(name); - if (interp == nullptr) + if (interp == nullptr) { throw BoutException("Could not find interpolation method '%s'", name.c_str()); + } return interp(localmesh); } diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index eb62510406..1d37b1fe7a 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -53,8 +53,9 @@ inline BoutReal sgn(BoutReal val) { return (BoutReal(0) < val) - (val < BoutReal // Calculate all the coefficients needed for the spline interpolation // dir MUST be either +1 or -1 -FCIMap::FCIMap(Mesh &mesh, int dir, bool zperiodic) - : dir(dir), boundary_mask(mesh), corner_boundary_mask(mesh), y_prime(&mesh) { +FCIMap::FCIMap(Mesh &mesh_in, int dir, bool zperiodic) + : mesh(mesh_in), dir(dir), boundary_mask(mesh_in), corner_boundary_mask(mesh_in), + y_prime(&mesh_in) { interp = InterpolationFactory::getInstance()->create(&mesh); interp->setYOffset(dir); @@ -249,10 +250,10 @@ const Field3D FCIMap::integrate(Field3D &f) const { result.allocate(); result.setLocation(f.getLocation()); - int nz = mesh->LocalNz; + int nz = mesh.LocalNz; - for(int x = mesh->xstart; x <= mesh->xend; x++) { - for(int y = mesh->ystart; y <= mesh->yend; y++) { + for(int x = mesh.xstart; x <= mesh.xend; x++) { + for(int y = mesh.ystart; y <= mesh.yend; y++) { int ynext = y+dir; @@ -269,7 +270,7 @@ const Field3D FCIMap::integrate(Field3D &f) const { if (corner_boundary_mask(x, y, z) || corner_boundary_mask(x - 1, y, z) || corner_boundary_mask(x, y, zm) || corner_boundary_mask(x - 1, y, zm) || - (x == mesh->xstart)) { + (x == mesh.xstart)) { // One of the corners leaves the domain. // Use the cell centre value, since boundary conditions are not // currently applied to corners. diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 293113b91d..d9cfe5af0c 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -40,11 +40,13 @@ class FCIMap { Interpolation *interp; // Cell centre Interpolation *interp_corner; // Cell corner at (x+1, z+1) + Mesh& mesh; + /// Private constructor - must be initialised with mesh FCIMap(); public: /// dir MUST be either +1 or -1 - FCIMap(Mesh& mesh, int dir, bool zperiodic); + FCIMap(Mesh& mesh_in, int dir, bool zperiodic); int dir; /**< Direction of map */ diff --git a/src/mesh/parallel_boundary_op.cxx b/src/mesh/parallel_boundary_op.cxx index 1fdd82bcaa..1b13325682 100644 --- a/src/mesh/parallel_boundary_op.cxx +++ b/src/mesh/parallel_boundary_op.cxx @@ -7,6 +7,8 @@ BoutReal BoundaryOpPar::getValue(int x, int y, int z, BoutReal t) { + Mesh* mesh = bndry->localmesh; + BoutReal xnorm; BoutReal ynorm; BoutReal znorm; @@ -35,6 +37,8 @@ BoutReal BoundaryOpPar::getValue(int x, int y, int z, BoutReal t) { BoutReal BoundaryOpPar::getValue(const BoundaryRegionPar &bndry, BoutReal t) { + Mesh* mesh = bndry.localmesh; + BoutReal xnorm; BoutReal ynorm; BoutReal znorm; diff --git a/src/physics/gyro_average.cxx b/src/physics/gyro_average.cxx index e67f11663a..b8f1a0cd83 100644 --- a/src/physics/gyro_average.cxx +++ b/src/physics/gyro_average.cxx @@ -27,6 +27,7 @@ * **************************************************************/ +#include #include #include #include @@ -89,7 +90,7 @@ const Field2D gyroPade1(const Field2D &f, const Field2D &rho, int flags) { const Field3D gyroPade2(const Field3D &f, BoutReal rho, int flags) { Field3D result = gyroPade1(gyroPade1(f, rho, flags), rho, flags); - mesh->communicate(result); + result.getMesh()->communicate(result); result = 0.5*rho*rho*Delp2( result ); result.applyBoundary("dirichlet"); return result; @@ -97,7 +98,7 @@ const Field3D gyroPade2(const Field3D &f, BoutReal rho, int flags) { const Field3D gyroPade2(const Field3D &f, const Field2D &rho, int flags) { Field3D result = gyroPade1(gyroPade1(f, rho, flags), rho, flags); - mesh->communicate(result); + result.getMesh()->communicate(result); result = 0.5*rho*rho*Delp2( result ); result.applyBoundary("dirichlet"); return result; diff --git a/src/physics/physicsmodel.cxx b/src/physics/physicsmodel.cxx index 98c56a329b..6235b1be4c 100644 --- a/src/physics/physicsmodel.cxx +++ b/src/physics/physicsmodel.cxx @@ -28,7 +28,11 @@ * **************************************************************************/ +#define BOUT_NO_USING_NAMESPACE_BOUTGLOBALS #include +#undef BOUT_NO_USING_NAMESPACE_BOUTGLOBALS + +#include PhysicsModel::PhysicsModel() : solver(nullptr), modelMonitor(this), splitop(false), userprecon(nullptr), @@ -127,7 +131,7 @@ int PhysicsModel::postInit(bool restarting) { // Add mesh information to restart file // Note this is done after reading, so mesh variables // are not overwritten. - mesh->outputVars(restart); + bout::globals::mesh->outputVars(restart); // Version expected by collect routine restart.addOnce(const_cast(BOUT_VERSION), "BOUT_VERSION"); diff --git a/src/physics/smoothing.cxx b/src/physics/smoothing.cxx index 4a3bc507b4..77326e2cb8 100644 --- a/src/physics/smoothing.cxx +++ b/src/physics/smoothing.cxx @@ -32,6 +32,7 @@ #include +#include #include #include #include @@ -461,6 +462,6 @@ const Field3D nl_filter(const Field3D &f, BoutReal w) { /// Perform filtering in Z, Y then X Field3D result = nl_filter_x(nl_filter_y(nl_filter_z(f, w), w), w); /// Communicate boundaries - mesh->communicate(result); + f.getMesh()->communicate(result); return result; } diff --git a/src/physics/sourcex.cxx b/src/physics/sourcex.cxx index 4094251960..1a549e50dc 100644 --- a/src/physics/sourcex.cxx +++ b/src/physics/sourcex.cxx @@ -5,6 +5,8 @@ #include #include +#include +#include #include #include @@ -16,38 +18,42 @@ BoutReal TanH(BoutReal a) { } // create radial buffer zones to set jpar zero near radial boundaries -const Field2D source_tanhx(const Field2D &UNUSED(f), BoutReal swidth, BoutReal slength) { - Field2D result; +const Field2D source_tanhx(const Field2D &f, BoutReal swidth, BoutReal slength) { + Mesh* localmesh = f.getMesh(); + + Field2D result(localmesh); result.allocate(); // create a radial buffer zone to set jpar zero near radial boundary BOUT_FOR(i, result.getRegion("RGN_ALL")) { - BoutReal lx = mesh->GlobalX(i.x()) - slength; + BoutReal lx = localmesh->GlobalX(i.x()) - slength; BoutReal dampl = TanH(lx / swidth); result[i] = 0.5 * (1.0 - dampl); } // Need to communicate boundaries - mesh->communicate(result); + localmesh->communicate(result); return result; } // create radial buffer zones to set jpar zero near radial boundaries -const Field2D source_expx2(const Field2D &UNUSED(f), BoutReal swidth, BoutReal slength) { - Field2D result; +const Field2D source_expx2(const Field2D &f, BoutReal swidth, BoutReal slength) { + Mesh* localmesh = f.getMesh(); + + Field2D result(localmesh); result.allocate(); // create a radial buffer zone to set jpar zero near radial boundary BOUT_FOR(i, result.getRegion("RGN_ALL")) { - BoutReal lx = mesh->GlobalX(i.x()) - slength; + BoutReal lx = localmesh->GlobalX(i.x()) - slength; BoutReal dampl = exp(-lx * lx / swidth / swidth); result[i] = dampl; } // Need to communicate boundaries - mesh->communicate(result); + localmesh->communicate(result); return result; } @@ -55,18 +61,20 @@ const Field2D source_expx2(const Field2D &UNUSED(f), BoutReal swidth, BoutReal s // create radial buffer zones to set jpar zero near radial boundaries const Field3D sink_tanhx(const Field2D &UNUSED(f0), const Field3D &f, BoutReal swidth, BoutReal slength, bool UNUSED(BoutRealspace)) { - Field3D result; + Mesh* localmesh = f.getMesh(); + + Field3D result(localmesh); result.allocate(); // create a radial buffer zone to set jpar zero near radial boundary BOUT_FOR(i, result.getRegion("RGN_ALL")) { - BoutReal rlx = 1. - mesh->GlobalX(i.x()) - slength; + BoutReal rlx = 1. - localmesh->GlobalX(i.x()) - slength; BoutReal dampr = TanH(rlx / swidth); result[i] = 0.5 * (1.0 - dampr) * f[i]; } // Need to communicate boundaries - mesh->communicate(result); + localmesh->communicate(result); return result; } @@ -75,12 +83,14 @@ const Field3D sink_tanhx(const Field2D &UNUSED(f0), const Field3D &f, BoutReal s const Field3D mask_x(const Field3D &f, bool UNUSED(BoutRealspace)) { TRACE("mask_x"); - Field3D result; + Mesh* localmesh = f.getMesh(); + + Field3D result(localmesh); result.allocate(); // create a radial buffer zone to set jpar zero near radial boundary BOUT_FOR(i, result.getRegion("RGN_ALL")) { - BoutReal lx = mesh->GlobalX(i.x()); + BoutReal lx = localmesh->GlobalX(i.x()); BoutReal dampl = TanH(lx / 40.0); BoutReal dampr = TanH((1. - lx) / 40.0); @@ -88,7 +98,7 @@ const Field3D mask_x(const Field3D &f, bool UNUSED(BoutRealspace)) { } // Need to communicate boundaries - mesh->communicate(result); + localmesh->communicate(result); return result; } @@ -98,19 +108,21 @@ const Field3D sink_tanhxl(const Field2D &UNUSED(f0), const Field3D &f, BoutReal BoutReal slength, bool UNUSED(BoutRealspace)) { TRACE("sink_tanhx"); - Field3D result; + Mesh* localmesh = f.getMesh(); + + Field3D result(localmesh); result.allocate(); BOUT_FOR(i, result.getRegion("RGN_ALL")) { - BoutReal lx = mesh->GlobalX(i.x()) - slength; + BoutReal lx = localmesh->GlobalX(i.x()) - slength; BoutReal dampl = TanH(lx / swidth); result[i] = 0.5 * (1.0 - dampl) * f[i]; } // Need to communicate boundaries - mesh->communicate(result); + localmesh->communicate(result); return result; } @@ -119,18 +131,21 @@ const Field3D sink_tanhxl(const Field2D &UNUSED(f0), const Field3D &f, BoutReal const Field3D sink_tanhxr(const Field2D &UNUSED(f0), const Field3D &f, BoutReal swidth, BoutReal slength, bool UNUSED(BoutRealspace)) { TRACE("sink_tanhxr"); - Field3D result; + + Mesh* localmesh = f.getMesh(); + + Field3D result(localmesh); result.allocate(); BOUT_FOR(i, result.getRegion("RGN_ALL")) { - BoutReal rlx = 1. - mesh->GlobalX(i.x()) - slength; + BoutReal rlx = 1. - localmesh->GlobalX(i.x()) - slength; BoutReal dampr = TanH(rlx / swidth); result[i] = 0.5 * (1.0 - dampr) * f[i]; } // Need to communicate boundaries - mesh->communicate(result); + localmesh->communicate(result); return result; } @@ -139,7 +154,9 @@ const Field3D sink_tanhxr(const Field2D &UNUSED(f0), const Field3D &f, BoutReal const Field3D buff_x(const Field3D &f, bool UNUSED(BoutRealspace)) { TRACE("buff_x"); - Field3D result; + Mesh* localmesh = f.getMesh(); + + Field3D result(localmesh); result.allocate(); const BoutReal dampl = 1.e0; @@ -148,7 +165,7 @@ const Field3D buff_x(const Field3D &f, bool UNUSED(BoutRealspace)) { const BoutReal deltar = 0.05; BOUT_FOR(i, result.getRegion("RGN_ALL")) { - BoutReal lx = mesh->GlobalX(i.x()); + BoutReal lx = localmesh->GlobalX(i.x()); BoutReal rlx = 1. - lx; result[i] = (dampl * exp(-(lx * lx) / (deltal * deltal)) + @@ -157,7 +174,7 @@ const Field3D buff_x(const Field3D &f, bool UNUSED(BoutRealspace)) { } // Need to communicate boundaries - mesh->communicate(result); + localmesh->communicate(result); return result; } diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx index 20316aac44..884f64fc4e 100644 --- a/src/solver/impls/arkode/arkode.cxx +++ b/src/solver/impls/arkode/arkode.cxx @@ -135,15 +135,26 @@ int ArkodeSolver::init(int nout, BoutReal tstep) { bool use_precon, use_jacobian, use_vector_abstol,set_linear; BoutReal start_timestep, max_timestep, min_timestep,fixed_timestep; bool imex,expl,impl; // Time-integration method - int MXSUB = mesh->xend - mesh->xstart + 1; + + // Compute band_width_default from actually added fields, to allow for multiple Mesh objects + // + // Previous implementation was equivalent to: + // int MXSUB = mesh->xend - mesh->xstart + 1; + // int band_width_default = n3Dvars()*(MXSUB+2); + int band_width_default = 0; + for (auto fvar : f3d) { + Mesh* localmesh = fvar.var->getMesh(); + band_width_default += localmesh->xend - localmesh->xstart + 3; + } + BoutReal cfl_frac; bool fixed_step; int mxsteps; // Maximum number of steps to take between outputs int mxorder; // Maximum lmm order to be used by the solver {TRACE("Getting options"); - options->get("mudq", mudq, n3Dvars()*(MXSUB+2)); - options->get("mldq", mldq, n3Dvars()*(MXSUB+2)); + options->get("mudq", mudq, band_width_default); + options->get("mldq", mldq, band_width_default); options->get("mukeep", mukeep, n3Dvars()+n2Dvars()); options->get("mlkeep", mlkeep, n3Dvars()+n2Dvars()); options->get("ATOL", abstol, 1.0e-12); @@ -712,11 +723,11 @@ void ArkodeSolver::set_abstol_values(BoutReal *abstolvec_data, int p = 0; // Counter for location in abstolvec_data array // All boundaries - for (const auto &i2d : mesh->getRegion2D("RGN_BNDRY")) { + for (const auto &i2d : bout::globals::mesh->getRegion2D("RGN_BNDRY")) { loop_abstol_values_op(i2d, abstolvec_data, p, f2dtols, f3dtols, true); } // Bulk of points - for (const auto &i2d : mesh->getRegion2D("RGN_NOBNDRY")) { + for (const auto &i2d : bout::globals::mesh->getRegion2D("RGN_NOBNDRY")) { loop_abstol_values_op(i2d, abstolvec_data, p, f2dtols, f3dtols, false); } } @@ -733,7 +744,7 @@ void ArkodeSolver::loop_abstol_values_op(Ind2D UNUSED(i2d), p++; } - for (int jz=0; jz < mesh->LocalNz; jz++) { + for (int jz=0; jz < bout::globals::mesh->LocalNz; jz++) { // Loop over 3D variables for(std::vector::size_type i=0; ixend - mesh->xstart + 1; + + // Compute band_width_default from actually added fields, to allow for multiple Mesh objects + // + // Previous implementation was equivalent to: + // int MXSUB = mesh->xend - mesh->xstart + 1; + // int band_width_default = n3Dvars()*(MXSUB+2); + int band_width_default = 0; + for (auto fvar : f3d) { + Mesh* localmesh = fvar.var->getMesh(); + band_width_default += localmesh->xend - localmesh->xstart + 3; + } + int mxsteps; // Maximum number of steps to take between outputs int mxorder; // Maximum lmm order to be used by the solver int lmm = CV_BDF; int iter = CV_NEWTON; {TRACE("Getting options"); - options->get("mudq", mudq, n3Dvars()*(MXSUB+2)); - options->get("mldq", mldq, n3Dvars()*(MXSUB+2)); + options->get("mudq", mudq, band_width_default); + options->get("mldq", mldq, band_width_default); options->get("mukeep", mukeep, n3Dvars()+n2Dvars()); options->get("mlkeep", mlkeep, n3Dvars()+n2Dvars()); options->get("ATOL", abstol, 1.0e-12); @@ -581,11 +592,11 @@ void CvodeSolver::set_abstol_values(BoutReal* abstolvec_data, std::vectorgetRegion2D("RGN_BNDRY")) { + for (const auto &i2d : bout::globals::mesh->getRegion2D("RGN_BNDRY")) { loop_abstol_values_op(i2d, abstolvec_data, p, f2dtols, f3dtols, true); } // Bulk of points - for (const auto &i2d : mesh->getRegion2D("RGN_NOBNDRY")) { + for (const auto &i2d : bout::globals::mesh->getRegion2D("RGN_NOBNDRY")) { loop_abstol_values_op(i2d, abstolvec_data, p, f2dtols, f3dtols, false); } } @@ -603,7 +614,7 @@ void CvodeSolver::loop_abstol_values_op(Ind2D UNUSED(i2d), p++; } - for (int jz=0; jz < mesh->LocalNz; jz++) { + for (int jz=0; jz < bout::globals::mesh->LocalNz; jz++) { // Loop over 3D variables for(std::vector::size_type i=0; ixend - mesh->xstart + 1; + // Compute band_width_default from actually added fields, to allow for multiple Mesh objects + // + // Previous implementation was equivalent to: + // int MXSUB = mesh->xend - mesh->xstart + 1; + // int band_width_default = n3Dvars()*(MXSUB+2); + int band_width_default = 0; + for (auto fvar : f3d) { + Mesh* localmesh = fvar.var->getMesh(); + band_width_default += localmesh->xend - localmesh->xstart + 3; + } BoutReal abstol, reltol; int maxl; @@ -134,8 +143,8 @@ IdaSolver::~IdaSolver() { } bool use_precon; bool correct_start; - OPTION(options, mudq, n3d*(MXSUB+2)); - OPTION(options, mldq, n3d*(MXSUB+2)); + OPTION(options, mudq, band_width_default); + OPTION(options, mldq, band_width_default); OPTION(options, mukeep, n3d); OPTION(options, mlkeep, n3d); options->get("ATOL", abstol, 1.0e-12); diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.cxx b/src/solver/impls/imex-bdf2/imex-bdf2.cxx index 341bf6a474..b3dc3cf311 100644 --- a/src/solver/impls/imex-bdf2/imex-bdf2.cxx +++ b/src/solver/impls/imex-bdf2/imex-bdf2.cxx @@ -3,6 +3,7 @@ #include "imex-bdf2.hxx" +#include #include #include #include @@ -205,6 +206,9 @@ int IMEXBDF2::init(int nout, BoutReal tstep) { //Set up a snes object stored at the specified location void IMEXBDF2::constructSNES(SNES *snesIn){ + // Use global mesh for now + Mesh* mesh = bout::globals::mesh; + // Nonlinear solver interface (SNES) SNESCreate(BoutComm::get(),snesIn); @@ -1310,6 +1314,9 @@ PetscErrorCode IMEXBDF2::precon(Vec x, Vec f) { */ template< class Op > void IMEXBDF2::loopVars(BoutReal *u) { + // Use global mesh for now + Mesh* mesh = bout::globals::mesh; + // Loop over 2D variables for(auto it = f2d.begin(); it != f2d.end(); ++it) { Op op(it->var, it->F_var); // Initialise the operator diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 90edc53398..c1ae7c8c1d 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -173,10 +173,19 @@ int PetscSolver::init(int NOUT, BoutReal TIMESTEP) { ierr = VecDestroy(&rhs_vec); ///////////// GET OPTIONS ///////////// - int MXSUB = mesh->xend - mesh->xstart + 1; + // Compute band_width_default from actually added fields, to allow for multiple Mesh objects + // + // Previous implementation was equivalent to: + // int MXSUB = mesh->xend - mesh->xstart + 1; + // int band_width_default = n3Dvars()*(MXSUB+2); + int band_width_default = 0; + for (auto fvar : f3d) { + Mesh* localmesh = fvar.var->getMesh(); + band_width_default += localmesh->xend - localmesh->xstart + 3; + } - OPTION(options, mudq, n3d*(MXSUB+2)); - OPTION(options, mldq, n3d*(MXSUB+2)); + OPTION(options, mudq, band_width_default); + OPTION(options, mldq, band_width_default); OPTION(options, mukeep, 0); OPTION(options, mlkeep, 0); OPTION(options, use_precon, false); @@ -376,14 +385,14 @@ int PetscSolver::init(int NOUT, BoutReal TIMESTEP) { PetscInt dof = n3Dvars(); // Maximum allowable size of stencil in x is the number of guard cells - PetscInt stencil_width = mesh->xstart; + PetscInt stencil_width_estimate = options->operator[]("stencil_width_estimate").withDefault(bout::globals::mesh->xstart); // This is the stencil in each direction (*2) along each dimension // (*3), plus the point itself. Not sure if this is correct // though, on several levels: // 1. Ignores corner points used in e.g. brackets // 2. Could have different stencil widths in each dimension // 3. FFTs couple every single point together - PetscInt cols = stencil_width*2*3+1; + PetscInt cols = stencil_width_estimate*2*3+1; PetscInt prealloc; // = cols*dof; ierr = MatCreate(comm,&J);CHKERRQ(ierr); diff --git a/src/solver/impls/power/power.cxx b/src/solver/impls/power/power.cxx index 04a1dd3b18..7f2f777899 100644 --- a/src/solver/impls/power/power.cxx +++ b/src/solver/impls/power/power.cxx @@ -37,8 +37,6 @@ int PowerSolver::init(int nout, BoutReal tstep) { // Allocate memory f0 = Array(nlocal); - // Save the eigenvalue to the output - dump.add(eigenvalue, "eigenvalue", true); eigenvalue = 0.0; // Put starting values into f0 diff --git a/src/solver/impls/power/power.hxx b/src/solver/impls/power/power.hxx index a753167820..913bd583f0 100644 --- a/src/solver/impls/power/power.hxx +++ b/src/solver/impls/power/power.hxx @@ -48,6 +48,14 @@ class PowerSolver : public Solver { int init(int nout, BoutReal tstep) override; int run() override; + + void outputVars(Datafile &outputfile, bool save_repeat=true) override { + // Include base class functionality + this->Solver::outputVars(outputfile, save_repeat); + + // Save the eigenvalue to the output + outputfile.add(eigenvalue, "eigenvalue", true); + } private: BoutReal curtime; // Current simulation time (fixed) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index 6881e92110..be9d4cbfa9 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -27,6 +27,7 @@ #ifdef BOUT_HAS_PVODE +#include #include #include #include @@ -121,10 +122,20 @@ int PvodeSolver::init(int nout, BoutReal tstep) { ///////////// GET OPTIONS ///////////// int pvode_mxstep; - int MXSUB = mesh->xend - mesh->xstart + 1; + // Compute band_width_default from actually added fields, to allow for multiple Mesh objects + // + // Previous implementation was equivalent to: + // int MXSUB = mesh->xend - mesh->xstart + 1; + // int band_width_default = n3Dvars()*(MXSUB+2); + int band_width_default = 0; + for (auto fvar : f3d) { + Mesh* localmesh = fvar.var->getMesh(); + band_width_default += localmesh->xend - localmesh->xstart + 3; + } + - options->get("mudq", mudq, n3d*(MXSUB+2)); - options->get("mldq", mldq, n3d*(MXSUB+2)); + options->get("mudq", mudq, band_width_default); + options->get("mldq", mldq, band_width_default); options->get("mukeep", mukeep, 0); options->get("mlkeep", mlkeep, 0); options->get("ATOL", abstol, 1.0e-12); diff --git a/src/solver/impls/slepc/slepc.cxx b/src/solver/impls/slepc/slepc.cxx index f327137062..aaf5d420c4 100644 --- a/src/solver/impls/slepc/slepc.cxx +++ b/src/solver/impls/slepc/slepc.cxx @@ -606,10 +606,10 @@ void SlepcSolver::monitor(PetscInt its, PetscInt nconv, PetscScalar eigr[], output<create2D("solution", Options::getRoot()->getSection(name), mesh); + v = fact->create2D("solution", Options::getRoot()->getSection(name), v.getMesh()); } else { initial_profile(name, v); } @@ -184,6 +184,8 @@ void Solver::add(Field2D &v, const std::string name) { void Solver::add(Field3D &v, const std::string name) { TRACE("Adding 3D field: Solver::add(%s)", name.c_str()); + Mesh* mesh = v.getMesh(); + #if CHECK > 0 if (varAdded(name)) throw BoutException("Variable '%s' already added to Solver", name.c_str()); @@ -772,6 +774,9 @@ int Solver::call_timestep_monitors(BoutReal simtime, BoutReal lastdt) { int Solver::getLocalN() { + // Use global mesh: FIX THIS! + Mesh* mesh = bout::globals::mesh; + /// Cache the value, so this is not repeatedly called. /// This value should not change after initialisation static int cacheLocalN = -1; @@ -828,6 +833,9 @@ Solver* Solver::create(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; switch(op) { @@ -962,6 +970,9 @@ void Solver::loop_vars_op(Ind2D i2d, BoutReal *udata, int &p, SOLVER_VAR_OP op, /// Loop over variables and domain. Used for all data operations for consistency void Solver::loop_vars(BoutReal *udata, SOLVER_VAR_OP op) { + // Use global mesh: FIX THIS! + Mesh* mesh = bout::globals::mesh; + int p = 0; // Counter for location in udata array // All boundaries @@ -1075,6 +1086,9 @@ void Solver::set_id(BoutReal *udata) { * */ const Field3D Solver::globalIndex(int localStart) { + // Use global mesh: FIX THIS! + Mesh* mesh = bout::globals::mesh; + Field3D index(-1, mesh); // Set to -1, indicating out of domain int n2d = f2d.size(); @@ -1358,11 +1372,11 @@ void Solver::add_mms_sources(BoutReal t) { // Iterate over 2D variables for(const auto& f : f2d) { - *f.F_var += fact->create2D("source", Options::getRoot()->getSection(f.name), mesh, (f.var)->getLocation(), t); + *f.F_var += fact->create2D("source", Options::getRoot()->getSection(f.name), f.var->getMesh(), (f.var)->getLocation(), t); } for(const auto& f : f3d) { - *f.F_var += fact->create3D("source", Options::getRoot()->getSection(f.name), mesh, (f.var)->getLocation(), t); + *f.F_var += fact->create3D("source", Options::getRoot()->getSection(f.name), f.var->getMesh(), (f.var)->getLocation(), t); } } @@ -1371,7 +1385,7 @@ void Solver::calculate_mms_error(BoutReal t) { FieldFactory *fact = FieldFactory::get(); for(const auto& f : f3d) { - Field3D solution = fact->create3D("solution", Options::getRoot()->getSection(f.name), mesh, (f.var)->getLocation(), t); + Field3D solution = fact->create3D("solution", Options::getRoot()->getSection(f.name), f.var->getMesh(), (f.var)->getLocation(), t); *(f.MMS_err) = *(f.var) - solution; } diff --git a/tests/integrated/test-petsc_laplace/data/BOUT.inp b/tests/integrated/test-petsc_laplace/data/BOUT.inp index 2fa9cca300..f152111b41 100644 --- a/tests/integrated/test-petsc_laplace/data/BOUT.inp +++ b/tests/integrated/test-petsc_laplace/data/BOUT.inp @@ -1,8 +1,11 @@ mz = 129 -#MYG = 0 -grid = "grids/flat_grid.nc" +MYG = 0 dump_format = "nc" +[mesh] +nx = 132 +ny = 1 + [mesh:ddx] first=C4 second=C4 diff --git a/tests/integrated/test-petsc_laplace/grids/flat_grid.nc b/tests/integrated/test-petsc_laplace/grids/flat_grid.nc deleted file mode 100644 index 1d09c4a7b4..0000000000 Binary files a/tests/integrated/test-petsc_laplace/grids/flat_grid.nc and /dev/null differ diff --git a/tests/unit/field/test_field2d.cxx b/tests/unit/field/test_field2d.cxx index 5174afba11..d396e2d9c3 100644 --- a/tests/unit/field/test_field2d.cxx +++ b/tests/unit/field/test_field2d.cxx @@ -18,7 +18,14 @@ #include /// Global mesh +namespace bout{ +namespace globals{ extern Mesh *mesh; +} // namespace globals +} // namespace bout + +// The unit tests use the global mesh +using namespace bout::globals; // Reuse the "standard" fixture for FakeMesh using Field2DTest = FakeMeshFixture; diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 3fa1c709d1..fde483a370 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -18,7 +18,14 @@ #include /// Global mesh +namespace bout{ +namespace globals{ extern Mesh *mesh; +} // namespace globals +} // namespace bout + +// The unit tests use the global mesh +using namespace bout::globals; // Reuse the "standard" fixture for FakeMesh using Field3DTest = FakeMeshFixture; diff --git a/tests/unit/field/test_fieldperp.cxx b/tests/unit/field/test_fieldperp.cxx index f986a9968b..1ebb8ff3a8 100644 --- a/tests/unit/field/test_fieldperp.cxx +++ b/tests/unit/field/test_fieldperp.cxx @@ -17,7 +17,14 @@ #include /// Global mesh +namespace bout{ +namespace globals{ extern Mesh *mesh; +} // namespace globals +} // namespace bout + +// The unit tests use the global mesh +using namespace bout::globals; /// Test fixture to make sure the global mesh is our fake one using FieldPerpTest = FakeMeshFixture; diff --git a/tests/unit/field/test_vector2d.cxx b/tests/unit/field/test_vector2d.cxx index a6ec890c34..9dfcbc5c95 100644 --- a/tests/unit/field/test_vector2d.cxx +++ b/tests/unit/field/test_vector2d.cxx @@ -10,7 +10,14 @@ #include "vector3d.hxx" /// Global mesh +namespace bout{ +namespace globals{ extern Mesh *mesh; +} // namespace globals +} // namespace bout + +// The unit tests use the global mesh +using namespace bout::globals; /// Test fixture to make sure the global mesh is our fake one class Vector2DTest : public ::testing::Test { diff --git a/tests/unit/field/test_vector3d.cxx b/tests/unit/field/test_vector3d.cxx index eae149755d..965e28f9ae 100644 --- a/tests/unit/field/test_vector3d.cxx +++ b/tests/unit/field/test_vector3d.cxx @@ -9,7 +9,14 @@ #include "vector3d.hxx" /// Global mesh +namespace bout{ +namespace globals{ extern Mesh *mesh; +} // namespace globals +} // namespace bout + +// The unit tests use the global mesh +using namespace bout::globals; /// Test fixture to make sure the global mesh is our fake one class Vector3DTest : public ::testing::Test { diff --git a/tests/unit/include/bout/test_region.cxx b/tests/unit/include/bout/test_region.cxx index a8496e232c..4e7d59be02 100644 --- a/tests/unit/include/bout/test_region.cxx +++ b/tests/unit/include/bout/test_region.cxx @@ -15,7 +15,14 @@ #include /// Global mesh +namespace bout{ +namespace globals{ extern Mesh *mesh; +} // namespace globals +} // namespace bout + +// The unit tests use the global mesh +using namespace bout::globals; /// Test fixture to make sure the global mesh is our fake one using RegionTest = FakeMeshFixture; diff --git a/tests/unit/include/test_derivs.cxx b/tests/unit/include/test_derivs.cxx index dadd9e8910..eef95ee9f2 100644 --- a/tests/unit/include/test_derivs.cxx +++ b/tests/unit/include/test_derivs.cxx @@ -13,6 +13,9 @@ #include #include +// The unit tests use the global mesh +using namespace bout::globals; + // Some basic sanity checks for the derivative kernels. Checks the // derivatives of sin(R) where R = {X, Y, Z} for each R // individually. To make this as fast as possible, we use only a diff --git a/tests/unit/mesh/parallel/test_shiftedmetric.cxx b/tests/unit/mesh/parallel/test_shiftedmetric.cxx index 4d8c2885a7..6d36c776cc 100644 --- a/tests/unit/mesh/parallel/test_shiftedmetric.cxx +++ b/tests/unit/mesh/parallel/test_shiftedmetric.cxx @@ -3,7 +3,14 @@ #include "bout/paralleltransform.hxx" #include "test_extras.hxx" -extern Mesh* mesh; +// The unit tests use the global mesh +using namespace bout::globals; + +namespace bout{ +namespace globals{ +extern Mesh *mesh; +} // namespace globals +} // namespace bout class ShiftedMetricTest : public ::testing::Test { public: diff --git a/tests/unit/mesh/test_boundary_factory.cxx b/tests/unit/mesh/test_boundary_factory.cxx index c21a37c7a4..15f69686f6 100644 --- a/tests/unit/mesh/test_boundary_factory.cxx +++ b/tests/unit/mesh/test_boundary_factory.cxx @@ -6,7 +6,14 @@ #include "test_extras.hxx" /// Global mesh +namespace bout{ +namespace globals{ extern Mesh *mesh; +} // namespace globals +} // namespace bout + +// The unit tests use the global mesh +using namespace bout::globals; class TestBoundary : public BoundaryOp { public: diff --git a/tests/unit/mesh/test_coordinates.cxx b/tests/unit/mesh/test_coordinates.cxx index 91781334b9..bb62537011 100644 --- a/tests/unit/mesh/test_coordinates.cxx +++ b/tests/unit/mesh/test_coordinates.cxx @@ -7,7 +7,13 @@ #include "test_extras.hxx" /// Global mesh -extern Mesh *mesh; +namespace bout { +namespace globals { +extern Mesh* mesh; +} +} // namespace bout + +using bout::globals::mesh; using CoordinatesTest = FakeMeshFixture; diff --git a/tests/unit/mesh/test_interpolation.cxx b/tests/unit/mesh/test_interpolation.cxx index f48b5c5c52..a4ad6ffe9d 100644 --- a/tests/unit/mesh/test_interpolation.cxx +++ b/tests/unit/mesh/test_interpolation.cxx @@ -19,7 +19,14 @@ /////// /// Global mesh +namespace bout{ +namespace globals{ extern Mesh *mesh; +} // namespace globals +} // namespace bout + +// The unit tests use the global mesh +using namespace bout::globals; /// Test fixture to make sure the global mesh is our fake one class Field3DInterpToTest : public ::testing::Test { diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx index 0dd4a91b1f..02309ce750 100644 --- a/tests/unit/test_extras.hxx +++ b/tests/unit/test_extras.hxx @@ -206,19 +206,19 @@ class FakeMeshFixture : public ::testing::Test { public: FakeMeshFixture() { // Delete any existing mesh - if (mesh != nullptr) { - delete mesh; - mesh = nullptr; + if (bout::globals::mesh != nullptr) { + delete bout::globals::mesh; + bout::globals::mesh = nullptr; } - mesh = new FakeMesh(nx, ny, nz); + bout::globals::mesh = new FakeMesh(nx, ny, nz); output_info.disable(); - mesh->createDefaultRegions(); + bout::globals::mesh->createDefaultRegions(); output_info.enable(); } ~FakeMeshFixture() { - delete mesh; - mesh = nullptr; + delete bout::globals::mesh; + bout::globals::mesh = nullptr; } static constexpr int nx = 3; diff --git a/tools/pylib/_boutcore_build/helper.cxx.in b/tools/pylib/_boutcore_build/helper.cxx.in index 9e6ffc1076..2007c06e65 100644 --- a/tools/pylib/_boutcore_build/helper.cxx.in +++ b/tools/pylib/_boutcore_build/helper.cxx.in @@ -6,6 +6,7 @@ cat < #include #include +#include #include EOF for ftype in $fields @@ -163,7 +164,7 @@ void c_set_f2d_part_(Field2D * f2d, const double data,int xs,int xe, int dx,int } Mesh * c_get_global_mesh(){ - return mesh; + return bout::globals::mesh; } void c_mesh_normalise(Mesh * msh, double norm){ diff --git a/tools/pylib/_boutcore_build/helper.h.in b/tools/pylib/_boutcore_build/helper.h.in index ec86ee7e82..24c26bdf04 100644 --- a/tools/pylib/_boutcore_build/helper.h.in +++ b/tools/pylib/_boutcore_build/helper.h.in @@ -93,7 +93,7 @@ public: solver->setModel(this); bout_monitor = new BoutMonitor(); solver->addMonitor(bout_monitor, Solver::BACK); - solver->outputVars(dump); + solver->outputVars(bout::globals::dump); }; void free(){ delete solver; @@ -123,7 +123,7 @@ void throw_BoutException(std::string err){ } Datafile * c_get_global_datafile(){ - return &dump; + return &bout::globals::dump; } EOF