Pvcode + Cvode improvements

include/bout/field.hxx

@@ -683,4 +683,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 f;
+}
+
 #endif /* FIELD_H */

include/bout/field3d.hxx

@@ -295,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, Options& tracking);
+
+  /// Disable tracking
+  Field3D& disableTracking() {
+    tracking = nullptr;
+    tracking_state = 0;
+    return *this;
+  }
+
   /////////////////////////////////////////////////////////
   // Data access
 
@@ -499,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};
@@ -514,6 +527,13 @@ private:
 
   /// RegionID over which the field is valid
   std::optional<size_t> regionID;
+
+  int tracking_state{0};
+  Options* tracking{nullptr};
+  std::string selfname;
+  template <class T>
+  Options* track(const T& change, std::string operation);
+  Options* track(const BoutReal& change, std::string operation);
 };
 
 // Non-member overloaded operators

include/bout/options.hxx

@@ -949,6 +949,9 @@ Tensor<BoutReal> Options::as<Tensor<BoutReal>>(const Tensor<BoutReal>& 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

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

@@ -245,6 +245,7 @@ Field3D& Field3D::operator=(const Field3D& rhs) {
   }
 
   TRACE("Field3D: Assignment from Field3D");
+  track(rhs, "operator=");
 
   // Copy base slice
   Field::operator=(rhs);
@@ -265,6 +266,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);
@@ -285,6 +287,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());
@@ -329,6 +332,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.
@@ -833,3 +837,47 @@ 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 <class T>
+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", changename},
+#endif
+    });
+    return &(*tracking)[outname];
+  }
+  return nullptr;
+}
+
+template Options* Field3D::track<Field3D>(const Field3D&, std::string);
+template Options* Field3D::track<Field2D>(const Field2D&, std::string);
+template Options* Field3D::track<FieldPerp>(const FieldPerp&, std::string);
+
+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", operation},
+        {"rhs.name", "BR"},
+    });
+    return &(*tracking)[outname];
+  }
+  return nullptr;
+}

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;