Skip to content

Add cell volume and areas #3326

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ZedThree merged 31 commits into from
Jun 19, 2026
Merged

Add cell volume and areas #3326

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
31 commits
Select commit Hold shift + click to select a range
3f521d0
Add metric info on cell faces for FCI
dschwoerer Jan 20, 2026
3065f88
Do not used preserved names
dschwoerer Jan 30, 2026
1c0e385
Add non-const versions
dschwoerer Jan 30, 2026
5c28542
Add cell areas and cell volumes
dschwoerer Mar 12, 2026
078322d
Fix formatting
dschwoerer Mar 12, 2026
7547d9e
Fix return statements
dschwoerer Mar 12, 2026
01d776b
Remove duplicate line
dschwoerer Mar 12, 2026
13fdb8a
Remove non-used Jg()
dschwoerer Mar 13, 2026
db66e35
add asField3DParallel stub to Field2D
dschwoerer Mar 13, 2026
cca09ac
Compile time fixes
dschwoerer Mar 13, 2026
a0d4954
Remove no-op check
dschwoerer Mar 13, 2026
aae5b5d
Fixup compilation
dschwoerer Mar 13, 2026
63ab8b4
Switch to ASSERT3
dschwoerer Mar 13, 2026
6b66188
Cleanup code
dschwoerer Mar 16, 2026
4c0c495
Add unit tests for cell areas and volume
dschwoerer Mar 16, 2026
fd12c54
Add missing headers
dschwoerer Mar 16, 2026
1b83ae8
Fixup unit test
dschwoerer Mar 17, 2026
f9cf8d2
Remove non-working tests
dschwoerer Mar 17, 2026
7ad4d05
Fix get region call
dschwoerer Mar 18, 2026
6751afd
Fix comments
dschwoerer Mar 18, 2026
042cf65
Switch test to fake_mesh_fixture
dschwoerer Mar 18, 2026
cc172bc
Remove Jxz from public API
dschwoerer Mar 18, 2026
1258b6f
Merge branch 'next' into cell_area_volume
dschwoerer Apr 27, 2026
46536c6
Fixup bad merge
dschwoerer Apr 27, 2026
48a4a68
Merge branch 'next' into cell_area_volume
bendudson Jun 12, 2026
afc091b
Prefer const int
dschwoerer Jun 15, 2026
90c3d49
Ensure optional values are set before using them
dschwoerer Jun 15, 2026
1334f5b
Include used headers
dschwoerer Jun 15, 2026
99a212e
Add missing header
dschwoerer Jun 15, 2026
68b1dc8
Remove Jxz from coordinates
dschwoerer Jun 16, 2026
22ab659
Also fill the boundaries to the extend possible
dschwoerer Jun 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
120 changes: 120 additions & 0 deletions include/bout/coordinates.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down Expand Up @@ -107,6 +109,124 @@ 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();
Comment on lines +115 to +116

@ZedThree ZedThree Mar 13, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably don't need these overloads?

@dschwoerer dschwoerer Mar 16, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We probably want to add the cell length, both in the cell centre, as well as at the faces, i.e. the distance between the cell centres.

That would replace g_22 then, right now they are still needed.

@ZedThree ZedThree Mar 16, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I meant the non-const ones!

@dschwoerer dschwoerer Mar 17, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We need them for normalization. But I am starting to be worried, that that might be a bad idea.

Right now:

coords->g_22 *= norm;
coords->g_22_ylow() *= norm;
coords->g_22_yhigh() *= norm;

gives the correct thing for FCI, but wrong for FA.

coords->g_22_ylow() *= norm;
coords->g_22_yhigh() *= norm;
coords->g_22 *= norm;

or

coords->g_22 *= norm;
if (coord->g_22.isFci()) {
  coords->g_22_ylow() *= norm;
  coords->g_22_yhigh() *= norm;
}

works for both
and

coords->g_22 *= norm;
coords->g_22_ylow();
coords->g_22_yhigh();

works only for FA.

So we probably want to remove write access to the fields, and let BOUT++ handle all of the normalisation. I think that is part of github.com//pull/2873 or #3046 - but that is 900 commits behind next and with its 500 comits has probably more merge conflicts that what I can deal with right now :-(

@ZedThree ZedThree Mar 17, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

#3046 is into #2873 by the way.

We (@tomc271 and myself) can take care of bringing those PRs up-to-date -- can you review them? Starting with #3046

@dschwoerer dschwoerer Mar 18, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was waiting for the wrappers that make coords->dy[i] working again, so not all the code needs to be changed. I thought you still wanted to do that?

