Add cell volume and areas

include/bout/coordinates.hxx

@@ -33,12 +33,14 @@
 #ifndef BOUT_COORDINATES_H
 #define BOUT_COORDINATES_H
 
+#include "bout/assert.hxx"
 #include "bout/field_data.hxx"
 #include <bout/bout_types.hxx>
 #include <bout/build_defines.hxx>
 #include <bout/field2d.hxx>
 #include <bout/field3d.hxx>
 #include <bout/paralleltransform.hxx>
+#include <optional>
 
 #include <array>
 #include <memory>
@@ -107,6 +109,126 @@ public:
   /// Covariant metric tensor
   FieldMetric g_11, g_22, g_33, g_12, g_13, g_23;
 
+  /// get g_22 at the cell faces;
+  const FieldMetric& g_22_ylow() const;
+  const FieldMetric& g_22_yhigh() const;
+  FieldMetric& g_22_ylow();
+  FieldMetric& g_22_yhigh();
+  // Cell Areas
+  const FieldMetric& cell_area_xlow() const {
+    if (!_cell_area_xlow.has_value()) {
+      _compute_cell_area_x();
+    }
+    ASSERT2(_cell_area_xlow.has_value());
+    return *_cell_area_xlow;
+  }
+  const FieldMetric& cell_area_xhigh() const {
+    if (!_cell_area_xhigh.has_value()) {
+      _compute_cell_area_x();
+    }
+    ASSERT2(_cell_area_xhigh.has_value());
+    return *_cell_area_xhigh;
+  }
+  const FieldMetric& cell_area_ylow() const {
+    if (!_cell_area_ylow.has_value()) {
+      _compute_cell_area_y();
+    }
+    ASSERT2(_cell_area_ylow.has_value());
+    return *_cell_area_ylow;
+  }
+  const FieldMetric& cell_area_yhigh() const {
+    if (!_cell_area_yhigh.has_value()) {
+      _compute_cell_area_y();
+    }
+    ASSERT2(_cell_area_yhigh.has_value());
+    return *_cell_area_yhigh;
+  }
+  const FieldMetric& cell_area_zlow() const {
+    if (!_cell_area_zlow.has_value()) {
+      _compute_cell_area_z();
+    }
+    ASSERT2(_cell_area_zlow.has_value());
+    return *_cell_area_zlow;
+  }
+  const FieldMetric& cell_area_zhigh() const {
+    if (!_cell_area_zhigh.has_value()) {
+      _compute_cell_area_z();
+    }
+    ASSERT2(_cell_area_zhigh.has_value());
+    return *_cell_area_zhigh;
+  }
+  FieldMetric& cell_area_xlow() {
+    if (!_cell_area_xlow.has_value()) {
+      _compute_cell_area_x();
+    }
+    ASSERT2(_cell_area_xlow.has_value());
+    return *_cell_area_xlow;
+  }
+  FieldMetric& cell_area_xhigh() {
+    if (!_cell_area_xhigh.has_value()) {
+      _compute_cell_area_x();
+    }
+    ASSERT2(_cell_area_xhigh.has_value());
+    return *_cell_area_xhigh;
+  }
+  FieldMetric& cell_area_ylow() {
+    if (!_cell_area_ylow.has_value()) {
+      _compute_cell_area_y();
+    }
+    ASSERT2(_cell_area_ylow.has_value());
+    return *_cell_area_ylow;
+  }
+  FieldMetric& cell_area_yhigh() {
+    if (!_cell_area_yhigh.has_value()) {
+      _compute_cell_area_y();
+    }
+    ASSERT2(_cell_area_yhigh.has_value());
+    return *_cell_area_yhigh;
+  }
+  FieldMetric& cell_area_zlow() {
+    if (!_cell_area_zlow.has_value()) {
+      _compute_cell_area_z();
+    }
+    ASSERT2(_cell_area_zlow.has_value());
+    return *_cell_area_zlow;
+  }
+  FieldMetric& cell_area_zhigh() {
+    if (!_cell_area_zhigh.has_value()) {
+      _compute_cell_area_z();
+    }
+    ASSERT2(_cell_area_zhigh.has_value());
+    return *_cell_area_zhigh;
+  }
+  // Cell Volume
+  const FieldMetric& cell_volume() const {
+    if (!_cell_volume.has_value()) {
+      _compute_cell_volume();
+    }
+    ASSERT2(_cell_volume.has_value());
+    return *_cell_volume;
+  }
+  FieldMetric& cell_volume() {
+    if (!_cell_volume.has_value()) {
+      _compute_cell_volume();
+    }
+    ASSERT2(_cell_volume.has_value());
+    return *_cell_volume;
+  }
+
+private:
+  mutable std::optional<FieldMetric> _g_22_ylow, _g_22_yhigh;
+  mutable std::optional<FieldMetric> _jxz_ylow, _jxz_yhigh, _jxz_centre;
+  mutable std::optional<FieldMetric> _cell_area_xlow, _cell_area_xhigh;
+  mutable std::optional<FieldMetric> _cell_area_ylow, _cell_area_yhigh;
+  mutable std::optional<FieldMetric> _cell_area_zlow, _cell_area_zhigh;
+  mutable std::optional<FieldMetric> _cell_volume;
+  void _compute_Jxz_cell_faces() const;
+  void _compute_cell_area_x() const;
+  void _compute_cell_area_y() const;
+  void _compute_cell_area_z() const;
+  void _compute_cell_volume() const;
+
+public:
   /// Christoffel symbol of the second kind (connection coefficients)
   FieldMetric G1_11, G1_22, G1_33, G1_12, G1_13, G1_23;
   FieldMetric G2_11, G2_22, G2_33, G2_12, G2_13, G2_23;

