Merge pull request #2889 from boutproject/pvcode-cvode-improvements

Track field operations for debugging

Author: dschwoerer

.git-blame-ignore-revs

@@ -5,6 +5,8 @@ ed2117e6d6826a98b6988e2f18c0c34e408563b6
 17ac13c28aa3b34a0e46dbe87bb3874f6b25e706
 # Added by the bot
 4b010b7634aee1045743be80c268d4644522cd29
+52301380586fdbf890f620c04f689b08d89a6c34
 a71cad2dd6ace5741a754e2ca7daacd4bb094e0e
+83cf77923a4c72e44303354923021acf932b4fd2
 2c2402ed59c91164eaff46dee0f79386b7347e9e
 05b7c571544c3bcb153fce67d12b9ac48947fc2d

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 <memory>
 #include <optional>
+#include <string>
 #include <vector>
 
 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<Options> 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<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};
@@ -508,6 +525,22 @@ private:
 
   /// RegionID over which the field is valid
   std::optional<size_t> regionID;
+
+  /// counter for tracking, to assign unique names to the variable names
+  int tracking_state{0};
+  std::weak_ptr<Options> 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 <typename T>
+  void track(const T& change, const std::string& operation) {
+    if (tracking_state != 0) {
+      _track(change, operation);
+    }
+  }
+  template <typename T, typename = bout::utils::EnableIfField<T>>
+  void _track(const T& change, std::string operation);
+  void _track(const BoutReal& change, std::string operation);
 };
 
 // Non-member overloaded operators

include/bout/options.hxx

@@ -1012,6 +1012,9 @@ auto Options::as(const Tensor<int>& similar_to) const -> Tensor<int>;
 /// 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

include/bout/options_io.hxx

@@ -5,29 +5,36 @@
 ///
 /// 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
 ///    `createOutput`.
 ///
 /// 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 <memory>
@@ -77,6 +85,11 @@ public:
   /// This uses the default file type and default options.
   static std::unique_ptr<OptionsIO> 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"

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

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

src/field/field3d.cxx

@@ -31,6 +31,9 @@
 #include <bout/globals.hxx>
 
 #include <cmath>
+#include <cpptrace/cpptrace.hpp>
+#include <memory>
+#include <utility>
 
 #include "bout/parallel_boundary_op.hxx"
 #include "bout/parallel_boundary_region.hxx"
@@ -46,6 +49,8 @@
 #include <bout/output.hxx>
 #include <bout/utils.hxx>
 
+#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<Options> _tracking) {
+  tracking = std::move(_tracking);
+  tracking_state = 1;
+  selfname = name;
+  return *this;
+}
+
+template <typename T, typename>
+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<Field3D, bout::utils::EnableIfField<Field3D>>(const Field3D&,
+                                                              std::string);
+template void Field3D::_track<Field2D>(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},
+  });
+}

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;