Refactor `Coordinates` with more encapsulation by ZedThree · Pull Request #3425 · boutproject/BOUT-dev

ZedThree · 2026-07-08T17:08:02Z

Replaces #2873

Refactor Coordinates with a more encapsulated interface. Derived quantities
are now automatically updated when the metric components are changed.

The most impactful change is the conversion of the metric components and grid
sizes to const& methods:

- coords.dx
+ coords.dx()

Among other things, this change allows us to defer computation of many of the
derived quantities (such as Christoffel symbols) and remove the
Coordinates::geometry method.

We can now also move all of the derivative operators out of Coordinates into free functions (we could've done most of these already, but we needed Coordinates::DD* for the metric derivatives, which we now defer computing)

github-actions

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 141. Check the log or trigger a new build to see more.

github-actions · 2026-07-08T17:37:53Z

    } else {
      // Use full expression with all terms
-      result = b0xGrad_dot_Grad(p, f) / coord->Bxy;
+      result = b0xGrad_dot_Grad(p, f) / coord->Bxy();


warning: no header providing "b0xGrad_dot_Grad" is directly included [misc-include-cleaner]

result = b0xGrad_dot_Grad(p, f) / coord->Bxy(); ^

github-actions · 2026-07-08T17:37:53Z

@@ -331,7 +331,7 @@ class CWM : public PhysicsModel {
      result = VDDZ(-DDX(p), f);


warning: no matching function for call to 'VDDZ' [clang-diagnostic-error]

result = VDDZ(-DDX(p), f); ^

Additional context

include/bout/derivs.hxx:472: candidate function not viable: no known conversion from 'const Field3D' to 'const Field2D' for 2nd argument

^

include/bout/derivs.hxx:455: candidate function not viable: no known conversion from 'BinaryExpr<Field2D, Constant, Field2D, bout::op::Mul>' (aka 'BinaryExpr<Field2D, Constant, Field2D, bout::op::Mul>') to 'const Field3D' for 1st argument

^

include/bout/derivs.hxx:489: candidate function not viable: no known conversion from 'BinaryExpr<Field2D, Constant, Field2D, bout::op::Mul>' (aka 'BinaryExpr<Field2D, Constant, Field2D, bout::op::Mul>') to 'const Field3D' for 1st argument

^

github-actions · 2026-07-08T17:37:53Z

    } else {
      // Use full expression with all terms
-      result = b0xGrad_dot_Grad(p, f) / coord->Bxy;
+      result = b0xGrad_dot_Grad(p, f) / coord->Bxy();


warning: no header providing "b0xGrad_dot_Grad" is directly included [misc-include-cleaner]

result = b0xGrad_dot_Grad(p, f) / coord->Bxy(); ^

github-actions · 2026-07-08T17:37:53Z

    } else {
      // Use full expression with all terms
-      result = b0xGrad_dot_Grad(p, f) / coord->Bxy;
+      result = b0xGrad_dot_Grad(p, f) / coord->Bxy();


warning: no header providing "b0xGrad_dot_Grad" is directly included [misc-include-cleaner]

result = b0xGrad_dot_Grad(p, f) / coord->Bxy(); ^

github-actions · 2026-07-08T17:37:53Z


-      hyper_mu_x = hyperviscos * metric->g_11 * SQ(metric->dx)
-                   * abs(metric->g11 * D2DX2(U)) / (abs(U) + 1e-3);
+      hyper_mu_x = hyperviscos * metric->g_11() * SQ(metric->dx())


warning: no header providing "SQ" is directly included [misc-include-cleaner]

hyper_mu_x = hyperviscos * metric->g_11() * SQ(metric->dx()) ^

github-actions

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 116. Check the log or trigger a new build to see more.

dschwoerer

LGTM

Refactor `Coordinates` with a more encapsulated interface. Derived quantities are now automatically updated when the metric components are changed. The most impactful change is the conversion of the metric components and grid sizes to `const&` methods: ```diff - coords.dx + coords.dx() ``` Among other things, this change allows us to defer computation of many of the derived quantities (such as Christoffel symbols) and remove the `Coordinates::geometry` method.

Covariant metrics have `_`. The new code is consistent with e.g. metric reads around line 450.

`new_coordinates` should be moved into `coords_map` after it has been recalculated, not before.

`CoordinatesAccessor` and non-uniform geometry quantities depend on cell spacing. Reset cache so they are updated when needed. Bxy is not calculated on demand, so is calculated in `set*MetricTensor` methods.

z index should be in range. Indexed accessors return BoutReal, not reference to field.

ZedThree · 2026-07-09T15:27:21Z

Thanks for the fixes @bendudson!

github-actions

clang-tidy made some suggestions

There were too many comments to post at once. Showing the first 25 out of 32. Check the log or trigger a new build to see more.

github-actions · 2026-07-09T15:51:11Z

 #include <bout/fft.hxx>
+#include <bout/field3d.hxx>
+#include <bout/output.hxx>
 #include <bout/sys/timer.hxx>


warning: included header output.hxx is not used directly [misc-include-cleaner]

Suggested change

#include <bout/sys/timer.hxx>

#include <bout/sys/timer.hxx>

github-actions · 2026-07-09T15:51:12Z


  // NOTE: For now the X-Z terms are omitted, so check that they are small
-  ASSERT2(max(abs(coord->g13)) < 1e-5);
+  ASSERT2(max(abs(coord->g13())) < 1e-5);


warning: no header providing "ASSERT2" is directly included [misc-include-cleaner]

ASSERT2(max(abs(coord->g13())) < 1e-5); ^

github-actions · 2026-07-09T15:51:12Z


  // NOTE: For now the X-Z terms are omitted, so check that they are small
-  ASSERT2(max(abs(coord->g13)) < 1e-5);
+  ASSERT2(max(abs(coord->g13())) < 1e-5);


warning: no header providing "abs" is directly included [misc-include-cleaner]

src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.cxx:2:

+ #include <cstdlib>

github-actions · 2026-07-09T15:51:12Z

-  const Field2D dy = coord->dy;
-  const Field2D J = coord->J;
-  const Field2D g_22 = coord->g_22;
+  const Field2D dy = coord->dy();


warning: no header providing "Field2D" is directly included [misc-include-cleaner]

const Field2D dy = coord->dy(); ^

github-actions · 2026-07-09T15:51:12Z

-#define COPY_STRIPE1(symbol)        \
-  if (coords->symbol.isAllocated()) \
-    data[stripe_size * ind.ind + static_cast<int>(Offset::symbol)] = coords->symbol[ind];
+#define COPY_STRIPE1(symbol)                                         \


warning: function-like macro 'COPY_STRIPE1' used; consider a 'constexpr' template function [cppcoreguidelines-macro-usage]

#define COPY_STRIPE1(symbol) \ ^

github-actions · 2026-07-09T15:51:12Z

-  dx4 = SQ(SQ(coords->dx));
-  dy4 = SQ(SQ(coords->dy));
-  dz4 = SQ(SQ(coords->dz));
+  dx4 = SQ(SQ(coords->dx()));


warning: no header providing "SQ" is directly included [misc-include-cleaner]

dx4 = SQ(SQ(coords->dx())); ^

github-actions · 2026-07-09T15:51:13Z

-  dz4 = SQ(SQ(coords->dz));
+  dx4 = SQ(SQ(coords->dx()));
+  dy4 = SQ(SQ(coords->dy()));
+  dz4 = SQ(SQ(coords->dz()));


warning: assigning to 'BoutReal' (aka 'double') from incompatible type 'BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Square>, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Square>, bout::op::Square>' [clang-diagnostic-error]

dz4 = SQ(SQ(coords->dz())); ^

github-actions · 2026-07-09T15:51:13Z


  SAVE_REPEAT(Ve);

  output.write("dx = {:e}, dy = {:e}, dz = {:e}\n", coords->dx(2, 2), coords->dy(2, 2),


warning: call to consteval function 'fmt::fstring<const double &, const double &, const Field2D &>::fstring<33UL>' is not a constant expression [clang-diagnostic-error]

output.write("dx = {:e}, dy = {:e}, dz = {:e}\n", coords->dx(2, 2), coords->dy(2, 2), ^

Additional context

externalpackages/fmt/include/fmt/base.h:1766: non-constexpr function 'report_error' cannot be used in a constant expression

report_error(message); ^

externalpackages/fmt/include/fmt/base.h:1650: in call to 'handler.on_error(&"unknown format specifier"[0])'

return handler.on_error("unknown format specifier"), end; ^

externalpackages/fmt/include/fmt/base.h:1663: in call to 'parse_replacement_field<char, fmt::detail::format_string_checker<char, 3, 0, false> &>(&"dx = {:e}, dy = {:e}, dz = {:e}\n"[29], &"dx = {:e}, dy = {:e}, dz = {:e}\n"[32], checker(s, arg_pack()))'

begin = p = parse_replacement_field(p - 1, end, handler); ^

externalpackages/fmt/include/fmt/base.h:2758: in call to 'parse_format_string<char, fmt::detail::format_string_checker<char, 3, 0, false>>({&"dx = {:e}, dy = {:e}, dz = {:e}\n"[0], 32}, checker(s, arg_pack()))'

if (FMT_USE_CONSTEVAL) parse_format_string<char>(s, checker(s, arg_pack())); ^

tests/MMS/GBS/gbs.cxx:301: in call to 'fstring<33UL>("dx = {:e}, dy = {:e}, dz = {:e}\n")'

output.write("dx = {:e}, dy = {:e}, dz = {:e}\n", coords->dx(2, 2), coords->dy(2, 2), ^

externalpackages/fmt/include/fmt/base.h:662: declared here

FMT_NORETURN FMT_API void report_error(const char* message); ^

github-actions · 2026-07-09T15:51:13Z

  if (ionvis) {
    Field3D tau_i = Omega_ci * tau_i0 * pow(Ti, 1.5) / Ne;
-    Gi = -(0.96 * Ti * Ne * tau_i) * (2. * Grad_par(Vi) + C(phi) / coords->Bxy);
+    Gi = -(0.96 * Ti * Ne * tau_i) * (2. * Grad_par(Vi) + C(phi) / coords->Bxy());


warning: no header providing "Grad_par" is directly included [misc-include-cleaner]

Gi = -(0.96 * Ti * Ne * tau_i) * (2. * Grad_par(Vi) + C(phi) / coords->Bxy()); ^

github-actions · 2026-07-09T15:51:13Z

    const auto& I = tokamak_coords.I_unnormalised;

-    B0 = tokamak_coords.Bxy;
+    B0 = tokamak_coords.Bxy();


warning: no matching function for call to object of type 'const FieldMetric' (aka 'const Field2D') [clang-diagnostic-error]

B0 = tokamak_coords.Bxy(); ^

Additional context

include/bout/field2d.hxx:351: candidate function not viable: requires single argument 'i', but no arguments were provided

BOUT_DEVICE inline BoutReal operator()(int i) { return View()(i); } ^

include/bout/field2d.hxx:352: candidate function not viable: requires single argument 'i', but no arguments were provided

BOUT_DEVICE inline BoutReal operator()(int i) const { return View()(i); } ^

include/bout/field2d.hxx:235: candidate function not viable: requires 2 arguments, but 0 were provided

inline BoutReal& operator()(int jx, int jy) { ^

include/bout/field2d.hxx:249: candidate function not viable: requires 2 arguments, but 0 were provided

inline const BoutReal& operator()(int jx, int jy) const { ^

include/bout/field2d.hxx:268: candidate function not viable: requires 3 arguments, but 0 were provided

BoutReal& operator()(int jx, int jy, [[maybe_unused]] int jz) { ^

include/bout/field2d.hxx:271: candidate function not viable: requires 3 arguments, but 0 were provided

const BoutReal& operator()(int jx, int jy, [[maybe_unused]] int jz) const { ^

Convert input to Field3DParallel before dividing by Bxy, to ensure that yup/ydown fields are retained.

github-actions

clang-tidy made some suggestions

github-actions · 2026-07-09T19:24:23Z

    }

-    ddt(U) -= 10 * (SQ(SQ(coords->dx)) * D4DX4(U) + SQ(SQ(coords->dz)) * D4DZ4(U));
+    ddt(U) -= 10 * (SQ(SQ(coords->dx())) * D4DX4(U) + SQ(SQ(coords->dz())) * D4DZ4(U));


warning: no header providing "SQ" is directly included [misc-include-cleaner]

ddt(U) -= 10 * (SQ(SQ(coords->dx())) * D4DX4(U) + SQ(SQ(coords->dz())) * D4DZ4(U)); ^

github-actions · 2026-07-09T19:24:23Z

    // Central differencing
-    ddt(f) = DDX(g, CELL_CENTRE); // + 20*SQ(coord->dx)*D2DX2(f);
-    ddt(g) = DDX(f, CELL_XLOW);   // + 20*SQ(coord->dx)*D2DX2(g);
+    ddt(f) = DDX(g, CELL_CENTRE); // + 20*SQ(coord->dx())*D2DX2(f);


warning: no header providing "CELL_CENTRE" is directly included [misc-include-cleaner]

ddt(f) = DDX(g, CELL_CENTRE); // + 20*SQ(coord->dx())*D2DX2(f); ^

github-actions · 2026-07-09T19:24:23Z


 // just define a macro for V_E dot Grad
-#define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy)
+#define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy())


warning: function-like macro 'vE_Grad' used; consider a 'constexpr' template function [cppcoreguidelines-macro-usage]

#define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy()) ^

github-actions · 2026-07-09T19:24:23Z


 // just define a macro for V_E dot Grad
-#define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy)
+#define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy())


warning: no header providing "b0xGrad_dot_Grad" is directly included [misc-include-cleaner]

#define vE_Grad(f, p) (b0xGrad_dot_Grad(p, f) / coord->Bxy()) ^

github-actions · 2026-07-09T19:24:23Z


    // DENSITY EQUATION
-    ddt(Ni) = -b0xGrad_dot_Grad(phi, Ni0) / coord->Bxy;
+    ddt(Ni) = -b0xGrad_dot_Grad(phi, Ni0) / coord->Bxy();


warning: no header providing "b0xGrad_dot_Grad" is directly included [misc-include-cleaner]

tests/integrated/test-interchange-instability/2fluid.cxx:5:

- #include <bout/bout.hxx> + #include "bout/difops.hxx" + #include <bout/bout.hxx>

github-actions · 2026-07-09T19:24:23Z


    // DENSITY EQUATION
-    ddt(Ni) = -b0xGrad_dot_Grad(phi, Ni0) / coord->Bxy;
+    ddt(Ni) = -b0xGrad_dot_Grad(phi, Ni0) / coord->Bxy();


warning: no header providing "ddt" is directly included [misc-include-cleaner]

tests/integrated/test-interchange-instability/2fluid.cxx:5:

- #include <bout/bout.hxx> + #include "bout/field3d.hxx" + #include <bout/bout.hxx>

github-actions · 2026-07-09T19:24:23Z


    // VORTICITY
-    ddt(rho) = 2.0 * coord->Bxy * b0xcv * Grad(pei);
+    ddt(rho) = 2.0 * coord->Bxy() * b0xcv * Grad(pei);


warning: no header providing "Grad" is directly included [misc-include-cleaner]

tests/integrated/test-interchange-instability/2fluid.cxx:5:

- #include <bout/bout.hxx> + #include "bout/vecops.hxx" + #include <bout/bout.hxx>

Should now be consistent with `next` branch. Avoids calculation of parallel slices for derived metric quantities.

Makes lazy getters thread-safe by adding named critical sections. Named so that nested sections (code in critical section calling a function that also contains a critical section) doesn't lead to deadlock.

bendudson · 2026-07-09T20:42:47Z

Don't know what's going on here. Tests pass except with OpenMP.

ZedThree · 2026-07-24T17:36:37Z

I can't reproduce this crash locally and because it's coming from inside an OpenMP region, the stacktrace is not super helpful:

#3 in XZLagrange4pt::interpolate(Field3D const&, std::string const&) const [clone ._omp_fn.0]
   at /home/runner/work/BOUT-dev/BOUT-dev/src/mesh/interpolation/lagrange_4pt_xz.cxx:137:5
     135:     // Then in X
     136:     f_interp(x, y_next, z) = lagrange_4pt(xvals, t_x(x, y, z));
   > 137:     ASSERT2(std::isfinite(f_interp(x, y_next, z)));
              ^
     138:   }
     139:   const auto region2 = y_offset != 0 ? fmt::format("RGN_YPAR_{:+d}", y_offset) : region;
#4 at /lib/x86_64-linux-gnu/libgomp.so.1
#5 in start_thread
   at ./nptl/pthread_create.c:447:8
#6 in clone3
   at ./misc/../sysdeps/unix/sysv/linux/x86_64/clone3.S:78:0

Some testing in CI gives some hints:

One run had:

====== Exception thrown ======
not finite at (2, 2, 9)
xvals: 1.1244546177789345e+257, 0.004739247964904323, 0.005235504377539103, -nan
t_x: 1

Another run:

====== Exception thrown ======
not finite at (2, 0, 0)
xvals: -nan, 0.004739247964904324, 0.005235504377539103, 3.8836107342926553e+166
t_x: 1

====== Exception thrown ======
not finite at (2, 5, 0)
xvals: -nan, 0.004739247964904324, 0.005235504377539103, -5.520902839380176e+165
t_x: 1

Another run, with more info:

====== Exception thrown ======
not finite at (2, 4, 0)
xvals: 1.0875872532615936e+257, 0.004739247964904324, 0.005235504377539103, -nan
t_x: 1
jx vals: (0, 1, 2, 3)
jz vals: (15, 0, 1, 2)
f vals (x->, z^): [
  -7.0929115670299516e-161, 0.004739247964904324, 0.005235504377539103, 1.6244820466831905e+189
  -2.3081895383550503e-149, 0.004739247964904324, 0.005235504377539103, 4.013589778025809e+189
  -8.37630765434779e-169, 0.004739247964904324, 0.005235504377539103, 3.224079296911227e-85
  -1.7740686685036805e+258, 0.004739247964904324, 0.005235504377539103, -nan
]

It seems to be that all the values in x=0 and x=3, i.e. the guard cells (the grid has nx=4 and MXG=1), are pretty suspicious.

I think this is happening in mesh->communicate(K), K = cos(y - z) so this should be constant in x and known, but I guess the communicate leaves the X guards unspecified? Do we need to apply a BC on X for this test?

Still not sure why this is only materialising now -- I don't really see how this PR should touch this test. Probably need to compile with ASAN and see what that says, but that's for next week

Also make G-value computation more like main branch

ZedThree requested review from bendudson and dschwoerer July 8, 2026 17:08

github-actions Bot reviewed Jul 8, 2026

View reviewed changes

bendudson reviewed Jul 8, 2026

View reviewed changes

Comment thread src/mesh/coordinates_accessor.cxx

bendudson reviewed Jul 8, 2026

View reviewed changes

Comment thread tests/unit/mesh/test_coordinates.cxx Outdated

bendudson reviewed Jul 8, 2026

View reviewed changes

Comment thread tests/unit/mesh/test_coordinates.cxx Outdated

bendudson reviewed Jul 9, 2026

View reviewed changes

Comment thread include/bout/coordinates.hxx

github-actions Bot reviewed Jul 9, 2026

View reviewed changes

dschwoerer previously approved these changes Jul 9, 2026

View reviewed changes

ZedThree dismissed dschwoerer’s stale review via b966661 July 9, 2026 15:15

ZedThree force-pushed the refactor-coords-rebase-rebase branch from 603ae3a to b966661 Compare July 9, 2026 15:15

ZedThree and others added 10 commits July 9, 2026 16:21

Move derivative operators out of Coordinates

c761c01

Clang tidy+format

adec652

Coordinates::readParallelMetricComponents swap co/contra-variant

448fadb

Covariant metrics have `_`. The new code is consistent with e.g. metric reads around line 450.

Mesh::recalculateStaggeredCoordinates move after reset

b7b7f12

`new_coordinates` should be moved into `coords_map` after it has been recalculated, not before.

Coordinates: Invalidate caches when cell spacing set

3580dc5

`CoordinatesAccessor` and non-uniform geometry quantities depend on cell spacing. Reset cache so they are updated when needed. Bxy is not calculated on demand, so is calculated in `set*MetricTensor` methods.

test_coordinates: Fix IndexedAccessors

7414796

z index should be in range. Indexed accessors return BoutReal, not reference to field.

Fix setting location of derivatives

740f40f

Fix broken examples

69950d1

Fix some clang-tidy warnings

7da06af

ZedThree force-pushed the refactor-coords-rebase-rebase branch from b966661 to 7da06af Compare July 9, 2026 15:26

github-actions Bot reviewed Jul 9, 2026

View reviewed changes

Div_par: Use Field3DParallel

d546bec

Convert input to Field3DParallel before dividing by Bxy, to ensure that yup/ydown fields are retained.

github-actions Bot reviewed Jul 9, 2026

View reviewed changes

bendudson added 2 commits July 9, 2026 13:25

Coordinates: Use communicate_no_slices

0870329

Should now be consistent with `next` branch. Avoids calculation of parallel slices for derived metric quantities.

Coordinates: Getters use BOUT_OMP_SAFE

312bfc3

Makes lazy getters thread-safe by adding named critical sections. Named so that nested sections (code in critical section calling a function that also contains a critical section) doesn't lead to deadlock.

Don't compute G-values unnecessarily or when using FCI

e4b126d

Also make G-value computation more like main branch

		@@ -331,7 +331,7 @@ class CWM : public PhysicsModel {
		result = VDDZ(-DDX(p), f);


		SAVE_REPEAT(Ve);

		output.write("dx = {:e}, dy = {:e}, dz = {:e}\n", coords->dx(2, 2), coords->dy(2, 2),

Uh oh!

Conversation

ZedThree commented Jul 8, 2026

Uh oh!

github-actions Bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 8, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 8, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 8, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 8, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 8, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

github-actions Bot left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

dschwoerer left a comment

Choose a reason for hiding this comment

Uh oh!

ZedThree commented Jul 9, 2026

Uh oh!

github-actions Bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 9, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 9, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 9, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 9, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 9, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 9, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 9, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 9, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 9, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 9, 2026