diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index 010100588c..2f273fc94b 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -5,6 +5,8 @@ ed2117e6d6826a98b6988e2f18c0c34e408563b6 17ac13c28aa3b34a0e46dbe87bb3874f6b25e706 # Added by the bot 4b010b7634aee1045743be80c268d4644522cd29 +52301380586fdbf890f620c04f689b08d89a6c34 a71cad2dd6ace5741a754e2ca7daacd4bb094e0e +83cf77923a4c72e44303354923021acf932b4fd2 2c2402ed59c91164eaff46dee0f79386b7347e9e 05b7c571544c3bcb153fce67d12b9ac48947fc2d 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/include/bout/options_io.hxx b/include/bout/options_io.hxx index 509305aad0..00dd0c4cf3 100644 --- a/include/bout/options_io.hxx +++ b/include/bout/options_io.hxx @@ -5,8 +5,8 @@ /// /// 1. Dump files, containing time history: /// -/// auto dump = OptionsIOFactory::getInstance().createOutput(); -/// dump->write(data); +/// auto dump = OptionsIOFactory::getInstance().createOutput(); +/// dump->write(data); /// /// where data is an Options tree. By default dump files are configured /// with the root `output` section, or an Option tree can be passed to @@ -14,20 +14,27 @@ /// /// 2. Restart files: /// -/// auto restart = OptionsIOFactory::getInstance().createOutput(); -/// restart->write(data); +/// auto restart = OptionsIOFactory::getInstance().createRestart(); +/// restart->write(data); /// /// where data is an Options tree. By default restart files are configured /// with the root `restart_files` section, or an Option tree can be passed to /// `createRestart`. /// /// 3. Ad-hoc single files -/// Note: The caller should consider how multiple processors interact with the file. +/// Note: The caller should consider how multiple processors interact with the file. /// -/// auto file = OptionsIOFactory::getInstance().createFile("some_file.nc"); -/// or -/// auto file = OptionsIO::create("some_file.nc"); +/// auto file = OptionsIOFactory::getInstance().createFile("some_file.nc"); +/// or +/// auto file = OptionsIO::create("some_file.nc"); /// +/// 4. Ad-hoc parallel files +/// This adds also metric information, such that the file can be read with +/// boutdata or xBOUT +/// +/// OptionIO::write("some_file", data, mesh); +/// +/// if mesh is omitted, no grid information is added. /// #pragma once @@ -37,6 +44,7 @@ #include "bout/build_defines.hxx" #include "bout/generic_factory.hxx" +#include "bout/mesh.hxx" #include "bout/options.hxx" #include @@ -77,6 +85,11 @@ public: /// This uses the default file type and default options. static std::unique_ptr create(const std::string& file); + /// Write some data to a file with a given name prefix + /// This will be done in parallel. If Mesh is given, also mesh data will be + /// added, which is needed for xBOUT or boutdata to read the files. + static void write(const std::string& prefix, Options& data, Mesh* mesh = nullptr); + /// Create an OptionsIO for I/O to the given file. /// The file will be configured using the given `config` options: /// - "type" : string The file type e.g. "netcdf" or "adios" diff --git a/manual/sphinx/developer_docs/debugging.rst b/manual/sphinx/developer_docs/debugging.rst index f16d0b6a3b..4f47cd63e3 100644 --- a/manual/sphinx/developer_docs/debugging.rst +++ b/manual/sphinx/developer_docs/debugging.rst @@ -158,3 +158,30 @@ If you need to capture runtime information in the message, you can use the ``fmt`` syntax also used by the loggers:: TRACE("Value of i={}, some arbitrary {}", i, "string"); + + +Time evolution +============== + +It can be useful to know what happened when the simulation failed. The pvode +solver can dump the state of the simulation, at the time the solver +failed. This information includes the individual terms in the derivative. This +allows to identify which term is causing the issue. Additionally, the +residuum is dumped. This identifies not only which term is causing the issue, +but also where in the domain the solver is struggling. This can be enabled +with:: + + solver:type=pvode solver:debug_on_failure=true + +It can be also useful for understanding why the solver is slow. Forcing a +higher min_timestep, the solver will fail to evolve the system as it +encounters the situation, and provides information where it is happening. +This can be done with:: + + solver:type=pvode solver:debug_on_failure=true solver:min_timestep=1e2 + +It is also possible to dump at a specific time using the euler solver. +This can be useful for tracking down what is causing differences between two +different versions. It can be used with:: + + solver:type=euler solver:dump_at_time=0 input:error_on_unused_options=false diff --git a/manual/sphinx/user_docs/bout_options.rst b/manual/sphinx/user_docs/bout_options.rst index 926c201c16..4a9c568e89 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 b3448238f5..2cd55b0957 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -31,6 +31,9 @@ #include #include +#include +#include +#include #include "bout/parallel_boundary_op.hxx" #include "bout/parallel_boundary_region.hxx" @@ -46,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) { @@ -235,6 +240,8 @@ Field3D& Field3D::operator=(const Field3D& rhs) { return (*this); // skip this assignment } + track(rhs, "operator="); + // Copy base slice Field::operator=(rhs); @@ -254,6 +261,7 @@ Field3D& Field3D::operator=(const Field3D& rhs) { } Field3D& Field3D::operator=(Field3D&& rhs) { + track(rhs, "operator="); // Move parallel slices or delete existing ones. yup_fields = std::move(rhs.yup_fields); @@ -274,6 +282,7 @@ Field3D& Field3D::operator=(Field3D&& rhs) { } Field3D& Field3D::operator=(const Field2D& rhs) { + track(rhs, "operator="); /// Check that the data is allocated ASSERT1(rhs.isAllocated()); @@ -318,6 +327,7 @@ void Field3D::operator=(const FieldPerp& rhs) { } Field3D& Field3D::operator=(const BoutReal val) { + track(val, "operator="); // Delete existing parallel slices. We don't copy parallel slices, so any // that currently exist will be incorrect. @@ -831,3 +841,63 @@ 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"); + + const std::string trace = cpptrace::generate_trace().to_string(); + + // 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 + {"trace", trace}, + }); +} + +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 trace = cpptrace::generate_trace().to_string(); + 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"}, + {"trace", trace}, + }); +} diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja index ecd4e628cc..5476242d50 100644 --- a/src/field/gen_fieldops.jinja +++ b/src/field/gen_fieldops.jinja @@ -129,9 +129,16 @@ } {% endif %} + {% if lhs == "Field3D" %} + track(rhs, "operator{{operator}}="); + {% 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..24bfde425c 100644 --- a/src/field/generated_fieldops.cxx +++ b/src/field/generated_fieldops.cxx @@ -42,9 +42,12 @@ Field3D& Field3D::operator*=(const Field3D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs[index]; } + track(rhs, "operator*="); + checkData(*this); } else { + track(rhs, "operator*="); (*this) = (*this) * rhs; } return *this; @@ -86,9 +89,12 @@ Field3D& Field3D::operator/=(const Field3D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] /= rhs[index]; } + track(rhs, "operator/="); + checkData(*this); } else { + track(rhs, "operator/="); (*this) = (*this) / rhs; } return *this; @@ -130,9 +136,12 @@ Field3D& Field3D::operator+=(const Field3D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs[index]; } + track(rhs, "operator+="); + checkData(*this); } else { + track(rhs, "operator+="); (*this) = (*this) + rhs; } return *this; @@ -174,9 +183,12 @@ Field3D& Field3D::operator-=(const Field3D& rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs[index]; } + track(rhs, "operator-="); + checkData(*this); } else { + track(rhs, "operator-="); (*this) = (*this) - rhs; } return *this; @@ -226,9 +238,12 @@ Field3D& Field3D::operator*=(const Field2D& rhs) { } } + track(rhs, "operator*="); + checkData(*this); } else { + track(rhs, "operator*="); (*this) = (*this) * rhs; } return *this; @@ -280,9 +295,12 @@ Field3D& Field3D::operator/=(const Field2D& rhs) { } } + track(rhs, "operator/="); + checkData(*this); } else { + track(rhs, "operator/="); (*this) = (*this) / rhs; } return *this; @@ -332,9 +350,12 @@ Field3D& Field3D::operator+=(const Field2D& rhs) { } } + track(rhs, "operator+="); + checkData(*this); } else { + track(rhs, "operator+="); (*this) = (*this) + rhs; } return *this; @@ -384,9 +405,12 @@ Field3D& Field3D::operator-=(const Field2D& rhs) { } } + track(rhs, "operator-="); + checkData(*this); } else { + track(rhs, "operator-="); (*this) = (*this) - rhs; } return *this; @@ -504,9 +528,12 @@ Field3D& Field3D::operator*=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] *= rhs; } + track(rhs, "operator*="); + checkData(*this); } else { + track(rhs, "operator*="); (*this) = (*this) * rhs; } return *this; @@ -546,9 +573,12 @@ 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/="); + checkData(*this); } else { + track(rhs, "operator/="); (*this) = (*this) / rhs; } return *this; @@ -586,9 +616,12 @@ Field3D& Field3D::operator+=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] += rhs; } + track(rhs, "operator+="); + checkData(*this); } else { + track(rhs, "operator+="); (*this) = (*this) + rhs; } return *this; @@ -626,9 +659,12 @@ Field3D& Field3D::operator-=(const BoutReal rhs) { BOUT_FOR(index, this->getRegion("RGN_ALL")) { (*this)[index] -= rhs; } + track(rhs, "operator-="); + checkData(*this); } else { + track(rhs, "operator-="); (*this) = (*this) - rhs; } return *this; diff --git a/src/physics/smoothing.cxx b/src/physics/smoothing.cxx index 309fff1542..1b437b4352 100644 --- a/src/physics/smoothing.cxx +++ b/src/physics/smoothing.cxx @@ -340,7 +340,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 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 a3911e78d6..9754a68826 100644 --- a/src/solver/impls/euler/euler.cxx +++ b/src/solver/impls/euler/euler.cxx @@ -1,11 +1,16 @@ - #include "euler.hxx" #include #include +#include #include +#include #include +#include +#include +#include + EulerSolver::EulerSolver(Options* options) : Solver(options), mxstep((*options)["mxstep"] .doc("Maximum number of steps between outputs") @@ -15,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) { @@ -134,7 +142,37 @@ 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); + 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(); + } + + const std::string outnumber = + dump_at_time < -3 ? fmt::format(".{}", debug_counter++) : ""; + const std::string outname = fmt::format("BOUT.debug{}", outnumber); + bout::OptionsIO::write(outname, debug, mesh); + 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 cb6e81d80a..8ccd5f134a 100644 --- a/src/solver/impls/pvode/pvode.cxx +++ b/src/solver/impls/pvode/pvode.cxx @@ -31,23 +31,59 @@ #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; + const 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]; @@ -184,10 +220,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"); @@ -291,7 +364,43 @@ 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) { + const 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.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); + } + run_rhs(simtime); + + for (auto& f : f3d) { + saveParallel(debug, f.name, *f.var); + } + bout::OptionsIO::write("BOUT.debug", debug, mesh); + MPI_Barrier(BoutComm::get()); + } + return -1.0; } return simtime; @@ -301,7 +410,8 @@ BoutReal PvodeSolver::run(BoutReal tout) { * RHS function **************************************************************************/ -void PvodeSolver::rhs(int UNUSED(N), BoutReal t, BoutReal* udata, BoutReal* dudata) { +void PvodeSolver::rhs([[maybe_unused]] int N, BoutReal t, BoutReal* udata, + BoutReal* dudata) { TRACE("Running RHS: PvodeSolver::rhs({})", t); // Get current timestep @@ -317,7 +427,8 @@ void PvodeSolver::rhs(int UNUSED(N), BoutReal t, BoutReal* udata, BoutReal* duda save_derivs(dudata); } -void PvodeSolver::gloc(int UNUSED(N), BoutReal t, BoutReal* udata, BoutReal* dudata) { +void PvodeSolver::gloc([[maybe_unused]] int N, BoutReal t, BoutReal* udata, + BoutReal* dudata) { TRACE("Running RHS: PvodeSolver::gloc({})", t); Timer timer("rhs"); @@ -358,8 +469,8 @@ void solver_gloc(integer N, BoutReal t, BoutReal* u, BoutReal* udot, void* f_dat } // Preconditioner communication function -void solver_cfn(integer UNUSED(N), BoutReal UNUSED(t), N_Vector UNUSED(u), - void* UNUSED(f_data)) { +void solver_cfn([[maybe_unused]] integer N, [[maybe_unused]] BoutReal t, + [[maybe_unused]] N_Vector u, [[maybe_unused]] void* f_data) { // doesn't do anything at the moment } 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 bb4c920b90..0a68bab124 100644 --- a/src/sys/options.cxx +++ b/src/sys/options.cxx @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -354,6 +355,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, diff --git a/src/sys/options/options_io.cxx b/src/sys/options/options_io.cxx index 6717b6b07d..256ec56b28 100644 --- a/src/sys/options/options_io.cxx +++ b/src/sys/options/options_io.cxx @@ -2,6 +2,7 @@ #include "bout/bout.hxx" #include "bout/globals.hxx" #include "bout/mesh.hxx" +#include "bout/version.hxx" #include "options_adios.hxx" #include "options_netcdf.hxx" @@ -55,4 +56,12 @@ void writeDefaultOutputFile(Options& data) { OptionsIOFactory::getInstance().createOutput()->write(data); } +void OptionsIO::write(const std::string& prefix, Options& data, Mesh* mesh) { + Options file_options = {{"prefix", prefix}}; + data["BOUT_VERSION"].force(bout::version::as_double); + if (mesh != nullptr) { + mesh->outputVars(data); + } + OptionsIOFactory::getInstance().createOutput(&file_options)->write(data); +} } // namespace bout