Add optional boolean parameter communicate to Coordinates setters and MetricTensor inverse() method

Add method Coordinates::communicateMetricTensor() and call after constructing Coordinates

to avoid circular dependency (between Mesh and Coordinates) as a result of the setter methods calling communicate()

Author: tomchapman

include/bout/coordinates.hxx

@@ -108,9 +108,9 @@ public:
   const BoutReal& J(int x, int y) const { return J()(x, y); }
 #endif
 
-  void setDx(FieldMetric dx);
-  void setDy(FieldMetric dy);
-  void setDz(FieldMetric dz);
+  void setDx(FieldMetric dx, const bool communicate = true);
+  void setDy(FieldMetric dy, const bool communicate = true);
+  void setDz(FieldMetric dz, const bool communicate = true);
 
   void setD1_dx(FieldMetric d1_dx) { d1_dx_ = std::move(d1_dx); }
   void setD1_dy(FieldMetric d1_dy) { d1_dy_ = std::move(d1_dy); }
@@ -235,15 +235,17 @@ public:
   void setMetricTensor(const ContravariantMetricTensor& contravariant_metric_tensor,
                        const CovariantMetricTensor& covariant_metric_tensor);
 
+  void communicateMetricTensor();
+
   ///< Coordinate system Jacobian, so volume of cell is J*dx*dy*dz
   FieldMetric& J() const;
 
   ///< Magnitude of B = nabla z times nabla x
   const FieldMetric& Bxy() const { return Bxy_; }
 
-  void setJ(const FieldMetric& J);
+  void setJ(const FieldMetric& J, const bool communicate = true);
 
-  void setBxy(FieldMetric Bxy);
+  void setBxy(FieldMetric Bxy, const bool communicate = true);
 
   /// d pitch angle / dx. Needed for vector differentials (Curl)
   const FieldMetric& ShiftTorsion() const { return ShiftTorsion_; }

include/bout/mesh.hxx

@@ -653,6 +653,8 @@ public:
     inserted.first->second->recalculateAndReset(recalculate_staggered,
                                                 force_interpolate_from_centre);
 
+    inserted.first->second->communicateMetricTensor();
+
     return inserted.first->second;
   }
 

include/bout/metric_tensor.hxx

@@ -61,7 +61,7 @@ public:
     g23_m = metric_tensor.g23();
   }
 
-  MetricTensor inverse(const std::string& region = "RGN_ALL");
+  MetricTensor inverse(const std::string& region = "RGN_ALL", const bool communicate = true);
 
   // Transforms the MetricTensor by applying the given function to every component
   template <class F>
@@ -74,8 +74,11 @@ public:
     g23_m = function(g23_m);
   }
 
+  void communicate() const;
+
 private:
   FieldMetric g11_m, g22_m, g33_m, g12_m, g13_m, g23_m;
+
 };
 
 class CovariantMetricTensor : public MetricTensor {

src/mesh/coordinates.cxx

@@ -501,9 +501,9 @@ void Coordinates::interpolateFromCoordinates(Options* mesh_options,
   checkCovariant();
 
   setJ(interpolateAndExtrapolate(coords_in->J(), location, true, true, false,
-                                 transform.get()));
+                                 transform.get()), false);
   setBxy(interpolateAndExtrapolate(coords_in->Bxy(), location, true, true, false,
-                                   transform.get()));
+                                   transform.get()), false);
 
   bout::checkFinite(J(), "The Jacobian", "RGN_NOCORNERS");
   bout::checkPositive(J(), "The Jacobian", "RGN_NOCORNERS");
@@ -567,9 +567,9 @@ void Coordinates::readFromMesh(Options* mesh_options, const std::string& suffix)
   dy_ = interpolateAndExtrapolate(dy_, location, extrapolate_x, extrapolate_y, false,
                                   transform.get());
 
-  setDx(dx_);
-  setDy(dy_);
-  setDz(dz_);
+  setDx(dx_, false);
+  setDy(dy_, false);
+  setDz(dz_, false);
 
   // grid data source has staggered fields, so read instead of interpolating
   // Diagonal components of metric tensor g^{ij} (default to 1)
@@ -618,7 +618,7 @@ void Coordinates::readFromMesh(Options* mesh_options, const std::string& suffix)
     output_warn.write("\tWARNING! Covariant components of metric tensor set manually. "
                       "Contravariant components NOT recalculated\n");
   } else {
-    covariantMetricTensor.setMetricTensor(contravariantMetricTensor.inverse());
+    covariantMetricTensor.setMetricTensor(contravariantMetricTensor.inverse("RGN_ALL", false));
     output_warn.write("Not all covariant components of metric tensor found. "
                       "Calculating all from the contravariant tensor\n");
   }
@@ -639,15 +639,15 @@ void Coordinates::readFromMesh(Options* mesh_options, const std::string& suffix)
     const auto J_from_file = getAtLoc(localmesh, "J", suffix, location);
     // Compare calculated and loaded values
     output_warn.write("\tMaximum difference in J is {:e}\n", max(abs(J() - J_from_file)));
-    setJ(J_from_file);
+    setJ(J_from_file, false);
 
     communicate(J());
   }
 
   // More robust to extrapolate derived quantities directly, rather than
   // deriving from extrapolated covariant metric components
   setJ(interpolateAndExtrapolate(J(), location, extrapolate_x, extrapolate_y, false,
-                                 transform.get()));
+                                 transform.get()), false);
 
   // Check jacobian
   bout::checkFinite(J(), "J" + suffix, "RGN_NOCORNERS");
