From 3f521d0f30eefbbc2038db2b3acdbea530e78c32 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 20 Jan 2026 09:25:47 +0100 Subject: [PATCH 01/22] Add metric info on cell faces for FCI Information like perpendicular J, Bxyz and g_22 are only known to the grid generator, but are needed for more correct parallel derivatives --- include/bout/coordinates.hxx | 16 ++++ src/mesh/coordinates.cxx | 140 ++++++++++++++++++++++++++++++++++- 2 files changed, 154 insertions(+), 2 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index e7ead42ee5..0b1ebb51e6 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -102,6 +102,20 @@ 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; + /// get Jxz at the cell faces or cell centre + const FieldMetric& Jxz_ylow() const; + const FieldMetric& Jxz_yhigh() const; + const FieldMetric& Jxz() const; + +private: + mutable std::optional _g_22_ylow, _g_22_yhigh; + mutable std::optional _Jxz_ylow, _Jxz_yhigh, _Jxz_centre; + void _compute_Jxz_cell_faces() 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; @@ -236,12 +250,14 @@ private: /// Cache variable for Grad2_par2 mutable std::map> Grad2_par2_DDY_invSgCache; mutable std::unique_ptr invSgCache{nullptr}; + mutable std::optional JgCache; /// Set the parallel (y) transform from the options file. /// Used in the constructor to create the transform object. void setParallelTransform(Options* options); const FieldMetric& invSg() const; + const FieldMetric& Jg() const; const FieldMetric& Grad2_par2_DDY_invSg(CELL_LOC outloc, const std::string& method) const; diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index d2634d3f88..d8ba0a02ee 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1259,6 +1259,31 @@ int Coordinates::calcCovariant(const std::string& region) { output_info.write("\tLocal maximum error in off-diagonal inversion is {:e}\n", maxerr); + if (Bxy.isFci()) { + BoutReal maxError = 0; + Options* options = Options::getRoot(); + auto BJg = Bxy.asField3DParallel() * J / sqrt(g_22.asField3DParallel()); + auto* mesh = localmesh; + for (int p = -mesh->ystart; p <= mesh->ystart; p++) { + if (p == 0) { + continue; + } + BOUT_FOR(i, BJg.getRegion("RGN_NOBNDRY")) { + auto local = BJg[i] / BJg.ynext(p)[i.yp(p)]; + maxError = std::max(std::abs(local - 1), maxError); + } + } + const BoutReal allowedError = (*options)["allowedFluxError"].withDefault(1e-6); + if (maxError < allowedError / 100) { + output_info.write("\tInfo: The maximum flux conservation error is {:e}", maxError); + } else if (maxError < allowedError) { + output_warn.write("\tWarning: The maximum flux conservation error is {:e}", + maxError); + } else { + throw BoutException("Error: The maximum flux conservation error is {:e}", maxError); + } + } + return 0; } @@ -1878,13 +1903,31 @@ Field2D Coordinates::Laplace_perpXY([[maybe_unused]] const Field2D& A, const Coordinates::FieldMetric& Coordinates::invSg() const { if (invSgCache == nullptr) { - auto ptr = std::make_unique(); + auto ptr = std::make_unique(); (*ptr) = 1.0 / sqrt(g_22); invSgCache = std::move(ptr); } return *invSgCache; } +const Coordinates::FieldMetric& Coordinates::Jg() const { + if (not JgCache.has_value()) { + auto* coords = this; // + // Need to modify yup and ydown fields + Field3D Jg = coords->Jxz(); + Jg.splitParallelSlicesAndAllocate(); + Field3DParallel B = coords->Bxy; + for (size_t j = 0; j < B.numberParallelSlices(); ++j) { + BOUT_FOR(i, B.getRegion("RGN_NOBNDRY")) { + Jg.yup(j)[i.yp(j + 1)] = Jg[i] * B.yup(j)[i.yp(j + 1)] / B[i]; + Jg.ydown(j)[i.ym(j + 1)] = Jg[i] * B.ydown(j)[i.ym(j + 1)] / B[i]; + } + } + JgCache = Jg; + } + return *JgCache; +} + const Coordinates::FieldMetric& Coordinates::Grad2_par2_DDY_invSg(CELL_LOC outloc, const std::string& method) const { if (auto search = Grad2_par2_DDY_invSgCache.find(method); @@ -1898,7 +1941,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(); + auto ptr = std::make_unique(); *ptr = DDY(*invSgCache, outloc, method) * invSg(); Grad2_par2_DDY_invSgCache[method] = std::move(ptr); return *Grad2_par2_DDY_invSgCache[method]; @@ -2007,3 +2050,96 @@ 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)) { //, CELL_YLOW)) { + 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 * (sqrt(g_22[i]) + 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)) { //, CELL_YHIGH)) { + 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 * (sqrt(g_22[i]) + sqrt(g_22[i.yp()]))); + } + } + return g_22_yhigh(); +} + +const Coordinates::FieldMetric& Coordinates::Jxz_ylow() const { + if (!_Jxz_ylow.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_Jxz_ylow; +} +const Coordinates::FieldMetric& Coordinates::Jxz_yhigh() const { + if (!_Jxz_yhigh.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_Jxz_yhigh; +} +const Coordinates::FieldMetric& Coordinates::Jxz() const { + if (!_Jxz_centre.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_Jxz_centre; +} + +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)) { + throw BoutException("The grid file does not contain `By`."); + } + if (mesh->get(By_l, "By_cell_ylow", 0.0, false)) { //, CELL_YLOW)) { + throw BoutException("The grid file does not contain `By_cell_ylow`."); + } + if (mesh->get(By_h, "By_cell_yhigh", 0.0, false)) { //, CELL_YHIGH)) { + 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()]); + } + } +} From 3065f88847c5bc9a5a897c4aa699b7dd612549c1 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 30 Jan 2026 10:08:19 +0100 Subject: [PATCH 02/22] Do not used preserved names _[A-Z].* is preserved for compilers, thus use _[a-z].* --- include/bout/coordinates.hxx | 2 +- src/mesh/coordinates.cxx | 35 ++++++++++++++++++----------------- 2 files changed, 19 insertions(+), 18 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 0b1ebb51e6..0d405a0123 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -112,7 +112,7 @@ public: private: mutable std::optional _g_22_ylow, _g_22_yhigh; - mutable std::optional _Jxz_ylow, _Jxz_yhigh, _Jxz_centre; + mutable std::optional _jxz_ylow, _jxz_yhigh, _jxz_centre; void _compute_Jxz_cell_faces() const; public: diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index d8ba0a02ee..a2f1a56222 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2093,31 +2093,32 @@ const Coordinates::FieldMetric& Coordinates::g_22_yhigh() const { } const Coordinates::FieldMetric& Coordinates::Jxz_ylow() const { - if (!_Jxz_ylow.has_value()) { + if (!_jxz_ylow.has_value()) { _compute_Jxz_cell_faces(); } - return *_Jxz_ylow; + return *_jxz_ylow; } const Coordinates::FieldMetric& Coordinates::Jxz_yhigh() const { - if (!_Jxz_yhigh.has_value()) { + if (!_jxz_yhigh.has_value()) { _compute_Jxz_cell_faces(); } - return *_Jxz_yhigh; + return *_jxz_yhigh; } const Coordinates::FieldMetric& Coordinates::Jxz() const { - if (!_Jxz_centre.has_value()) { + if (!_jxz_centre.has_value()) { _compute_Jxz_cell_faces(); } - return *_Jxz_centre; + return *_jxz_centre; +} } 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(); + _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; @@ -2132,14 +2133,14 @@ void Coordinates::_compute_Jxz_cell_faces() const { 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]; + (*_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()]); + 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()]); } } } From 1c0e385ca51301fee213d26254e38d02e0e8d828 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 30 Jan 2026 10:08:33 +0100 Subject: [PATCH 03/22] Add non-const versions --- include/bout/coordinates.hxx | 5 +++++ src/mesh/coordinates.cxx | 28 ++++++++++++++++++++++++++++ 2 files changed, 33 insertions(+) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 0d405a0123..cf8e7cf39f 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -105,10 +105,15 @@ public: /// 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(); /// get Jxz at the cell faces or cell centre const FieldMetric& Jxz_ylow() const; const FieldMetric& Jxz_yhigh() const; const FieldMetric& Jxz() const; + FieldMetric& Jxz_ylow(); + FieldMetric& Jxz_yhigh(); + FieldMetric& Jxz(); private: mutable std::optional _g_22_ylow, _g_22_yhigh; diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index a2f1a56222..23fa77bbce 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2092,6 +2092,16 @@ const Coordinates::FieldMetric& Coordinates::g_22_yhigh() const { return g_22_yhigh(); } +Coordinates::FieldMetric& Coordinates::g_22_yhigh() { + return const_cast( + const_cast(this)->g_22_yhigh()); +} + +Coordinates::FieldMetric& Coordinates::g_22_ylow() { + return const_cast( + const_cast(this)->g_22_ylow()); +} + const Coordinates::FieldMetric& Coordinates::Jxz_ylow() const { if (!_jxz_ylow.has_value()) { _compute_Jxz_cell_faces(); @@ -2110,6 +2120,24 @@ const Coordinates::FieldMetric& Coordinates::Jxz() const { } return *_jxz_centre; } + +Coordinates::FieldMetric& Coordinates::Jxz_ylow() { + if (!_jxz_ylow.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_jxz_ylow; +} +Coordinates::FieldMetric& Coordinates::Jxz_yhigh() { + if (!_jxz_yhigh.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_jxz_yhigh; +} +Coordinates::FieldMetric& Coordinates::Jxz() { + if (!_jxz_centre.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_jxz_centre; } void Coordinates::_compute_Jxz_cell_faces() const { From 5c28542f9400fcd8593c8281efae10371d674045 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 12 Mar 2026 13:12:01 +0100 Subject: [PATCH 04/22] Add cell areas and cell volumes --- include/bout/coordinates.hxx | 137 +++++++++++++++++++++++++++++++++-- src/mesh/coordinates.cxx | 100 +++++++++++++------------ 2 files changed, 183 insertions(+), 54 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index cf8e7cf39f..402efecd20 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -108,17 +108,142 @@ public: FieldMetric& g_22_ylow(); FieldMetric& g_22_yhigh(); /// get Jxz at the cell faces or cell centre - const FieldMetric& Jxz_ylow() const; - const FieldMetric& Jxz_yhigh() const; - const FieldMetric& Jxz() const; - FieldMetric& Jxz_ylow(); - FieldMetric& Jxz_yhigh(); - FieldMetric& Jxz(); + const FieldMetric& Jxz_ylow() const { + if (!_jxz_ylow.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_jxz_ylow; + } + const FieldMetric& Jxz_yhigh() const { + if (!_jxz_yhigh.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_jxz_yhigh; + } + const FieldMetric& Jxz() const { + if (!_jxz_centre.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_jxz_centre; + } + FieldMetric& Jxz_ylow() { + if (!_jxz_ylow.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_jxz_ylow; + } + FieldMetric& Jxz_yhigh() { + if (!_jxz_yhigh.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_jxz_yhigh; + } + FieldMetric& Jxz() { + if (!_jxz_centre.has_value()) { + _compute_Jxz_cell_faces(); + } + return *_jxz_centre; + } + // Cell Areas + const FieldMetric& cell_area_xlow() const { + if (!_cell_area_xlow.has_value()) { + _compute_cell_area_x(); + } + return _cell_area_xlow; + } + const FieldMetric& cell_area_xhigh() const { + if (!_cell_area_xhigh.has_value()) { + _compute_cell_area_x(); + } + return _cell_area_xhigh; + } + const FieldMetric& cell_area_ylow() const { + if (!_cell_area_ylow.has_value()) { + _compute_cell_area_y(); + } + return _cell_area_ylow; + } + const FieldMetric& cell_area_yhigh() const { + if (!_cell_area_yhigh.has_value()) { + _compute_cell_area_y(); + } + return _cell_area_yhigh; + } + const FieldMetric& cell_area_zlow() const { + if (!_cell_area_zlow.has_value()) { + _compute_cell_area_z(); + } + return _cell_area_zlow; + } + const FieldMetric& cell_area_zhigh() const { + if (!_cell_area_zhigh.has_value()) { + _compute_cell_area_z(); + } + return _cell_area_zhigh; + } + FieldMetric& cell_area_xlow() { + if (!_cell_area_xlow.has_value()) { + _compute_cell_area_x(); + } + return _cell_area_xlow; + } + FieldMetric& cell_area_xhigh() { + if (!_cell_area_xhigh.has_value()) { + _compute_cell_area_x(); + } + return _cell_area_xhigh; + } + FieldMetric& cell_area_ylow() { + if (!_cell_area_ylow.has_value()) { + _compute_cell_area_y(); + } + return _cell_area_ylow; + } + FieldMetric& cell_area_yhigh() { + if (!_cell_area_yhigh.has_value()) { + _compute_cell_area_y(); + } + return _cell_area_yhigh; + } + FieldMetric& cell_area_zlow() { + if (!_cell_area_zlow.has_value()) { + _compute_cell_area_z(); + } + return _cell_area_zlow; + } + FieldMetric& cell_area_zhigh() { + if (!_cell_area_zhigh.has_value()) { + _compute_cell_area_z(); + } + return _cell_area_zhigh; + } + // Cell Volume + const FieldMetric& cell_volume() const { + if (!_cell_volume.has_value()) { + _compute_cell_volume(); + } + return _cell_volume; + } + FieldMetric& cell_volume() { + if (!_cell_volume.has_value()) { + _compute_cell_volume(); + } + return _cell_volume; + } private: mutable std::optional _g_22_ylow, _g_22_yhigh; mutable std::optional _jxz_ylow, _jxz_yhigh, _jxz_centre; + mutable std::optional _cell_area_xlow, _cell_area_xhigh; + mutable std::optional _cell_area_ylow, _cell_area_yhigh; + mutable std::optional _cell_area_zlow, _cell_area_zhigh; + mutable std::optional _cell_volume; + mutable std::optional _jxz_ylow, _jxz_yhigh, _jxz_centre; 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) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 23fa77bbce..0e1d610410 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2092,54 +2092,6 @@ const Coordinates::FieldMetric& Coordinates::g_22_yhigh() const { return g_22_yhigh(); } -Coordinates::FieldMetric& Coordinates::g_22_yhigh() { - return const_cast( - const_cast(this)->g_22_yhigh()); -} - -Coordinates::FieldMetric& Coordinates::g_22_ylow() { - return const_cast( - const_cast(this)->g_22_ylow()); -} - -const Coordinates::FieldMetric& Coordinates::Jxz_ylow() const { - if (!_jxz_ylow.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_ylow; -} -const Coordinates::FieldMetric& Coordinates::Jxz_yhigh() const { - if (!_jxz_yhigh.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_yhigh; -} -const Coordinates::FieldMetric& Coordinates::Jxz() const { - if (!_jxz_centre.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_centre; -} - -Coordinates::FieldMetric& Coordinates::Jxz_ylow() { - if (!_jxz_ylow.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_ylow; -} -Coordinates::FieldMetric& Coordinates::Jxz_yhigh() { - if (!_jxz_yhigh.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_yhigh; -} -Coordinates::FieldMetric& Coordinates::Jxz() { - if (!_jxz_centre.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_centre; -} - 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())); @@ -2172,3 +2124,55 @@ void Coordinates::_compute_Jxz_cell_faces() const { } } } + +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. + ASSERT0(mesh->xstart > 0); + BOUT_FOR(i, _jxz_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 { + const auto Jxz = Jxz(); + if (Jxz.isFci()) { + ASSERT4(isUniform(dx, true, "RGN_ALL")); + ASSERT2(isUniform(dx, false, "RGN_ALL")); + ASSERT4(isUniform(dz, true, "RGN_ALL")); + ASSERT2(isUniform(dz, false, "RGN_ALL")); + _cell_area_ylow = _Jxz_ylow * dx * dz; + _cell_area_yhigh = _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. + ASSERT0(mesh->ystart > 0); + BOUT_FOR(i, _jxz_centre->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, _jxz_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); } From 078322d06b08fc5cd7c9985b3d92acdfee6652e6 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 12 Mar 2026 13:15:12 +0100 Subject: [PATCH 05/22] Fix formatting --- include/bout/coordinates.hxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 402efecd20..13cf6331bc 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -3,16 +3,16 @@ * * ChangeLog * ========= - * + * * 2014-11-10 Ben Dudson * * Created by separating metric from Mesh * - * + * ************************************************************************** * Copyright 2014-2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -402,10 +402,10 @@ private: class TokamakCoordinates : public Coordinates { public: TokamakCoordinates(Mesh *mesh) : Coordinates(mesh) { - + } private: - + }; */ From 7547d9ec7df1616c7b261a0f316207769c40e8c8 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 12 Mar 2026 13:20:34 +0100 Subject: [PATCH 06/22] Fix return statements --- include/bout/coordinates.hxx | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 13cf6331bc..a7a09cfd9d 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -149,86 +149,86 @@ public: if (!_cell_area_xlow.has_value()) { _compute_cell_area_x(); } - return _cell_area_xlow; + return *_cell_area_xlow; } const FieldMetric& cell_area_xhigh() const { if (!_cell_area_xhigh.has_value()) { _compute_cell_area_x(); } - return _cell_area_xhigh; + return *_cell_area_xhigh; } const FieldMetric& cell_area_ylow() const { if (!_cell_area_ylow.has_value()) { _compute_cell_area_y(); } - return _cell_area_ylow; + return *_cell_area_ylow; } const FieldMetric& cell_area_yhigh() const { if (!_cell_area_yhigh.has_value()) { _compute_cell_area_y(); } - return _cell_area_yhigh; + return *_cell_area_yhigh; } const FieldMetric& cell_area_zlow() const { if (!_cell_area_zlow.has_value()) { _compute_cell_area_z(); } - return _cell_area_zlow; + return *_cell_area_zlow; } const FieldMetric& cell_area_zhigh() const { if (!_cell_area_zhigh.has_value()) { _compute_cell_area_z(); } - return _cell_area_zhigh; + return *_cell_area_zhigh; } FieldMetric& cell_area_xlow() { if (!_cell_area_xlow.has_value()) { _compute_cell_area_x(); } - return _cell_area_xlow; + return *_cell_area_xlow; } FieldMetric& cell_area_xhigh() { if (!_cell_area_xhigh.has_value()) { _compute_cell_area_x(); } - return _cell_area_xhigh; + return *_cell_area_xhigh; } FieldMetric& cell_area_ylow() { if (!_cell_area_ylow.has_value()) { _compute_cell_area_y(); } - return _cell_area_ylow; + return *_cell_area_ylow; } FieldMetric& cell_area_yhigh() { if (!_cell_area_yhigh.has_value()) { _compute_cell_area_y(); } - return _cell_area_yhigh; + return *_cell_area_yhigh; } FieldMetric& cell_area_zlow() { if (!_cell_area_zlow.has_value()) { _compute_cell_area_z(); } - return _cell_area_zlow; + return *_cell_area_zlow; } FieldMetric& cell_area_zhigh() { if (!_cell_area_zhigh.has_value()) { _compute_cell_area_z(); } - return _cell_area_zhigh; + return *_cell_area_zhigh; } // Cell Volume const FieldMetric& cell_volume() const { if (!_cell_volume.has_value()) { _compute_cell_volume(); } - return _cell_volume; + return *_cell_volume; } FieldMetric& cell_volume() { if (!_cell_volume.has_value()) { _compute_cell_volume(); } - return _cell_volume; + return *_cell_volume; } private: From 01d776ba2a4621a01481b409fb120e29e2c99459 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 12 Mar 2026 13:20:55 +0100 Subject: [PATCH 07/22] Remove duplicate line --- include/bout/coordinates.hxx | 1 - 1 file changed, 1 deletion(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index a7a09cfd9d..c7b8b83a99 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -238,7 +238,6 @@ private: mutable std::optional _cell_area_ylow, _cell_area_yhigh; mutable std::optional _cell_area_zlow, _cell_area_zhigh; mutable std::optional _cell_volume; - mutable std::optional _jxz_ylow, _jxz_yhigh, _jxz_centre; void _compute_Jxz_cell_faces() const; void _compute_cell_area_x() const; void _compute_cell_area_y() const; From 13fdb8ad50804394c498fe183dff645dabddf4b4 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 13 Mar 2026 12:34:44 +0100 Subject: [PATCH 08/22] Remove non-used Jg() --- include/bout/coordinates.hxx | 2 -- src/mesh/coordinates.cxx | 18 ------------------ 2 files changed, 20 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index c7b8b83a99..85aac3de6c 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -379,14 +379,12 @@ private: /// Cache variable for Grad2_par2 mutable std::map> Grad2_par2_DDY_invSgCache; mutable std::unique_ptr invSgCache{nullptr}; - mutable std::optional JgCache; /// Set the parallel (y) transform from the options file. /// Used in the constructor to create the transform object. void setParallelTransform(Options* options); const FieldMetric& invSg() const; - const FieldMetric& Jg() const; const FieldMetric& Grad2_par2_DDY_invSg(CELL_LOC outloc, const std::string& method) const; diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 0e1d610410..b14bef12be 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1910,24 +1910,6 @@ const Coordinates::FieldMetric& Coordinates::invSg() const { return *invSgCache; } -const Coordinates::FieldMetric& Coordinates::Jg() const { - if (not JgCache.has_value()) { - auto* coords = this; // - // Need to modify yup and ydown fields - Field3D Jg = coords->Jxz(); - Jg.splitParallelSlicesAndAllocate(); - Field3DParallel B = coords->Bxy; - for (size_t j = 0; j < B.numberParallelSlices(); ++j) { - BOUT_FOR(i, B.getRegion("RGN_NOBNDRY")) { - Jg.yup(j)[i.yp(j + 1)] = Jg[i] * B.yup(j)[i.yp(j + 1)] / B[i]; - Jg.ydown(j)[i.ym(j + 1)] = Jg[i] * B.ydown(j)[i.ym(j + 1)] / B[i]; - } - } - JgCache = Jg; - } - return *JgCache; -} - const Coordinates::FieldMetric& Coordinates::Grad2_par2_DDY_invSg(CELL_LOC outloc, const std::string& method) const { if (auto search = Grad2_par2_DDY_invSgCache.find(method); From db66e35f5991db3f8f7892836b6f2a1762ff26de Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 13 Mar 2026 12:35:08 +0100 Subject: [PATCH 09/22] add asField3DParallel stub to Field2D --- include/bout/field2d.hxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx index e068747ee3..e003c4b316 100644 --- a/include/bout/field2d.hxx +++ b/include/bout/field2d.hxx @@ -275,6 +275,9 @@ public: int size() const override { return nx * ny; } + Field2D& asField3DParallel() { return *this; } + const Field2D& asField3DParallel() const { return *this; } + private: /// Internal data array. Handles allocation/freeing of memory Array data; From cca09aceeb63c8949529e907255918e41e368d4d Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 13 Mar 2026 12:35:34 +0100 Subject: [PATCH 10/22] Compile time fixes --- src/mesh/coordinates.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index b14bef12be..4e7a9e1c24 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2121,14 +2121,14 @@ void Coordinates::_compute_cell_area_x() const { } void Coordinates::_compute_cell_area_y() const { - const auto Jxz = Jxz(); - if (Jxz.isFci()) { + Jxz(); + if (_jxz_centre->isFci()) { ASSERT4(isUniform(dx, true, "RGN_ALL")); ASSERT2(isUniform(dx, false, "RGN_ALL")); ASSERT4(isUniform(dz, true, "RGN_ALL")); ASSERT2(isUniform(dz, false, "RGN_ALL")); - _cell_area_ylow = _Jxz_ylow * dx * dz; - _cell_area_yhigh = _Jxz_yhigh * dx * dz; + _cell_area_ylow = _jxz_ylow * dx * dz; + _cell_area_yhigh = _jxz_yhigh * dx * dz; } else { // Field aligned const auto area_centre = sqrt(g_11 * g_33 - SQ(g_13)) * dx * dz; From a0d495459105e2f292ee2e9514cc0ca6f395ed5f Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 13 Mar 2026 12:36:13 +0100 Subject: [PATCH 11/22] Remove no-op check At this time the parallel transform is not yet set, so isFci() return false. --- src/mesh/coordinates.cxx | 25 ------------------------- 1 file changed, 25 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 4e7a9e1c24..4f8183cff0 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1259,31 +1259,6 @@ int Coordinates::calcCovariant(const std::string& region) { output_info.write("\tLocal maximum error in off-diagonal inversion is {:e}\n", maxerr); - if (Bxy.isFci()) { - BoutReal maxError = 0; - Options* options = Options::getRoot(); - auto BJg = Bxy.asField3DParallel() * J / sqrt(g_22.asField3DParallel()); - auto* mesh = localmesh; - for (int p = -mesh->ystart; p <= mesh->ystart; p++) { - if (p == 0) { - continue; - } - BOUT_FOR(i, BJg.getRegion("RGN_NOBNDRY")) { - auto local = BJg[i] / BJg.ynext(p)[i.yp(p)]; - maxError = std::max(std::abs(local - 1), maxError); - } - } - const BoutReal allowedError = (*options)["allowedFluxError"].withDefault(1e-6); - if (maxError < allowedError / 100) { - output_info.write("\tInfo: The maximum flux conservation error is {:e}", maxError); - } else if (maxError < allowedError) { - output_warn.write("\tWarning: The maximum flux conservation error is {:e}", - maxError); - } else { - throw BoutException("Error: The maximum flux conservation error is {:e}", maxError); - } - } - return 0; } From aae5b5dcb3e1dec9b9975e320ab5e6254ee2c511 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 13 Mar 2026 12:51:12 +0100 Subject: [PATCH 12/22] Fixup compilation --- src/mesh/coordinates.cxx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 4f8183cff0..c3a4b727d1 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2088,6 +2088,7 @@ void Coordinates::_compute_cell_area_x() const { _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, _jxz_centre->getRegion("RGN_NOX")) { (*_cell_area_xlow)[i] = 0.5 * (area_centre[i] + area_centre[i.xm()]); @@ -2102,8 +2103,8 @@ void Coordinates::_compute_cell_area_y() const { ASSERT2(isUniform(dx, false, "RGN_ALL")); ASSERT4(isUniform(dz, true, "RGN_ALL")); ASSERT2(isUniform(dz, false, "RGN_ALL")); - _cell_area_ylow = _jxz_ylow * dx * dz; - _cell_area_yhigh = _jxz_yhigh * dx * dz; + _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; @@ -2111,6 +2112,7 @@ void Coordinates::_compute_cell_area_y() const { _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, _jxz_centre->getRegion("RGN_NOY")) { (*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]); From 63ab8b40245ac8b4ebf03acbe46e390ffcfe9c53 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 13 Mar 2026 13:00:09 +0100 Subject: [PATCH 13/22] Switch to ASSERT3 --- src/mesh/coordinates.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index c3a4b727d1..8dd62e11b6 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2099,9 +2099,9 @@ void Coordinates::_compute_cell_area_x() const { void Coordinates::_compute_cell_area_y() const { Jxz(); if (_jxz_centre->isFci()) { - ASSERT4(isUniform(dx, true, "RGN_ALL")); + ASSERT3(isUniform(dx, true, "RGN_ALL")); ASSERT2(isUniform(dx, false, "RGN_ALL")); - ASSERT4(isUniform(dz, true, "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); From 6b66188168a953c0be07df7b67a9015e0ff64d1b Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 16 Mar 2026 09:07:47 +0100 Subject: [PATCH 14/22] Cleanup code --- include/bout/coordinates.hxx | 1 + src/mesh/coordinates.cxx | 19 +++++++++---------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 85aac3de6c..da22fdd126 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -38,6 +38,7 @@ #include #include #include +#include class Mesh; diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 8dd62e11b6..643ff5c6c0 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2014,9 +2014,9 @@ const Coordinates::FieldMetric& Coordinates::g_22_ylow() const { } _g_22_ylow.emplace(emptyFrom(g_22)); //_g_22_ylow->setLocation(CELL_YLOW); - auto mesh = Bxy.getMesh(); + auto* mesh = Bxy.getMesh(); if (Bxy.isFci()) { - if (mesh->get(_g_22_ylow.value(), "g_22_cell_ylow", 0.0, false)) { //, CELL_YLOW)) { + 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 { @@ -2034,10 +2034,9 @@ const Coordinates::FieldMetric& Coordinates::g_22_yhigh() const { } _g_22_yhigh.emplace(emptyFrom(g_22)); //_g_22_yhigh->setLocation(CELL_YHIGH); - auto mesh = Bxy.getMesh(); + auto* mesh = Bxy.getMesh(); if (Bxy.isFci()) { - if (mesh->get(_g_22_yhigh.value(), "g_22_cell_yhigh", 0.0, - false)) { //, CELL_YHIGH)) { + 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 { @@ -2060,13 +2059,13 @@ void Coordinates::_compute_Jxz_cell_faces() const { Coordinates::FieldMetric By_c; Coordinates::FieldMetric By_h; Coordinates::FieldMetric By_l; - if (mesh->get(By_c, "By", 0.0, false, CELL_CENTRE)) { + 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)) { //, CELL_YLOW)) { + 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)) { //, CELL_YHIGH)) { + 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")) { @@ -2088,7 +2087,7 @@ void Coordinates::_compute_cell_area_x() const { _cell_area_xhigh.emplace(emptyFrom(area_centre)); // We cannot setLocation, as that would trigger the computation of staggered // metrics. - auto mesh = Bxy.getMesh(); + auto* mesh = Bxy.getMesh(); ASSERT0(mesh->xstart > 0); BOUT_FOR(i, _jxz_centre->getRegion("RGN_NOX")) { (*_cell_area_xlow)[i] = 0.5 * (area_centre[i] + area_centre[i.xm()]); @@ -2112,7 +2111,7 @@ void Coordinates::_compute_cell_area_y() const { _cell_area_yhigh.emplace(emptyFrom(area_centre)); // We cannot setLocation, as that would trigger the computation of staggered // metrics. - auto mesh = Bxy.getMesh(); + auto* mesh = Bxy.getMesh(); ASSERT0(mesh->ystart > 0); BOUT_FOR(i, _jxz_centre->getRegion("RGN_NOY")) { (*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]); From 4c0c495954667bc071ae43bf7534a19e9bdce561 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 16 Mar 2026 10:48:23 +0100 Subject: [PATCH 15/22] Add unit tests for cell areas and volume --- tests/unit/mesh/test_coordinates.cxx | 71 ++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/tests/unit/mesh/test_coordinates.cxx b/tests/unit/mesh/test_coordinates.cxx index fb400e593e..181282a064 100644 --- a/tests/unit/mesh/test_coordinates.cxx +++ b/tests/unit/mesh/test_coordinates.cxx @@ -271,3 +271,74 @@ TEST_F(CoordinatesTest, NegativeB) { EXPECT_THROW(Coordinates coords(mesh), BoutException); } + +TEST_F(CoordinatesTest, CellAreas) { + + static_cast(bout::globals::mesh) + ->setGridDataSource( + new FakeGridDataSource({{"g_11", 4.0}, {"g_22", 1.0}, {"g_33", 9}})); + + Coordinates coords(mesh); + + EXPECT_TRUE(IsFieldEqual(coords.dx, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.dy, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.dz, 1.0)); + + EXPECT_TRUE(IsFieldEqual(coords.g11, 4.0)); + EXPECT_TRUE(IsFieldEqual(coords.g22, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.g33, 9.0)); + EXPECT_TRUE(IsFieldEqual(coords.g12, 0.0)); + EXPECT_TRUE(IsFieldEqual(coords.g13, 0.0)); + EXPECT_TRUE(IsFieldEqual(coords.g23, 0.0)); + + EXPECT_TRUE(IsFieldEqual(coords.J, 6.0)); + EXPECT_TRUE(IsFieldEqual(coords.Bxy, 1.0)); + + EXPECT_TRUE(IsFieldEqual(coords.cell_area_xlow(), 3.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_xhigh(), 3.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_ylow(), 6.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_yhigh(), 6.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zlow(), 2.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zhigh(), 2.0)); + + EXPECT_TRUE(IsFieldEqual(coords.cell_volume(), 6.0)); +} + +TEST_F(CoordinatesTest, CellAreasUpdate) { + + static_cast(bout::globals::mesh) + ->setGridDataSource(new FakeGridDataSource()); + + Coordinates coords(mesh); + + EXPECT_TRUE(IsFieldEqual(coords.dx, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.dy, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.dz, 1.0)); + + EXPECT_TRUE(IsFieldEqual(coords.g11, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.g22, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.g33, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.g12, 0.0)); + EXPECT_TRUE(IsFieldEqual(coords.g13, 0.0)); + EXPECT_TRUE(IsFieldEqual(coords.g23, 0.0)); + + EXPECT_TRUE(IsFieldEqual(coords.J, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.Bxy, 1.0)); + + coords.cell_area_xlow() *= 2; + coords.cell_area_xhigh() *= 3; + coords.cell_area_ylow() *= 4; + coords.cell_area_yhigh() *= 5; + coords.cell_area_zlow() *= 6; + coords.cell_area_zhigh() *= 7; + coords.cell_volume() *= 8; + + EXPECT_TRUE(IsFieldEqual(coords.cell_area_xlow(), 2.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_xhigh(), 3.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_ylow(), 4.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_yhigh(), 5.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zlow(), 6.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zhigh(), 7.0)); + + EXPECT_TRUE(IsFieldEqual(coords.cell_volume(), 8.0)); +} From fd12c5428e8e393d7cc83745f0a492b3bada55a1 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 16 Mar 2026 10:59:03 +0100 Subject: [PATCH 16/22] Add missing headers --- src/mesh/coordinates.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 643ff5c6c0..8ff6958dcb 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -8,7 +8,9 @@ #include #include #include +#include #include +#include #include #include From 1b83ae843e5d289a7bc9400fc998b862ce0fd907 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 17 Mar 2026 08:51:25 +0100 Subject: [PATCH 17/22] Fixup unit test --- tests/unit/mesh/test_coordinates.cxx | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/tests/unit/mesh/test_coordinates.cxx b/tests/unit/mesh/test_coordinates.cxx index 181282a064..c013d5a261 100644 --- a/tests/unit/mesh/test_coordinates.cxx +++ b/tests/unit/mesh/test_coordinates.cxx @@ -276,7 +276,7 @@ TEST_F(CoordinatesTest, CellAreas) { static_cast(bout::globals::mesh) ->setGridDataSource( - new FakeGridDataSource({{"g_11", 4.0}, {"g_22", 1.0}, {"g_33", 9}})); + new FakeGridDataSource({{"g_11", 4.0}, {"g_22", 1.0}, {"g_33", 9}, {"dz", 1}})); Coordinates coords(mesh); @@ -284,12 +284,9 @@ TEST_F(CoordinatesTest, CellAreas) { EXPECT_TRUE(IsFieldEqual(coords.dy, 1.0)); EXPECT_TRUE(IsFieldEqual(coords.dz, 1.0)); - EXPECT_TRUE(IsFieldEqual(coords.g11, 4.0)); - EXPECT_TRUE(IsFieldEqual(coords.g22, 1.0)); - EXPECT_TRUE(IsFieldEqual(coords.g33, 9.0)); - EXPECT_TRUE(IsFieldEqual(coords.g12, 0.0)); - EXPECT_TRUE(IsFieldEqual(coords.g13, 0.0)); - EXPECT_TRUE(IsFieldEqual(coords.g23, 0.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_11, 4.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_22, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_33, 9.0)); EXPECT_TRUE(IsFieldEqual(coords.J, 6.0)); EXPECT_TRUE(IsFieldEqual(coords.Bxy, 1.0)); From f9cf8d23750ed8f91c714c554c58ac1e92e8bbc0 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 17 Mar 2026 10:12:47 +0100 Subject: [PATCH 18/22] Remove non-working tests --- tests/unit/mesh/test_coordinates.cxx | 9 --------- 1 file changed, 9 deletions(-) diff --git a/tests/unit/mesh/test_coordinates.cxx b/tests/unit/mesh/test_coordinates.cxx index c013d5a261..e6a9e53028 100644 --- a/tests/unit/mesh/test_coordinates.cxx +++ b/tests/unit/mesh/test_coordinates.cxx @@ -280,15 +280,6 @@ TEST_F(CoordinatesTest, CellAreas) { Coordinates coords(mesh); - EXPECT_TRUE(IsFieldEqual(coords.dx, 1.0)); - EXPECT_TRUE(IsFieldEqual(coords.dy, 1.0)); - EXPECT_TRUE(IsFieldEqual(coords.dz, 1.0)); - - EXPECT_TRUE(IsFieldEqual(coords.g_11, 4.0)); - EXPECT_TRUE(IsFieldEqual(coords.g_22, 1.0)); - EXPECT_TRUE(IsFieldEqual(coords.g_33, 9.0)); - - EXPECT_TRUE(IsFieldEqual(coords.J, 6.0)); EXPECT_TRUE(IsFieldEqual(coords.Bxy, 1.0)); EXPECT_TRUE(IsFieldEqual(coords.cell_area_xlow(), 3.0)); From 7ad4d05ca07fbaa97a4e787d036f58f1bc780ae9 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 18 Mar 2026 12:15:27 +0100 Subject: [PATCH 19/22] Fix get region call --- src/mesh/coordinates.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 8ff6958dcb..ec3483d452 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2091,7 +2091,7 @@ void Coordinates::_compute_cell_area_x() const { // metrics. auto* mesh = Bxy.getMesh(); ASSERT0(mesh->xstart > 0); - BOUT_FOR(i, _jxz_centre->getRegion("RGN_NOX")) { + 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()]); } @@ -2115,7 +2115,7 @@ void Coordinates::_compute_cell_area_y() const { // metrics. auto* mesh = Bxy.getMesh(); ASSERT0(mesh->ystart > 0); - BOUT_FOR(i, _jxz_centre->getRegion("RGN_NOY")) { + 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()]); } @@ -2129,7 +2129,7 @@ void Coordinates::_compute_cell_area_z() const { // We cannot setLocation, as that would trigger the computation of staggered // metrics. //ASSERT0(mesh->zstart > 0); - BOUT_FOR(i, _jxz_centre->getRegion("RGN_NOZ")) { + 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()]); } From 6751afda756618998ef8428be004e9f94b4f8d32 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 18 Mar 2026 12:15:39 +0100 Subject: [PATCH 20/22] Fix comments --- tests/unit/mesh/test_coordinates.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/unit/mesh/test_coordinates.cxx b/tests/unit/mesh/test_coordinates.cxx index e6a9e53028..efe051fe6d 100644 --- a/tests/unit/mesh/test_coordinates.cxx +++ b/tests/unit/mesh/test_coordinates.cxx @@ -36,7 +36,7 @@ TEST_F(CoordinatesTest, ZLength) { FieldMetric{0.0}, // g23 FieldMetric{1.0}, // g_11 FieldMetric{1.0}, // g_22 - FieldMetric{1.0}, // g_23 + FieldMetric{1.0}, // g_33 FieldMetric{0.0}, // g_12 FieldMetric{0.0}, // g_13 FieldMetric{0.0}, // g_23 @@ -69,7 +69,7 @@ TEST_F(CoordinatesTest, ZLength3D) { FieldMetric{0.0}, // g23 FieldMetric{1.0}, // g_11 FieldMetric{1.0}, // g_22 - FieldMetric{1.0}, // g_23 + FieldMetric{1.0}, // g_33 FieldMetric{0.0}, // g_12 FieldMetric{0.0}, // g_13 FieldMetric{0.0}, // g_23 @@ -96,7 +96,7 @@ TEST_F(CoordinatesTest, Jacobian) { FieldMetric{0.0}, // g23 FieldMetric{1.0}, // g_11 FieldMetric{1.0}, // g_22 - FieldMetric{1.0}, // g_23 + FieldMetric{1.0}, // g_33 FieldMetric{0.0}, // g_12 FieldMetric{0.0}, // g_13 FieldMetric{0.0}, // g_23 @@ -127,7 +127,7 @@ TEST_F(CoordinatesTest, CalcContravariant) { FieldMetric{0.0}, // g23 FieldMetric{0.0}, // g_11 FieldMetric{0.0}, // g_22 - FieldMetric{0.0}, // g_23 + FieldMetric{0.0}, // g_33 FieldMetric{0.0}, // g_12 FieldMetric{0.0}, // g_13 FieldMetric{0.0}, // g_23 @@ -160,7 +160,7 @@ TEST_F(CoordinatesTest, CalcCovariant) { FieldMetric{0.0}, // g23 FieldMetric{1.0}, // g_11 FieldMetric{1.0}, // g_22 - FieldMetric{1.0}, // g_23 + FieldMetric{1.0}, // g_33 FieldMetric{0.0}, // g_12 FieldMetric{0.0}, // g_13 FieldMetric{0.0}, // g_23 From 042cf65e824b06eabcc5ed91e4c44959c2794ed4 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 18 Mar 2026 12:16:15 +0100 Subject: [PATCH 21/22] Switch test to fake_mesh_fixture Also specify area for testing --- tests/unit/mesh/test_coordinates.cxx | 98 ++++++++++++++++++++-------- 1 file changed, 71 insertions(+), 27 deletions(-) diff --git a/tests/unit/mesh/test_coordinates.cxx b/tests/unit/mesh/test_coordinates.cxx index efe051fe6d..ea233cd9ca 100644 --- a/tests/unit/mesh/test_coordinates.cxx +++ b/tests/unit/mesh/test_coordinates.cxx @@ -274,41 +274,85 @@ TEST_F(CoordinatesTest, NegativeB) { TEST_F(CoordinatesTest, CellAreas) { - static_cast(bout::globals::mesh) - ->setGridDataSource( - new FakeGridDataSource({{"g_11", 4.0}, {"g_22", 1.0}, {"g_33", 9}, {"dz", 1}})); + Coordinates coords{mesh, + FieldMetric{1.0}, // dx + FieldMetric{1.0}, // dy + FieldMetric{1.0}, // dz + FieldMetric{6.0}, // J + FieldMetric{1.0}, // Bxy + FieldMetric{1.0}, // g11 + FieldMetric{1.0}, // g22 + FieldMetric{1.0}, // g33 + FieldMetric{0.0}, // g12 + FieldMetric{0.0}, // g13 + FieldMetric{0.0}, // g23 + FieldMetric{4.0}, // g_11 + FieldMetric{1.0}, // g_22 + FieldMetric{9.0}, // g_33 + FieldMetric{0.0}, // g_12 + FieldMetric{0.0}, // g_13 + FieldMetric{0.0}, // g_23 + FieldMetric{0.0}, // ShiftTorsion + FieldMetric{0.0}}; // IntShiftTorsion + // No call to Coordinates::geometry() needed here - Coordinates coords(mesh); + EXPECT_TRUE(IsFieldEqual(coords.dx, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.dy, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.dz, 1.0)); + + EXPECT_TRUE(IsFieldEqual(coords.g_11, 4.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_22, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_33, 9.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_12, 0.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_13, 0.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_23, 0.0)); + EXPECT_TRUE(IsFieldEqual(coords.J, 6.0)); EXPECT_TRUE(IsFieldEqual(coords.Bxy, 1.0)); - EXPECT_TRUE(IsFieldEqual(coords.cell_area_xlow(), 3.0)); - EXPECT_TRUE(IsFieldEqual(coords.cell_area_xhigh(), 3.0)); - EXPECT_TRUE(IsFieldEqual(coords.cell_area_ylow(), 6.0)); - EXPECT_TRUE(IsFieldEqual(coords.cell_area_yhigh(), 6.0)); - EXPECT_TRUE(IsFieldEqual(coords.cell_area_zlow(), 2.0)); - EXPECT_TRUE(IsFieldEqual(coords.cell_area_zhigh(), 2.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_xlow(), 3.0, "RGN_NOX")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_xhigh(), 3.0, "RGN_NOX")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_ylow(), 6.0, "RGN_NOY")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_yhigh(), 6.0, "RGN_NOY")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zlow(), 2.0, "RGN_NOZ")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zhigh(), 2.0, "RGN_NOZ")); EXPECT_TRUE(IsFieldEqual(coords.cell_volume(), 6.0)); } TEST_F(CoordinatesTest, CellAreasUpdate) { - - static_cast(bout::globals::mesh) - ->setGridDataSource(new FakeGridDataSource()); - - Coordinates coords(mesh); + Coordinates coords{mesh, + FieldMetric{1.0}, // dx + FieldMetric{1.0}, // dy + FieldMetric{1.0}, // dz + FieldMetric{1.0}, // J + FieldMetric{1.0}, // Bxy + FieldMetric{1.0}, // g11 + FieldMetric{1.0}, // g22 + FieldMetric{1.0}, // g33 + FieldMetric{0.0}, // g12 + FieldMetric{0.0}, // g13 + FieldMetric{0.0}, // g23 + FieldMetric{1.0}, // g_11 + FieldMetric{1.0}, // g_22 + FieldMetric{1.0}, // g_33 + FieldMetric{0.0}, // g_12 + FieldMetric{0.0}, // g_13 + FieldMetric{0.0}, // g_23 + FieldMetric{0.0}, // ShiftTorsion + FieldMetric{0.0}}; // IntShiftTorsion + // No call to Coordinates::geometry() needed here EXPECT_TRUE(IsFieldEqual(coords.dx, 1.0)); EXPECT_TRUE(IsFieldEqual(coords.dy, 1.0)); EXPECT_TRUE(IsFieldEqual(coords.dz, 1.0)); - EXPECT_TRUE(IsFieldEqual(coords.g11, 1.0)); - EXPECT_TRUE(IsFieldEqual(coords.g22, 1.0)); - EXPECT_TRUE(IsFieldEqual(coords.g33, 1.0)); - EXPECT_TRUE(IsFieldEqual(coords.g12, 0.0)); - EXPECT_TRUE(IsFieldEqual(coords.g13, 0.0)); - EXPECT_TRUE(IsFieldEqual(coords.g23, 0.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_11, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_22, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_33, 1.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_12, 0.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_13, 0.0)); + EXPECT_TRUE(IsFieldEqual(coords.g_23, 0.0)); EXPECT_TRUE(IsFieldEqual(coords.J, 1.0)); EXPECT_TRUE(IsFieldEqual(coords.Bxy, 1.0)); @@ -321,12 +365,12 @@ TEST_F(CoordinatesTest, CellAreasUpdate) { coords.cell_area_zhigh() *= 7; coords.cell_volume() *= 8; - EXPECT_TRUE(IsFieldEqual(coords.cell_area_xlow(), 2.0)); - EXPECT_TRUE(IsFieldEqual(coords.cell_area_xhigh(), 3.0)); - EXPECT_TRUE(IsFieldEqual(coords.cell_area_ylow(), 4.0)); - EXPECT_TRUE(IsFieldEqual(coords.cell_area_yhigh(), 5.0)); - EXPECT_TRUE(IsFieldEqual(coords.cell_area_zlow(), 6.0)); - EXPECT_TRUE(IsFieldEqual(coords.cell_area_zhigh(), 7.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_xlow(), 2.0, "RGN_NOX")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_xhigh(), 3.0, "RGN_NOX")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_ylow(), 4.0, "RGN_NOY")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_yhigh(), 5.0, "RGN_NOY")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zlow(), 6.0, "RGN_NOZ")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zhigh(), 7.0, "RGN_NOZ")); EXPECT_TRUE(IsFieldEqual(coords.cell_volume(), 8.0)); } From cc172bc524e068f291f6e4be2dbae3b02fa474ce Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 18 Mar 2026 14:12:45 +0100 Subject: [PATCH 22/22] Remove Jxz from public API --- include/bout/coordinates.hxx | 37 ------------------------------------ src/mesh/coordinates.cxx | 2 +- 2 files changed, 1 insertion(+), 38 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index da22fdd126..22a2de7ac7 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -108,43 +108,6 @@ public: const FieldMetric& g_22_yhigh() const; FieldMetric& g_22_ylow(); FieldMetric& g_22_yhigh(); - /// get Jxz at the cell faces or cell centre - const FieldMetric& Jxz_ylow() const { - if (!_jxz_ylow.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_ylow; - } - const FieldMetric& Jxz_yhigh() const { - if (!_jxz_yhigh.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_yhigh; - } - const FieldMetric& Jxz() const { - if (!_jxz_centre.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_centre; - } - FieldMetric& Jxz_ylow() { - if (!_jxz_ylow.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_ylow; - } - FieldMetric& Jxz_yhigh() { - if (!_jxz_yhigh.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_yhigh; - } - FieldMetric& Jxz() { - if (!_jxz_centre.has_value()) { - _compute_Jxz_cell_faces(); - } - return *_jxz_centre; - } // Cell Areas const FieldMetric& cell_area_xlow() const { if (!_cell_area_xlow.has_value()) { diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index ec3483d452..511edfdbe6 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2098,7 +2098,7 @@ void Coordinates::_compute_cell_area_x() const { } void Coordinates::_compute_cell_area_y() const { - Jxz(); + _compute_Jxz_cell_faces(); if (_jxz_centre->isFci()) { ASSERT3(isUniform(dx, true, "RGN_ALL")); ASSERT2(isUniform(dx, false, "RGN_ALL"));