Merge pull request #3056 from boutproject/fix-after-reactoring-commun…

…icate Fixes to address failing tests after changes to refactor-coordinates
boutproject · Jan 17, 2025 · e117fa5 · e117fa5
2 parents 47cff9e + 834b090
commit e117fa5
Show file tree

Hide file tree

Showing 7 changed files with 86 additions and 45 deletions.
diff --git a/include/bout/coordinates.hxx b/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,19 @@ public:
   void setMetricTensor(const ContravariantMetricTensor& contravariant_metric_tensor,
                        const CovariantMetricTensor& covariant_metric_tensor);
 
+  void communicateMetricTensor();
+
+  void communicateDz();
+
   ///< 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_; }

diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
@@ -653,6 +653,9 @@ public:
     inserted.first->second->recalculateAndReset(recalculate_staggered,
                                                 force_interpolate_from_centre);
 
+    inserted.first->second->communicateMetricTensor();
+    inserted.first->second->communicateDz();
+
     return inserted.first->second;
   }
 

diff --git a/include/bout/metric_tensor.hxx b/include/bout/metric_tensor.hxx
@@ -61,7 +61,8 @@ 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,6 +75,8 @@ public:
     g23_m = function(g23_m);
   }
 
+  void communicate() const;
+
 private:
   FieldMetric g11_m, g22_m, g33_m, g12_m, g13_m, g23_m;
 };

diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx
@@ -501,9 +501,11 @@ 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 +569,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 +620,8 @@ 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 +642,16 @@ 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 +666,16 @@ 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 +767,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 +1499,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 +1541,10 @@ void Coordinates::setMetricTensor(
   contravariantMetricTensor.setMetricTensor(contravariant_metric_tensor);
   covariantMetricTensor.setMetricTensor(covariant_metric_tensor);
 }
+
+void Coordinates::communicateMetricTensor() {
+  contravariantMetricTensor.communicate();
+  covariantMetricTensor.communicate();
+}
+
+void Coordinates::communicateDz() { localmesh->communicate(dz_); }
diff --git a/src/mesh/metric_tensor.cxx b/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);
+}
diff --git a/tests/integrated/test-snb/test_snb.cxx b/tests/integrated/test-snb/test_snb.cxx
@@ -209,20 +209,18 @@ int main(int argc, char** argv) {
     // Change the mesh spacing and cell volume (Jdy)
     Coordinates* coord = Te.getCoordinates();
 
-    {
-      auto dy = emptyFrom(coord->dx());
-      auto J = emptyFrom(coord->J());
-      for (int x = mesh->xstart; x <= mesh->xend; x++) {
-        for (int y = mesh->ystart; y <= mesh->yend; y++) {
-          const double y_n = (double(y) + 0.5) / double(mesh->yend + 1);
-
-          dy(x, y) = 1. - 0.9 * y_n;
-          J(x, y) = 1. + y_n * y_n;
-        }
+    auto dy_copy = coord->dy();
+    auto J_copy = coord->J();
+    for (int x = mesh->xstart; x <= mesh->xend; x++) {
+      for (int y = mesh->ystart; y <= mesh->yend; y++) {
+        const double y_n = (double(y) + 0.5) / double(mesh->yend + 1);
+
+        dy_copy(x, y) = 1. - 0.9 * y_n;
+        J_copy(x, y) = 1. + y_n * y_n;
       }
-      coord->setDy(dy);
-      coord->setJ(J);
     }
+    coord->setDy(dy_copy);
+    coord->setJ(J_copy);
 
     HeatFluxSNB snb;
 

diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx
@@ -440,6 +440,7 @@ public:
         Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{0.0},
         Field2D{0.0}, Field2D{0.0}, Field2D{1.0}, Field2D{1.0}, Field2D{1.0},
         Field2D{0.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0});
+    static_cast<FakeMesh*>(bout::globals::mesh)->setCoordinates(test_coords);
 
     // Set nonuniform corrections
     test_coords->setNon_uniform(true);
@@ -488,6 +489,7 @@ public:
         Field2D{0.0, mesh_staggered}, Field2D{0.0, mesh_staggered},
         Field2D{0.0, mesh_staggered}, Field2D{0.0, mesh_staggered},
         Field2D{0.0, mesh_staggered});
+    static_cast<FakeMesh*>(mesh_staggered)->setCoordinates(test_coords_staggered);
 
     // Set nonuniform corrections
     test_coords_staggered->setNon_uniform(true);