From f32587acf9c6c2a540c481c86d135926cdbca29d Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 13 Aug 2026 14:01:35 +0200 Subject: [PATCH 01/10] Ensure we do not break B for FCI --- src/mesh/coordinates.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 5e461602d0..34a9dbc993 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -823,6 +823,7 @@ Coordinates::FieldMetric Coordinates::recalculateJacobian() const { } Coordinates::FieldMetric Coordinates::recalculateBxy() const { + ASSERT2(not J().isFci()); return sqrt(g_22()) / J(); } From 25b79ade8a145a59ec79e2235a69bdda3c41bab4 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 13 Aug 2026 14:06:22 +0200 Subject: [PATCH 02/10] Do not add parallel slices to parallel slices --- src/mesh/coordinates.cxx | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 34a9dbc993..fca933621e 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1278,6 +1278,9 @@ void Coordinates::communicateMetricTensor() { void Coordinates::communicateDz() { localmesh->communicate(dz_); } void Coordinates::splitBxyParallelSlices() { - Bxy_.splitParallelSlices(); - Bxy_.yup() = Bxy_.ydown() = Bxy_; + if (not Bxy_.hasParallelSlices()) { + auto copy = Bxy_; + Bxy_.splitParallelSlices(); + Bxy_.yup() = Bxy_.ydown() = copy; + } } From d1d09e8672b7a80a1beb1eeb2403156e98c6a3e2 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 13 Aug 2026 14:29:36 +0200 Subject: [PATCH 03/10] Add normalisation for metrics --- include/bout/coordinates.hxx | 7 +++++ src/mesh/coordinates.cxx | 55 ++++++++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index f01e3883e0..dc88d20b3c 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -372,6 +372,13 @@ public: void communicateMetricTensor(); +private: + void normaliseMetricTokamak(BoutReal rho_s0, BoutReal Bnorm); + void normaliseMetricFCI(BoutReal rho_s0, BoutReal Bnorm); + +public: + void normaliseMetric(BoutReal rho_s0, BoutReal Bnorm); + void communicateDz(); ///< Coordinate system Jacobian, so volume of cell is J*dx*dy*dz diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index fca933621e..df9daea0d0 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1284,3 +1284,58 @@ void Coordinates::splitBxyParallelSlices() { Bxy_.yup() = Bxy_.ydown() = copy; } } + +void Coordinates::normaliseMetricTokamak(const BoutReal rho_s0, const BoutReal Bnorm) { + contravariantMetricTensor = ContravariantMetricTensor( + contravariantMetricTensor.g11() / SQ(Bnorm * rho_s0), + contravariantMetricTensor.g22() * SQ(rho_s0), + contravariantMetricTensor.g33() * SQ(rho_s0), + contravariantMetricTensor.g12() / Bnorm, contravariantMetricTensor.g13() / Bnorm, + contravariantMetricTensor.g23() * SQ(rho_s0)); + + covariantMetricTensor = CovariantMetricTensor( + covariantMetricTensor.g11() * SQ(Bnorm * rho_s0), + covariantMetricTensor.g22() / SQ(rho_s0), covariantMetricTensor.g33() / SQ(rho_s0), + covariantMetricTensor.g12() * Bnorm, covariantMetricTensor.g13() * Bnorm, + covariantMetricTensor.g23() / SQ(rho_s0)); + + setBxy(Bxy() / Bnorm); + setJ(J() * Bnorm / rho_s0); + invalidateMetricCaches(); +} + +void Coordinates::normaliseMetricFCI(const BoutReal rho_s0, const BoutReal Bnorm) { + BoutReal rhoSQ = SQ(rho_s0); + + // coord->g11.asField3DParallel() *= rhoSQ; + contravariantMetricTensor.map([&](auto f) -> FieldMetric { + if (f.hasParallelSlices()) { + return f.asField3DParallel() * rhoSQ; + }; + return f * rhoSQ; + }); + + //coord->g_11.asField3DParallel() /= rhoSQ; + covariantMetricTensor.map([&](auto f) -> FieldMetric { + if (f.hasParallelSlices()) { + return f.asField3DParallel() / rhoSQ; + }; + return f / rhoSQ; + }); + // coord->J.asField3DParallel() /= rho_s0 * rho_s0 * rho_s0; + if (J().hasParallelSlices()) { + setJ(Field3DParallel{J() * (rho_s0 * rho_s0 * rho_s0)}); + } else { + setJ(J() * (rho_s0 * rho_s0 * rho_s0)); + } + setBxy(Field3DParallel(Bxy() / Bnorm)); + invalidateMetricCaches(); +} + +void Coordinates::normaliseMetric(const BoutReal rho_s0, const BoutReal Bnorm) { + ASSERT2(hasParallelTransform()); + if (Bxy().isFci()) { + return normaliseMetricFCI(rho_s0, Bnorm); + } + return normaliseMetricTokamak(rho_s0, Bnorm); +} From 1dd158300224c8509da5f31f7b1ca1cdb7b30018 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 13 Aug 2026 14:57:26 +0200 Subject: [PATCH 04/10] Fix compilation for 2D metrics --- src/mesh/coordinates.cxx | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index df9daea0d0..7f34758107 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1323,12 +1323,17 @@ void Coordinates::normaliseMetricFCI(const BoutReal rho_s0, const BoutReal Bnorm return f / rhoSQ; }); // coord->J.asField3DParallel() /= rho_s0 * rho_s0 * rho_s0; - if (J().hasParallelSlices()) { - setJ(Field3DParallel{J() * (rho_s0 * rho_s0 * rho_s0)}); + if constexpr (bout::build::use_metric_3d) { + if (J().hasParallelSlices()) { + setJ(Field3DParallel{J() * (rho_s0 * rho_s0 * rho_s0)}); + } else { + setJ(J() * (rho_s0 * rho_s0 * rho_s0)); + } + setBxy(Field3DParallel(Bxy() / Bnorm)); } else { setJ(J() * (rho_s0 * rho_s0 * rho_s0)); + setBxy(Bxy() / Bnorm); } - setBxy(Field3DParallel(Bxy() / Bnorm)); invalidateMetricCaches(); } From 2735dcb4a1f6547e720f9cf5c30ba4cacc7ea45b Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 13 Aug 2026 14:57:47 +0200 Subject: [PATCH 05/10] Suppress unused var warnings --- src/mesh/difops.cxx | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 3c2a15b184..3d72763bef 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -762,10 +762,13 @@ Field3D Laplace(const Field3D& f, CELL_LOC outloc, * Inverse of Laplacian operator in LaplaceXY solver *******************************************************************************/ -Field2D Laplace_perpXY(const Field2D& A, const Field2D& f) { #if BOUT_USE_METRIC_3D +Field2D Laplace_perpXY([[maybe_unused]] const Field2D& A, + [[maybe_unused]] const Field2D& f) { throw BoutException("Coordinates::Laplace_perpXY for 3D metric not implemented"); +} #else +Field2D Laplace_perpXY(const Field2D& A, const Field2D& f) { const auto& coords = *f.getCoordinates(); Field2D result; @@ -821,8 +824,8 @@ Field2D Laplace_perpXY(const Field2D& A, const Field2D& f) { } return result; -#endif } +#endif /******************************************************************************* * b0xGrad_dot_Grad From 6282f56f7d35b091f65ff211893411d8fb825cfb Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 13 Aug 2026 15:35:52 +0200 Subject: [PATCH 06/10] Switch to preprocessor #if --- src/mesh/coordinates.cxx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 7f34758107..b4e937f559 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1323,17 +1323,17 @@ void Coordinates::normaliseMetricFCI(const BoutReal rho_s0, const BoutReal Bnorm return f / rhoSQ; }); // coord->J.asField3DParallel() /= rho_s0 * rho_s0 * rho_s0; - if constexpr (bout::build::use_metric_3d) { - if (J().hasParallelSlices()) { - setJ(Field3DParallel{J() * (rho_s0 * rho_s0 * rho_s0)}); - } else { - setJ(J() * (rho_s0 * rho_s0 * rho_s0)); - } - setBxy(Field3DParallel(Bxy() / Bnorm)); +#if BOUT_USE_METRIC_3D + if (J().hasParallelSlices()) { + setJ(Field3DParallel{J() * (rho_s0 * rho_s0 * rho_s0)}); } else { setJ(J() * (rho_s0 * rho_s0 * rho_s0)); - setBxy(Bxy() / Bnorm); } + setBxy(Field3DParallel(Bxy() / Bnorm)); +#else + setJ(J() * (rho_s0 * rho_s0 * rho_s0)); + setBxy(Bxy() / Bnorm); +#endif invalidateMetricCaches(); } From 75d4fb9fa4e14682cd88d978be04c134f0689184 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 14 Aug 2026 10:43:19 +0200 Subject: [PATCH 07/10] Add dx normalisation for Tokamak geometry --- src/mesh/coordinates.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index b4e937f559..a54fbaa5aa 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1299,6 +1299,7 @@ void Coordinates::normaliseMetricTokamak(const BoutReal rho_s0, const BoutReal B covariantMetricTensor.g12() * Bnorm, covariantMetricTensor.g13() * Bnorm, covariantMetricTensor.g23() / SQ(rho_s0)); + setDx(dx() / (rho_s0 * rho_s0 * Bnorm)); setBxy(Bxy() / Bnorm); setJ(J() * Bnorm / rho_s0); invalidateMetricCaches(); From 993e633f9a4ebf20d1d145bc73fedd80bd75200b Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 14 Aug 2026 10:45:18 +0200 Subject: [PATCH 08/10] Improve comments --- 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 a54fbaa5aa..c8abceaa16 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1308,7 +1308,7 @@ void Coordinates::normaliseMetricTokamak(const BoutReal rho_s0, const BoutReal B void Coordinates::normaliseMetricFCI(const BoutReal rho_s0, const BoutReal Bnorm) { BoutReal rhoSQ = SQ(rho_s0); - // coord->g11.asField3DParallel() *= rhoSQ; + // g??.asField3DParallel() *= rhoSQ; contravariantMetricTensor.map([&](auto f) -> FieldMetric { if (f.hasParallelSlices()) { return f.asField3DParallel() * rhoSQ; @@ -1316,14 +1316,14 @@ void Coordinates::normaliseMetricFCI(const BoutReal rho_s0, const BoutReal Bnorm return f * rhoSQ; }); - //coord->g_11.asField3DParallel() /= rhoSQ; + // g_??.asField3DParallel() /= rhoSQ; covariantMetricTensor.map([&](auto f) -> FieldMetric { if (f.hasParallelSlices()) { return f.asField3DParallel() / rhoSQ; }; return f / rhoSQ; }); - // coord->J.asField3DParallel() /= rho_s0 * rho_s0 * rho_s0; + // J.asField3DParallel() /= rho_s0 * rho_s0 * rho_s0; #if BOUT_USE_METRIC_3D if (J().hasParallelSlices()) { setJ(Field3DParallel{J() * (rho_s0 * rho_s0 * rho_s0)}); From bb4964e993659960848117f308346767c7f3a57f Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 14 Aug 2026 10:45:31 +0200 Subject: [PATCH 09/10] Fix J normalisation --- 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 c8abceaa16..ee80827094 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1326,13 +1326,13 @@ void Coordinates::normaliseMetricFCI(const BoutReal rho_s0, const BoutReal Bnorm // J.asField3DParallel() /= rho_s0 * rho_s0 * rho_s0; #if BOUT_USE_METRIC_3D if (J().hasParallelSlices()) { - setJ(Field3DParallel{J() * (rho_s0 * rho_s0 * rho_s0)}); + setJ(Field3DParallel{J() / (rho_s0 * rho_s0 * rho_s0)}); } else { - setJ(J() * (rho_s0 * rho_s0 * rho_s0)); + setJ(J() / (rho_s0 * rho_s0 * rho_s0)); } setBxy(Field3DParallel(Bxy() / Bnorm)); #else - setJ(J() * (rho_s0 * rho_s0 * rho_s0)); + setJ(J() / (rho_s0 * rho_s0 * rho_s0)); setBxy(Bxy() / Bnorm); #endif invalidateMetricCaches(); From eddc7e1d36ac7c7e4e71a9ab39e8235e1079bf6a Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 17 Aug 2026 12:43:17 +0200 Subject: [PATCH 10/10] Avoid returning void from function call --- src/mesh/coordinates.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index ee80827094..bacbbe81c4 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -1341,7 +1341,8 @@ void Coordinates::normaliseMetricFCI(const BoutReal rho_s0, const BoutReal Bnorm void Coordinates::normaliseMetric(const BoutReal rho_s0, const BoutReal Bnorm) { ASSERT2(hasParallelTransform()); if (Bxy().isFci()) { - return normaliseMetricFCI(rho_s0, Bnorm); + normaliseMetricFCI(rho_s0, Bnorm); + return; } - return normaliseMetricTokamak(rho_s0, Bnorm); + normaliseMetricTokamak(rho_s0, Bnorm); }