// Cell Areas
const FieldMetric& cell_area_xlow() const {
if (!_cell_area_xlow.has_value()) {
_compute_cell_area_x();
}
ASSERT2(_cell_area_xlow.has_value());
Comment thread
dschwoerer marked this conversation as resolved.
return *_cell_area_xlow;
Comment thread
dschwoerer marked this conversation as resolved.
}
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;
Comment thread
dschwoerer marked this conversation as resolved.
}
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;
Comment thread
dschwoerer marked this conversation as resolved.
}
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;
Comment thread
dschwoerer marked this conversation as resolved.
}
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;
Comment thread
dschwoerer marked this conversation as resolved.
}
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;
Comment thread
dschwoerer marked this conversation as resolved.
}
FieldMetric& cell_area_xlow() {
if (!_cell_area_xlow.has_value()) {
_compute_cell_area_x();
}
ASSERT2(_cell_area_xlow.has_value());
return *_cell_area_xlow;
Comment thread
dschwoerer marked this conversation as resolved.
}
FieldMetric& cell_area_xhigh() {
if (!_cell_area_xhigh.has_value()) {
_compute_cell_area_x();
}
ASSERT2(_cell_area_xhigh.has_value());
return *_cell_area_xhigh;
Comment thread
dschwoerer marked this conversation as resolved.
}
FieldMetric& cell_area_ylow() {
if (!_cell_area_ylow.has_value()) {
_compute_cell_area_y();
}
ASSERT2(_cell_area_ylow.has_value());
return *_cell_area_ylow;
Comment thread
dschwoerer marked this conversation as resolved.
}
FieldMetric& cell_area_yhigh() {
if (!_cell_area_yhigh.has_value()) {
_compute_cell_area_y();
}
ASSERT2(_cell_area_yhigh.has_value());
return *_cell_area_yhigh;
Comment thread
dschwoerer marked this conversation as resolved.
}
FieldMetric& cell_area_zlow() {
if (!_cell_area_zlow.has_value()) {
_compute_cell_area_z();
}
ASSERT2(_cell_area_zlow.has_value());
return *_cell_area_zlow;
Comment thread
dschwoerer marked this conversation as resolved.
}
FieldMetric& cell_area_zhigh() {
if (!_cell_area_zhigh.has_value()) {
_compute_cell_area_z();
}
ASSERT2(_cell_area_zhigh.has_value());
return *_cell_area_zhigh;
Comment thread
dschwoerer marked this conversation as resolved.
}
// Cell Volume
const FieldMetric& cell_volume() const {
if (!_cell_volume.has_value()) {
_compute_cell_volume();
}
ASSERT2(_cell_volume.has_value());
return *_cell_volume;
Comment thread
dschwoerer marked this conversation as resolved.
}
FieldMetric& cell_volume() {
if (!_cell_volume.has_value()) {
_compute_cell_volume();
}
ASSERT2(_cell_volume.has_value());
return *_cell_volume;
Comment thread
dschwoerer marked this conversation as resolved.
}

private:
mutable std::optional<FieldMetric> _g_22_ylow, _g_22_yhigh;
Comment thread
dschwoerer marked this conversation as resolved.
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_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;
Expand Down
146 changes: 136 additions & 10 deletions src/mesh/coordinates.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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;
}
}
Expand All @@ -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.;
}
Expand All @@ -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);
Expand Down Expand Up @@ -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");
}
}
Expand Down Expand Up @@ -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>();
Comment thread
dschwoerer marked this conversation as resolved.
(*ptr) = 1.0 / sqrt(g_22);
invSgCache = std::move(ptr);
}
Expand All @@ -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];
Expand Down Expand Up @@ -2003,6 +2009,126 @@ 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));
Comment thread
dschwoerer marked this conversation as resolved.
//_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")) {
Comment thread
dschwoerer marked this conversation as resolved.
_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_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 {
auto* mesh = Bxy.getMesh();
if (g_11.isFci()) {
const auto jxz_centre = sqrt(g_11 * g_33 - SQ(g_13));
auto jxz_ylow = emptyFrom(jxz_centre);
auto jxz_yhigh = emptyFrom(jxz_centre);

auto By_c = emptyFrom(jxz_centre);
auto By_h = emptyFrom(jxz_yhigh);
auto By_l = emptyFrom(jxz_ylow);
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];
}
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.
BOUT_FOR(i, mesh->getRegion("RGN_ALL")) {
if (i.y() > 0) {
(*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]);
} else {
(*_cell_area_ylow)[i] = BoutNaN;
}
if (i.y() < mesh->LocalNy - 1) {
(*_cell_area_yhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.yp()]);
} else {
(*_cell_area_yhigh)[i] = BoutNaN;
}
}
}
}

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);
}
Loading
Loading