diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 7b874ae7aa..e260e40456 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -39,6 +39,7 @@ #include #include #include +#include #include #include @@ -107,6 +108,112 @@ 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(); + } + 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; + 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; diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 836f0e8ce3..a1c4c92480 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -7,20 +7,26 @@ #include "bout/field3d.hxx" #include "bout/field_data.hxx" #include +#include +#include #include #include #include #include #include +#include #include #include #include #include +#include +#include #include #include #include #include +#include #include "invert3x3.hxx" #include "parallel/fci.hxx" @@ -63,9 +69,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 +89,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 +115,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; + int xi = bndry->x + (i * bndry->bx); + 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 +506,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 +1878,7 @@ 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); } @@ -1893,7 +1898,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]; @@ -2003,6 +2008,133 @@ 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 * (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) != 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 * (sqrt(g_22[i]) + 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(); + 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 Coordinates::makeYBoundary(YBndryType type) const { return std::make_shared(type, localoptions, *localmesh); } diff --git a/tests/unit/mesh/test_coordinates.cxx b/tests/unit/mesh/test_coordinates.cxx index fb400e593e..ea233cd9ca 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 @@ -271,3 +271,106 @@ TEST_F(CoordinatesTest, NegativeB) { EXPECT_THROW(Coordinates coords(mesh), BoutException); } + +TEST_F(CoordinatesTest, CellAreas) { + + 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 + + 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, "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) { + 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.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)); + + 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, "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)); +}