Pvcode + Cvode improvements

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 <typename T, typename = bout::utils::EnableIfField<T>, 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 <typename T, typename = bout::utils::EnableIfField<T>, 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<T>(f);
+}
+
 #endif /* FIELD_H */

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

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,8 @@
 #include <bout/globals.hxx>
 
 #include <cmath>
+#include <memory>
+#include <utility>
 
 #include "bout/parallel_boundary_op.hxx"
 #include "bout/parallel_boundary_region.hxx"
@@ -47,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) {
@@ -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<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");
+  // 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<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 outname{fmt::format("track_{:s}_{:d}", selfname, tracking_state++)};
+  locked->set(outname, change, "tracking");
+  (*locked)[outname].setAttributes({
+      {"operation", operation},
+      {"rhs.name", "BoutReal"},
+  });
+}

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;

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 <bout/field2d.hxx>
-#include <bout/field3d.hxx>
-#include <bout/globals.hxx>
-#include <bout/interpolation.hxx>
-#include <bout/mesh.hxx>
-#include <bout/region.hxx>
+#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"
 """