From 639808c795879e6948a305c6b99631f00c2ba2ff Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 16 Jun 2026 14:08:32 +0200 Subject: [PATCH 1/9] Use cell area and volumes for div_par_mod --- include/bout/fv_ops.hxx | 59 +++++++++++------------------------------ 1 file changed, 15 insertions(+), 44 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 3a1a4eaee8..9ffc235918 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -710,12 +710,10 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const auto iym = i.ym(); result[i] = (0.25 * (f_in[i] + f_up[iyp]) * (v_in[i] + v_up[iyp]) - * (coord->J[i] + coord->J.yup()[iyp]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.yup()[iyp])) + * coord->cell_area_yhigh()[i] - 0.25 * (f_in[i] + f_down[iym]) * (v_in[i] + v_down[iym]) - * (coord->J[i] + coord->J.ydown()[iym]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.ydown()[iym]))) - / (coord->dy[i] * coord->J[i]); + * coord->cell_area_ylow()[i]) + / coord->cell_volume()[i]; } return result; } @@ -756,57 +754,30 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, // Pre-calculate factors which multiply fluxes #if not(BOUT_USE_METRIC_3D) // For right cell boundaries - const BoutReal common_factor_r = - (coord->J(i, j) + coord->J(i, j + 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); + const BoutReal area_rp = coord->cell_area_yhigh()(i, j); - const BoutReal flux_factor_rc = - common_factor_r / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_rp = - common_factor_r / (coord->dy(i, j + 1) * coord->J(i, j + 1)); - - const BoutReal area_rp = - common_factor_r * coord->dx(i, j + 1) * coord->dz(i, j + 1); + const BoutReal flux_factor_rc = area_rp / coord->cell_volume()(i, j); + const BoutReal flux_factor_rp = area_rp / coord->cell_volume()(i, j + 1); // For left cell boundaries - const BoutReal common_factor_l = - (coord->J(i, j) + coord->J(i, j - 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); + const BoutReal area_lc = coord->cell_area_ylow()(i, j); - const BoutReal flux_factor_lc = - common_factor_l / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_lm = - common_factor_l / (coord->dy(i, j - 1) * coord->J(i, j - 1)); - - const BoutReal area_lc = common_factor_l * coord->dx(i, j) * coord->dz(i, j); + const BoutReal flux_factor_lc = area_lc / coord->cell_volume()(i, j); + const BoutReal flux_factor_lm = area_lc / coord->cell_volume()(i, j - 1); #endif for (int k = 0; k < mesh->LocalNz; k++) { #if BOUT_USE_METRIC_3D // For right cell boundaries - const BoutReal common_factor_r = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); - - const BoutReal flux_factor_rc = - common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_rp = - common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + const BoutReal area_rp = coord->cell_area_yhigh()(i, j, k); - const BoutReal area_rp = - common_factor_r * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k); + const BoutReal flux_factor_rc = area_rp / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_rp = area_rp / coord->cell_volume()(i, j + 1, k); // For left cell boundaries - const BoutReal common_factor_l = - (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - const BoutReal flux_factor_lc = - common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_lm = - common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + const BoutReal area_lc = coord->cell_area_ylow()(i, j, k); - const BoutReal area_lc = - common_factor_l * coord->dx(i, j, k) * coord->dz(i, j, k); + const BoutReal flux_factor_lc = area_lc / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_lm = area_lc / coord->cell_volume()(i, j - 1, k); #endif //////////////////////////////////////////// From 5a2a3dc0e5e19a16c3e3e652717cb917f09a13d7 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 16 Jun 2026 14:11:43 +0200 Subject: [PATCH 2/9] Use cell area and volumes for Div_par_fvv --- include/bout/fv_ops.hxx | 22 +++++++--------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 9ffc235918..6198ca81d6 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -974,24 +974,16 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, for (int k = 0; k < mesh->LocalNz; k++) { // For right cell boundaries - const BoutReal common_factor_r = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + const BoutReal area_r = coord->cell_area_yhigh()(i, j, k); - const BoutReal flux_factor_rc = - common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_rp = - common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1, k); // For left cell boundaries - const BoutReal common_factor_l = - (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - const BoutReal flux_factor_lc = - common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_lm = - common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + const BoutReal area_l = coord->cell_area_ylow()(i, j, k); + + const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1, k); //////////////////////////////////////////// // Reconstruct f at the cell faces From 820e2219e365ee8250eda82d5ae403e712e6abc6 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 16 Jun 2026 14:17:06 +0200 Subject: [PATCH 3/9] Use cell area and volumes for Div_par --- include/bout/fv_ops.hxx | 41 ++++++++++++----------------------------- 1 file changed, 12 insertions(+), 29 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 6198ca81d6..e96dfca1d7 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -382,41 +382,24 @@ Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_sp // Pre-calculate factors which multiply fluxes #if not(BOUT_USE_METRIC_3D) // For right cell boundaries - BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - - const BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); - + const BoutReal area_r = coord->cell_area_yhigh()(i, j); + const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j); + const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1); // For left cell boundaries - common_factor = (coord->J(i, j) + coord->J(i, j - 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - - const BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); + const BoutReal area_l = coord->cell_area_ylow()(i, j); + const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j); + const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1); #endif for (int k = mesh->zstart; k <= mesh->zend; k++) { #if BOUT_USE_METRIC_3D // For right cell boundaries - BoutReal common_factor = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); - - BoutReal flux_factor_rc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); - + const BoutReal area_r = coord->cell_area_yhigh()(i, j, k); + const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1, k); // For left cell boundaries - common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - BoutReal flux_factor_lc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + const BoutReal area_l = coord->cell_area_ylow()(i, j, k); + const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1, k); #endif //////////////////////////////////////////// From 3b0117ea11652b576fbcf8e2dcb2be9163227ce3 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 17 Jun 2026 10:46:58 +0200 Subject: [PATCH 4/9] Also fill the boundaries to the extend possible --- src/mesh/coordinates.cxx | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 39bcf16936..2e1c4ad389 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2099,10 +2099,17 @@ 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. - 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()]); + 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; + } } } } From bc980b1366fc6d3feec0c033e21f0c9464de6334 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 12:07:02 +0200 Subject: [PATCH 5/9] Ensure calls to _compute_... are guarded --- include/bout/coordinates.hxx | 140 ++++++++++++++++++++++++----------- 1 file changed, 98 insertions(+), 42 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 6d803ab28a..a8f4c09c5c 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -116,103 +116,159 @@ public: FieldMetric& g_22_yhigh(); // Cell Areas const FieldMetric& cell_area_xlow() const { - if (!_cell_area_xlow.has_value()) { + if (_cell_area_xlow.has_value()) { + return *_cell_area_xlow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_x(); + ASSERT2(_cell_area_xlow.has_value()); + return *_cell_area_xlow; } - ASSERT2(_cell_area_xlow.has_value()); - return *_cell_area_xlow; } const FieldMetric& cell_area_xhigh() const { - if (!_cell_area_xhigh.has_value()) { + if (_cell_area_xhigh.has_value()) { + return *_cell_area_xhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_x(); + ASSERT2(_cell_area_xhigh.has_value()); + return *_cell_area_xhigh; } - ASSERT2(_cell_area_xhigh.has_value()); - return *_cell_area_xhigh; } const FieldMetric& cell_area_ylow() const { - if (!_cell_area_ylow.has_value()) { + if (_cell_area_ylow.has_value()) { + return *_cell_area_ylow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_y(); + ASSERT2(_cell_area_ylow.has_value()); + return *_cell_area_ylow; } - ASSERT2(_cell_area_ylow.has_value()); - return *_cell_area_ylow; } const FieldMetric& cell_area_yhigh() const { - if (!_cell_area_yhigh.has_value()) { + if (_cell_area_yhigh.has_value()) { + return *_cell_area_yhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_y(); + ASSERT2(_cell_area_yhigh.has_value()); + return *_cell_area_yhigh; } - ASSERT2(_cell_area_yhigh.has_value()); - return *_cell_area_yhigh; } const FieldMetric& cell_area_zlow() const { - if (!_cell_area_zlow.has_value()) { + if (_cell_area_zlow.has_value()) { + return *_cell_area_zlow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_z(); + ASSERT2(_cell_area_zlow.has_value()); + return *_cell_area_zlow; } - ASSERT2(_cell_area_zlow.has_value()); - return *_cell_area_zlow; } const FieldMetric& cell_area_zhigh() const { - if (!_cell_area_zhigh.has_value()) { + if (_cell_area_zhigh.has_value()) { + return *_cell_area_zhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_z(); + ASSERT2(_cell_area_zhigh.has_value()); + return *_cell_area_zhigh; } - ASSERT2(_cell_area_zhigh.has_value()); - return *_cell_area_zhigh; } FieldMetric& cell_area_xlow() { - if (!_cell_area_xlow.has_value()) { + if (_cell_area_xlow.has_value()) { + return *_cell_area_xlow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_x(); + ASSERT2(_cell_area_xlow.has_value()); + return *_cell_area_xlow; } - ASSERT2(_cell_area_xlow.has_value()); - return *_cell_area_xlow; } FieldMetric& cell_area_xhigh() { - if (!_cell_area_xhigh.has_value()) { + if (_cell_area_xhigh.has_value()) { + return *_cell_area_xhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_x(); + ASSERT2(_cell_area_xhigh.has_value()); + return *_cell_area_xhigh; } - ASSERT2(_cell_area_xhigh.has_value()); - return *_cell_area_xhigh; } FieldMetric& cell_area_ylow() { - if (!_cell_area_ylow.has_value()) { + if (_cell_area_ylow.has_value()) { + return *_cell_area_ylow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_y(); + ASSERT2(_cell_area_ylow.has_value()); + return *_cell_area_ylow; } - ASSERT2(_cell_area_ylow.has_value()); - return *_cell_area_ylow; } FieldMetric& cell_area_yhigh() { - if (!_cell_area_yhigh.has_value()) { + if (_cell_area_yhigh.has_value()) { + return *_cell_area_yhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_y(); + ASSERT2(_cell_area_yhigh.has_value()); + return *_cell_area_yhigh; } - ASSERT2(_cell_area_yhigh.has_value()); - return *_cell_area_yhigh; } FieldMetric& cell_area_zlow() { - if (!_cell_area_zlow.has_value()) { + if (_cell_area_zlow.has_value()) { + return *_cell_area_zlow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_z(); + ASSERT2(_cell_area_zlow.has_value()); + return *_cell_area_zlow; } - ASSERT2(_cell_area_zlow.has_value()); - return *_cell_area_zlow; } FieldMetric& cell_area_zhigh() { - if (!_cell_area_zhigh.has_value()) { + if (_cell_area_zhigh.has_value()) { + return *_cell_area_zhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_z(); + ASSERT2(_cell_area_zhigh.has_value()); + return *_cell_area_zhigh; } - ASSERT2(_cell_area_zhigh.has_value()); - return *_cell_area_zhigh; } // Cell Volume const FieldMetric& cell_volume() const { - if (!_cell_volume.has_value()) { + if (_cell_volume.has_value()) { + return *_cell_volume; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_volume(); + ASSERT2(_cell_volume.has_value()); + return *_cell_volume; } - ASSERT2(_cell_volume.has_value()); - return *_cell_volume; } FieldMetric& cell_volume() { - if (!_cell_volume.has_value()) { + if (_cell_volume.has_value()) { + return *_cell_volume; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_volume(); + ASSERT2(_cell_volume.has_value()); + return *_cell_volume; } - ASSERT2(_cell_volume.has_value()); - return *_cell_volume; } private: From e68d1659aca232c6e01278325817dd9c01dff56a Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 12:53:49 +0200 Subject: [PATCH 6/9] Do not return from critical areas --- 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 a8f4c09c5c..2917e39388 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -123,8 +123,8 @@ public: { _compute_cell_area_x(); ASSERT2(_cell_area_xlow.has_value()); - return *_cell_area_xlow; } + return *_cell_area_xlow; } const FieldMetric& cell_area_xhigh() const { if (_cell_area_xhigh.has_value()) { @@ -134,8 +134,8 @@ public: { _compute_cell_area_x(); ASSERT2(_cell_area_xhigh.has_value()); - return *_cell_area_xhigh; } + return *_cell_area_xhigh; } const FieldMetric& cell_area_ylow() const { if (_cell_area_ylow.has_value()) { @@ -145,8 +145,8 @@ public: { _compute_cell_area_y(); ASSERT2(_cell_area_ylow.has_value()); - return *_cell_area_ylow; } + return *_cell_area_ylow; } const FieldMetric& cell_area_yhigh() const { if (_cell_area_yhigh.has_value()) { @@ -156,8 +156,8 @@ public: { _compute_cell_area_y(); ASSERT2(_cell_area_yhigh.has_value()); - return *_cell_area_yhigh; } + return *_cell_area_yhigh; } const FieldMetric& cell_area_zlow() const { if (_cell_area_zlow.has_value()) { @@ -167,8 +167,8 @@ public: { _compute_cell_area_z(); ASSERT2(_cell_area_zlow.has_value()); - return *_cell_area_zlow; } + return *_cell_area_zlow; } const FieldMetric& cell_area_zhigh() const { if (_cell_area_zhigh.has_value()) { @@ -178,8 +178,8 @@ public: { _compute_cell_area_z(); ASSERT2(_cell_area_zhigh.has_value()); - return *_cell_area_zhigh; } + return *_cell_area_zhigh; } FieldMetric& cell_area_xlow() { if (_cell_area_xlow.has_value()) { @@ -189,8 +189,8 @@ public: { _compute_cell_area_x(); ASSERT2(_cell_area_xlow.has_value()); - return *_cell_area_xlow; } + return *_cell_area_xlow; } FieldMetric& cell_area_xhigh() { if (_cell_area_xhigh.has_value()) { @@ -200,8 +200,8 @@ public: { _compute_cell_area_x(); ASSERT2(_cell_area_xhigh.has_value()); - return *_cell_area_xhigh; } + return *_cell_area_xhigh; } FieldMetric& cell_area_ylow() { if (_cell_area_ylow.has_value()) { @@ -211,8 +211,8 @@ public: { _compute_cell_area_y(); ASSERT2(_cell_area_ylow.has_value()); - return *_cell_area_ylow; } + return *_cell_area_ylow; } FieldMetric& cell_area_yhigh() { if (_cell_area_yhigh.has_value()) { @@ -222,8 +222,8 @@ public: { _compute_cell_area_y(); ASSERT2(_cell_area_yhigh.has_value()); - return *_cell_area_yhigh; } + return *_cell_area_yhigh; } FieldMetric& cell_area_zlow() { if (_cell_area_zlow.has_value()) { @@ -233,8 +233,8 @@ public: { _compute_cell_area_z(); ASSERT2(_cell_area_zlow.has_value()); - return *_cell_area_zlow; } + return *_cell_area_zlow; } FieldMetric& cell_area_zhigh() { if (_cell_area_zhigh.has_value()) { @@ -244,8 +244,8 @@ public: { _compute_cell_area_z(); ASSERT2(_cell_area_zhigh.has_value()); - return *_cell_area_zhigh; } + return *_cell_area_zhigh; } // Cell Volume const FieldMetric& cell_volume() const { @@ -256,8 +256,8 @@ public: { _compute_cell_volume(); ASSERT2(_cell_volume.has_value()); - return *_cell_volume; } + return *_cell_volume; } FieldMetric& cell_volume() { if (_cell_volume.has_value()) { @@ -267,8 +267,8 @@ public: { _compute_cell_volume(); ASSERT2(_cell_volume.has_value()); - return *_cell_volume; } + return *_cell_volume; } private: From 1672208242c2bf2c6428129c6945efa49cb24f3d Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 14:10:57 +0200 Subject: [PATCH 7/9] Move critical section to compute functions --- include/bout/coordinates.hxx | 98 +++++------------ src/mesh/coordinates.cxx | 202 ++++++++++++++++++++--------------- 2 files changed, 146 insertions(+), 154 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 2917e39388..a21f62cd4d 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -119,132 +119,96 @@ public: if (_cell_area_xlow.has_value()) { return *_cell_area_xlow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_x(); - ASSERT2(_cell_area_xlow.has_value()); - } + _compute_cell_area_x(); + ASSERT2(_cell_area_xlow.has_value()); return *_cell_area_xlow; } const FieldMetric& cell_area_xhigh() const { if (_cell_area_xhigh.has_value()) { return *_cell_area_xhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_x(); - ASSERT2(_cell_area_xhigh.has_value()); - } + _compute_cell_area_x(); + ASSERT2(_cell_area_xhigh.has_value()); return *_cell_area_xhigh; } const FieldMetric& cell_area_ylow() const { if (_cell_area_ylow.has_value()) { return *_cell_area_ylow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_y(); - ASSERT2(_cell_area_ylow.has_value()); - } + _compute_cell_area_y(); + ASSERT2(_cell_area_ylow.has_value()); return *_cell_area_ylow; } const FieldMetric& cell_area_yhigh() const { if (_cell_area_yhigh.has_value()) { return *_cell_area_yhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_y(); - ASSERT2(_cell_area_yhigh.has_value()); - } + _compute_cell_area_y(); + ASSERT2(_cell_area_yhigh.has_value()); return *_cell_area_yhigh; } const FieldMetric& cell_area_zlow() const { if (_cell_area_zlow.has_value()) { return *_cell_area_zlow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_z(); - ASSERT2(_cell_area_zlow.has_value()); - } + _compute_cell_area_z(); + ASSERT2(_cell_area_zlow.has_value()); return *_cell_area_zlow; } const FieldMetric& cell_area_zhigh() const { if (_cell_area_zhigh.has_value()) { return *_cell_area_zhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_z(); - ASSERT2(_cell_area_zhigh.has_value()); - } + _compute_cell_area_z(); + ASSERT2(_cell_area_zhigh.has_value()); return *_cell_area_zhigh; } FieldMetric& cell_area_xlow() { if (_cell_area_xlow.has_value()) { return *_cell_area_xlow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_x(); - ASSERT2(_cell_area_xlow.has_value()); - } + _compute_cell_area_x(); + ASSERT2(_cell_area_xlow.has_value()); return *_cell_area_xlow; } FieldMetric& cell_area_xhigh() { if (_cell_area_xhigh.has_value()) { return *_cell_area_xhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_x(); - ASSERT2(_cell_area_xhigh.has_value()); - } + _compute_cell_area_x(); + ASSERT2(_cell_area_xhigh.has_value()); return *_cell_area_xhigh; } FieldMetric& cell_area_ylow() { if (_cell_area_ylow.has_value()) { return *_cell_area_ylow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_y(); - ASSERT2(_cell_area_ylow.has_value()); - } + _compute_cell_area_y(); + ASSERT2(_cell_area_ylow.has_value()); return *_cell_area_ylow; } FieldMetric& cell_area_yhigh() { if (_cell_area_yhigh.has_value()) { return *_cell_area_yhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_y(); - ASSERT2(_cell_area_yhigh.has_value()); - } + _compute_cell_area_y(); + ASSERT2(_cell_area_yhigh.has_value()); return *_cell_area_yhigh; } FieldMetric& cell_area_zlow() { if (_cell_area_zlow.has_value()) { return *_cell_area_zlow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_z(); - ASSERT2(_cell_area_zlow.has_value()); - } + _compute_cell_area_z(); + ASSERT2(_cell_area_zlow.has_value()); return *_cell_area_zlow; } FieldMetric& cell_area_zhigh() { if (_cell_area_zhigh.has_value()) { return *_cell_area_zhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_z(); - ASSERT2(_cell_area_zhigh.has_value()); - } + _compute_cell_area_z(); + ASSERT2(_cell_area_zhigh.has_value()); return *_cell_area_zhigh; } // Cell Volume @@ -252,22 +216,16 @@ public: if (_cell_volume.has_value()) { return *_cell_volume; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_volume(); - ASSERT2(_cell_volume.has_value()); - } + _compute_cell_volume(); + ASSERT2(_cell_volume.has_value()); return *_cell_volume; } FieldMetric& cell_volume() { if (_cell_volume.has_value()) { return *_cell_volume; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_volume(); - ASSERT2(_cell_volume.has_value()); - } + _compute_cell_volume(); + ASSERT2(_cell_volume.has_value()); return *_cell_volume; } diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 057ffa65b5..afb63f981d 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2015,17 +2015,23 @@ const Coordinates::FieldMetric& Coordinates::g_22_ylow() const { if (_g_22_ylow.has_value()) { return *_g_22_ylow; } - _g_22_ylow.emplace(emptyFrom(g_22)); - //_g_22_ylow->setLocation(CELL_YLOW); - auto* mesh = Bxy.getMesh(); - if (Bxy.isFci()) { - if (mesh->get(_g_22_ylow.value(), "g_22_cell_ylow", 0.0, false) != 0) { - throw BoutException("The grid file does not contain `g_22_cell_ylow`."); - } - } else { - ASSERT0(mesh->ystart > 0); - BOUT_FOR(i, g_22.getRegion("RGN_NOY")) { - _g_22_ylow.value()[i] = SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.ym()]))); + BOUT_OMP_SAFE(critical) + { + if (!_g_22_ylow.has_value()) { + _g_22_ylow.emplace(emptyFrom(g_22)); + //_g_22_ylow->setLocation(CELL_YLOW); + auto* mesh = Bxy.getMesh(); + if (Bxy.isFci()) { + if (mesh->get(_g_22_ylow.value(), "g_22_cell_ylow", 0.0, false) != 0) { + throw BoutException("The grid file does not contain `g_22_cell_ylow`."); + } + } else { + ASSERT0(mesh->ystart > 0); + BOUT_FOR(i, g_22.getRegion("RGN_NOY")) { + _g_22_ylow.value()[i] = + SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.ym()]))); + } + } } } return g_22_ylow(); @@ -2035,99 +2041,127 @@ 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)); - 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()]))); + BOUT_OMP_SAFE(critical) + { + if (!_g_22_yhigh.has_value()) { + _g_22_yhigh.emplace(emptyFrom(g_22)); + 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 FieldMetric 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()]); + BOUT_OMP_SAFE(critical) + { + if (!_cell_area_xlow.has_value()) { + const FieldMetric 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 FieldMetric 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 FieldMetric 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()]); + BOUT_OMP_SAFE(critical) + { + if (!_cell_area_ylow.has_value()) { + auto* mesh = Bxy.getMesh(); + if (g_11.isFci()) { + const FieldMetric 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 { - (*_cell_area_yhigh)[i] = BoutNaN; + // Field aligned + const FieldMetric 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 FieldMetric 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. - 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()]); + BOUT_OMP_SAFE(critical) + { + if (!_cell_volume.has_value()) { + const FieldMetric 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. + 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); } +void Coordinates::_compute_cell_volume() const { + BOUT_OMP_SAFE(critical) + { + if (!_cell_volume.has_value()) { + _cell_volume.emplace(J * dx * dy * dz); + } + } +} std::shared_ptr Coordinates::makeYBoundary(YBndryType type) const { return std::make_shared(type, localoptions, *localmesh); From b89867dd3a54d17d1c7b4174a79baabebec44c54 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 14:11:38 +0200 Subject: [PATCH 8/9] Add DDY template for expression --- include/bout/derivs.hxx | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/include/bout/derivs.hxx b/include/bout/derivs.hxx index a8d9279378..92f094f550 100644 --- a/include/bout/derivs.hxx +++ b/include/bout/derivs.hxx @@ -86,6 +86,13 @@ Field3D DDY(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY"); +template +std::enable_if_t && !bout::utils::is_Field3D_v, const Field3D> +DDY(const E& expr, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", + const std::string& region = "RGN_NOBNDRY") { + return DDY(Field3D{expr}, outloc, method, region); +} + /// Calculate first partial derivative in Y /// /// \f$\partial / \partial y\f$ From 26876b7d6bf786b6c47dcc9ebd06b009939db8fc Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 26 Jun 2026 14:09:49 +0200 Subject: [PATCH 9/9] Add missing used header --- include/bout/derivs.hxx | 1 + 1 file changed, 1 insertion(+) diff --git a/include/bout/derivs.hxx b/include/bout/derivs.hxx index 92f094f550..14bc5c2824 100644 --- a/include/bout/derivs.hxx +++ b/include/bout/derivs.hxx @@ -35,6 +35,7 @@ #include "bout/vector3d.hxx" #include "bout/bout_types.hxx" +#include ////////// FIRST DERIVATIVES //////////