src/mesh/coordinates.cxx

@@ -7,20 +7,27 @@
 #include "bout/field3d.hxx"
 #include "bout/field_data.hxx"
 #include <bout/assert.hxx>
+#include <bout/bout_types.hxx>
+#include <bout/boutexception.hxx>
 #include <bout/build_defines.hxx>
 #include <bout/constants.hxx>
 #include <bout/coordinates.hxx>
 #include <bout/derivs.hxx>
 #include <bout/fft.hxx>
+#include <bout/field.hxx>
 #include <bout/globals.hxx>
 #include <bout/interpolation.hxx>
 #include <bout/output.hxx>
 #include <bout/output_bout_types.hxx>
+#include <bout/paralleltransform.hxx>
+#include <bout/region.hxx>
 #include <bout/sys/timer.hxx>
 #include <bout/utils.hxx>
 #include <bout/yboundary_regions.hxx>
 
+#include <cmath>
 #include <memory>
+#include <string>
 
 #include "invert3x3.hxx"
 #include "parallel/fci.hxx"
@@ -63,9 +70,8 @@ Field2D interpolateAndExtrapolate(const Field2D& f, CELL_LOC location, bool extr
         // Can use no_extra_interpolate argument to skip the extra interpolation when we
         // want to extrapolate the Christoffel symbol terms which come from derivatives so
         // don't have the extra point set already
-        if ((location == CELL_XLOW) && (bndry->bx > 0)) {
-          extrap_start = 1;
-        } else if ((location == CELL_YLOW) && (bndry->by > 0)) {
+        if (((location == CELL_XLOW) && (bndry->bx > 0))
+            || ((location == CELL_YLOW) && (bndry->by > 0))) {
           extrap_start = 1;
         }
       }
@@ -84,7 +90,7 @@ Field2D interpolateAndExtrapolate(const Field2D& f, CELL_LOC location, bool extr
               (9.
                    * (f(bndry->x - bndry->bx, bndry->y - bndry->by)
                       + f(bndry->x, bndry->y))
-               - f(bndry->x - 2 * bndry->bx, bndry->y - 2 * bndry->by)
+               - f(bndry->x - (2 * bndry->bx), bndry->y - (2 * bndry->by))
                - f(bndry->x + bndry->bx, bndry->y + bndry->by))
               / 16.;
         }
@@ -110,8 +116,8 @@ Field2D interpolateAndExtrapolate(const Field2D& f, CELL_LOC location, bool extr
           }
           // extrapolate into boundary guard cells if there are enough grid points
           for (int i = extrap_start; i < bndry->width; i++) {
-            int xi = bndry->x + i * bndry->bx;
-            int yi = bndry->y + i * bndry->by;
+            const int xi = bndry->x + (i * bndry->bx);
+            const int yi = bndry->y + (i * bndry->by);
             result(xi, yi) = 3.0 * result(xi - bndry->bx, yi - bndry->by)
                              - 3.0 * result(xi - 2 * bndry->bx, yi - 2 * bndry->by)
                              + result(xi - 3 * bndry->bx, yi - 3 * bndry->by);
@@ -501,13 +507,13 @@ Coordinates::Coordinates(Mesh* mesh, Options* options)
       output_warn.write("Not all covariant components of metric tensor found. "
                         "Calculating all from the contravariant tensor\n");
       /// Calculate contravariant metric components if not found
-      if (calcCovariant("RGN_NOCORNERS")) {
+      if (calcCovariant("RGN_NOCORNERS") != 0) {
         throw BoutException("Error in calcCovariant call");
       }
     }
   } else {
     /// Calculate contravariant metric components if not found
-    if (calcCovariant("RGN_NOCORNERS")) {
+    if (calcCovariant("RGN_NOCORNERS") != 0) {
       throw BoutException("Error in calcCovariant call");
     }
   }
@@ -1873,7 +1879,7 @@ Field2D Coordinates::Laplace_perpXY([[maybe_unused]] const Field2D& A,
 
 const Coordinates::FieldMetric& Coordinates::invSg() const {
   if (invSgCache == nullptr) {
-    auto ptr = std::make_unique<FieldMetric>();
+    auto ptr = std::make_unique<Coordinates::FieldMetric>();
     (*ptr) = 1.0 / sqrt(g_22);
     invSgCache = std::move(ptr);
   }
@@ -1893,7 +1899,7 @@ Coordinates::Grad2_par2_DDY_invSg(CELL_LOC outloc, const std::string& method) co
   invSgCache->applyParallelBoundary("parallel_neumann_o2");
 
   // cache
-  auto ptr = std::make_unique<FieldMetric>();
+  auto ptr = std::make_unique<Coordinates::FieldMetric>();
   *ptr = DDY(*invSgCache, outloc, method) * invSg();
   Grad2_par2_DDY_invSgCache[method] = std::move(ptr);
   return *Grad2_par2_DDY_invSgCache[method];
@@ -2003,6 +2009,136 @@ void Coordinates::checkContravariant() {
   }
 }
 
+const Coordinates::FieldMetric& Coordinates::g_22_ylow() const {
+  if (_g_22_ylow.has_value()) {
+    return *_g_22_ylow;
+  }
+  _g_22_ylow.emplace(emptyFrom(g_22));
+  //_g_22_ylow->setLocation(CELL_YLOW);
+  auto* mesh = Bxy.getMesh();
+  if (Bxy.isFci()) {
+    if (mesh->get(_g_22_ylow.value(), "g_22_cell_ylow", 0.0, false) != 0) {
+      throw BoutException("The grid file does not contain `g_22_cell_ylow`.");
+    }
+  } else {
+    ASSERT0(mesh->ystart > 0);
+    BOUT_FOR(i, g_22.getRegion("RGN_NOY")) {
+      _g_22_ylow.value()[i] = SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.ym()])));
+    }
+  }
+  return g_22_ylow();
+}
+
+const Coordinates::FieldMetric& Coordinates::g_22_yhigh() const {
+  if (_g_22_yhigh.has_value()) {
+    return *_g_22_yhigh;
+  }
+  _g_22_yhigh.emplace(emptyFrom(g_22));
+  //_g_22_yhigh->setLocation(CELL_YHIGH);
+  auto* mesh = Bxy.getMesh();
+  if (Bxy.isFci()) {
+    if (mesh->get(_g_22_yhigh.value(), "g_22_cell_yhigh", 0.0, false) != 0) {
+      throw BoutException("The grid file does not contain `g_22_cell_yhigh`.");
+    }
+  } else {
+    ASSERT0(mesh->ystart > 0);
+    BOUT_FOR(i, g_22.getRegion("RGN_NOY")) {
+      _g_22_yhigh.value()[i] = SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.yp()])));
+    }
+  }
+  return g_22_yhigh();
+}
+
+void Coordinates::_compute_Jxz_cell_faces() const {
+  _jxz_centre.emplace(sqrt(g_11 * g_33 - SQ(g_13)));
+  _jxz_ylow.emplace(emptyFrom(_jxz_centre.value()));
+  //_jxz_ylow->setLocation(CELL_YLOW);
+  _jxz_yhigh.emplace(emptyFrom(_jxz_centre.value()));
+  //_jxz_yhigh->setLocation(CELL_YHIGH);
+  auto* mesh = _jxz_centre->getMesh();
+  if (Bxy.isFci()) {
+    Coordinates::FieldMetric By_c;
+    Coordinates::FieldMetric By_h;
+    Coordinates::FieldMetric By_l;
+    if (mesh->get(By_c, "By", 0.0, false, CELL_CENTRE) != 0) {
+      throw BoutException("The grid file does not contain `By`.");
+    }
+    if (mesh->get(By_l, "By_cell_ylow", 0.0, false) != 0) {
+      throw BoutException("The grid file does not contain `By_cell_ylow`.");
+    }
+    if (mesh->get(By_h, "By_cell_yhigh", 0.0, false) != 0) {
+      throw BoutException("The grid file does not contain `By_cell_yhigh`.");
+    }
+    BOUT_FOR(i, By_c.getRegion("RGN_NOY")) {
+      (*_jxz_ylow)[i] = By_c[i] / By_l[i] * (*_jxz_centre)[i];
+      (*_jxz_yhigh)[i] = By_c[i] / By_h[i] * (*_jxz_centre)[i];
+    }
+  } else {
+    ASSERT0(mesh->ystart > 0);
+    BOUT_FOR(i, _jxz_centre->getRegion("RGN_NOY")) {
+      (*_jxz_ylow)[i] = 0.5 * ((*_jxz_centre)[i] + (*_jxz_centre)[i.ym()]);
+      (*_jxz_yhigh)[i] = 0.5 * ((*_jxz_centre)[i] + (*_jxz_centre)[i.yp()]);
+    }
+  }
+}
+
+void Coordinates::_compute_cell_area_x() const {
+  const auto area_centre = sqrt(g_22 * g_33 - SQ(g_23)) * dy * dz;
+  _cell_area_xlow.emplace(emptyFrom(area_centre));
+  _cell_area_xhigh.emplace(emptyFrom(area_centre));
+  // We cannot setLocation, as that would trigger the computation of staggered
+  // metrics.
+  auto* mesh = Bxy.getMesh();
+  ASSERT0(mesh->xstart > 0);
+  BOUT_FOR(i, area_centre.getRegion("RGN_NOX")) {
+    (*_cell_area_xlow)[i] = 0.5 * (area_centre[i] + area_centre[i.xm()]);
+    (*_cell_area_xhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.xp()]);
+  }
+}
+
+void Coordinates::_compute_cell_area_y() const {
+  _compute_Jxz_cell_faces();
+  ASSERT2(_jxz_centre.has_value());
+  ASSERT2(_jxz_ylow.has_value());
+  ASSERT2(_jxz_yhigh.has_value());
+  if (_jxz_centre->isFci()) {
+    ASSERT3(isUniform(dx, true, "RGN_ALL"));
+    ASSERT2(isUniform(dx, false, "RGN_ALL"));
+    ASSERT3(isUniform(dz, true, "RGN_ALL"));
+    ASSERT2(isUniform(dz, false, "RGN_ALL"));
+    _cell_area_ylow.emplace(*_jxz_ylow * dx * dz);
+    _cell_area_yhigh.emplace(*_jxz_yhigh * dx * dz);
+  } else {
+    // Field aligned
+    const auto area_centre = sqrt(g_11 * g_33 - SQ(g_13)) * dx * dz;
+    _cell_area_ylow.emplace(emptyFrom(area_centre));
+    _cell_area_yhigh.emplace(emptyFrom(area_centre));
+    // We cannot setLocation, as that would trigger the computation of staggered
+    // metrics.
+    auto* mesh = Bxy.getMesh();
+    ASSERT0(mesh->ystart > 0);
+    BOUT_FOR(i, mesh->getRegion("RGN_NOY")) {
+      (*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]);
+      (*_cell_area_yhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.yp()]);
+    }
+  }
+}
+
+void Coordinates::_compute_cell_area_z() const {
+  const auto area_centre = sqrt(g_11 * g_22 - SQ(g_12)) * dx * dy;
+  _cell_area_zlow.emplace(emptyFrom(area_centre));
+  _cell_area_zhigh.emplace(emptyFrom(area_centre));
+  // We cannot setLocation, as that would trigger the computation of staggered
+  // metrics.
+  //ASSERT0(mesh->zstart > 0);
+  BOUT_FOR(i, area_centre.getRegion("RGN_NOZ")) {
+    (*_cell_area_zlow)[i] = 0.5 * (area_centre[i] + area_centre[i.zm()]);
+    (*_cell_area_zhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.zp()]);
+  }
+}
+
+void Coordinates::_compute_cell_volume() const { _cell_volume.emplace(J * dx * dy * dz); }
+
 std::shared_ptr<YBoundary> Coordinates::makeYBoundary(YBndryType type) const {
   return std::make_shared<YBoundary>(type, localoptions, *localmesh);
 }