diff --git a/include/bout/field.hxx b/include/bout/field.hxx index 61edff6723..2d813bb275 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -36,6 +36,7 @@ class Field; #include "bout/bout_types.hxx" #include "bout/boutcomm.hxx" #include "bout/boutexception.hxx" +#include "bout/build_defines.hxx" #include "bout/field_data.hxx" #include "bout/msg_stack.hxx" #include "bout/region.hxx" @@ -678,4 +679,19 @@ 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 + f.name = fmt::format(name, args...); +#endif + return std::forward(f); +} + #endif /* FIELD_H */ diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index a75e38df36..fad4815b92 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -33,11 +33,15 @@ class Field3D; #include "bout/field2d.hxx" #include "bout/fieldperp.hxx" #include "bout/region.hxx" +#include "bout/traits.hxx" +#include #include +#include #include class Mesh; +class Options; /// Class for 3D X-Y-Z scalar fields /*! @@ -291,6 +295,17 @@ 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, std::weak_ptr tracking); + + /// Disable tracking + Field3D& disableTracking() { + tracking.reset(); + tracking_state = 0; + return *this; + } + ///////////////////////////////////////////////////////// // Data access @@ -493,6 +508,8 @@ public: int size() const override { return nx * ny * nz; }; + std::weak_ptr getTracking() { return tracking; }; + private: /// Array sizes (from fieldmesh). These are valid only if fieldmesh is not null int nx{-1}, ny{-1}, nz{-1}; @@ -508,6 +525,22 @@ private: /// RegionID over which the field is valid std::optional regionID; + + /// counter for tracking, to assign unique names to the variable names + int tracking_state{0}; + std::weak_ptr tracking; + // name is changed if we assign to the variable, while selfname is a + // non-changing copy that is used for the variable names in the dump files + std::string selfname; + template + void track(const T& change, const std::string& operation) { + if (tracking_state != 0) { + _track(change, operation); + } + } + template > + void _track(const T& change, std::string operation); + void _track(const BoutReal& change, std::string operation); }; // Non-member overloaded operators diff --git a/include/bout/options.hxx b/include/bout/options.hxx index 1c552cd4f9..a731f73692 100644 --- a/include/bout/options.hxx +++ b/include/bout/options.hxx @@ -1012,6 +1012,9 @@ auto Options::as(const Tensor& similar_to) const -> Tensor; /// 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/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 diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 4c5ffa9748..7ece06a454 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -31,6 +31,8 @@ #include #include +#include +#include #include "bout/parallel_boundary_op.hxx" #include "bout/parallel_boundary_region.hxx" @@ -47,6 +49,8 @@ #include #include +#include "fmt/format.h" + /// Constructor Field3D::Field3D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in) : Field(localmesh, location_in, directions_in) { @@ -246,6 +250,7 @@ Field3D& Field3D::operator=(const Field3D& rhs) { } TRACE("Field3D: Assignment from Field3D"); + track(rhs, "operator="); // Copy base slice Field::operator=(rhs); @@ -267,6 +272,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); @@ -288,6 +294,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()); @@ -334,6 +341,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. @@ -864,3 +872,57 @@ 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, + std::weak_ptr _tracking) { + tracking = std::move(_tracking); + tracking_state = 1; + selfname = name; + return *this; +} + +template +void Field3D::_track(const T& change, std::string operation) { + if (tracking_state == 0) { + return; + } + auto locked = tracking.lock(); + if (locked == nullptr) { + return; + } + const std::string outname{fmt::format("track_{:s}_{:d}", selfname, tracking_state++)}; + + locked->set(outname, change, "tracking"); + // Workaround for bug in gcc9.4 +#if BOUT_USE_TRACK + const std::string changename = change.name; +#endif + (*locked)[outname].setAttributes({ + {"operation", operation}, +#if BOUT_USE_TRACK + {"rhs.name", changename}, +#endif + }); +} + +template void +Field3D::_track>(const Field3D&, + std::string); +template void Field3D::_track(const Field2D&, std::string); +template void Field3D::_track<>(const FieldPerp&, std::string); + +void Field3D::_track(const BoutReal& change, std::string operation) { + if (tracking_state == 0) { + return; + } + auto locked = tracking.lock(); + if (locked == nullptr) { + return; + } + const std::string outname{fmt::format("track_{:s}_{:d}", selfname, tracking_state++)}; + locked->set(outname, change, "tracking"); + (*locked)[outname].setAttributes({ + {"operation", operation}, + {"rhs.name", "BoutReal"}, + }); +} 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/gen_fieldops.py b/src/field/gen_fieldops.py index 29631ff7aa..b6f99549d2 100755 --- a/src/field/gen_fieldops.py +++ b/src/field/gen_fieldops.py @@ -65,12 +65,15 @@ def smart_open(filename, mode="r"): ) header = """// This file is autogenerated - see gen_fieldops.py -#include -#include -#include -#include -#include -#include +#include "bout/build_defines.hxx" +#include "bout/field2d.hxx" +#include "bout/field3d.hxx" +#include "bout/globals.hxx" +#include "bout/interpolation.hxx" +#include "bout/mesh.hxx" +#include "bout/region.hxx" + +#include "fmt/format.h" """ diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx index 6b778acee3..f37a870c74 100644 --- a/src/field/generated_fieldops.cxx +++ b/src/field/generated_fieldops.cxx @@ -1,10 +1,13 @@ // This file is autogenerated - see gen_fieldops.py -#include -#include -#include -#include -#include -#include +#include "bout/build_defines.hxx" +#include "bout/field2d.hxx" +#include "bout/field3d.hxx" +#include "bout/globals.hxx" +#include "bout/interpolation.hxx" +#include "bout/mesh.hxx" +#include "bout/region.hxx" + +#include "fmt/format.h" // Provide the C++ wrapper for multiplication of Field3D and Field3D Field3D operator*(const Field3D& lhs, const Field3D& rhs) { @@ -20,6 +23,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 +48,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 +76,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 +101,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 +129,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 +154,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 +182,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 +207,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 +240,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 +268,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 +302,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 +331,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 +364,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 +392,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 +425,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 +453,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 +483,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 +506,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 +529,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 +552,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 +572,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 +594,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 +622,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 +645,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 +672,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 +694,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 +721,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 +743,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 +776,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 +802,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 +828,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 +854,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 +873,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 +892,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 +916,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 +935,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 +959,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 +978,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 +1002,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 +1021,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 +1049,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 +1072,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 +1095,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 +1118,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 +1136,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 +1154,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 +1178,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 +1197,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 +1220,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 +1238,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 +1261,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 +1279,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 +1307,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 +1332,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 +1360,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 +1385,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 +1413,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 +1438,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 +1466,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 +1491,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 +1519,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 +1544,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 +1572,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 +1597,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 +1625,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 +1650,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 +1678,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 +1703,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 +1727,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 +1746,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 +1770,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 +1789,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 +1813,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 +1832,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 +1856,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 +1875,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 +1898,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 +1916,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 +1940,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 +1958,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 +1981,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 +1999,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 +2022,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 +2040,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 +2065,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 +2085,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 +2105,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 +2125,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 +2143,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 +2161,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 +2179,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 +2197,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 +2215,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 +2233,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 +2251,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 +2269,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; } diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 34c524d1e7..e805df00f2 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -1544,7 +1545,7 @@ 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(); + return setName(::DDY(var, outloc, method) * invSg(), "Grad_par({:s})", var.name); } ///////////////////////////////////////////////////////// @@ -1603,7 +1604,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), "C:Div_par({:s})", f.name); } ///////////////////////////////////////////////////////// diff --git a/src/physics/smoothing.cxx b/src/physics/smoothing.cxx index 0a1391907f..036c475cd7 100644 --- a/src/physics/smoothing.cxx +++ b/src/physics/smoothing.cxx @@ -350,7 +350,7 @@ BoutReal Average_XY(const Field2D& var) { return Vol_Glb; } -BoutReal Vol_Integral(const Field2D& var) { +BoutReal Vol_Integral([[maybe_unused]] const Field2D& var) { #if BOUT_USE_METRIC_3D AUTO_TRACE(); throw BoutException("Vol_Intregral currently incompatible with 3D metrics"); diff --git a/src/solver/impls/euler/euler.cxx b/src/solver/impls/euler/euler.cxx index 3976f4402c..a4804590dc 100644 --- a/src/solver/impls/euler/euler.cxx +++ b/src/solver/impls/euler/euler.cxx @@ -1,13 +1,20 @@ #include "euler.hxx" +#include "bout/field.hxx" +#include "bout/version.hxx" #include #include +#include #include #include +#include #include +#include + #include +#include #include @@ -20,7 +27,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 +151,47 @@ 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) || dump_at_time < -3; + std::shared_ptr debug_ptr; + if (dump_now) { + debug_ptr = std::make_shared(); + Options& debug = *debug_ptr; + for (auto& f : f3d) { + f.F_var->enableTracking(fmt::format("ddt_{:s}", f.name), debug_ptr); + setName(*f.var, f.name); + debug[fmt::format("pre_{:s}", f.name)] = *f.var; + f.var->allocate(); + } + } + run_rhs(curtime); + if (dump_now) { + Options& debug = *debug_ptr; + Mesh* mesh{nullptr}; + for (auto& f : f3d) { + saveParallel(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 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()); + for (auto& f : f3d) { + f.F_var->disableTracking(); + } + } 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..fc9b7f53bb 100644 --- a/src/solver/impls/euler/euler.hxx +++ b/src/solver/impls/euler/euler.hxx @@ -64,6 +64,9 @@ private: /// Take a single step to calculate f1 void take_step(BoutReal curtime, BoutReal dt, Array& start, Array& result); + + BoutReal dump_at_time{-1.}; + int debug_counter{0}; }; #endif // BOUT_KARNIADAKIS_SOLVER_H diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx index 63d3c11753..2b83d4f155 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -31,23 +31,60 @@ #include #include +#include +#include #include #include +#include #include #include +#include +#include -#include "bout/unused.hxx" - +#include #include // use CVSPGMR linear solver each internal step #include // contains the enum for types of preconditioning #include // band preconditioner function prototypes +#include "fmt/format.h" + +#include +#include +#include +#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, const 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 + auto evolve_bndry = evolve_bndrys.cbegin(); + for (auto 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]; @@ -185,10 +222,47 @@ 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; - cvode_mem = CVodeMalloc(neq, solver_f, simtime, u, BDF, NEWTON, SS, &reltol, &abstol, - this, nullptr, optIn, iopt, ropt, machEnv); + { + /* 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)); + + 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); if (cvode_mem == nullptr) { throw BoutException("\tError: CVodeMalloc failed.\n"); @@ -293,7 +367,55 @@ BoutReal PvodeSolver::run(BoutReal tout) { // Check return flag if (flag != SUCCESS) { output_error.write("ERROR CVODE step failed, flag = {:d}\n", flag); - return (-1.0); + if (debug_on_failure) { + CVodeMemRec* cv_mem = static_cast(cvode_mem); + if (!(f2d.empty() and v2d.empty() and v3d.empty())) { + output_warn.write("debug_on_failure is currently only supported for Field3Ds"); + return -1.0; + } + auto debug_ptr = std::make_shared(); + Options& debug = *debug_ptr; + 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)); + } + + for (auto& f : f3d) { + f.F_var->enableTracking(fmt::format("ddt_{:s}", f.name), debug_ptr); + setName(*f.var, f.name); + } + run_rhs(simtime); + + for (auto& f : f3d) { + saveParallel(debug, f.name, *f.var); + } + + 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()); + } + return -1.0; } return simtime; diff --git a/src/solver/impls/pvode/pvode.hxx b/src/solver/impls/pvode/pvode.hxx index 6425fc1868..cf68444b6c 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 diff --git a/src/sys/options.cxx b/src/sys/options.cxx index 077b2a6d85..ada63bd526 100644 --- a/src/sys/options.cxx +++ b/src/sys/options.cxx @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -358,6 +359,23 @@ Options& Options::assign<>(Tensor val, std::string source) { return *this; } +void saveParallel(Options& opt, const std::string& name, const Field3D& tosave) { + opt[name] = tosave; + 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(); + 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,