Skip to content

Refactor Coordinates with more encapsulation - #3425

Merged
ZedThree merged 20 commits into
nextfrom
refactor-coords-rebase-rebase
Aug 13, 2026
Merged

Refactor Coordinates with more encapsulation#3425
ZedThree merged 20 commits into
nextfrom
refactor-coords-rebase-rebase

Conversation

@ZedThree

@ZedThree ZedThree commented Jul 8, 2026

Copy link
Copy Markdown
Member

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)

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

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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.

Comment thread examples/conducting-wall-mode/cwm.cxx
Comment thread examples/conducting-wall-mode/cwm.cxx Outdated
Comment thread examples/conducting-wall-mode/cwm.cxx
Comment thread examples/conducting-wall-mode/cwm.cxx
Comment thread examples/elm-pb-outerloop/elm_pb_outerloop.cxx
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/mesh/coordinates_accessor.cxx
Comment thread tests/unit/mesh/test_coordinates.cxx Outdated
Comment thread tests/unit/mesh/test_coordinates.cxx Outdated
Comment thread include/bout/coordinates.hxx

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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.

Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector2d.cxx Outdated
Comment thread src/field/vector3d.cxx Outdated
dschwoerer
dschwoerer previously approved these changes Jul 9, 2026

@dschwoerer dschwoerer left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

LGTM

ZedThree and others added 10 commits July 9, 2026 16:21
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
ZedThree force-pushed the refactor-coords-rebase-rebase branch from b966661 to 7da06af Compare July 9, 2026 15:26
@ZedThree

ZedThree commented Jul 9, 2026

Copy link
Copy Markdown
Member Author

Thanks for the fixes @bendudson!

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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.

Comment thread src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.cxx
Comment thread src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.cxx
Comment thread src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.cxx
Comment thread src/invert/pardiv/impls/cyclic/pardiv_cyclic.cxx
#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) \

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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

#define COPY_STRIPE1(symbol)                                         \
        ^

Comment thread tests/MMS/GBS/gbs.cxx
dx4 = SQ(SQ(coords->dx));
dy4 = SQ(SQ(coords->dy));
dz4 = SQ(SQ(coords->dz));
dx4 = SQ(SQ(coords->dx()));

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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

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

Comment thread tests/MMS/GBS/gbs.cxx
dz4 = SQ(SQ(coords->dz));
dx4 = SQ(SQ(coords->dx()));
dy4 = SQ(SQ(coords->dy()));
dz4 = SQ(SQ(coords->dz()));

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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()));
        ^

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Maybe we should delete it?
It does not compile, but even with -DBOUT_RUN_ALL_TESTS=ON I do not get it to compile.
So this has not been working for a while.

Somewhat related, most MMS tests that can be enabled with -DBOUT_RUN_ALL_TESTS=ON are failing at runtime.

Comment thread tests/MMS/GBS/gbs.cxx

SAVE_REPEAT(Ve);

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

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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);
                          ^

Comment thread tests/MMS/GBS/gbs.cxx
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());

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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());
                                           ^

const auto& I = tokamak_coords.I_unnormalised;

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

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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 github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

clang-tidy made some suggestions

}

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

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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));
                    ^

// 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);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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

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


// 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())

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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())
        ^


// 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())

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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())
                       ^


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

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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>


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

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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>


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

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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>

bendudson added 2 commits July 9, 2026 13:25
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

Copy link
Copy Markdown
Contributor

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

@ZedThree

Copy link
Copy Markdown
Member Author

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

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

clang-tidy made some suggestions

Comment thread tests/MMS/GBS/gbs.cxx
dz4 = SQ(SQ(coords->dz));
dx4 = SQ(SQ(coords->dx()));
dy4 = SQ(SQ(coords->dy()));
dz4 = SQ(SQ(coords->dz()));

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

warning: assigning to 'BoutReal' (aka 'double') from incompatible type 'BinaryExpr<UnaryResT, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Square>, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Square>, bout::op::Square>' (aka '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()));
        ^

const auto& I = tokamak_coords.I_unnormalised;

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

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

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:361: 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:362: 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:238: candidate function not viable: requires 2 arguments, but 0 were provided

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

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

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

include/bout/field2d.hxx:271: 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:274: candidate function not viable: requires 3 arguments, but 0 were provided

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

@github-actions github-actions Bot left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

clang-tidy made some suggestions

Comment thread include/bout/derivs.hxx
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);
return DDY(Field3D{Field3DParallel{expr}}, outloc, method, region);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

warning: slicing object from type 'Field3DParallel' to 'Field3D' discards 1 bytes of state [cppcoreguidelines-slicing]

              ^

dschwoerer
dschwoerer previously approved these changes Aug 12, 2026
dschwoerer
dschwoerer previously approved these changes Aug 13, 2026
@ZedThree
ZedThree merged commit 5ea320a into next Aug 13, 2026
23 checks passed
@ZedThree
ZedThree deleted the refactor-coords-rebase-rebase branch August 13, 2026 09:39
@ZedThree

Copy link
Copy Markdown
Member Author

Thanks @tomc271 and @dschwoerer! A big change a long time coming!

@dschwoerer

Copy link
Copy Markdown
Contributor

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

- coords.dx
+ coords.dx()

Is there a tool that can do that somewhat automatically?

@ZedThree

Copy link
Copy Markdown
Member Author

bout-upgrader v6 v6_upgrader <files>...

@ZedThree ZedThree mentioned this pull request Aug 13, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants