From fee27c95cbf1ac53eb2e7612293e90716548f761 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 3 Nov 2023 11:51:03 +0100 Subject: [PATCH 01/30] Dump field before and after rhs evaluation for debugging --- src/solver/impls/pvode/pvode.cxx | 43 ++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index b2cfd233a9..8a51a9cb19 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -293,6 +293,49 @@ BoutReal PvodeSolver::run(BoutReal tout) { // Check return flag if (flag != SUCCESS) { output_error.write("ERROR CVODE step failed, flag = {:d}\n", flag); + CVodeMemRec* cv_mem = (CVodeMem)cvode_mem; + if (f2d.empty() and v2d.empty() and v3d.empty()) { + Options debug{}; + using namespace std::string_literals; + Mesh* mesh{}; + for (const auto& prefix : {"pre_"s, "residuum_"s}) { + std::vector ffs{}; + std::vector evolve_bndrys{}; + for (const auto& f : f3d) { + Field3D ff{0.}; + ff.allocate(); + ff.setLocation(f.location); + mesh = ff.getMesh(); + debug[fmt::format("{:s}{:s}", prefix, f.name)] = ff; + ffs.push_back(ff); + evolve_bndrys.push_back(f.evolve_bndry); + } + pvode_load_data_f3d(evolve_bndrys, ffs, + prefix == "pre_"s ? udata : N_VDATA(cv_mem->cv_acor)); + } + + for (auto& f : f3d) { + f.F_var->enableTracking(fmt::format("ddt_{:s}", f.name), debug); + setName(f.var, f.name); + } + run_rhs(simtime); + + for (auto& f : f3d) { + debug[f.name] = *f.var; + } + + if (mesh) { + mesh->outputVars(debug); + debug["BOUT_VERSION"].force(bout::version::as_double); + } + + std::string outname = fmt::format( + "{}/BOUT.debug.{}.nc", + Options::root()["datadir"].withDefault("data"), BoutComm::rank()); + + bout::OptionsNetCDF(outname).write(debug); + MPI_Barrier(BoutComm::get()); + } return (-1.0); } From 96da2e9f88e313f4dd388b2c83d701046aedf96c Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 3 Nov 2023 11:49:32 +0100 Subject: [PATCH 02/30] Add setName function --- include/bout/field.hxx | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index c0693ec0fb..0867560c3b 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -683,4 +683,12 @@ inline T floor(const T& var, BoutReal f, const std::string& rgn = "RGN_ALL") { #undef FIELD_FUNC +template , class... Types> +inline T setName(T&& f, const std::string& name, Types... args) { +#if BOUT_USE_TRACK + f.name = fmt::format(name, args...); +#endif + return f; +} + #endif /* FIELD_H */ From e56981ceeb0409d7d48d81e1225e7332c0346519 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 19 Jun 2023 09:33:00 +0200 Subject: [PATCH 03/30] Set div_par and grad_par names --- src/mesh/coordinates.cxx | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 01f0fe46ca..3948c75b94 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1542,7 +1542,11 @@ Field3D Coordinates::Grad_par(const Field3D& var, CELL_LOC outloc, TRACE("Coordinates::Grad_par( Field3D )"); ASSERT1(location == outloc || outloc == CELL_DEFAULT); - return ::DDY(var, outloc, method) * invSg(); + if (invSg == nullptr) { + invSg = std::make_unique(); + (*invSg) = 1.0 / sqrt(g_22); + } + return setName(::DDY(var, outloc, method) * invSg(), "Grad_par({:s})", var.name); } ///////////////////////////////////////////////////////// @@ -1601,7 +1605,7 @@ Field3D Coordinates::Div_par(const Field3D& f, CELL_LOC outloc, f_B.yup(i) = f.yup(i) / Bxy_floc.yup(i); f_B.ydown(i) = f.ydown(i) / Bxy_floc.ydown(i); } - return Bxy * Grad_par(f_B, outloc, method); + return setName(Bxy * Grad_par(f_B, outloc, method), "Div_par({:s})", f.name); } ///////////////////////////////////////////////////////// From 6b2c132e3db61fa4f453e43a6988e8534f6c0a43 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 19 Jun 2023 09:32:25 +0200 Subject: [PATCH 04/30] Dump debug file if PVODE fails Use the new track feature (better name required) to dump the different components of the ddt() as well as the residuum for the evolved fields. --- src/solver/impls/pvode/pvode.cxx | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index 8a51a9cb19..7283b7d0eb 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -42,12 +42,39 @@ #include // contains the enum for types of preconditioning #include // band preconditioner function prototypes +#include + using namespace pvode; void solver_f(integer N, BoutReal t, N_Vector u, N_Vector udot, void* f_data); void solver_gloc(integer N, BoutReal t, BoutReal* u, BoutReal* udot, void* f_data); void solver_cfn(integer N, BoutReal t, N_Vector u, void* f_data); +namespace { +// local only +void pvode_load_data_f3d(const std::vector& evolve_bndrys, + std::vector& ffs, BoutReal* udata) { + int p = 0; + Mesh* mesh = ffs[0].getMesh(); + const int nz = mesh->LocalNz; + for (const auto& bndry : {true, false}) { + for (const auto& i2d : mesh->getRegion2D(bndry ? "RGN_BNDRY" : "RGN_NOBNDRY")) { + for (int jz = 0; jz < nz; jz++) { + // Loop over 3D variables + std::vector::const_iterator evolve_bndry = evolve_bndrys.begin(); + for (std::vector::iterator ff = ffs.begin(); ff != ffs.end(); ++ff) { + if (bndry && !*evolve_bndry) + continue; + (*ff)[mesh->ind2Dto3D(i2d, jz)] = udata[p]; + p++; + } + ++evolve_bndry; + } + } + } +} +} // namespace + const BoutReal ZERO = 0.0; long int iopt[OPT_SIZE]; From df2d66189f4d2e87cc024521ad287936b9c0ba85 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 19 Jun 2023 09:27:19 +0200 Subject: [PATCH 05/30] Add tracking to Field3D This keeps track of all the changes done to the field and stores them to a OptionsObject. --- include/bout/field3d.hxx | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index ba8c8e879e..8992575be6 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -514,6 +514,13 @@ private: /// RegionID over which the field is valid std::optional regionID; + + int tracking_state{0}; + Options* tracking{nullptr}; + std::string selfname{""}; + template + Options* track(const T& change, std::string op); + Options* track(const BoutReal& change, std::string op); }; // Non-member overloaded operators From 263f9fedaad854832bfdbf6a5ac9322eb492c71e Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 19 Jun 2023 09:27:19 +0200 Subject: [PATCH 06/30] Add tracking to Field3D This keeps track of all the changes done to the field and stores them to a OptionsObject. --- include/bout/field3d.hxx | 4 + src/field/field3d.cxx | 44 ++++ src/field/gen_fieldops.jinja | 14 ++ src/field/generated_fieldops.cxx | 348 +++++++++++++++++++++++++++++++ 4 files changed, 410 insertions(+) diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 8992575be6..bf7a9cc180 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -295,6 +295,10 @@ public: /// cuts on closed field lines? bool requiresTwistShift(bool twist_shift_enabled); + /// Enable a special tracking mode for debugging + /// Save all changes that, are done to the field, to tracking + Field3D& enableTracking(const std::string& name, Options& tracking); + ///////////////////////////////////////////////////////// // Data access diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 011353f34a..f0f088b656 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -243,6 +243,7 @@ Field3D& Field3D::operator=(const Field3D& rhs) { } TRACE("Field3D: Assignment from Field3D"); + track(rhs, "operator="); // Copy base slice Field::operator=(rhs); @@ -263,6 +264,7 @@ Field3D& Field3D::operator=(const Field3D& rhs) { Field3D& Field3D::operator=(Field3D&& rhs) { TRACE("Field3D: Assignment from Field3D"); + track(rhs, "operator="); // Move parallel slices or delete existing ones. yup_fields = std::move(rhs.yup_fields); @@ -283,6 +285,7 @@ Field3D& Field3D::operator=(Field3D&& rhs) { Field3D& Field3D::operator=(const Field2D& rhs) { TRACE("Field3D = Field2D"); + track(rhs, "operator="); /// Check that the data is allocated ASSERT1(rhs.isAllocated()); @@ -327,6 +330,7 @@ void Field3D::operator=(const FieldPerp& rhs) { Field3D& Field3D::operator=(const BoutReal val) { TRACE("Field3D = BoutReal"); + track(val, "operator="); // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. @@ -831,3 +835,43 @@ Field3D::getValidRegionWithDefault(const std::string& region_name) const { void Field3D::setRegion(const std::string& region_name) { regionID = fieldmesh->getRegionID(region_name); } + +Field3D& Field3D::enableTracking(const std::string& name, Options& _tracking) { + tracking = &_tracking; + tracking_state = 1; + selfname = name; + return *this; +} + +template +Options* Field3D::track(const T& change, std::string op) { + if (tracking and tracking_state) { + const std::string outname{fmt::format("track_{:s}_{:d}", selfname, tracking_state++)}; + tracking->set(outname, change, "tracking"); + (*tracking)[outname].setAttributes({ + {"operation", op}, +#if BOUT_USE_TRACK + {"rhs.name", change.name}, +#endif + }); + return &(*tracking)[outname]; + } + return nullptr; +} + +template Options* Field3D::track(const Field3D&, std::string); +template Options* Field3D::track(const Field2D&, std::string); +template Options* Field3D::track(const FieldPerp&, std::string); + +Options* Field3D::track(const BoutReal& change, std::string op) { + if (tracking and tracking_state) { + const std::string outname{fmt::format("track_{:s}_{:d}", selfname, tracking_state++)}; + tracking->set(outname, change, "tracking"); + (*tracking)[outname].setAttributes({ + {"operation", op}, + {"rhs.name", "BR"}, + }); + return &(*tracking)[outname]; + } + return nullptr; +} diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja index ecd4e628cc..58b1ae28ba 100644 --- a/src/field/gen_fieldops.jinja +++ b/src/field/gen_fieldops.jinja @@ -61,6 +61,10 @@ } {% endif %} +#if BOUT_USE_TRACK + {{out.name}}.name = fmt::format("{:s} {{operator}} {:s}", {{'"BR"' if lhs == "BoutReal" else lhs.name + ".name"}} + , {{'"BR"' if rhs == "BoutReal" else rhs.name + ".name"}}); +#endif checkData({{out.name}}); return {{out.name}}; } @@ -129,9 +133,19 @@ } {% endif %} + {% if lhs == "Field3D" %} + track(rhs, "operator{{operator}}="); + {% endif %} +#if BOUT_USE_TRACK + name = fmt::format("{:s} {{operator}}= {:s}", this->name, {{'"BR"' if rhs == "BoutReal" else rhs.name + ".name"}}); +#endif + checkData(*this); } else { + {% if lhs == "Field3D" %} + track(rhs, "operator{{operator}}="); + {% endif %} (*this) = (*this) {{operator}} {{rhs.name}}; } return *this; diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx index 6b778acee3..3495d87dbc 100644 --- a/src/field/generated_fieldops.cxx +++ b/src/field/generated_fieldops.cxx @@ -20,6 +20,9 @@ Field3D operator*(const Field3D& lhs, const Field3D& rhs) { result[index] = lhs[index] * rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -42,9 +45,15 @@ Field3D& Field3D::operator*=(const Field3D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs[index]; } + track(rhs, "operator*="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { + track(rhs, "operator*="); (*this) = (*this) * rhs; } return *this; @@ -64,6 +73,9 @@ Field3D operator/(const Field3D& lhs, const Field3D& rhs) { result[index] = lhs[index] / rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -86,9 +98,15 @@ Field3D& Field3D::operator/=(const Field3D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] /= rhs[index]; } + track(rhs, "operator/="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { + track(rhs, "operator/="); (*this) = (*this) / rhs; } return *this; @@ -108,6 +126,9 @@ Field3D operator+(const Field3D& lhs, const Field3D& rhs) { result[index] = lhs[index] + rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -130,9 +151,15 @@ Field3D& Field3D::operator+=(const Field3D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs[index]; } + track(rhs, "operator+="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { + track(rhs, "operator+="); (*this) = (*this) + rhs; } return *this; @@ -152,6 +179,9 @@ Field3D operator-(const Field3D& lhs, const Field3D& rhs) { result[index] = lhs[index] - rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -174,9 +204,15 @@ Field3D& Field3D::operator-=(const Field3D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs[index]; } + track(rhs, "operator-="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { + track(rhs, "operator-="); (*this) = (*this) - rhs; } return *this; @@ -201,6 +237,9 @@ Field3D operator*(const Field3D& lhs, const Field2D& rhs) { } } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -226,9 +265,15 @@ Field3D& Field3D::operator*=(const Field2D& rhs) { } } + track(rhs, "operator*="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { + track(rhs, "operator*="); (*this) = (*this) * rhs; } return *this; @@ -254,6 +299,9 @@ Field3D operator/(const Field3D& lhs, const Field2D& rhs) { } } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -280,9 +328,15 @@ Field3D& Field3D::operator/=(const Field2D& rhs) { } } + track(rhs, "operator/="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { + track(rhs, "operator/="); (*this) = (*this) / rhs; } return *this; @@ -307,6 +361,9 @@ Field3D operator+(const Field3D& lhs, const Field2D& rhs) { } } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -332,9 +389,15 @@ Field3D& Field3D::operator+=(const Field2D& rhs) { } } + track(rhs, "operator+="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { + track(rhs, "operator+="); (*this) = (*this) + rhs; } return *this; @@ -359,6 +422,9 @@ Field3D operator-(const Field3D& lhs, const Field2D& rhs) { } } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -384,9 +450,15 @@ Field3D& Field3D::operator-=(const Field2D& rhs) { } } + track(rhs, "operator-="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { + track(rhs, "operator-="); (*this) = (*this) - rhs; } return *this; @@ -408,6 +480,9 @@ FieldPerp operator*(const Field3D& lhs, const FieldPerp& rhs) { result[index] = lhs[base_ind] * rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -428,6 +503,9 @@ FieldPerp operator/(const Field3D& lhs, const FieldPerp& rhs) { result[index] = lhs[base_ind] / rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -448,6 +526,9 @@ FieldPerp operator+(const Field3D& lhs, const FieldPerp& rhs) { result[index] = lhs[base_ind] + rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -468,6 +549,9 @@ FieldPerp operator-(const Field3D& lhs, const FieldPerp& rhs) { result[index] = lhs[base_ind] - rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -485,6 +569,9 @@ Field3D operator*(const Field3D& lhs, const BoutReal rhs) { result[index] = lhs[index] * rhs; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -504,9 +591,15 @@ Field3D& Field3D::operator*=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs; } + track(rhs, "operator*="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, "BR"); +#endif + checkData(*this); } else { + track(rhs, "operator*="); (*this) = (*this) * rhs; } return *this; @@ -526,6 +619,9 @@ Field3D operator/(const Field3D& lhs, const BoutReal rhs) { result[index] = lhs[index] * tmp; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -546,9 +642,15 @@ Field3D& Field3D::operator/=(const BoutReal rhs) { const auto tmp = 1.0 / rhs; BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= tmp; } + track(rhs, "operator/="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, "BR"); +#endif + checkData(*this); } else { + track(rhs, "operator/="); (*this) = (*this) / rhs; } return *this; @@ -567,6 +669,9 @@ Field3D operator+(const Field3D& lhs, const BoutReal rhs) { result[index] = lhs[index] + rhs; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -586,9 +691,15 @@ Field3D& Field3D::operator+=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs; } + track(rhs, "operator+="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, "BR"); +#endif + checkData(*this); } else { + track(rhs, "operator+="); (*this) = (*this) + rhs; } return *this; @@ -607,6 +718,9 @@ Field3D operator-(const Field3D& lhs, const BoutReal rhs) { result[index] = lhs[index] - rhs; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -626,9 +740,15 @@ Field3D& Field3D::operator-=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs; } + track(rhs, "operator-="); +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, "BR"); +#endif + checkData(*this); } else { + track(rhs, "operator-="); (*this) = (*this) - rhs; } return *this; @@ -653,6 +773,9 @@ Field3D operator*(const Field2D& lhs, const Field3D& rhs) { } } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -676,6 +799,9 @@ Field3D operator/(const Field2D& lhs, const Field3D& rhs) { } } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -699,6 +825,9 @@ Field3D operator+(const Field2D& lhs, const Field3D& rhs) { } } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -722,6 +851,9 @@ Field3D operator-(const Field2D& lhs, const Field3D& rhs) { } } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -738,6 +870,9 @@ Field2D operator*(const Field2D& lhs, const Field2D& rhs) { result[index] = lhs[index] * rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -754,6 +889,10 @@ Field2D& Field2D::operator*=(const Field2D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs[index]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -774,6 +913,9 @@ Field2D operator/(const Field2D& lhs, const Field2D& rhs) { result[index] = lhs[index] / rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -790,6 +932,10 @@ Field2D& Field2D::operator/=(const Field2D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] /= rhs[index]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -810,6 +956,9 @@ Field2D operator+(const Field2D& lhs, const Field2D& rhs) { result[index] = lhs[index] + rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -826,6 +975,10 @@ Field2D& Field2D::operator+=(const Field2D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs[index]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -846,6 +999,9 @@ Field2D operator-(const Field2D& lhs, const Field2D& rhs) { result[index] = lhs[index] - rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -862,6 +1018,10 @@ Field2D& Field2D::operator-=(const Field2D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs[index]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -886,6 +1046,9 @@ FieldPerp operator*(const Field2D& lhs, const FieldPerp& rhs) { result[index] = lhs[base_ind] * rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -906,6 +1069,9 @@ FieldPerp operator/(const Field2D& lhs, const FieldPerp& rhs) { result[index] = lhs[base_ind] / rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -926,6 +1092,9 @@ FieldPerp operator+(const Field2D& lhs, const FieldPerp& rhs) { result[index] = lhs[base_ind] + rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -946,6 +1115,9 @@ FieldPerp operator-(const Field2D& lhs, const FieldPerp& rhs) { result[index] = lhs[base_ind] - rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -961,6 +1133,9 @@ Field2D operator*(const Field2D& lhs, const BoutReal rhs) { result[index] = lhs[index] * rhs; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -976,6 +1151,10 @@ Field2D& Field2D::operator*=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, "BR"); +#endif + checkData(*this); } else { @@ -996,6 +1175,9 @@ Field2D operator/(const Field2D& lhs, const BoutReal rhs) { result[index] = lhs[index] * tmp; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -1012,6 +1194,10 @@ Field2D& Field2D::operator/=(const BoutReal rhs) { const auto tmp = 1.0 / rhs; BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= tmp; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, "BR"); +#endif + checkData(*this); } else { @@ -1031,6 +1217,9 @@ Field2D operator+(const Field2D& lhs, const BoutReal rhs) { result[index] = lhs[index] + rhs; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -1046,6 +1235,10 @@ Field2D& Field2D::operator+=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, "BR"); +#endif + checkData(*this); } else { @@ -1065,6 +1258,9 @@ Field2D operator-(const Field2D& lhs, const BoutReal rhs) { result[index] = lhs[index] - rhs; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -1080,6 +1276,10 @@ Field2D& Field2D::operator-=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, "BR"); +#endif + checkData(*this); } else { @@ -1104,6 +1304,9 @@ FieldPerp operator*(const FieldPerp& lhs, const Field3D& rhs) { result[index] = lhs[index] * rhs[base_ind]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1126,6 +1329,10 @@ FieldPerp& FieldPerp::operator*=(const Field3D& rhs) { (*this)[index] *= rhs[base_ind]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1150,6 +1357,9 @@ FieldPerp operator/(const FieldPerp& lhs, const Field3D& rhs) { result[index] = lhs[index] / rhs[base_ind]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1172,6 +1382,10 @@ FieldPerp& FieldPerp::operator/=(const Field3D& rhs) { (*this)[index] /= rhs[base_ind]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1196,6 +1410,9 @@ FieldPerp operator+(const FieldPerp& lhs, const Field3D& rhs) { result[index] = lhs[index] + rhs[base_ind]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1218,6 +1435,10 @@ FieldPerp& FieldPerp::operator+=(const Field3D& rhs) { (*this)[index] += rhs[base_ind]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1242,6 +1463,9 @@ FieldPerp operator-(const FieldPerp& lhs, const Field3D& rhs) { result[index] = lhs[index] - rhs[base_ind]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1264,6 +1488,10 @@ FieldPerp& FieldPerp::operator-=(const Field3D& rhs) { (*this)[index] -= rhs[base_ind]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1288,6 +1516,9 @@ FieldPerp operator*(const FieldPerp& lhs, const Field2D& rhs) { result[index] = lhs[index] * rhs[base_ind]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1310,6 +1541,10 @@ FieldPerp& FieldPerp::operator*=(const Field2D& rhs) { (*this)[index] *= rhs[base_ind]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1334,6 +1569,9 @@ FieldPerp operator/(const FieldPerp& lhs, const Field2D& rhs) { result[index] = lhs[index] / rhs[base_ind]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1356,6 +1594,10 @@ FieldPerp& FieldPerp::operator/=(const Field2D& rhs) { (*this)[index] /= rhs[base_ind]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1380,6 +1622,9 @@ FieldPerp operator+(const FieldPerp& lhs, const Field2D& rhs) { result[index] = lhs[index] + rhs[base_ind]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1402,6 +1647,10 @@ FieldPerp& FieldPerp::operator+=(const Field2D& rhs) { (*this)[index] += rhs[base_ind]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1426,6 +1675,9 @@ FieldPerp operator-(const FieldPerp& lhs, const Field2D& rhs) { result[index] = lhs[index] - rhs[base_ind]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1448,6 +1700,10 @@ FieldPerp& FieldPerp::operator-=(const Field2D& rhs) { (*this)[index] -= rhs[base_ind]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1468,6 +1724,9 @@ FieldPerp operator*(const FieldPerp& lhs, const FieldPerp& rhs) { result[index] = lhs[index] * rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1484,6 +1743,10 @@ FieldPerp& FieldPerp::operator*=(const FieldPerp& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs[index]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1504,6 +1767,9 @@ FieldPerp operator/(const FieldPerp& lhs, const FieldPerp& rhs) { result[index] = lhs[index] / rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1520,6 +1786,10 @@ FieldPerp& FieldPerp::operator/=(const FieldPerp& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] /= rhs[index]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1540,6 +1810,9 @@ FieldPerp operator+(const FieldPerp& lhs, const FieldPerp& rhs) { result[index] = lhs[index] + rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1556,6 +1829,10 @@ FieldPerp& FieldPerp::operator+=(const FieldPerp& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs[index]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1576,6 +1853,9 @@ FieldPerp operator-(const FieldPerp& lhs, const FieldPerp& rhs) { result[index] = lhs[index] - rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, rhs.name); +#endif checkData(result); return result; } @@ -1592,6 +1872,10 @@ FieldPerp& FieldPerp::operator-=(const FieldPerp& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs[index]; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, rhs.name); +#endif + checkData(*this); } else { @@ -1611,6 +1895,9 @@ FieldPerp operator*(const FieldPerp& lhs, const BoutReal rhs) { result[index] = lhs[index] * rhs; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -1626,6 +1913,10 @@ FieldPerp& FieldPerp::operator*=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} *= {:s}", this->name, "BR"); +#endif + checkData(*this); } else { @@ -1646,6 +1937,9 @@ FieldPerp operator/(const FieldPerp& lhs, const BoutReal rhs) { result[index] = lhs[index] * tmp; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -1661,6 +1955,10 @@ FieldPerp& FieldPerp::operator/=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] /= rhs; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} /= {:s}", this->name, "BR"); +#endif + checkData(*this); } else { @@ -1680,6 +1978,9 @@ FieldPerp operator+(const FieldPerp& lhs, const BoutReal rhs) { result[index] = lhs[index] + rhs; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -1695,6 +1996,10 @@ FieldPerp& FieldPerp::operator+=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} += {:s}", this->name, "BR"); +#endif + checkData(*this); } else { @@ -1714,6 +2019,9 @@ FieldPerp operator-(const FieldPerp& lhs, const BoutReal rhs) { result[index] = lhs[index] - rhs; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", lhs.name, "BR"); +#endif checkData(result); return result; } @@ -1729,6 +2037,10 @@ FieldPerp& FieldPerp::operator-=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs; } +#if BOUT_USE_TRACK + name = fmt::format("{:s} -= {:s}", this->name, "BR"); +#endif + checkData(*this); } else { @@ -1750,6 +2062,9 @@ Field3D operator*(const BoutReal lhs, const Field3D& rhs) { result[index] = lhs * rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", "BR", rhs.name); +#endif checkData(result); return result; } @@ -1767,6 +2082,9 @@ Field3D operator/(const BoutReal lhs, const Field3D& rhs) { result[index] = lhs / rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", "BR", rhs.name); +#endif checkData(result); return result; } @@ -1784,6 +2102,9 @@ Field3D operator+(const BoutReal lhs, const Field3D& rhs) { result[index] = lhs + rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", "BR", rhs.name); +#endif checkData(result); return result; } @@ -1801,6 +2122,9 @@ Field3D operator-(const BoutReal lhs, const Field3D& rhs) { result[index] = lhs - rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", "BR", rhs.name); +#endif checkData(result); return result; } @@ -1816,6 +2140,9 @@ Field2D operator*(const BoutReal lhs, const Field2D& rhs) { result[index] = lhs * rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", "BR", rhs.name); +#endif checkData(result); return result; } @@ -1831,6 +2158,9 @@ Field2D operator/(const BoutReal lhs, const Field2D& rhs) { result[index] = lhs / rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", "BR", rhs.name); +#endif checkData(result); return result; } @@ -1846,6 +2176,9 @@ Field2D operator+(const BoutReal lhs, const Field2D& rhs) { result[index] = lhs + rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", "BR", rhs.name); +#endif checkData(result); return result; } @@ -1861,6 +2194,9 @@ Field2D operator-(const BoutReal lhs, const Field2D& rhs) { result[index] = lhs - rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", "BR", rhs.name); +#endif checkData(result); return result; } @@ -1876,6 +2212,9 @@ FieldPerp operator*(const BoutReal lhs, const FieldPerp& rhs) { result[index] = lhs * rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} * {:s}", "BR", rhs.name); +#endif checkData(result); return result; } @@ -1891,6 +2230,9 @@ FieldPerp operator/(const BoutReal lhs, const FieldPerp& rhs) { result[index] = lhs / rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} / {:s}", "BR", rhs.name); +#endif checkData(result); return result; } @@ -1906,6 +2248,9 @@ FieldPerp operator+(const BoutReal lhs, const FieldPerp& rhs) { result[index] = lhs + rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} + {:s}", "BR", rhs.name); +#endif checkData(result); return result; } @@ -1921,6 +2266,9 @@ FieldPerp operator-(const BoutReal lhs, const FieldPerp& rhs) { result[index] = lhs - rhs[index]; } +#if BOUT_USE_TRACK + result.name = fmt::format("{:s} - {:s}", "BR", rhs.name); +#endif checkData(result); return result; } From bac4ca92a31b60c40e87dd47b2ef764f571a58be Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 25 Apr 2023 09:16:38 +0200 Subject: [PATCH 07/30] cvode: Add option to use Adams Moulton solver instead of BDF --- src/solver/impls/pvode/pvode.cxx | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index 7283b7d0eb..ae5cd783a8 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -214,7 +214,34 @@ int PvodeSolver::init() { } iopt[MXSTEP] = pvode_mxstep; - cvode_mem = CVodeMalloc(neq, solver_f, simtime, u, BDF, NEWTON, SS, &reltol, &abstol, + { + /* ropt[H0] : initial step size. Optional input. */ + + /* ropt[HMAX] : maximum absolute value of step size allowed. * + * Optional input. (Default is infinity). */ + const BoutReal hmax( + (*options)["max_timestep"].doc("Maximum internal timestep").withDefault(-1.)); + if (hmax > 0) { + ropt[HMAX] = hmax; + } + /* ropt[HMIN] : minimum absolute value of step size allowed. * + * Optional input. (Default is 0.0). */ + const BoutReal hmin( + (*options)["min_timestep"].doc("Minimum internal timestep").withDefault(-1.)); + if (hmin > 0) { + ropt[HMIN] = hmin; + } + /* iopt[MAXORD] : maximum lmm order to be used by the solver. * + * Optional input. (Default = 12 for ADAMS, 5 for * + * BDF). */ + const int maxOrder((*options)["max_order"].doc("Maximum order").withDefault(-1)); + if (maxOrder > 0) { + iopt[MAXORD] = maxOrder; + } + } + const bool use_adam((*options)["adams_moulton"].doc("Use Adams Moulton solver instead of BDF").withDefault(false)); + + cvode_mem = CVodeMalloc(neq, solver_f, simtime, u, use_adam ? ADAMS : BDF, NEWTON, SS, &reltol, &abstol, this, nullptr, optIn, iopt, ropt, machEnv); if (cvode_mem == nullptr) { From 708bdcb2ff0a5c346edb43f7e609949b7c7afbd9 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 29 Mar 2023 12:59:22 +0200 Subject: [PATCH 08/30] Expose more pvode option to user --- src/solver/impls/pvode/pvode.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index ae5cd783a8..ac6e981b50 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -212,6 +212,9 @@ int PvodeSolver::init() { for (i = 0; i < OPT_SIZE; i++) { ropt[i] = ZERO; } + /* iopt[MXSTEP] : maximum number of internal steps to be taken by * + * the solver in its attempt to reach tout. * + * Optional input. (Default = 500). */ iopt[MXSTEP] = pvode_mxstep; { From d88b454aa2b2924a26a180440f956911c05a0abc Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 19 Mar 2024 16:04:48 +0100 Subject: [PATCH 09/30] Fix bad cherry-pick --- src/mesh/coordinates.cxx | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 3948c75b94..32774d6229 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1542,10 +1542,6 @@ Field3D Coordinates::Grad_par(const Field3D& var, CELL_LOC outloc, TRACE("Coordinates::Grad_par( Field3D )"); ASSERT1(location == outloc || outloc == CELL_DEFAULT); - if (invSg == nullptr) { - invSg = std::make_unique(); - (*invSg) = 1.0 / sqrt(g_22); - } return setName(::DDY(var, outloc, method) * invSg(), "Grad_par({:s})", var.name); } From 023bc41730de39040a50ae245363945d2447d63b Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 19 Mar 2024 16:35:43 +0100 Subject: [PATCH 10/30] Update to new API --- include/bout/field.hxx | 7 +++++++ src/solver/impls/pvode/pvode.cxx | 4 ++-- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index 0867560c3b..04035f5b76 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -683,6 +683,13 @@ inline T floor(const T& var, BoutReal f, const std::string& rgn = "RGN_ALL") { #undef FIELD_FUNC +template , class... Types> +inline void setName(T& f, const std::string& name, Types... args) { +#if BOUT_USE_TRACK + f.name = fmt::format(name, args...); +#endif +} + template , class... Types> inline T setName(T&& f, const std::string& name, Types... args) { #if BOUT_USE_TRACK diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index ac6e981b50..762fba32d1 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -373,7 +373,7 @@ BoutReal PvodeSolver::run(BoutReal tout) { for (auto& f : f3d) { f.F_var->enableTracking(fmt::format("ddt_{:s}", f.name), debug); - setName(f.var, f.name); + setName(*f.var, f.name); } run_rhs(simtime); @@ -390,7 +390,7 @@ BoutReal PvodeSolver::run(BoutReal tout) { "{}/BOUT.debug.{}.nc", Options::root()["datadir"].withDefault("data"), BoutComm::rank()); - bout::OptionsNetCDF(outname).write(debug); + bout::OptionsIO::create(outname)->write(debug); MPI_Barrier(BoutComm::get()); } return (-1.0); From affc995c4ba482c79677c890cc44ccb47d45b648 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 19 Mar 2024 16:35:53 +0100 Subject: [PATCH 11/30] Fix documentation --- manual/sphinx/user_docs/bout_options.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/manual/sphinx/user_docs/bout_options.rst b/manual/sphinx/user_docs/bout_options.rst index 85a8a17d59..330a0dad7e 100644 --- a/manual/sphinx/user_docs/bout_options.rst +++ b/manual/sphinx/user_docs/bout_options.rst @@ -889,7 +889,7 @@ Fields can also be stored and written:: Options fields; fields["f2d"] = Field2D(1.0); fields["f3d"] = Field3D(2.0); - bout::OptionsIO::create("fields.nc").write(fields); + bout::OptionsIO::create("fields.nc")->write(fields); This allows the input settings and evolving variables to be combined into a single tree (see above on joining trees) and written From 71f5b6adb6a8ad7b8941ba783773897906d870d2 Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Tue, 19 Mar 2024 15:50:33 +0000 Subject: [PATCH 12/30] Apply clang-format changes --- src/solver/impls/pvode/pvode.cxx | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index 762fba32d1..a12f330964 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -242,10 +242,12 @@ int PvodeSolver::init() { iopt[MAXORD] = maxOrder; } } - const bool use_adam((*options)["adams_moulton"].doc("Use Adams Moulton solver instead of BDF").withDefault(false)); + const bool use_adam((*options)["adams_moulton"] + .doc("Use Adams Moulton solver instead of BDF") + .withDefault(false)); - cvode_mem = CVodeMalloc(neq, solver_f, simtime, u, use_adam ? ADAMS : BDF, NEWTON, SS, &reltol, &abstol, - this, nullptr, optIn, iopt, ropt, machEnv); + cvode_mem = CVodeMalloc(neq, solver_f, simtime, u, use_adam ? ADAMS : BDF, NEWTON, SS, + &reltol, &abstol, this, nullptr, optIn, iopt, ropt, machEnv); if (cvode_mem == nullptr) { throw BoutException("\tError: CVodeMalloc failed.\n"); @@ -373,12 +375,12 @@ BoutReal PvodeSolver::run(BoutReal tout) { for (auto& f : f3d) { f.F_var->enableTracking(fmt::format("ddt_{:s}", f.name), debug); - setName(*f.var, f.name); + setName(*f.var, f.name); } run_rhs(simtime); for (auto& f : f3d) { - debug[f.name] = *f.var; + debug[f.name] = *f.var; } if (mesh) { From 31fd46153fad6977524394179d0b83ac51f26b9e Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 20 Mar 2024 09:59:13 +0100 Subject: [PATCH 13/30] Apply recomendations from code-review --- include/bout/field3d.hxx | 6 +++--- src/field/field3d.cxx | 10 +++++----- src/solver/impls/pvode/pvode.cxx | 5 +++-- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index bf7a9cc180..cfde9e5328 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -521,10 +521,10 @@ private: int tracking_state{0}; Options* tracking{nullptr}; - std::string selfname{""}; + std::string selfname; template - Options* track(const T& change, std::string op); - Options* track(const BoutReal& change, std::string op); + Options* track(const T& change, std::string operation); + Options* track(const BoutReal& change, std::string operation); }; // Non-member overloaded operators diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index f0f088b656..2196d6eea4 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -844,12 +844,12 @@ Field3D& Field3D::enableTracking(const std::string& name, Options& _tracking) { } template -Options* Field3D::track(const T& change, std::string op) { - if (tracking and tracking_state) { +Options* Field3D::track(const T& change, std::string operation) { + if (tracking != nullptr and tracking_state != 0) { const std::string outname{fmt::format("track_{:s}_{:d}", selfname, tracking_state++)}; tracking->set(outname, change, "tracking"); (*tracking)[outname].setAttributes({ - {"operation", op}, + {"operation", operation}, #if BOUT_USE_TRACK {"rhs.name", change.name}, #endif @@ -863,12 +863,12 @@ template Options* Field3D::track(const Field3D&, std::string); template Options* Field3D::track(const Field2D&, std::string); template Options* Field3D::track(const FieldPerp&, std::string); -Options* Field3D::track(const BoutReal& change, std::string op) { +Options* Field3D::track(const BoutReal& change, std::string operation) { if (tracking and tracking_state) { const std::string outname{fmt::format("track_{:s}_{:d}", selfname, tracking_state++)}; tracking->set(outname, change, "tracking"); (*tracking)[outname].setAttributes({ - {"operation", op}, + {"operation", operation}, {"rhs.name", "BR"}, }); return &(*tracking)[outname]; diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index a12f330964..f3a96b03af 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -53,7 +53,7 @@ void solver_cfn(integer N, BoutReal t, N_Vector u, void* f_data); namespace { // local only void pvode_load_data_f3d(const std::vector& evolve_bndrys, - std::vector& ffs, BoutReal* udata) { + std::vector& ffs, const BoutReal* udata) { int p = 0; Mesh* mesh = ffs[0].getMesh(); const int nz = mesh->LocalNz; @@ -63,8 +63,9 @@ void pvode_load_data_f3d(const std::vector& evolve_bndrys, // Loop over 3D variables std::vector::const_iterator evolve_bndry = evolve_bndrys.begin(); for (std::vector::iterator ff = ffs.begin(); ff != ffs.end(); ++ff) { - if (bndry && !*evolve_bndry) + if (bndry && !*evolve_bndry) { continue; + } (*ff)[mesh->ind2Dto3D(i2d, jz)] = udata[p]; p++; } From 17e46cfc1c0a835fea474ac68f64a9addbf4379f Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 20 Mar 2024 10:02:28 +0100 Subject: [PATCH 14/30] Use more meaningful names --- src/solver/impls/pvode/pvode.cxx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index f3a96b03af..c389a3c0d1 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -359,18 +359,18 @@ BoutReal PvodeSolver::run(BoutReal tout) { using namespace std::string_literals; Mesh* mesh{}; for (const auto& prefix : {"pre_"s, "residuum_"s}) { - std::vector ffs{}; + std::vector list_of_fields{}; std::vector evolve_bndrys{}; for (const auto& f : f3d) { - Field3D ff{0.}; - ff.allocate(); - ff.setLocation(f.location); - mesh = ff.getMesh(); - debug[fmt::format("{:s}{:s}", prefix, f.name)] = ff; - ffs.push_back(ff); + mesh = f.var->getMesh(); + Field3D to_load{0., mesh}; + to_load.allocate(); + to_load.setLocation(f.location); + debug[fmt::format("{:s}{:s}", prefix, f.name)] = to_load; + list_of_fields.push_back(to_load); evolve_bndrys.push_back(f.evolve_bndry); } - pvode_load_data_f3d(evolve_bndrys, ffs, + pvode_load_data_f3d(evolve_bndrys, list_of_fields, prefix == "pre_"s ? udata : N_VDATA(cv_mem->cv_acor)); } From 9c0ae16ed905588b50f3e4fe634dcedf47de22b5 Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Wed, 20 Mar 2024 09:03:10 +0000 Subject: [PATCH 15/30] Apply clang-format changes --- src/solver/impls/pvode/pvode.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index c389a3c0d1..fe231e1086 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -65,7 +65,7 @@ void pvode_load_data_f3d(const std::vector& evolve_bndrys, for (std::vector::iterator ff = ffs.begin(); ff != ffs.end(); ++ff) { if (bndry && !*evolve_bndry) { continue; - } + } (*ff)[mesh->ind2Dto3D(i2d, jz)] = udata[p]; p++; } From 4a17b4982df9788fc26407db70f14d0ce16098e3 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 20 Mar 2024 10:04:52 +0100 Subject: [PATCH 16/30] Apply suggestions from code review --- src/solver/impls/pvode/pvode.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index fe231e1086..db28f64d86 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -384,12 +384,12 @@ BoutReal PvodeSolver::run(BoutReal tout) { debug[f.name] = *f.var; } - if (mesh) { + if (mesh != nullptr) { mesh->outputVars(debug); debug["BOUT_VERSION"].force(bout::version::as_double); } - std::string outname = fmt::format( + const std::string outname = fmt::format( "{}/BOUT.debug.{}.nc", Options::root()["datadir"].withDefault("data"), BoutComm::rank()); From 2f7c3c0664c954c016b949ea8c199f6f35ac289f Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 20 Mar 2024 16:02:22 +0100 Subject: [PATCH 17/30] Workaround for gcc 9.4 gcc 9.4 is unable to correctly parse the construction for the function argument, if name.change is used directly. First making a copy seems to work around that issue. --- src/field/field3d.cxx | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 2196d6eea4..e0d4dda01d 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -848,10 +848,14 @@ Options* Field3D::track(const T& change, std::string operation) { if (tracking != nullptr and tracking_state != 0) { const std::string outname{fmt::format("track_{:s}_{:d}", selfname, tracking_state++)}; tracking->set(outname, change, "tracking"); + // Workaround for bug in gcc9.4 +#if BOUT_USE_TRACK + const std::string changename = change.name; +#endif (*tracking)[outname].setAttributes({ {"operation", operation}, #if BOUT_USE_TRACK - {"rhs.name", change.name}, + {"rhs.name", changename}, #endif }); return &(*tracking)[outname]; From d611758ae6b02a93e51eeeaebe025b628a23c9b3 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 9 Apr 2024 15:15:30 +0200 Subject: [PATCH 18/30] Add option to debug on failure --- src/solver/impls/pvode/pvode.cxx | 80 ++++++++++++++++++-------------- src/solver/impls/pvode/pvode.hxx | 9 +++- 2 files changed, 52 insertions(+), 37 deletions(-) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index db28f64d86..9dce5d357f 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -247,6 +247,11 @@ int PvodeSolver::init() { .doc("Use Adams Moulton solver instead of BDF") .withDefault(false)); + debug_on_failure = + (*options)["debug_on_failure"] + .doc("Run an aditional rhs if the solver fails with extra tracking") + .withDefault(false); + cvode_mem = CVodeMalloc(neq, solver_f, simtime, u, use_adam ? ADAMS : BDF, NEWTON, SS, &reltol, &abstol, this, nullptr, optIn, iopt, ropt, machEnv); @@ -354,47 +359,52 @@ BoutReal PvodeSolver::run(BoutReal tout) { if (flag != SUCCESS) { output_error.write("ERROR CVODE step failed, flag = {:d}\n", flag); CVodeMemRec* cv_mem = (CVodeMem)cvode_mem; - if (f2d.empty() and v2d.empty() and v3d.empty()) { - Options debug{}; - using namespace std::string_literals; - Mesh* mesh{}; - for (const auto& prefix : {"pre_"s, "residuum_"s}) { - std::vector list_of_fields{}; - std::vector evolve_bndrys{}; - for (const auto& f : f3d) { - mesh = f.var->getMesh(); - Field3D to_load{0., mesh}; - to_load.allocate(); - to_load.setLocation(f.location); - debug[fmt::format("{:s}{:s}", prefix, f.name)] = to_load; - list_of_fields.push_back(to_load); - evolve_bndrys.push_back(f.evolve_bndry); + if (debug_on_failure) { + if (f2d.empty() and v2d.empty() and v3d.empty()) { + Options debug{}; + using namespace std::string_literals; + Mesh* mesh{}; + for (const auto& prefix : {"pre_"s, "residuum_"s}) { + std::vector list_of_fields{}; + std::vector evolve_bndrys{}; + for (const auto& f : f3d) { + mesh = f.var->getMesh(); + Field3D to_load{0., mesh}; + to_load.allocate(); + to_load.setLocation(f.location); + debug[fmt::format("{:s}{:s}", prefix, f.name)] = to_load; + list_of_fields.push_back(to_load); + evolve_bndrys.push_back(f.evolve_bndry); + } + pvode_load_data_f3d(evolve_bndrys, list_of_fields, + prefix == "pre_"s ? udata : N_VDATA(cv_mem->cv_acor)); } - pvode_load_data_f3d(evolve_bndrys, list_of_fields, - prefix == "pre_"s ? udata : N_VDATA(cv_mem->cv_acor)); - } - for (auto& f : f3d) { - f.F_var->enableTracking(fmt::format("ddt_{:s}", f.name), debug); - setName(*f.var, f.name); - } - run_rhs(simtime); + for (auto& f : f3d) { + f.F_var->enableTracking(fmt::format("ddt_{:s}", f.name), debug); + setName(*f.var, f.name); + } + run_rhs(simtime); - for (auto& f : f3d) { - debug[f.name] = *f.var; - } + for (auto& f : f3d) { + debug[f.name] = *f.var; + } - if (mesh != nullptr) { - mesh->outputVars(debug); - debug["BOUT_VERSION"].force(bout::version::as_double); - } + if (mesh != nullptr) { + mesh->outputVars(debug); + debug["BOUT_VERSION"].force(bout::version::as_double); + } - const std::string outname = fmt::format( - "{}/BOUT.debug.{}.nc", - Options::root()["datadir"].withDefault("data"), BoutComm::rank()); + const std::string outname = + fmt::format("{}/BOUT.debug.{}.nc", + Options::root()["datadir"].withDefault("data"), + BoutComm::rank()); - bout::OptionsIO::create(outname)->write(debug); - MPI_Barrier(BoutComm::get()); + bout::OptionsIO::create(outname)->write(debug); + MPI_Barrier(BoutComm::get()); + } else { + output_warn.write("debug_on_failure is currently only supported for Field3Ds"); + } } return (-1.0); } diff --git a/src/solver/impls/pvode/pvode.hxx b/src/solver/impls/pvode/pvode.hxx index d29135d02e..3b385af99d 100644 --- a/src/solver/impls/pvode/pvode.hxx +++ b/src/solver/impls/pvode/pvode.hxx @@ -75,10 +75,15 @@ private: pvode::machEnvType machEnv{nullptr}; void* cvode_mem{nullptr}; - BoutReal abstol, reltol; // addresses passed in init must be preserved + BoutReal abstol, reltol; + // addresses passed in init must be preserved pvode::PVBBDData pdata{nullptr}; - bool pvode_initialised = false; + /// is pvode already initialised? + bool pvode_initialised{false}; + + /// Add debugging data if solver fails + bool debug_on_failure{false}; }; #endif // BOUT_PVODE_SOLVER_H From db96b7ecc3df8613a80203d59a421bf4186f25d1 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 9 Apr 2024 15:56:25 +0200 Subject: [PATCH 19/30] Add option to euler solver to dump debug info --- src/solver/impls/euler/euler.cxx | 36 +++++++++++++++++++++++++++++++- src/solver/impls/euler/euler.hxx | 2 ++ 2 files changed, 37 insertions(+), 1 deletion(-) diff --git a/src/solver/impls/euler/euler.cxx b/src/solver/impls/euler/euler.cxx index 3976f4402c..45ba5ccdbf 100644 --- a/src/solver/impls/euler/euler.cxx +++ b/src/solver/impls/euler/euler.cxx @@ -20,7 +20,10 @@ EulerSolver::EulerSolver(Options* options) .withDefault(2.)), timestep((*options)["timestep"] .doc("Internal timestep (defaults to output timestep)") - .withDefault(getOutputTimestep())) {} + .withDefault(getOutputTimestep())), + dump_at_time((*options)["dump_at_time"] + .doc("Dump debug info about the simulation") + .withDefault(-1)) {} void EulerSolver::setMaxTimestep(BoutReal dt) { if (dt >= cfl_factor * timestep) { @@ -141,7 +144,38 @@ void EulerSolver::take_step(BoutReal curtime, BoutReal dt, Array& star Array& result) { load_vars(std::begin(start)); + const bool dump_now = dump_at_time > 0 && std::abs(dump_at_time - curtime) < dt; + std::unique_ptr debug_ptr; + if (dump_now) { + debug_ptr = std::make_unique(); + Options& debug = *debug_ptr; + for (auto& f : f3d) { + f.F_var->enableTracking(fmt::format("ddt_{:s}", f.name), debug); + setName(*f.var, f.name); + } + } + run_rhs(curtime); + if (dump_now) { + Options& debug = *debug_ptr; + Mesh* mesh; + for (auto& f : f3d) { + debug[f.name] = *f.var; + mesh = f.var->getMesh(); + } + + if (mesh != nullptr) { + mesh->outputVars(debug); + debug["BOUT_VERSION"].force(bout::version::as_double); + } + + const std::string outname = fmt::format( + "{}/BOUT.debug.{}.nc", + Options::root()["datadir"].withDefault("data"), BoutComm::rank()); + + bout::OptionsIO::create(outname)->write(debug); + MPI_Barrier(BoutComm::get()); + } save_derivs(std::begin(result)); BOUT_OMP_PERF(parallel for) diff --git a/src/solver/impls/euler/euler.hxx b/src/solver/impls/euler/euler.hxx index 0ee81a3d33..4b6dc60a62 100644 --- a/src/solver/impls/euler/euler.hxx +++ b/src/solver/impls/euler/euler.hxx @@ -64,6 +64,8 @@ private: /// Take a single step to calculate f1 void take_step(BoutReal curtime, BoutReal dt, Array& start, Array& result); + + BoutReal dump_at_time{-1.}; }; #endif // BOUT_KARNIADAKIS_SOLVER_H From 84bfcef2d10dc2dabf280c792a548519158414df Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 15 Apr 2024 11:25:59 +0200 Subject: [PATCH 20/30] Disable tracking once we are done --- include/bout/field3d.hxx | 7 +++++++ src/solver/impls/euler/euler.cxx | 5 ++++- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index cfde9e5328..97118cb70a 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -299,6 +299,13 @@ public: /// Save all changes that, are done to the field, to tracking Field3D& enableTracking(const std::string& name, Options& tracking); + /// Disable tracking + Field3D& disableTracking(){ + tracking = nullptr; + tracking_state = 0; + return *this; + } + ///////////////////////////////////////////////////////// // Data access diff --git a/src/solver/impls/euler/euler.cxx b/src/solver/impls/euler/euler.cxx index 45ba5ccdbf..227d7ba454 100644 --- a/src/solver/impls/euler/euler.cxx +++ b/src/solver/impls/euler/euler.cxx @@ -158,7 +158,7 @@ void EulerSolver::take_step(BoutReal curtime, BoutReal dt, Array& star run_rhs(curtime); if (dump_now) { Options& debug = *debug_ptr; - Mesh* mesh; + Mesh* mesh{nullptr}; for (auto& f : f3d) { debug[f.name] = *f.var; mesh = f.var->getMesh(); @@ -175,6 +175,9 @@ void EulerSolver::take_step(BoutReal curtime, BoutReal dt, Array& star bout::OptionsIO::create(outname)->write(debug); MPI_Barrier(BoutComm::get()); + for (auto& f : f3d) { + f.F_var->disableTracking(); + } } save_derivs(std::begin(result)); From 73265df7ebafd0a4d63763765977dd90cda2ce69 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 15 Apr 2024 11:27:07 +0200 Subject: [PATCH 21/30] Allow to dump every timestep with euler --- src/solver/impls/euler/euler.cxx | 14 ++++++++++---- src/solver/impls/euler/euler.hxx | 1 + 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/solver/impls/euler/euler.cxx b/src/solver/impls/euler/euler.cxx index 227d7ba454..ded28984aa 100644 --- a/src/solver/impls/euler/euler.cxx +++ b/src/solver/impls/euler/euler.cxx @@ -144,7 +144,8 @@ void EulerSolver::take_step(BoutReal curtime, BoutReal dt, Array& star Array& result) { load_vars(std::begin(start)); - const bool dump_now = dump_at_time > 0 && std::abs(dump_at_time - curtime) < dt; + const bool dump_now = + (dump_at_time >= 0 && std::abs(dump_at_time - curtime) < dt) || dump_at_time < -3; std::unique_ptr debug_ptr; if (dump_now) { debug_ptr = std::make_unique(); @@ -152,6 +153,8 @@ void EulerSolver::take_step(BoutReal curtime, BoutReal dt, Array& star for (auto& f : f3d) { f.F_var->enableTracking(fmt::format("ddt_{:s}", f.name), debug); setName(*f.var, f.name); + debug[fmt::format("pre_{:s}", f.name)] = *f.var; + f.var->allocate(); } } @@ -169,9 +172,12 @@ void EulerSolver::take_step(BoutReal curtime, BoutReal dt, Array& star debug["BOUT_VERSION"].force(bout::version::as_double); } - const std::string outname = fmt::format( - "{}/BOUT.debug.{}.nc", - Options::root()["datadir"].withDefault("data"), BoutComm::rank()); + const std::string outnumber = + dump_at_time < -3 ? fmt::format(".{}", debug_counter++) : ""; + const std::string outname = + fmt::format("{}/BOUT.debug{}.{}.nc", + Options::root()["datadir"].withDefault("data"), + outnumber, BoutComm::rank()); bout::OptionsIO::create(outname)->write(debug); MPI_Barrier(BoutComm::get()); diff --git a/src/solver/impls/euler/euler.hxx b/src/solver/impls/euler/euler.hxx index 4b6dc60a62..fc9b7f53bb 100644 --- a/src/solver/impls/euler/euler.hxx +++ b/src/solver/impls/euler/euler.hxx @@ -66,6 +66,7 @@ private: Array& result); BoutReal dump_at_time{-1.}; + int debug_counter{0}; }; #endif // BOUT_KARNIADAKIS_SOLVER_H From 8ff388aeb9f74c23fb6cc8a073565cc19276d268 Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Mon, 15 Apr 2024 09:30:12 +0000 Subject: [PATCH 22/30] Apply clang-format changes --- include/bout/field3d.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 97118cb70a..395e3201e4 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -300,7 +300,7 @@ public: Field3D& enableTracking(const std::string& name, Options& tracking); /// Disable tracking - Field3D& disableTracking(){ + Field3D& disableTracking() { tracking = nullptr; tracking_state = 0; return *this; From 6724e7548bf88a7824845a9caf1c77c49bef269b Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 9 Aug 2024 12:02:43 +0200 Subject: [PATCH 23/30] Dump also parallel fields by default --- src/solver/impls/pvode/pvode.cxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index 9dce5d357f..66344f7cde 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -358,8 +358,8 @@ BoutReal PvodeSolver::run(BoutReal tout) { // Check return flag if (flag != SUCCESS) { output_error.write("ERROR CVODE step failed, flag = {:d}\n", flag); - CVodeMemRec* cv_mem = (CVodeMem)cvode_mem; if (debug_on_failure) { + CVodeMemRec* cv_mem = (CVodeMem)cvode_mem; if (f2d.empty() and v2d.empty() and v3d.empty()) { Options debug{}; using namespace std::string_literals; @@ -388,6 +388,9 @@ BoutReal PvodeSolver::run(BoutReal tout) { for (auto& f : f3d) { debug[f.name] = *f.var; + if (f.var->hasParallelSlices()) { + saveParallel(debug, f.name, *f.var); + } } if (mesh != nullptr) { From 28212c2e97c725f20099fd302fb17efc86243c7b Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 9 Aug 2024 12:01:16 +0200 Subject: [PATCH 24/30] Also dump parallel fields --- src/solver/impls/euler/euler.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solver/impls/euler/euler.cxx b/src/solver/impls/euler/euler.cxx index ded28984aa..5477b5760b 100644 --- a/src/solver/impls/euler/euler.cxx +++ b/src/solver/impls/euler/euler.cxx @@ -163,7 +163,7 @@ void EulerSolver::take_step(BoutReal curtime, BoutReal dt, Array& star Options& debug = *debug_ptr; Mesh* mesh{nullptr}; for (auto& f : f3d) { - debug[f.name] = *f.var; + saveParallel(debug, f.name, *f.var); mesh = f.var->getMesh(); } From 9cf0fba6971c1c6367b3aee3c5d68f9bcf4121ed Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 9 Aug 2024 11:53:26 +0200 Subject: [PATCH 25/30] Clarify which div_par has been used --- src/mesh/coordinates.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 32774d6229..483b3e458c 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1601,7 +1601,7 @@ Field3D Coordinates::Div_par(const Field3D& f, CELL_LOC outloc, f_B.yup(i) = f.yup(i) / Bxy_floc.yup(i); f_B.ydown(i) = f.ydown(i) / Bxy_floc.ydown(i); } - return setName(Bxy * Grad_par(f_B, outloc, method), "Div_par({:s})", f.name); + return setName(Bxy * Grad_par(f_B, outloc, method), "C:Div_par({:s})", f.name); } ///////////////////////////////////////////////////////// From 97b67e187149d5c7e6624ffd1222ef8d242d5270 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 9 Aug 2024 10:30:32 +0200 Subject: [PATCH 26/30] Expose tracking Useful for physics module to record additional data if the solver is failing --- include/bout/field3d.hxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 395e3201e4..0f227a8da3 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -510,6 +510,8 @@ public: int size() const override { return nx * ny * nz; }; + Options* getTracking() { return tracking; }; + private: /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null int nx{-1}, ny{-1}, nz{-1}; From 77e08ec6f7e4e38978e1ec85218af8475991ce4b Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Tue, 22 Oct 2024 09:57:45 +0000 Subject: [PATCH 27/30] Apply clang-format changes --- src/solver/impls/pvode/pvode.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index 66344f7cde..d4ac14fef2 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -388,9 +388,9 @@ BoutReal PvodeSolver::run(BoutReal tout) { for (auto& f : f3d) { debug[f.name] = *f.var; - if (f.var->hasParallelSlices()) { - saveParallel(debug, f.name, *f.var); - } + if (f.var->hasParallelSlices()) { + saveParallel(debug, f.name, *f.var); + } } if (mesh != nullptr) { From ba6fc6cea816d0c2163dcbd3af4f64824a25f446 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 9 Aug 2024 10:32:17 +0200 Subject: [PATCH 28/30] Add simple interface to store parallel fields Just dumping the parallel slices does in general not work, as then collect discards that, especially if NYPE==ny --- include/bout/options.hxx | 3 +++ src/sys/options.cxx | 16 ++++++++++++++++ 2 files changed, 19 insertions(+) diff --git a/include/bout/options.hxx b/include/bout/options.hxx index 839c847289..e1f5ae68fa 100644 --- a/include/bout/options.hxx +++ b/include/bout/options.hxx @@ -946,6 +946,9 @@ Tensor Options::as>(const Tensor& similar_t /// Convert \p value to string std::string toString(const Options& value); +/// Save the parallel fields +void saveParallel(Options& opt, const std::string name, const Field3D& tosave); + /// Output a stringified \p value to a stream /// /// This is templated to avoid implict casting: anything is diff --git a/src/sys/options.cxx b/src/sys/options.cxx index 49a81cfa88..71339b6089 100644 --- a/src/sys/options.cxx +++ b/src/sys/options.cxx @@ -306,6 +306,22 @@ void Options::assign<>(Tensor val, std::string source) { _set_no_check(std::move(val), std::move(source)); } +void saveParallel(Options& opt, const std::string name, const Field3D& tosave){ + ASSERT2(tosave.hasParallelSlices()); + opt[name] = tosave; + for (size_t i0=1 ; i0 <= tosave.numberParallelSlices(); ++i0) { + for (int i: {i0, -i0} ) { + Field3D tmp; + tmp.allocate(); + const auto& fpar = tosave.ynext(i); + for (auto j: fpar.getValidRegionWithDefault("RGN_NO_BOUNDARY")){ + tmp[j.yp(-i)] = fpar[j]; + } + opt[fmt::format("{}_y{:+d}", name, i)] = tmp; + } + } +} + namespace { /// Use FieldFactory to evaluate expression double parseExpression(const Options::ValueType& value, const Options* options, From 3c0d6a8f6a4d20bdee9d140b9f90af9432a6810c Mon Sep 17 00:00:00 2001 From: dschwoerer Date: Tue, 22 Oct 2024 10:26:41 +0000 Subject: [PATCH 29/30] Apply clang-format changes --- src/sys/options.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/sys/options.cxx b/src/sys/options.cxx index 71339b6089..1193deb013 100644 --- a/src/sys/options.cxx +++ b/src/sys/options.cxx @@ -306,16 +306,16 @@ void Options::assign<>(Tensor val, std::string source) { _set_no_check(std::move(val), std::move(source)); } -void saveParallel(Options& opt, const std::string name, const Field3D& tosave){ +void saveParallel(Options& opt, const std::string name, const Field3D& tosave) { ASSERT2(tosave.hasParallelSlices()); opt[name] = tosave; - for (size_t i0=1 ; i0 <= tosave.numberParallelSlices(); ++i0) { - for (int i: {i0, -i0} ) { + for (size_t i0 = 1; i0 <= tosave.numberParallelSlices(); ++i0) { + for (int i : {i0, -i0}) { Field3D tmp; tmp.allocate(); const auto& fpar = tosave.ynext(i); - for (auto j: fpar.getValidRegionWithDefault("RGN_NO_BOUNDARY")){ - tmp[j.yp(-i)] = fpar[j]; + for (auto j : fpar.getValidRegionWithDefault("RGN_NO_BOUNDARY")) { + tmp[j.yp(-i)] = fpar[j]; } opt[fmt::format("{}_y{:+d}", name, i)] = tmp; } From d0669cdcafa76bec00116abb2d934e05405eae35 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 22 Oct 2024 15:37:56 +0200 Subject: [PATCH 30/30] Do not use numberParallelSlices() --- src/sys/options.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/sys/options.cxx b/src/sys/options.cxx index df1dee56f4..3f7a0c4071 100644 --- a/src/sys/options.cxx +++ b/src/sys/options.cxx @@ -338,9 +338,10 @@ void Options::assign<>(Tensor val, std::string source) { } void saveParallel(Options& opt, const std::string name, const Field3D& tosave) { - ASSERT2(tosave.hasParallelSlices()); opt[name] = tosave; - for (size_t i0 = 1; i0 <= tosave.numberParallelSlices(); ++i0) { + const size_t numberParallelSlices = + tosave.hasParallelSlices() ? 0 : tosave.getMesh()->ystart; + for (size_t i0 = 1; i0 <= numberParallelSlices; ++i0) { for (int i : {i0, -i0}) { Field3D tmp; tmp.allocate();