diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index e7ead42ee5..22a2de7ac7 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 @@ -38,6 +38,7 @@ #include #include #include +#include class Mesh; @@ -102,6 +103,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; @@ -256,10 +363,10 @@ private: class TokamakCoordinates : public Coordinates { public: TokamakCoordinates(Mesh *mesh) : Coordinates(mesh) { - + } private: - + }; */ 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; diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index d2634d3f88..511edfdbe6 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -8,7 +8,9 @@ #include #include #include +#include #include +#include #include #include @@ -1878,7 +1880,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); } @@ -1898,7 +1900,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 +2009,130 @@ 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); } 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)); +}