Improvements for hermitespline interpolation by dschwoerer · Pull Request #3298 · boutproject/BOUT-dev

dschwoerer · 2026-03-03T13:02:07Z

Includes X-splitting also for monotonichermitespline and more.

For FCI we need to be able to access "random" data from the adjacent slices. If they are split in x-direction, this requires some tricky communication pattern. It can be used like this: ``` // Create object GlobalField3DAccess fci_comm(thismesh); // let it know what data points will be required: // where IndG3D is an index in the global field, which would be the // normal Ind3D if there would be only one proc. fci_comm.get(IndG3D(i, ny, nz)); // If all index have been added, the communication pattern will be // established. This has to be called by all processors in parallel fci_comm.setup() // Once the data for a given field is needed, it needs to be // communicated: GlobalField3DAccessInstance global_data = fci_comm.communicate(f3d); // and can be accessed like this BoutReal data = global_data[IndG3D(i, ny, nz)]; // ny and nz in the IndG3D are always optional. ```

If they are two instances of the same template, this allows to have an if in the inner loop that can be optimised out.

lower_bound takes into account the data is sorted

Ensures we all ways check for monotonicity

Tags were different for sender and receiver

Otherwise mpi might wait for the wrong request.

to avoid confusion whether the offsets are for sending or receiving

Using a local set for each thread ensures we do not need a mutex for adding data, at the cost of having to merge the different sets later.

We want to skip sending if there is no data for this process ...

The result needs to be well understood, as the indices are a mix of local and global indices, as we need to access non-local data.

dschwoerer · 2026-03-10T11:16:54Z

Looks ok, but needs a lot more explanatory comments and docstrings. I can guess what a lot of things do, but I'd rather read a high-level overview first.

I have added some docs.

Lots of things need renaming too -- again, I can guess, but I'd really rather not. We have to think about new developers, are they going to be able to work out what xyzl or o_ids means?

I improved the naming

The template base class for the monotonic version is a good idea, but I still really don't like squashing together the PETSc, "new weights", and non-PETSc versions. It means three separate builds to test them all, rather than being able to test all of them in a single build. The interleaved preprocessor statements also make it much harder to understand what's going on.

I extended on the templating. Now all versions are templated.

Lastly, it would be good to make this already work for parallelised-Z so we have to do less work when that's ready. There's also already two other mechanisms for local-to-global mapping, GlobalIndexer and GlobalField (which I'm not sure is even used for anything seriously now), it would be good to understand if they could be used here at all.

I added an interface to the mesh, to do the mapping. If you think this should be reproduced here, I would need to know how you intend to do it.

All comments addressed, ready for re-review 👍

github-actions

clang-tidy made some suggestions

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

github-actions · 2026-03-10T11:20:05Z

 #endif

+namespace {
+enum class implementation_type { new_weights, petsc, legacy };


warning: enum 'implementation_type' uses a larger base type ('int', size: 4 bytes) than necessary for its value set, consider using 'std::uint8_t' (1 byte) as the base type to reduce its size [performance-enum-size]

enum class implementation_type { new_weights, petsc, legacy }; ^

github-actions · 2026-03-10T11:20:05Z

  PE_XIND = MYPE % NXPE;
  PE_ZIND = 0;

+  ASSERT2(MYPE == getProcIndex(PE_XIND, PE_YIND, PE_ZIND));


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

ASSERT2(MYPE == getProcIndex(PE_XIND, PE_YIND, PE_ZIND)); ^

github-actions · 2026-03-10T11:20:05Z

  }
 }

+int BoutMesh::getProcIndex(int X, int Y, int Z) { return Y * NXPE + X; }


warning: '*' has higher precedence than '+'; add parentheses to explicitly specify the order of operations [readability-math-missing-parentheses]

Suggested change

int BoutMesh::getProcIndex(int X, int Y, int Z) { return Y * NXPE + X; }

int BoutMesh::getProcIndex(int X, int Y, int Z) { return (Y * NXPE) + X; }

github-actions · 2026-03-10T11:20:05Z

  }
 }

+int BoutMesh::getProcIndex(int X, int Y, int Z) { return Y * NXPE + X; }


warning: method 'getProcIndex' can be made const [readability-make-member-function-const]

Suggested change

int BoutMesh::getProcIndex(int X, int Y, int Z) { return Y * NXPE + X; }

int BoutMesh::getProcIndex(int X, int Y, int Z) const { return Y * NXPE + X; }

