From c9a8bc00651107e8284f375868bc853309946071 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 17:50:55 +0000 Subject: [PATCH 01/42] Move global variables into namespace bout::globals --- include/bout/physicsmodel.hxx | 8 ++++ include/bout/solver.hxx | 2 + include/datafile.hxx | 84 +++++++++++++++++------------------ include/globals.hxx | 4 ++ src/physics/physicsmodel.cxx | 2 + 5 files changed, 58 insertions(+), 42 deletions(-) diff --git a/include/bout/physicsmodel.hxx b/include/bout/physicsmodel.hxx index 6e98ecebbf..13db3e8a8c 100644 --- a/include/bout/physicsmodel.hxx +++ b/include/bout/physicsmodel.hxx @@ -43,6 +43,14 @@ class PhysicsModel; #include "solver.hxx" #include "unused.hxx" #include "bout/macro_for_each.hxx" + +#ifndef BOUT_NO_USING_NAMESPACE_BOUTGLOBALS +// Include using statement by default in user code. +// Macro allows us to include physicsmodel.hxx without the using statement in +// library code. +using namespace bout::globals; +#endif // BOUT_NO_USING_NAMESPACE_BOUTGLOBALS + /*! Base class for physics models */ diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 411011f82e..24037e7f1d 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 diff --git a/include/datafile.hxx b/include/datafile.hxx index ecfd369312..a19246f34d 100644 --- a/include/datafile.hxx +++ b/include/datafile.hxx @@ -137,63 +137,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/globals.hxx b/include/globals.hxx index e98fe36345..5763ebe34d 100644 --- a/include/globals.hxx +++ b/include/globals.hxx @@ -31,6 +31,8 @@ #include "datafile.hxx" #include "bout/macro_for_each.hxx" +namespace bout { +namespace globals { #ifndef GLOBALORIGIN #define GLOBAL extern #define SETTING(name, val) extern name @@ -83,5 +85,7 @@ GLOBAL Datafile dump; #undef GLOBAL #undef SETTING +} // namespace globals +} // namespace bout #endif // __GLOBALS_H__ diff --git a/src/physics/physicsmodel.cxx b/src/physics/physicsmodel.cxx index 98c56a329b..cfd9cf1061 100644 --- a/src/physics/physicsmodel.cxx +++ b/src/physics/physicsmodel.cxx @@ -28,7 +28,9 @@ * **************************************************************************/ +#define BOUT_NO_USING_NAMESPACE_BOUTGLOBALS #include +#undef BOUT_NO_USING_NAMESPACE_BOUTGLOBALS PhysicsModel::PhysicsModel() : solver(nullptr), modelMonitor(this), splitop(false), userprecon(nullptr), From ffec7b607cf7e57ba97364ad0be117565d799652 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 18:42:39 +0000 Subject: [PATCH 02/42] Add namespaces to extern Mesh* mesh --- include/boundary_region.hxx | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) 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() {} From 759a6239c4c862494fcc30cc9efcfe7680bb644b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 18:43:32 +0000 Subject: [PATCH 03/42] Make Solver::outputVars virtual, overload in Power Allows Power solver to add "eigenvalue" to the output as part of outputVars, without needing the global 'dump' Datafile. --- include/bout/solver.hxx | 2 +- src/solver/impls/power/power.cxx | 2 -- src/solver/impls/power/power.hxx | 8 ++++++++ 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 24037e7f1d..5d0108f588 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -292,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/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) From de60d012857f41a228361a4cea52cd7a3b059ab3 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 18:59:28 +0000 Subject: [PATCH 04/42] Remove use of global mesh in Vector3D::operator* --- src/field/vector3d.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) 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) { From eda81e7691601c18dfdbb35db3c929772aac3c65 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 19:17:59 +0000 Subject: [PATCH 05/42] Use local Mesh* in Datafile and DataFormat Reduces dependence on bout::globals::mesh --- include/datafile.hxx | 3 ++- include/dataformat.hxx | 6 ++++++ src/bout++.cxx | 2 +- src/fileio/datafile.cxx | 15 +++++++++------ src/fileio/dataformat.cxx | 3 +++ src/fileio/formatfactory.cxx | 12 +++++++----- src/fileio/formatfactory.hxx | 3 ++- src/fileio/impls/hdf5/h5_format.cxx | 5 +++-- src/fileio/impls/hdf5/h5_format.hxx | 7 ++++--- src/fileio/impls/netcdf/nc_format.cxx | 4 ++-- src/fileio/impls/netcdf/nc_format.hxx | 7 ++++--- src/fileio/impls/netcdf4/ncxx4.cxx | 4 ++-- src/fileio/impls/netcdf4/ncxx4.hxx | 7 ++++--- src/fileio/impls/pnetcdf/pnetcdf.cxx | 4 ++-- src/fileio/impls/pnetcdf/pnetcdf.hxx | 7 ++++--- 15 files changed, 55 insertions(+), 34 deletions(-) diff --git a/include/datafile.hxx b/include/datafile.hxx index a19246f34d..97a5b8c299 100644 --- a/include/datafile.hxx +++ b/include/datafile.hxx @@ -37,7 +37,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 +78,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? 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/src/bout++.cxx b/src/bout++.cxx index 6e4ef0b5cd..ab4ae8cfb3 100644 --- a/src/bout++.cxx +++ b/src/bout++.cxx @@ -493,7 +493,7 @@ 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")); + dump = Datafile(options->getSection("output"), mesh); /// Open a file for the output if(append) { diff --git a/src/fileio/datafile.cxx b/src/fileio/datafile.cxx index ba6b259395..9f57c39944 100644 --- a/src/fileio/datafile.cxx +++ b/src/fileio/datafile.cxx @@ -46,10 +46,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 +93,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 +108,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 +192,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..2a9eb6069b 100644 --- a/src/fileio/dataformat.cxx +++ b/src/fileio/dataformat.cxx @@ -3,6 +3,9 @@ #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..2de3ca1661 100644 --- a/src/fileio/formatfactory.cxx +++ b/src/fileio/formatfactory.cxx @@ -25,27 +25,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..fbed69b0dc 100644 --- a/src/fileio/impls/hdf5/h5_format.cxx +++ b/src/fileio/impls/hdf5/h5_format.cxx @@ -34,7 +34,7 @@ #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 +69,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 9cc38a5d04..e06d564203 100644 --- a/src/fileio/impls/netcdf/nc_format.cxx +++ b/src/fileio/impls/netcdf/nc_format.cxx @@ -37,7 +37,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 +50,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..8591620496 100644 --- a/src/fileio/impls/netcdf4/ncxx4.cxx +++ b/src/fileio/impls/netcdf4/ncxx4.cxx @@ -38,7 +38,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 +51,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; From 9d056579bac9e7888901c0607c2f62b1d3af2bdb Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 23:33:25 +0000 Subject: [PATCH 06/42] Remove uses of global mesh pointer in BoutMesh --- src/mesh/impls/bout/boutmesh.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 0bf59de09a..e53f0d4454 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -2208,7 +2208,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, @@ -2225,7 +2225,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, From 41f8795c1873d550c9e7f2cb7ccd892a9af86559 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 23:37:15 +0000 Subject: [PATCH 07/42] Remove uses of global mesh pointer in FCITransform/FCIMap --- src/mesh/parallel/fci.cxx | 13 +++++++------ src/mesh/parallel/fci.hxx | 4 +++- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index eb62510406..3f5312b83b 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), corner_boundary_mask(mesh), + y_prime(&mesh) { 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 */ From 38191f6e1c18dcfea0b9df3b3e3f50757def5170 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 23:51:53 +0000 Subject: [PATCH 08/42] Replace some more uses of global mesh difops.cxx and interpolation_factory.cxx --- src/mesh/difops.cxx | 2 ++ src/mesh/interpolation/interpolation_factory.cxx | 11 +++++++---- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index c924ae11e0..f9ff129359 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -452,6 +452,7 @@ const Field3D Div_par_LtoC(const Field3D &var) { Field3D result; result.allocate(); + Mesh* mesh = var.getMesh(); Coordinates *metric = var.getCoordinates(CELL_CENTRE); // NOTE: Need to calculate one more point than centred vars @@ -481,6 +482,7 @@ const Field3D Div_par_CtoL(const Field3D &var) { Field3D result; result.allocate(); + Mesh* mesh = var.getMesh(); Coordinates *metric = var.getCoordinates(CELL_CENTRE); // NOTE: Need to calculate one more point than centred vars 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); } From b494e992449c78781d25d11510b11a1b2482c59b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 23:52:13 +0000 Subject: [PATCH 09/42] Replace uses of global mesh pointer in boundary_standard.cxx --- src/mesh/boundary_standard.cxx | 83 +++++++++++++++++++++++++++++++++- 1 file changed, 82 insertions(+), 1 deletion(-) diff --git a/src/mesh/boundary_standard.cxx b/src/mesh/boundary_standard.cxx index 5e5c0a1c16..6e9e6a250f 100644 --- a/src/mesh/boundary_standard.cxx +++ b/src/mesh/boundary_standard.cxx @@ -29,6 +29,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 +122,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 +311,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 +550,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 +580,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 +767,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 +968,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 +1001,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 +1200,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 +1406,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 +1436,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 +1453,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 +1475,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 +1517,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 +1584,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 +1783,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 +1988,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 +2012,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 +2071,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 +2128,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 +2160,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 +2180,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 +2207,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 +2264,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 +2289,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 +2339,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 +2431,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 +2529,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 +2599,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 +2710,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 +2809,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 +2929,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 +2952,8 @@ BoundaryOp* BoundaryFree_O3::clone(BoundaryRegion *region, const std::listlocalmesh; + ASSERT1(mesh = f.getMesh()); bndry->first(); // Check for staggered grids @@ -2987,6 +3052,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 +3181,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 +3234,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 +3310,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 +3332,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 +3358,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 +3379,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); From b150d10b554a793714c86cb891f419a720d278ea Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 00:04:21 +0000 Subject: [PATCH 10/42] Remove uses of global mesh pointer in FV operations --- include/bout/fv_ops.hxx | 7 +++++++ src/mesh/fv_ops.cxx | 6 ++++++ 2 files changed, 13 insertions(+) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index f1bbf8bd6e..797506c20d 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -182,6 +182,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 +345,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/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) From 057ddf7c7381fcee30f2469c9ba7393140c7f0ce Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 00:06:46 +0000 Subject: [PATCH 11/42] Remove uses of global mesh pointer in ParallelBoundaryOp --- src/mesh/parallel_boundary_op.cxx | 4 ++++ 1 file changed, 4 insertions(+) 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; From 17057e7aa1479094d71621fa607abf18c6489bd7 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 00:09:43 +0000 Subject: [PATCH 12/42] Remove reference to global mesh in smoothing.cxx --- src/physics/smoothing.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics/smoothing.cxx b/src/physics/smoothing.cxx index 4a3bc507b4..ee8888282b 100644 --- a/src/physics/smoothing.cxx +++ b/src/physics/smoothing.cxx @@ -461,6 +461,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; } From 7f5255fd86a00bebaa918cd7ca0da1631aa24c95 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 10:20:12 +0000 Subject: [PATCH 13/42] Remove uses of global mesh in sourcex.cxx --- src/physics/sourcex.cxx | 67 +++++++++++++++++++++++++---------------- 1 file changed, 41 insertions(+), 26 deletions(-) diff --git a/src/physics/sourcex.cxx b/src/physics/sourcex.cxx index 4094251960..21bc195486 100644 --- a/src/physics/sourcex.cxx +++ b/src/physics/sourcex.cxx @@ -16,57 +16,63 @@ 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; } // create radial buffer zones to set jpar zero near radial boundaries -const Field3D sink_tanhx(const Field2D &UNUSED(f0), const Field3D &f, BoutReal swidth, +const Field3D sink_tanhx(const Field2D &f0, const Field3D &f, BoutReal swidth, BoutReal slength, bool UNUSED(BoutRealspace)) { - Field3D result; + Mesh* localmesh = f0.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 +81,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,49 +96,54 @@ const Field3D mask_x(const Field3D &f, bool UNUSED(BoutRealspace)) { } // Need to communicate boundaries - mesh->communicate(result); + localmesh->communicate(result); return result; } // create radial buffer zones to set jpar zero near radial boundaries -const Field3D sink_tanhxl(const Field2D &UNUSED(f0), const Field3D &f, BoutReal swidth, +const Field3D sink_tanhxl(const Field2D &f0, const Field3D &f, BoutReal swidth, BoutReal slength, bool UNUSED(BoutRealspace)) { TRACE("sink_tanhx"); - Field3D result; + Mesh* localmesh = f0.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; } // create radial buffer zones to set jpar zero near radial boundaries -const Field3D sink_tanhxr(const Field2D &UNUSED(f0), const Field3D &f, BoutReal swidth, +const Field3D sink_tanhxr(const Field2D &f0, const Field3D &f, BoutReal swidth, BoutReal slength, bool UNUSED(BoutRealspace)) { TRACE("sink_tanhxr"); - Field3D result; + + Mesh* localmesh = f0.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 +152,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 +163,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 +172,7 @@ const Field3D buff_x(const Field3D &f, bool UNUSED(BoutRealspace)) { } // Need to communicate boundaries - mesh->communicate(result); + localmesh->communicate(result); return result; } From fb2238237978f59c94f66c7b61a2987dd9e1839b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 10:21:44 +0000 Subject: [PATCH 14/42] Replace uses of global mesh in gyro_average.cxx --- src/physics/gyro_average.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics/gyro_average.cxx b/src/physics/gyro_average.cxx index e67f11663a..a71fa23751 100644 --- a/src/physics/gyro_average.cxx +++ b/src/physics/gyro_average.cxx @@ -89,7 +89,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 +97,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; From ad932aaa47f29339254b28c176b05f769b75c5de Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 16:41:18 +0000 Subject: [PATCH 15/42] Initialize result with the mesh of the input field --- src/mesh/difops.cxx | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index f9ff129359..b76f892e22 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -449,10 +449,11 @@ 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(); - Mesh* mesh = var.getMesh(); Coordinates *metric = var.getCoordinates(CELL_CENTRE); // NOTE: Need to calculate one more point than centred vars @@ -479,10 +480,11 @@ 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(); - Mesh* mesh = var.getMesh(); Coordinates *metric = var.getCoordinates(CELL_CENTRE); // NOTE: Need to calculate one more point than centred vars From 1d6570dd295b71a00648ec35a402ff1aa106bc23 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 16:52:37 +0000 Subject: [PATCH 16/42] Get mesh from already-used argument, restore UNUSED where possible --- src/physics/sourcex.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/physics/sourcex.cxx b/src/physics/sourcex.cxx index 21bc195486..284bb0cb2b 100644 --- a/src/physics/sourcex.cxx +++ b/src/physics/sourcex.cxx @@ -57,9 +57,9 @@ const Field2D source_expx2(const Field2D &f, BoutReal swidth, BoutReal slength) } // create radial buffer zones to set jpar zero near radial boundaries -const Field3D sink_tanhx(const Field2D &f0, const Field3D &f, BoutReal swidth, +const Field3D sink_tanhx(const Field2D &UNUSED(f0), const Field3D &f, BoutReal swidth, BoutReal slength, bool UNUSED(BoutRealspace)) { - Mesh* localmesh = f0.getMesh(); + Mesh* localmesh = f.getMesh(); Field3D result(localmesh); result.allocate(); @@ -102,11 +102,11 @@ const Field3D mask_x(const Field3D &f, bool UNUSED(BoutRealspace)) { } // create radial buffer zones to set jpar zero near radial boundaries -const Field3D sink_tanhxl(const Field2D &f0, const Field3D &f, BoutReal swidth, +const Field3D sink_tanhxl(const Field2D &UNUSED(f0), const Field3D &f, BoutReal swidth, BoutReal slength, bool UNUSED(BoutRealspace)) { TRACE("sink_tanhx"); - Mesh* localmesh = f0.getMesh(); + Mesh* localmesh = f.getMesh(); Field3D result(localmesh); @@ -126,11 +126,11 @@ const Field3D sink_tanhxl(const Field2D &f0, const Field3D &f, BoutReal swidth, } // create radial buffer zones to set jpar zero near radial boundaries -const Field3D sink_tanhxr(const Field2D &f0, const Field3D &f, BoutReal swidth, +const Field3D sink_tanhxr(const Field2D &UNUSED(f0), const Field3D &f, BoutReal swidth, BoutReal slength, bool UNUSED(BoutRealspace)) { TRACE("sink_tanhxr"); - Mesh* localmesh = f0.getMesh(); + Mesh* localmesh = f.getMesh(); Field3D result(localmesh); result.allocate(); From 6ce32e66e4dce86f82af12a90c22fd4a16af73ca Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 18:58:26 +0000 Subject: [PATCH 17/42] Use bout::globals in bout++.cxx --- src/bout++.cxx | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/bout++.cxx b/src/bout++.cxx index ab4ae8cfb3..8dfd6c75f1 100644 --- a/src/bout++.cxx +++ b/src/bout++.cxx @@ -475,9 +475,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 +493,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"), mesh); + 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 +556,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 +621,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"); From a812c7c7229144b43e24bd7170ce181ec88cb575 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 18:59:42 +0000 Subject: [PATCH 18/42] Replace mesh with bout::globals::mesh Also replace default argument of 'mesh' with 'nullptr' and then handle 'nullptr' by setting local mesh to 'mesh' in a few places. --- include/field.hxx | 4 ++-- include/interpolation.hxx | 2 +- include/invert_laplace.hxx | 4 ++-- include/invert_parderiv.hxx | 6 +++--- include/mask.hxx | 2 +- src/field/field.cxx | 6 ++---- src/field/field2d.cxx | 2 +- src/field/field3d.cxx | 2 +- src/field/fieldperp.cxx | 2 +- 9 files changed, 14 insertions(+), 16 deletions(-) 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/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 17d1e48507..fa1d7d1ec7 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 = bout::globals::mesh); 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 = bout::globals::mesh); 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 ddca5c3366..cd5edaf7d3 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/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 02f2a2e480..e13f121570 100644 --- a/src/field/field2d.cxx +++ b/src/field/field2d.cxx @@ -93,7 +93,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; } diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 65511f428a..44141b0b01 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -127,7 +127,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; diff --git a/src/field/fieldperp.cxx b/src/field/fieldperp.cxx index 5787f0ce29..e184bfc29a 100644 --- a/src/field/fieldperp.cxx +++ b/src/field/fieldperp.cxx @@ -50,7 +50,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; } From 4427ce443720e89b0795487be86a7b2b29aaf4e8 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 18 Dec 2018 19:22:48 +0000 Subject: [PATCH 19/42] More replacements mesh->bout::globals::mesh in default arguments In Laplacian and InvertPar --- src/invert/laplace/impls/cyclic/cyclic_laplace.hxx | 2 +- src/invert/laplace/impls/multigrid/multigrid_laplace.hxx | 3 ++- src/invert/laplace/impls/mumps/mumps_laplace.hxx | 4 ++-- src/invert/laplace/impls/naulin/naulin_laplace.hxx | 2 +- src/invert/laplace/impls/pdd/pdd.hxx | 2 +- src/invert/laplace/impls/petsc/petsc_laplace.hxx | 4 ++-- src/invert/laplace/impls/serial_band/serial_band.hxx | 2 +- src/invert/laplace/impls/serial_tri/serial_tri.hxx | 2 +- src/invert/laplace/impls/shoot/shoot_laplace.hxx | 2 +- src/invert/laplace/impls/spt/spt.hxx | 2 +- src/invert/laplace/laplacefactory.hxx | 3 ++- src/invert/parderiv/impls/cyclic/cyclic.hxx | 3 ++- src/invert/parderiv/parderiv_factory.hxx | 6 +++--- 13 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx b/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx index b007542569..6d94cc8203 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 = bout::globals::mesh); ~LaplaceCyclic(); using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx index 0067e1bb51..5b2a6ffd18 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 = bout::globals::mesh); ~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..75ca899193 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) = bout::globals::mesh) { 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 = bout::globals::mesh); ~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..f6cd5b3ba2 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 = bout::globals::mesh); ~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..76e54e737e 100644 --- a/src/invert/laplace/impls/pdd/pdd.hxx +++ b/src/invert/laplace/impls/pdd/pdd.hxx @@ -40,7 +40,7 @@ class LaplacePDD; 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 = bout::globals::mesh) : 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.hxx b/src/invert/laplace/impls/petsc/petsc_laplace.hxx index 61b094ace2..f949ebd169 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) = bout::globals::mesh) { 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 = bout::globals::mesh); ~LaplacePetsc() { KSPDestroy( &ksp ); VecDestroy( &xs ); diff --git a/src/invert/laplace/impls/serial_band/serial_band.hxx b/src/invert/laplace/impls/serial_band/serial_band.hxx index 05a7649e32..98607d07a6 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 = bout::globals::mesh); ~LaplaceSerialBand(){}; using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/serial_tri/serial_tri.hxx b/src/invert/laplace/impls/serial_tri/serial_tri.hxx index 519dc9a361..52c025fca9 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 = bout::globals::mesh); ~LaplaceSerialTri(){}; using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/shoot/shoot_laplace.hxx b/src/invert/laplace/impls/shoot/shoot_laplace.hxx index db6e1bd617..5197781ce3 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 = bout::globals::mesh); ~LaplaceShoot(){}; using Laplacian::setCoefA; diff --git a/src/invert/laplace/impls/spt/spt.hxx b/src/invert/laplace/impls/spt/spt.hxx index 44b7d98eb8..52391df81e 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 = bout::globals::mesh); ~LaplaceSPT(); using Laplacian::setCoefA; diff --git a/src/invert/laplace/laplacefactory.hxx b/src/invert/laplace/laplacefactory.hxx index 215ee8790d..79a618ef00 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 = bout::globals::mesh); private: LaplaceFactory() {} // Prevent instantiation of this class 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) From e6ea1108abccd3326505725557546154aa36ea94 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 10:18:23 +0000 Subject: [PATCH 20/42] Add namespace to global mesh in PhysicsModel::postInit() --- src/physics/physicsmodel.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics/physicsmodel.cxx b/src/physics/physicsmodel.cxx index cfd9cf1061..6b408d8a0b 100644 --- a/src/physics/physicsmodel.cxx +++ b/src/physics/physicsmodel.cxx @@ -129,7 +129,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"); From 6deed69143df400414091dae3b7e0fdd4bd2379b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 11:03:07 +0000 Subject: [PATCH 21/42] Reduce uses of global mesh in solvers, add namespace to remaining Loops in solvers are hard-coded to use a single Mesh*: leave as the global 'bout::globals::mesh' for now, and leave removal of global mesh from solvers for later. --- src/solver/impls/arkode/arkode.cxx | 22 ++++++++++++++++------ src/solver/impls/cvode/cvode.cxx | 22 ++++++++++++++++------ src/solver/impls/ida/ida.cxx | 14 +++++++++++--- src/solver/impls/imex-bdf2/imex-bdf2.cxx | 6 ++++++ src/solver/impls/petsc/petsc.cxx | 18 +++++++++++++----- src/solver/impls/pvode/pvode.cxx | 15 ++++++++++++--- src/solver/solver.cxx | 22 ++++++++++++++++++---- 7 files changed, 92 insertions(+), 27 deletions(-) diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx index 20316aac44..d8ed659eeb 100644 --- a/src/solver/impls/arkode/arkode.cxx +++ b/src/solver/impls/arkode/arkode.cxx @@ -135,15 +135,25 @@ 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) { + band_width_default += fvar.var->getNx(); + } + 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 +722,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 +743,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) { + band_width_default += fvar.var->getNx(); + } + 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 +591,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 +613,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) { + band_width_default += fvar.var->getNx(); + } BoutReal abstol, reltol; int maxl; @@ -134,8 +142,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..e586326d27 100644 --- a/src/solver/impls/imex-bdf2/imex-bdf2.cxx +++ b/src/solver/impls/imex-bdf2/imex-bdf2.cxx @@ -205,6 +205,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 +1313,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 a4358118fb..15b3e36861 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -173,10 +173,18 @@ 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) { + band_width_default += fvar.var->getNx(); + } - 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 +384,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/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index 6881e92110..4214ba7281 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -121,10 +121,19 @@ 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) { + band_width_default += fvar.var->getNx(); + } + - 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/solver.cxx b/src/solver/solver.cxx index 4e6b57d289..d9f4e84686 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -158,7 +158,7 @@ void Solver::add(Field2D &v, const std::string name) { FieldFactory *fact = FieldFactory::get(); - v = fact->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; } From 7fff20fba85d7a2c76644b901a88bd9b55bbc4c3 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 17:02:03 +0000 Subject: [PATCH 22/42] Restore exactly original behaviour of mudq/mldq defaults --- src/solver/impls/arkode/arkode.cxx | 3 ++- src/solver/impls/cvode/cvode.cxx | 3 ++- src/solver/impls/ida/ida.cxx | 3 ++- src/solver/impls/petsc/petsc.cxx | 3 ++- src/solver/impls/pvode/pvode.cxx | 3 ++- 5 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx index d8ed659eeb..884f64fc4e 100644 --- a/src/solver/impls/arkode/arkode.cxx +++ b/src/solver/impls/arkode/arkode.cxx @@ -143,7 +143,8 @@ int ArkodeSolver::init(int nout, BoutReal tstep) { // int band_width_default = n3Dvars()*(MXSUB+2); int band_width_default = 0; for (auto fvar : f3d) { - band_width_default += fvar.var->getNx(); + Mesh* localmesh = fvar.var->getMesh(); + band_width_default += localmesh->xend - localmesh->xstart + 3; } BoutReal cfl_frac; diff --git a/src/solver/impls/cvode/cvode.cxx b/src/solver/impls/cvode/cvode.cxx index f1ee3f5b6a..58a3fc5edd 100644 --- a/src/solver/impls/cvode/cvode.cxx +++ b/src/solver/impls/cvode/cvode.cxx @@ -139,7 +139,8 @@ int CvodeSolver::init(int nout, BoutReal tstep) { // int band_width_default = n3Dvars()*(MXSUB+2); int band_width_default = 0; for (auto fvar : f3d) { - band_width_default += fvar.var->getNx(); + Mesh* localmesh = fvar.var->getMesh(); + band_width_default += localmesh->xend - localmesh->xstart + 3; } int mxsteps; // Maximum number of steps to take between outputs diff --git a/src/solver/impls/ida/ida.cxx b/src/solver/impls/ida/ida.cxx index ebcf665761..706ca47246 100644 --- a/src/solver/impls/ida/ida.cxx +++ b/src/solver/impls/ida/ida.cxx @@ -132,7 +132,8 @@ IdaSolver::~IdaSolver() { } // int band_width_default = n3Dvars()*(MXSUB+2); int band_width_default = 0; for (auto fvar : f3d) { - band_width_default += fvar.var->getNx(); + Mesh* localmesh = fvar.var->getMesh(); + band_width_default += localmesh->xend - localmesh->xstart + 3; } BoutReal abstol, reltol; diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 15b3e36861..66bc8ff7ae 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -180,7 +180,8 @@ int PetscSolver::init(int NOUT, BoutReal TIMESTEP) { // int band_width_default = n3Dvars()*(MXSUB+2); int band_width_default = 0; for (auto fvar : f3d) { - band_width_default += fvar.var->getNx(); + Mesh* localmesh = fvar.var->getMesh(); + band_width_default += localmesh->xend - localmesh->xstart + 3; } OPTION(options, mudq, band_width_default); diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index 4214ba7281..3ff6b40207 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -128,7 +128,8 @@ int PvodeSolver::init(int nout, BoutReal tstep) { // int band_width_default = n3Dvars()*(MXSUB+2); int band_width_default = 0; for (auto fvar : f3d) { - band_width_default += fvar.var->getNx(); + Mesh* localmesh = fvar.var->getMesh(); + band_width_default += localmesh->xend - localmesh->xstart + 3; } From aa1879e4117fa2d4e16753c420739375e4b6647e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 11:38:15 +0000 Subject: [PATCH 23/42] Use globals namespace in boutcore --- tools/pylib/_boutcore_build/helper.cxx.in | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/pylib/_boutcore_build/helper.cxx.in b/tools/pylib/_boutcore_build/helper.cxx.in index 73d2f46834..41e79c2385 100644 --- a/tools/pylib/_boutcore_build/helper.cxx.in +++ b/tools/pylib/_boutcore_build/helper.cxx.in @@ -160,7 +160,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){ From 95146f26e4bf1d8d750f86912f3e0ba7071fcd6e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 11:57:54 +0000 Subject: [PATCH 24/42] Use bout::globals namespace in InvertableOperator --- include/bout/invertable_operator.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/bout/invertable_operator.hxx b/include/bout/invertable_operator.hxx index 6886cb4a13..d46fa17278 100644 --- a/include/bout/invertable_operator.hxx +++ b/include/bout/invertable_operator.hxx @@ -133,7 +133,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(); }; From 566614f8014c2282054b08ea0bc1132d58fb4e10 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 12:46:19 +0000 Subject: [PATCH 25/42] Use bout::globals namespace in unit tests Unit tests use the global mesh pointer. --- tests/unit/field/test_field2d.cxx | 7 +++++++ tests/unit/field/test_field3d.cxx | 7 +++++++ tests/unit/field/test_fieldperp.cxx | 7 +++++++ tests/unit/field/test_vector2d.cxx | 7 +++++++ tests/unit/field/test_vector3d.cxx | 7 +++++++ tests/unit/include/bout/test_region.cxx | 7 +++++++ tests/unit/mesh/test_boundary_factory.cxx | 7 +++++++ tests/unit/mesh/test_interpolation.cxx | 7 +++++++ tests/unit/test_extras.hxx | 14 +++++++------- 9 files changed, 63 insertions(+), 7 deletions(-) diff --git a/tests/unit/field/test_field2d.cxx b/tests/unit/field/test_field2d.cxx index f9df75f162..9f5a043815 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 a8b4525b27..0471ff1d66 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 55757dc83e..3dfa1c599f 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/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_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 7bb41523e4..219eee7267 100644 --- a/tests/unit/test_extras.hxx +++ b/tests/unit/test_extras.hxx @@ -176,19 +176,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; From e2d18edfd67e63a392fe53360a564cbfbdf85566 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Thu, 20 Dec 2018 09:47:16 +0000 Subject: [PATCH 26/42] Add bout::globals:: in field_data.cxx FieldData uses the global mesh. Allow this to continue for now, pending a redesign of FieldData; add bout::globals:: namespace to uses of the global mesh pointer. --- src/field/field_data.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/field/field_data.cxx b/src/field/field_data.cxx index 4251ace1b4..9cb00a34a0 100644 --- a/src/field/field_data.cxx +++ b/src/field/field_data.cxx @@ -20,7 +20,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 +28,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 +43,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 +72,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 { From 721c1a8aa9b54f1c6f5fa463d8a6d7c13b4b10ba Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 2 Jan 2019 12:02:07 +0000 Subject: [PATCH 27/42] Replace mesh->bout::globals::mesh in laplacexy.cxx --- src/invert/laplacexy/laplacexy.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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) { From 653057f41e0e7cc521ec8244f420d74e71b1fb90 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 2 Jan 2019 12:29:30 +0000 Subject: [PATCH 28/42] Change dump->bout::globals::dump in boutcore --- tools/pylib/_boutcore_build/helper.h.in | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tools/pylib/_boutcore_build/helper.h.in b/tools/pylib/_boutcore_build/helper.h.in index 4e4e0384ce..f06ff0ab65 100644 --- a/tools/pylib/_boutcore_build/helper.h.in +++ b/tools/pylib/_boutcore_build/helper.h.in @@ -88,7 +88,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; @@ -101,7 +101,7 @@ public: _rhs=__rhs; }; void solve(){ - //solver->outputVars(dump); + //solver->outputVars(bout::globals::dump); solver->solve(); } Solver * getSolver(){ @@ -121,7 +121,7 @@ void throw_BoutException(std::string err){ } Datafile * c_get_global_datafile(){ - return &dump; + return &bout::globals::dump; } EOF From 8d762d735fe1a94090794cba97696831df0fed3a Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 12:01:25 +0000 Subject: [PATCH 29/42] Move 'using namespace bout::globals' to bout.hxx Tests/examples that define their own main function still import bout.hxx, so this is the appropriate place for the 'using' statement. Remove '#include ' from difops.cxx; it should not be included anywhere within the library except in bout++.cxx. --- include/bout.hxx | 7 +++++++ include/bout/physicsmodel.hxx | 7 ------- src/bout++.cxx | 5 ++++- src/mesh/difops.cxx | 2 +- 4 files changed, 12 insertions(+), 9 deletions(-) diff --git a/include/bout.hxx b/include/bout.hxx index 921d81bc7f..28fe47da09 100644 --- a/include/bout.hxx +++ b/include/bout.hxx @@ -65,6 +65,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/physicsmodel.hxx b/include/bout/physicsmodel.hxx index 13db3e8a8c..459ef68112 100644 --- a/include/bout/physicsmodel.hxx +++ b/include/bout/physicsmodel.hxx @@ -44,13 +44,6 @@ class PhysicsModel; #include "unused.hxx" #include "bout/macro_for_each.hxx" -#ifndef BOUT_NO_USING_NAMESPACE_BOUTGLOBALS -// Include using statement by default in user code. -// Macro allows us to include physicsmodel.hxx without the using statement in -// library code. -using namespace bout::globals; -#endif // BOUT_NO_USING_NAMESPACE_BOUTGLOBALS - /*! Base class for physics models */ diff --git a/src/bout++.cxx b/src/bout++.cxx index 8dfd6c75f1..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 diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index b76f892e22..b3ba1cc95f 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -24,7 +24,7 @@ **************************************************************************/ #include -#include +#include #include #include #include From de8552f0a953a6bae2a7c3a95a62783fec3246a5 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 00:06:25 +0000 Subject: [PATCH 30/42] Forward-declare Mesh in datafile.hxx --- include/datafile.hxx | 1 + 1 file changed, 1 insertion(+) diff --git a/include/datafile.hxx b/include/datafile.hxx index 97a5b8c299..008b140a85 100644 --- a/include/datafile.hxx +++ b/include/datafile.hxx @@ -26,6 +26,7 @@ class Datafile; #include #include +class Mesh; #include #include From 0e204f0fc5360cb57e1256fcb685dd0c1d80b055 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 09:50:20 +0000 Subject: [PATCH 31/42] Remove unused localmesh variables in Field2D and Field3D --- src/field/field2d.cxx | 2 -- src/field/field3d.cxx | 2 -- 2 files changed, 4 deletions(-) diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx index e13f121570..d8a6a08375 100644 --- a/src/field/field2d.cxx +++ b/src/field/field2d.cxx @@ -560,8 +560,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 44141b0b01..b41177e1ef 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -1094,8 +1094,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 From 146f60b983c5cc1bbe95cc467554243698e0b743 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 14:14:45 +0000 Subject: [PATCH 32/42] Convert test-petsc_laplace to have mesh in input file Also set ny=1, myg=1. --- .../integrated/test-petsc_laplace/data/BOUT.inp | 7 +++++-- .../test-petsc_laplace/grids/flat_grid.nc | Bin 13744 -> 0 bytes 2 files changed, 5 insertions(+), 2 deletions(-) delete mode 100644 tests/integrated/test-petsc_laplace/grids/flat_grid.nc 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 1d09c4a7b4c6dca7c491384863222714116407d7..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 13744 zcmeI3!AiqG5QfvF*rGjn@ZiD29`qDiITk3Gs2d7r=jS2>9UBQ3v*`K zZbp44qqRh8jXB3gAWoxs55iXy56t?rTzpkNmV3}QKp)Tt^a0-i9DoCG01m(bH~ Date: Wed, 19 Dec 2018 15:49:39 +0000 Subject: [PATCH 33/42] Don't include mesh.hxx in globals.hxx Forward declare class Mesh instead. Requires explicitly including mesh.hxx in several other files, plus a few other extra explicit includes. --- include/bout.hxx | 1 + include/bout/fv_ops.hxx | 1 + include/bout/invertable_operator.hxx | 1 + include/globals.hxx | 3 ++- src/field/field2d.cxx | 1 + src/field/field_data.cxx | 1 + src/field/fieldperp.cxx | 3 +-- src/field/initialprofiles.cxx | 3 +++ src/field/vecops.cxx | 2 ++ src/fileio/datafile.cxx | 2 ++ src/fileio/dataformat.cxx | 1 + src/fileio/formatfactory.cxx | 1 + src/fileio/impls/hdf5/h5_format.cxx | 1 + src/fileio/impls/netcdf/nc_format.cxx | 1 + src/fileio/impls/netcdf4/ncxx4.cxx | 1 + src/invert/lapack_routines.cxx | 1 + src/invert/laplace/impls/cyclic/cyclic_laplace.cxx | 1 + src/invert/laplace/impls/multigrid/multigrid_laplace.cxx | 1 + src/invert/laplace/impls/pdd/pdd.hxx | 1 + src/invert/laplace/impls/petsc/petsc_laplace.cxx | 1 + src/invert/laplace/impls/serial_band/serial_band.cxx | 1 + src/invert/laplace/impls/serial_tri/serial_tri.cxx | 1 + src/invert/laplace/impls/shoot/shoot_laplace.cxx | 1 + src/invert/laplace/impls/spt/spt.cxx | 1 + src/invert/laplace/invert_laplace.cxx | 1 + src/mesh/boundary_factory.cxx | 1 + src/mesh/boundary_region.cxx | 1 + src/mesh/boundary_standard.cxx | 1 + src/physics/gyro_average.cxx | 1 + src/physics/physicsmodel.cxx | 2 ++ src/physics/smoothing.cxx | 1 + src/physics/sourcex.cxx | 2 ++ src/solver/impls/imex-bdf2/imex-bdf2.cxx | 1 + src/solver/impls/pvode/pvode.cxx | 1 + tools/pylib/_boutcore_build/helper.cxx.in | 1 + 35 files changed, 42 insertions(+), 3 deletions(-) diff --git a/include/bout.hxx b/include/bout.hxx index 28fe47da09..278a2a7fa2 100644 --- a/include/bout.hxx +++ b/include/bout.hxx @@ -38,6 +38,7 @@ #include "boutcomm.hxx" +#include #include "globals.hxx" #include "field2d.hxx" diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 797506c20d..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 { /*! diff --git a/include/bout/invertable_operator.hxx b/include/bout/invertable_operator.hxx index d46fa17278..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 diff --git a/include/globals.hxx b/include/globals.hxx index 5763ebe34d..88bab32731 100644 --- a/include/globals.hxx +++ b/include/globals.hxx @@ -27,10 +27,11 @@ #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 diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx index d8a6a08375..3016957eee 100644 --- a/src/field/field2d.cxx +++ b/src/field/field2d.cxx @@ -39,6 +39,7 @@ #include #include +#include #include #include diff --git a/src/field/field_data.cxx b/src/field/field_data.cxx index 9cb00a34a0..8282f9cda4 100644 --- a/src/field/field_data.cxx +++ b/src/field/field_data.cxx @@ -1,4 +1,5 @@ +#include #include #include #include diff --git a/src/field/fieldperp.cxx b/src/field/fieldperp.cxx index e184bfc29a..6dc281e723 100644 --- a/src/field/fieldperp.cxx +++ b/src/field/fieldperp.cxx @@ -28,6 +28,7 @@ #include +#include #include #include #include @@ -542,8 +543,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 a14855a5ab..fcc3fb155e 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/fileio/datafile.cxx b/src/fileio/datafile.cxx index 9f57c39944..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 diff --git a/src/fileio/dataformat.cxx b/src/fileio/dataformat.cxx index 2a9eb6069b..cf09ea170b 100644 --- a/src/fileio/dataformat.cxx +++ b/src/fileio/dataformat.cxx @@ -1,4 +1,5 @@ +#include #include #include #include diff --git a/src/fileio/formatfactory.cxx b/src/fileio/formatfactory.cxx index 2de3ca1661..d2ee7fd8f8 100644 --- a/src/fileio/formatfactory.cxx +++ b/src/fileio/formatfactory.cxx @@ -12,6 +12,7 @@ #include #include +#include #include FormatFactory *FormatFactory::instance = nullptr; diff --git a/src/fileio/impls/hdf5/h5_format.cxx b/src/fileio/impls/hdf5/h5_format.cxx index fbed69b0dc..b938a52202 100644 --- a/src/fileio/impls/hdf5/h5_format.cxx +++ b/src/fileio/impls/hdf5/h5_format.cxx @@ -30,6 +30,7 @@ #include #include +#include #include #include #include diff --git a/src/fileio/impls/netcdf/nc_format.cxx b/src/fileio/impls/netcdf/nc_format.cxx index e06d564203..8e96547e44 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 diff --git a/src/fileio/impls/netcdf4/ncxx4.cxx b/src/fileio/impls/netcdf4/ncxx4.cxx index 8591620496..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 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/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/pdd/pdd.hxx b/src/invert/laplace/impls/pdd/pdd.hxx index 76e54e737e..e25354809e 100644 --- a/src/invert/laplace/impls/pdd/pdd.hxx +++ b/src/invert/laplace/impls/pdd/pdd.hxx @@ -34,6 +34,7 @@ class LaplacePDD; #ifndef __LAPLACE_PDD_H__ #define __LAPLACE_PDD_H__ +#include #include #include #include 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/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_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/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/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/invert_laplace.cxx b/src/invert/laplace/invert_laplace.cxx index c106a7a52e..ffa19a24e6 100644 --- a/src/invert/laplace/invert_laplace.cxx +++ b/src/invert/laplace/invert_laplace.cxx @@ -31,6 +31,7 @@ * */ +#include #include #include #include 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 6e9e6a250f..ff76562006 100644 --- a/src/mesh/boundary_standard.cxx +++ b/src/mesh/boundary_standard.cxx @@ -1,3 +1,4 @@ +#include #include #include #include diff --git a/src/physics/gyro_average.cxx b/src/physics/gyro_average.cxx index a71fa23751..b8f1a0cd83 100644 --- a/src/physics/gyro_average.cxx +++ b/src/physics/gyro_average.cxx @@ -27,6 +27,7 @@ * **************************************************************/ +#include #include #include #include diff --git a/src/physics/physicsmodel.cxx b/src/physics/physicsmodel.cxx index 6b408d8a0b..6235b1be4c 100644 --- a/src/physics/physicsmodel.cxx +++ b/src/physics/physicsmodel.cxx @@ -32,6 +32,8 @@ #include #undef BOUT_NO_USING_NAMESPACE_BOUTGLOBALS +#include + PhysicsModel::PhysicsModel() : solver(nullptr), modelMonitor(this), splitop(false), userprecon(nullptr), userjacobian(nullptr), initialised(false) { diff --git a/src/physics/smoothing.cxx b/src/physics/smoothing.cxx index ee8888282b..77326e2cb8 100644 --- a/src/physics/smoothing.cxx +++ b/src/physics/smoothing.cxx @@ -32,6 +32,7 @@ #include +#include #include #include #include diff --git a/src/physics/sourcex.cxx b/src/physics/sourcex.cxx index 284bb0cb2b..1a549e50dc 100644 --- a/src/physics/sourcex.cxx +++ b/src/physics/sourcex.cxx @@ -5,6 +5,8 @@ #include #include +#include +#include #include #include diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.cxx b/src/solver/impls/imex-bdf2/imex-bdf2.cxx index e586326d27..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 diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index 3ff6b40207..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 diff --git a/tools/pylib/_boutcore_build/helper.cxx.in b/tools/pylib/_boutcore_build/helper.cxx.in index 41e79c2385..fb6f9660b6 100644 --- a/tools/pylib/_boutcore_build/helper.cxx.in +++ b/tools/pylib/_boutcore_build/helper.cxx.in @@ -3,6 +3,7 @@ cat < #include #include +#include #include EOF for ftype in "Field3D" "Field2D" From 71d0a9eb64228a7eca1420e566050eb93649fc41 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 16:04:47 +0000 Subject: [PATCH 34/42] Use nullptr instead of bout::globals::mesh as default argument Handle the nullptr case in the initializer list instead. --- include/invert_laplace.hxx | 2 +- src/invert/laplace/invert_laplace.cxx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/include/invert_laplace.hxx b/include/invert_laplace.hxx index fa1d7d1ec7..9ec16ec4a0 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 = bout::globals::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 diff --git a/src/invert/laplace/invert_laplace.cxx b/src/invert/laplace/invert_laplace.cxx index ffa19a24e6..e2f4760100 100644 --- a/src/invert/laplace/invert_laplace.cxx +++ b/src/invert/laplace/invert_laplace.cxx @@ -53,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 From 32b172c029aae962ead6fa6a12eafc8b40cf0b50 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 19 Dec 2018 18:25:17 +0000 Subject: [PATCH 35/42] nullptr default Mesh* argument to Laplacian implementations constructors nullptr is handled by the base Laplacian class constructor and converted to bout::globals::mesh. --- src/invert/laplace/impls/cyclic/cyclic_laplace.hxx | 2 +- src/invert/laplace/impls/multigrid/multigrid_laplace.hxx | 2 +- src/invert/laplace/impls/mumps/mumps_laplace.hxx | 4 ++-- src/invert/laplace/impls/naulin/naulin_laplace.hxx | 2 +- src/invert/laplace/impls/pdd/pdd.hxx | 2 +- src/invert/laplace/impls/petsc/petsc_laplace.hxx | 4 ++-- src/invert/laplace/impls/serial_band/serial_band.hxx | 2 +- src/invert/laplace/impls/serial_tri/serial_tri.hxx | 2 +- src/invert/laplace/impls/shoot/shoot_laplace.hxx | 2 +- src/invert/laplace/impls/spt/spt.hxx | 2 +- 10 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx b/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx index 6d94cc8203..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 = bout::globals::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.hxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx index 5b2a6ffd18..e0e05328eb 100644 --- a/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx +++ b/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx @@ -133,7 +133,7 @@ private: class LaplaceMultigrid : public Laplacian { public: LaplaceMultigrid(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh *mesh_in = bout::globals::mesh); + 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 75ca899193..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) = bout::globals::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 = bout::globals::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 f6cd5b3ba2..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 = bout::globals::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 e25354809e..62c2031afb 100644 --- a/src/invert/laplace/impls/pdd/pdd.hxx +++ b/src/invert/laplace/impls/pdd/pdd.hxx @@ -41,7 +41,7 @@ class LaplacePDD; class LaplacePDD : public Laplacian { public: - LaplacePDD(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = bout::globals::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.hxx b/src/invert/laplace/impls/petsc/petsc_laplace.hxx index f949ebd169..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) = bout::globals::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 = bout::globals::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.hxx b/src/invert/laplace/impls/serial_band/serial_band.hxx index 98607d07a6..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 = bout::globals::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.hxx b/src/invert/laplace/impls/serial_tri/serial_tri.hxx index 52c025fca9..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 = bout::globals::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.hxx b/src/invert/laplace/impls/shoot/shoot_laplace.hxx index 5197781ce3..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 = bout::globals::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.hxx b/src/invert/laplace/impls/spt/spt.hxx index 52391df81e..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 = bout::globals::mesh); + LaplaceSPT(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE, Mesh *mesh_in = nullptr); ~LaplaceSPT(); using Laplacian::setCoefA; From caabb9c92d65dc60a2c070ffe027ccea77342273 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 23 Dec 2018 15:49:42 +0000 Subject: [PATCH 36/42] Fix some errors, warnings Use bout::globals::mesh in Coordinates, and test_coordinates --- include/bout/invert/laplacexz.hxx | 4 ++-- src/invert/laplacexz/laplacexz.cxx | 2 +- src/mesh/coordinates.cxx | 2 +- tests/unit/mesh/test_coordinates.cxx | 8 +++++++- 4 files changed, 11 insertions(+), 5 deletions(-) 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/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/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index d6ce7a38a0..e05dbdf7bb 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -788,7 +788,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/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; From d0401ee4e66680122045c352d57609be4e7027ab Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 2 Jan 2019 12:10:27 +0000 Subject: [PATCH 37/42] Use forward declares in datafile.hxx Forward declare Field2D, Field3D, Vector2D and Vector3D in datafile.hxx to avoid circular dependencies due to including datafile.hxx in globals.hxx. --- include/datafile.hxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/datafile.hxx b/include/datafile.hxx index 008b140a85..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" @@ -27,6 +23,10 @@ class Datafile; #include #include class Mesh; +class Field2D; +class Field3D; +class Vector2D; +class Vector3D; #include #include From 8886ff87fb642c2705f40f4d63397bfee8356b76 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Thu, 3 Jan 2019 10:14:27 +0000 Subject: [PATCH 38/42] Use bout::globals::dump in SlepcSolver Update SlepcSolver to use 'bout::globals::dump' instead of 'dump'. --- src/solver/impls/slepc/slepc.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) 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< Date: Thu, 3 Jan 2019 10:26:08 +0000 Subject: [PATCH 39/42] Use bout::globals::dump in BOUTMAIN macro BOUTMAIN is defined in physicsmodel.hxx, so a file using it will usually include 'using namespace bout::globals' from physicsmodel.hxx, but in case the 'using' statement is explicitly excluded using the BOUT_NO_USING_NAMESPACE_BOUTGLOBALS macro or if in future the 'using' statement is removed, use the full name 'bout::globals::dump' instead of 'dump'. --- include/bout/physicsmodel.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/bout/physicsmodel.hxx b/include/bout/physicsmodel.hxx index 459ef68112..8f4f3266d1 100644 --- a/include/bout/physicsmodel.hxx +++ b/include/bout/physicsmodel.hxx @@ -311,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; \ From 606c2013b9346661a55daa1270d5ec3c8adaaf44 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 4 Feb 2019 12:43:22 +0000 Subject: [PATCH 40/42] Pass mesh_in in initializer list of FCIMap::FCIMap constructor Avoids potential issues if order of FCIMap member declarations changes in the future. --- src/mesh/parallel/fci.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 3f5312b83b..1d37b1fe7a 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -54,8 +54,8 @@ 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_in, int dir, bool zperiodic) - : mesh(mesh_in), dir(dir), boundary_mask(mesh), corner_boundary_mask(mesh), - y_prime(&mesh) { + : 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); From fc2db927376bc733a09ea7e63ef9606eb116b15c Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 4 Feb 2019 15:36:08 +0000 Subject: [PATCH 41/42] Add 'using namespace bout::globals;' for new unit test files. --- tests/unit/include/test_derivs.cxx | 3 +++ tests/unit/mesh/parallel/test_shiftedmetric.cxx | 9 ++++++++- 2 files changed, 11 insertions(+), 1 deletion(-) 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: From 64cbef311b54db650b9af5b90d25b5c0c23f43f6 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 5 Feb 2019 16:04:13 +0000 Subject: [PATCH 42/42] Use nullptr as default Mesh* argument for Laplacian::create and also for LaplaceFactory::createLaplacian. LaplaceFactory uses bout::globals::mesh if the input 'Mesh* mesh_in' is nullptr. --- include/invert_laplace.hxx | 2 +- src/invert/laplace/laplacefactory.cxx | 4 ++++ src/invert/laplace/laplacefactory.hxx | 2 +- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/include/invert_laplace.hxx b/include/invert_laplace.hxx index a2cca2b2d0..e6a4417208 100644 --- a/include/invert_laplace.hxx +++ b/include/invert_laplace.hxx @@ -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 = bout::globals::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/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 79a618ef00..780e8c2410 100644 --- a/src/invert/laplace/laplacefactory.hxx +++ b/src/invert/laplace/laplacefactory.hxx @@ -12,7 +12,7 @@ class LaplaceFactory { static LaplaceFactory* getInstance(); Laplacian *createLaplacian(Options *options = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh *mesh_in = bout::globals::mesh); + Mesh *mesh_in = nullptr); private: LaplaceFactory() {} // Prevent instantiation of this class