@@ -662,15 +662,15 @@ void Coordinates::readFromMesh(Options* mesh_options, const std::string& suffix)
                       "Calculating from metric tensor\n",
                       suffix);
     // Re-evaluate Bxy using new J
-    setBxy(recalculateBxy());
+    setBxy(recalculateBxy(), false);
   } else {
     const auto Bcalc = getAtLoc(localmesh, "Bxy", suffix, location);
-    setBxy(Bcalc);
+    setBxy(Bcalc, false);
     output_warn.write("\tMaximum difference in Bxy is {:e}\n", max(abs(Bxy() - Bcalc)));
   }
 
   setBxy(interpolateAndExtrapolate(Bxy(), location, extrapolate_x, extrapolate_y, false,
-                                   transform.get()));
+                                   transform.get()), false);
 
   // Check Bxy
   bout::checkFinite(Bxy(), "Bxy" + suffix, "RGN_NOCORNERS");
@@ -762,31 +762,34 @@ const Field2D& Coordinates::zlength() const {
   return *zlength_cache;
 }
 
-void Coordinates::setDx(FieldMetric dx) {
+void Coordinates::setDx(FieldMetric dx, const bool communicate) {
   if (min(abs(dx)) < 1e-8) {
     throw BoutException("dx magnitude less than 1e-8");
   }
-
   dx_ = std::move(dx);
-  localmesh->communicate(dx_);
+  if (communicate) {
+    localmesh->communicate(dx_);
+  }
 }
 
-void Coordinates::setDy(FieldMetric dy) {
+void Coordinates::setDy(FieldMetric dy, const bool communicate) {
   if (min(abs(dy)) < 1e-8) {
     throw BoutException("dy magnitude less than 1e-8");
   }
-
   dy_ = std::move(dy);
-  localmesh->communicate(dy_);
+  if (communicate) {
+    localmesh->communicate(dy_);
+  }
 }
 
-void Coordinates::setDz(FieldMetric dz) {
+void Coordinates::setDz(FieldMetric dz, const bool communicate) {
   if (min(abs(dz)) < 1e-8) {
     throw BoutException("dz magnitude less than 1e-8");
   }
-
   dz_ = std::move(dz);
-  localmesh->communicate(dz_);
+  if (communicate) {
+    localmesh->communicate(dz_);
+  }
 }
 
 void Coordinates::recalculateAndReset(bool recalculate_staggered,
@@ -1491,19 +1494,23 @@ FieldMetric& Coordinates::J() const {
   return *jacobian_cache;
 }
 
-void Coordinates::setJ(const FieldMetric& J) {
+void Coordinates::setJ(const FieldMetric& J, const bool communicate) {
   bout::checkFinite(J, "J", "RGN_NOCORNERS");
   bout::checkPositive(J, "J", "RGN_NOCORNERS");
 
   //TODO: Calculate J and check value is close
   jacobian_cache = std::make_unique<FieldMetric>(J);
-  localmesh->communicate(*jacobian_cache);
+  if (communicate) {
+    localmesh->communicate(*jacobian_cache);
+  }
 }
 
-void Coordinates::setBxy(FieldMetric Bxy) {
+void Coordinates::setBxy(FieldMetric Bxy, const bool communicate) {
   //TODO: Calculate Bxy and check value is close
   Bxy_ = std::move(Bxy);
-  localmesh->communicate(Bxy_);
+  if (communicate) {
+    localmesh->communicate(Bxy_);
+  }
 }
 
 void Coordinates::setContravariantMetricTensor(
@@ -1529,3 +1536,8 @@ void Coordinates::setMetricTensor(
   contravariantMetricTensor.setMetricTensor(contravariant_metric_tensor);
   covariantMetricTensor.setMetricTensor(covariant_metric_tensor);
 }
+
+void Coordinates::communicateMetricTensor() {
+  contravariantMetricTensor.communicate();
+  covariantMetricTensor.communicate();
+}

src/mesh/metric_tensor.cxx

@@ -72,7 +72,7 @@ void MetricTensor::check(int ystart) {
   }
 }
 
-MetricTensor MetricTensor::inverse(const std::string& region) {
+MetricTensor MetricTensor::inverse(const std::string& region, const bool communicate) {
 
   TRACE("MetricTensor::inverse");
 
@@ -124,7 +124,19 @@ MetricTensor MetricTensor::inverse(const std::string& region) {
 
   output_info.write("\tMaximum error in off-diagonal inversion is {:e}\n",
                     off_diagonal_maxerr);
-
-  g_11.getMesh()->communicate(g_11, g_22, g_33, g_12, g_13, g_23);
+  if (communicate) {
+    MetricTensor::communicate();
+  }
   return MetricTensor(g_11, g_22, g_33, g_12, g_13, g_23);
 }
+
+void MetricTensor::communicate() const {
+  // copy to non-const variables (because communicate() cannot take const references)
+  auto tmp = g11_m;
+  auto tmp2 = g22_m;
+  auto tmp3 = g33_m;
+  auto tmp4 = g12_m;
+  auto tmp5 = g13_m;
+  auto tmp6 = g23_m;
+  g11_m.getMesh()->communicate(tmp, tmp2, tmp3, tmp4, tmp5, tmp6);
+}