github-actions · 2026-03-10T11:20:05Z

  }
 }

+int BoutMesh::getProcIndex(int X, int Y, int Z) { return Y * NXPE + X; }


warning: out-of-line definition of 'getProcIndex' does not match any declaration in 'BoutMesh' [clang-diagnostic-error]

int BoutMesh::getProcIndex(int X, int Y, int Z) { return Y * NXPE + X; } ^

Additional context

src/mesh/impls/bout/boutmesh.hxx:18: BoutMesh defined here

class BoutMesh : public Mesh { ^

src/mesh/impls/bout/boutmesh.hxx:65: member declaration does not match because it is const qualified

int getProcIndex(int X, int Y, int Z) const override; ^

github-actions · 2026-03-10T11:20:06Z

-#endif
+    if constexpr (monotonic) {
+      const auto gind =
+          gf3daccess->xyzg(i_corn, y + y_offset + y_global_offset, k_corner(x, y, z));


warning: no member named 'xyzg' in 'GlobalField3DAccess' [clang-diagnostic-error]

gf3daccess->xyzg(i_corn, y + y_offset + y_global_offset, k_corner(x, y, z)); ^

github-actions · 2026-03-10T11:20:06Z

+                                                           const BoutMask& mask,
+                                                           const std::string& region) {
  setMask(mask);
  calcWeights(delta_x, delta_z, region);


warning: no viable conversion from 'const std::string' (aka 'const basic_string') to 'const BoutMask' [clang-diagnostic-error]

calcWeights(delta_x, delta_z, region); ^

Additional context

src/mesh/interpolation/hermite_spline_xz.cxx:535: in instantiation of member function 'XZHermiteSplineBase<true, (anonymous namespace)::implementation_type::new_weights>::calcWeights' requested here

template class XZHermiteSplineBase<true, implementation_type::new_weights>; ^

include/bout/mask.hxx:47: candidate constructor (the implicit copy constructor) not viable: no known conversion from 'const std::string' (aka 'const basic_string') to 'const BoutMask &' for 1st argument

class BoutMask { ^

include/bout/mask.hxx:47: candidate constructor (the implicit move constructor) not viable: no known conversion from 'const std::string' (aka 'const basic_string') to 'BoutMask &&' for 1st argument

class BoutMask { ^

include/bout/mask.hxx:55: explicit constructor is not a candidate

explicit BoutMask(const Mesh& mesh, bool value = false) ^

include/bout/mask.hxx:57: explicit constructor is not a candidate

explicit BoutMask(const Mesh* mesh = nullptr, bool value = false) ^

/usr/include/c++/15/bits/basic_string.h:1035: candidate function

operator __sv_type() const noexcept ^

include/bout/interpolation_xz.hxx:209: passing argument to parameter 'mask' here

void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask, ^

github-actions · 2026-03-10T11:20:06Z

-      f_interp[iyp] += newWeights[w * 4 + 3][i] * f[ic.zp(2).xp(w - 1)];
+
+  if constexpr (imp_type == implementation_type::petsc) {
+    BoutReal* ptr;


warning: variable 'ptr' is not initialized [cppcoreguidelines-init-variables]

Suggested change

BoutReal* ptr;

BoutReal* ptr = nullptr;

github-actions · 2026-03-10T11:20:06Z

+
+  if constexpr (imp_type == implementation_type::petsc) {
+    BoutReal* ptr;
+    const BoutReal* cptr;


warning: variable 'cptr' is not initialized [cppcoreguidelines-init-variables]

Suggested change

const BoutReal* cptr;

const BoutReal* cptr = nullptr;

github-actions · 2026-03-10T11:20:06Z

+      if constexpr (monotonic) {
+        const auto iyp = i;
+        const auto i = iyp.ym(y_offset);
+        const auto corners = {(*gf)[IndG3D(g3dinds[i][0])], (*gf)[IndG3D(g3dinds[i][1])],


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

const auto corners = {(*gf)[IndG3D(g3dinds[i][0])], (*gf)[IndG3D(g3dinds[i][1])], ^

github-actions

clang-tidy made some suggestions

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

github-actions · 2026-03-10T11:59:32Z

  }
 }

+int BoutMesh::getProcIndex(int X, int Y, int Z) const { return Y * NXPE + X; }


warning: '*' has higher precedence than '+'; add parentheses to explicitly specify the order of operations [readability-math-missing-parentheses]

Suggested change

int BoutMesh::getProcIndex(int X, int Y, int Z) const { return Y * NXPE + X; }

int BoutMesh::getProcIndex(int X, int Y, int Z) const { return (Y * NXPE) + X; }

github-actions · 2026-03-10T11:59:32Z

  }
 }

+int BoutMesh::getProcIndex(int X, int Y, int Z) const { return Y * NXPE + X; }


warning: parameter 'Z' is unused [misc-unused-parameters]

Suggested change

int BoutMesh::getProcIndex(int X, int Y, int Z) const { return Y * NXPE + X; }

int BoutMesh::getProcIndex(int X, int Y, int /*Z*/) const { return Y * NXPE + X; }

github-actions · 2026-03-10T11:59:32Z

  }
 }

+int BoutMesh::getProcIndex(int X, int Y, int Z) const { return Y * NXPE + X; }


warning: unused parameter 'Z' [clang-diagnostic-unused-parameter]

int BoutMesh::getProcIndex(int X, int Y, int Z) const { return Y * NXPE + X; } ^

github-actions · 2026-03-10T11:59:32Z

-  if (localmesh->getNXPE() > 1) {
-    throw BoutException("Require PETSc for MPI splitting in X");
+  if constexpr (imp_type == implementation_type::petsc) {
+    petsclib = new PetscLib(


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

src/mesh/interpolation/hermite_spline_xz.cxx:31:

+ #include "bout/petsclib.hxx"

github-actions · 2026-03-10T11:59:32Z

+                                                           const BoutMask& mask,
+                                                           const std::string& region) {
  setMask(mask);
  calcWeights(delta_x, delta_z, region);


warning: no viable conversion from 'const std::string' (aka 'const basic_string') to 'const BoutMask' [clang-diagnostic-error]

calcWeights(delta_x, delta_z, region); ^

Additional context

src/mesh/interpolation/hermite_spline_xz.cxx:537: in instantiation of member function 'XZHermiteSplineBase<true, (anonymous namespace)::implementation_type::new_weights>::calcWeights' requested here

template class XZHermiteSplineBase<true, implementation_type::new_weights>; ^

include/bout/mask.hxx:47: candidate constructor (the implicit copy constructor) not viable: no known conversion from 'const std::string' (aka 'const basic_string') to 'const BoutMask &' for 1st argument

class BoutMask { ^

include/bout/mask.hxx:47: candidate constructor (the implicit move constructor) not viable: no known conversion from 'const std::string' (aka 'const basic_string') to 'BoutMask &&' for 1st argument

class BoutMask { ^

include/bout/mask.hxx:55: explicit constructor is not a candidate

explicit BoutMask(const Mesh& mesh, bool value = false) ^

include/bout/mask.hxx:57: explicit constructor is not a candidate

explicit BoutMask(const Mesh* mesh = nullptr, bool value = false) ^

/usr/include/c++/15/bits/basic_string.h:1035: candidate function

operator __sv_type() const noexcept ^

include/bout/interpolation_xz.hxx:209: passing argument to parameter 'mask' here

void calcWeights(const Field3D& delta_x, const Field3D& delta_z, const BoutMask& mask, ^

github-actions · 2026-03-10T11:59:33Z

+    const auto piy = global2local_y.convert(gind.y());
+    const auto piz = global2local_z.convert(gind.z());
+    ASSERT3(piz.proc == 0);
+    toGet[mesh.getProcIndex(pix.proc, piy.index, piz.index)].push_back(


warning: no member named 'index' in 'fci_comm::ProcLocal' [clang-diagnostic-error]

toGet[mesh.getProcIndex(pix.proc, piy.index, piz.index)].push_back( ^

github-actions · 2026-03-10T11:59:33Z

+    int offset = 0;
+    for (const auto& get : toGet) {
+      getOffsets.push_back(offset);
+      offset += get.size();


warning: narrowing conversion from 'size_type' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]

offset += get.size(); ^

github-actions · 2026-03-10T11:59:33Z

+    getOffsets.push_back(offset);
+  }
+  for (const auto id : ids) {
+    const IndG3D gind{id, global2local_y.getGlobalWith(), global2local_z.getGlobalWith()};


warning: use of undeclared identifier 'global2local_y' [clang-diagnostic-error]

const IndG3D gind{id, global2local_y.getGlobalWith(), global2local_z.getGlobalWith()}; ^

github-actions · 2026-03-10T11:59:33Z

+  for (const auto id : ids) {
+    const IndG3D gind{id, global2local_y.getGlobalWith(), global2local_z.getGlobalWith()};
+    const auto pix = global2local_x.convert(gind.x());
+    const auto piy = global2local_y.convert(gind.y());


warning: use of undeclared identifier 'global2local_y' [clang-diagnostic-error]

const auto piy = global2local_y.convert(gind.y()); ^

github-actions · 2026-03-10T11:59:33Z

+    const auto piy = global2local_y.convert(gind.y());
+    const auto piz = global2local_z.convert(gind.z());
+    ASSERT3(piz.proc == 0);
+    const auto proc = mesh.getProcIndex(pix.proc, piy.index, piz.index);


warning: member reference type 'Mesh *' is a pointer; did you mean to use '->'? [clang-diagnostic-error]

Suggested change

const auto proc = mesh.getProcIndex(pix.proc, piy.index, piz.index);

const auto proc = mesh->getProcIndex(pix.proc, piy.index, piz.index);

github-actions

clang-tidy made some suggestions

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

github-actions · 2026-03-10T12:17:49Z

+      continue;
+    }
+    auto ret = MPI_Irecv(data.data() + getOffsets[proc], toGet[proc].size(), MPI_DOUBLE,
+                         proc, 666, comm, reqs.data() + cnt1);


warning: narrowing conversion from 'size_t' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]

proc, 666, comm, reqs.data() + cnt1); ^

github-actions · 2026-03-10T12:17:49Z

+    for (auto i : toSend[proc]) {
+      sendBuffer[cnt++] = f[Ind3D(i)];
+    }
+    auto ret = MPI_Send(start, toSend[proc].size(), MPI_DOUBLE, proc, 666, comm);


warning: narrowing conversion from 'size_t' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]

auto ret = MPI_Send(start, toSend[proc].size(), MPI_DOUBLE, proc, 666, comm); ^

github-actions · 2026-03-10T12:17:49Z

+    for (auto i : toSend[proc]) {
+      sendBuffer[cnt++] = f[Ind3D(i)];
+    }
+    auto ret = MPI_Send(start, toSend[proc].size(), MPI_DOUBLE, proc, 666, comm);


warning: narrowing conversion from 'size_type' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]

auto ret = MPI_Send(start, toSend[proc].size(), MPI_DOUBLE, proc, 666, comm); ^

github-actions

clang-tidy made some suggestions

github-actions · 2026-03-10T12:41:57Z

+    auto it = std::lower_bound(vec.begin(), vec.end(), tofind);
+    ASSERT3(it != vec.end());
+    ASSERT3(*it == tofind);
+    mapping[id] = std::distance(vec.begin(), it) + getOffsets[proc];


warning: narrowing conversion from 'typename iterator_traits<__gnu_cxx::__normal_iterator<const int *, std::vector<int, std::allocator>>>::difference_type' (aka 'long') to signed type 'mapped_type' (aka 'int') is implementation-defined [bugprone-narrowing-conversions]

mapping[id] = std::distance(vec.begin(), it) + getOffsets[proc]; ^

github-actions · 2026-03-10T12:41:57Z

+    auto it = std::lower_bound(vec.begin(), vec.end(), tofind);
+    ASSERT3(it != vec.end());
+    ASSERT3(*it == tofind);
+    mapping[id] = std::distance(vec.begin(), it) + getOffsets[proc];


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

src/mesh/parallel/fci_comm.cxx:33:

- #include <mpi.h> + #include <iterator> + #include <mpi.h>

github-actions · 2026-03-10T12:41:57Z

+  fci_comm::GlobalToLocal1D global2local_z;
+
+public:
+  fci_comm::XYZ2Ind<Ind3D> xyzlocal;


warning: member variable 'xyzlocal' has public visibility [cppcoreguidelines-non-private-member-variables-in-classes]

fci_comm::XYZ2Ind<Ind3D> xyzlocal; ^

github-actions · 2026-03-10T12:41:57Z

+
+public:
+  fci_comm::XYZ2Ind<Ind3D> xyzlocal;
+  fci_comm::XYZ2Ind<IndG3D> xyzglobal;


warning: member variable 'xyzglobal' has public visibility [cppcoreguidelines-non-private-member-variables-in-classes]

fci_comm::XYZ2Ind<IndG3D> xyzglobal; ^

Also, the indices should all start at 0

dschwoerer added 30 commits March 3, 2026 12:45

Allow XZHermiteSpline also without y-offset

5502ce3

Unify XZMonotonicHermiteSpline and XZMonotonicHermiteSpline

afc68d9

If they are two instances of the same template, this allows to have an if in the inner loop that can be optimised out.

Use x-splitting for monotonichermitespline test

84237c8

Set region for lagrange4pt

0662449

Add monotonic check also to other code branches

f77849e

use lower_bound instead of find

820f0a5

lower_bound takes into account the data is sorted

Do not shadow mapping

1178243

Ensure setup has been called

3d6f587

Use pointer to data

265c8cd

Call setup before communicator is used

16a162e

Deduplicate code

e5878c1

Ensures we all ways check for monotonicity

Fix tags for comm

aee2713

Tags were different for sender and receiver

Use pointer instead of std::vector

b4f7be5

Do not reuse requests if the array is still in use

578ee8d

Otherwise mpi might wait for the wrong request.

rename offset to getOffsets

bf9c9eb

to avoid confusion whether the offsets are for sending or receiving

Fix: mixup of sending / receiving size

2155f06

Fix receive data offset

0662247

Add check to ensure the proc layout is as expected

28d4e28

clang-format

e4507d9

Use BOUT++ assert

d7ba297

Fix check

f55f8e4

Take periodicity into account

0e9f7d6

Sort a reference, not a copy

944adfb

Only communicate non-empty vectors

6daec20

Make check stricter

3444dd0

Fix: include global offset in monotonic spline

91bee68

make fci_comm openmp thread safe

d7aeae4

Using a local set for each thread ensures we do not need a mutex for adding data, at the cost of having to merge the different sets later.

Fix communication for fci

ac5cc70

We want to skip sending if there is no data for this process ...

Add check to reduce risk of bugs

885b465

The result needs to be well understood, as the indices are a mix of local and global indices, as we need to access non-local data.

dschwoerer added 10 commits March 10, 2026 09:59

Test serveral implementations

cb84b76

Cleanup interface and move more functions to cxx file

a81d3eb

Move processor ID calculation to mesh

4632df0

Use getNZPE()

4b0ed1c

respect mesh->periodicX

325dde3

Ensure we are periodic in Y

8dc7b4f

avoid static_cast<void*> for pointer

db53852

Merge remote-tracking branch 'origin/next' into hermitespline

29f9727

Add more docs

dcafb6c

Use better names

424136b

dschwoerer force-pushed the hermitespline branch from 06f3890 to 424136b Compare March 10, 2026 11:15

github-actions Bot reviewed Mar 10, 2026

View reviewed changes

dschwoerer added 2 commits March 10, 2026 12:43

Fixup getProcIndex

500989b

Add headers from clang-tidy

cabf678

dschwoerer force-pushed the hermitespline branch from 88827de to cabf678 Compare March 10, 2026 11:45

[bot] Apply format changes

01e2263

github-actions Bot reviewed Mar 10, 2026

View reviewed changes

Apply clang-tidy fixes

6bd15d2

dschwoerer force-pushed the hermitespline branch from fd92807 to 6bd15d2 Compare March 10, 2026 12:12

github-actions Bot reviewed Mar 10, 2026

View reviewed changes

Fixup compile issues

f2d860f

github-actions Bot reviewed Mar 10, 2026

View reviewed changes

dschwoerer added 4 commits March 11, 2026 09:49

Add getProcIndex to fake_mesh

159ce2f

Also, the indices should all start at 0

Add guards for petsc related code

1de6329

Fixup: always use proc for proc computation

47f93b1

Fixup rename

58db33a

bendudson approved these changes Mar 11, 2026

View reviewed changes

bendudson merged commit 568c70a into next Mar 11, 2026
32 of 33 checks passed

bendudson deleted the hermitespline branch March 11, 2026 16:31

	int BoutMesh::getProcIndex(int X, int Y, int Z) { return Y * NXPE + X; }
	int BoutMesh::getProcIndex(int X, int Y, int Z) { return (Y * NXPE) + X; }

	const auto proc = mesh.getProcIndex(pix.proc, piy.index, piz.index);
	const auto proc = mesh->getProcIndex(pix.proc, piy.index, piz.index);

Conversation

dschwoerer commented Mar 3, 2026

Uh oh!

dschwoerer commented Mar 10, 2026

Uh oh!

github-actions Bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 10, 2026

Choose a reason for hiding this comment

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!

github-actions Bot Mar 10, 2026