New Y-Boundaries for unified sheath code by dschwoerer · Pull Request #3308 · boutproject/BOUT-dev

dschwoerer · 2026-03-05T08:10:18Z

This allows to implement unified BCs for yup/ydown as well as FCI / FA.
Taking into acount things like whether the boundary is next to the domain, or between the parallel slices.

github-actions

clang-tidy made some suggestions

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

github-actions

clang-tidy made some suggestions

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

github-actions

clang-tidy made some suggestions

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

Makes it easier to reuse for other code

Mimicks the parallel case, to write FCI independent code

The boundary region does not match what is done for the parallel case, thus porting it to a iterator does not work.

hasBndryUpperY / hasBndryLowerY does not work for FCI and thus the request does not make sense / can be configured to throw. Thus it should not be checked if it is not needed.

BoundaryRegionIter has been delted

This allows to extend the boundary code to place the boundary further away from the boundary.

This is more general and takes the offset() into account, and thus works for cases where the boundary is between the first and second guard cell

Previously only the first boundary point was set

Taken from hermes-3, adopted for higher order

This allows to write code for FCI and non-FCI using templates.

Directly iterate over the points

github-actions

clang-tidy made some suggestions

github-actions

clang-tidy made some suggestions

github-actions · 2026-04-15T08:16:23Z

+std::shared_ptr<YBoundary> getYBoundary(Coordinates* coords,
+                                        YBndryType type = YBndryType::sheath) {
+  auto itype = static_cast<int>(type);
+  if (coords->ybndrys[itype] == nullptr) {


warning: do not use array subscript when the index is not an integer constant expression [cppcoreguidelines-pro-bounds-constant-array-index]

if (coords->ybndrys[itype] == nullptr) { ^

github-actions · 2026-04-15T08:16:23Z

+                                        YBndryType type = YBndryType::sheath) {
+  auto itype = static_cast<int>(type);
+  if (coords->ybndrys[itype] == nullptr) {
+    coords->ybndrys[itype] = coords->makeYBoundary(type);


warning: do not use array subscript when the index is not an integer constant expression [cppcoreguidelines-pro-bounds-constant-array-index]

coords->ybndrys[itype] = coords->makeYBoundary(type); ^

github-actions · 2026-04-15T08:16:23Z

+  if (coords->ybndrys[itype] == nullptr) {
+    coords->ybndrys[itype] = coords->makeYBoundary(type);
+  }
+  return coords->ybndrys[itype];


warning: do not use array subscript when the index is not an integer constant expression [cppcoreguidelines-pro-bounds-constant-array-index]

return coords->ybndrys[itype]; ^

Thanks to Mike Kryjak for noticing.

github-actions

clang-tidy made some suggestions

ZedThree

I've not finished going through this, but it looks like it's almost there. The thing that would make this really nice is to template the new boundary region stuff so we can handle X/Y boundaries without any runtime switches or checking bx/by etc.

There's lots of clang-tidy comments that should be fixed as well.

And more docstrings and explanatory comments please!

ZedThree · 2026-04-15T15:42:47Z

+#include <bout/assert.hxx>
+#include <bout/bout_types.hxx>
+
+namespace parallel_stencil {


Suggested change

namespace parallel_stencil {

namespace bout::parallel_stencil {

ZedThree · 2026-04-15T15:58:50Z

+inline BoutReal pow(BoutReal val, int exp) {
+  // constexpr int expval = exp;
+  // static_assert(expval == 2 or expval == 3, "This pow is only for exponent 2 or 3");
+  if (exp == 2) {
+    return val * val;
+  }
+  ASSERT3(exp == 3);
+  return val * val * val;
+}


It doesn't look like we call this with anything other than 2, and from what I can tell, this isn't actually faster than std::pow, so let's just use that?

I tested it with gcc 15.2.1:

g++ example.cxx -std=c++20 -ffast-math -O2 -march=native -isystem include -Lbuild/src -lbenchmark -lpthread -o mybenchmark && ./mybenchmark 2>&1 | tail -n2 BM_StdPow 0.234 ns 0.234 ns 3000470135 BM_multi 0.444 ns 0.444 ns 1573472940 [dave@raven03: benchmark]👍 g++ example.cxx -std=c++20 -O2 -march=native -isystem include -Lbuild/src -lbenchmark -lpthread -o mybenchmark && ./mybenchmark 2>&1 | tail -n2 BM_StdPow 20.5 ns 20.5 ns 34629261 BM_multi 0.441 ns 0.441 ns 1600645212

So std::pow is faster, but only with -ffast-math.

Correction: with -ffast-math they are the same speed, otherwise std::pow is a factor 100 slower.

ZedThree · 2026-04-15T16:06:51Z

-          bndry->setValid(0);
+      for (auto pnt : *bndry) {
+        for (auto pnt2 : *bndry2) {
+#warning this could likely be done faster


I think this is zip? Unfortunately, that's C++23, but we could add an implementation if necessary. Either way, we don't want to merge this with the warning!

ZedThree · 2026-04-15T16:51:42Z

+  BoutReal& ynext(Field3D& f) const { return f[ind().yp(by).xp(bx)]; }
+  const BoutReal& ynext(const Field3D& f) const { return f[ind().yp(by).xp(bx)]; }
+  BoutReal& yprev(Field3D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+  const BoutReal& yprev(const Field3D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+  BoutReal& ythis(Field3D& f) const { return f[ind()]; }


I'm a bit confused by these functions, they start with y, but the implementation has .xp in there, so I guess these are supposed to be generic for X and Y boundaries?

They are only used for y, but could also be used for x.

ZedThree · 2026-04-15T16:55:48Z

+  const BoutReal& yprev(const Field2D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+#endif
+
+  const int dir;


Please sort out all these clang-tidy comments, basically just add getters:

public: int dir() { return dir_; } private: int dir_;

plus docstrings

ZedThree · 2026-04-15T16:57:51Z

  // check that contravariant tensors are positive (if expected) and finite (always)
  void checkContravariant();
+
+  mutable std::array<std::shared_ptr<YBoundary>, 3> ybndrys;


Where does the 3 come from? How do we check that it is still the correct value?

ZedThree · 2026-04-15T17:01:36Z

+    return bndry_position != rhs.bndry_position;
+  }
+
+#define ITER() for (int i = 0; i < localmesh->ystart - abs_offset(); ++i)


We don't really need a macro for this

ZedThree · 2026-04-15T17:05:42Z

+#include <algorithm>
+#include <functional>
+
+class BoundaryRegionIter {


I guess ideally this would be templated on the direction (X or Y)

Does it need to be templated? It could also just be generic without templating?
Are you worried about performance, or what is your motivation for templating?

github-actions bot reviewed Mar 5, 2026

View reviewed changes

github-actions bot reviewed Mar 9, 2026

View reviewed changes

dschwoerer added 27 commits March 10, 2026 08:17

Move stencils to separte header

ae1f285

Makes it easier to reuse for other code

Add basic boundary region iterator

b9b42ff

Mimicks the parallel case, to write FCI independent code

Provide a new boundary iterator based on RangeIterator

9bdfef5

The boundary region does not match what is done for the parallel case, thus porting it to a iterator does not work.

Only check hasBndry*Y if they would be included

c51a85f

hasBndryUpperY / hasBndryLowerY does not work for FCI and thus the request does not make sense / can be configured to throw. Thus it should not be checked if it is not needed.

ensure we dont mix non-fci BCs with fci

06d581b

rename to include order

7f1d26f

Fix: remove broken code

84f8c83

BoundaryRegionIter has been delted

Convert test to new iterator scheeme

ea05bae

Add Field2D version for 2D metrics

17f8efa

Add setYPrevIfValid

0db3052

Add offset to parallel boundary region

d2f2cc8

This allows to extend the boundary code to place the boundary further away from the boundary.

Add function to check wehther point is in the boundary

1f48a59

Add setValid to BoundaryRegionParIterBase

4af5a4d

Reimplement ynext for the boundary region

e5fbcaf

This is more general and takes the offset() into account, and thus works for cases where the boundary is between the first and second guard cell

Fix extrapolaton / interpolation

ef7b03d

Fix parallel boundary to interpolate into the boundary

6dba751

Previously only the first boundary point was set

Ensure data is sorted

f105c3f

Calculate valid for the general case

f7c1aca

fix parallel neumann BC

55bdbcf

Add limitFree

6f40b4e

Taken from hermes-3, adopted for higher order

add setAll to parallel BC

36c1f94

Move yboundary iterator from hermes to BOUT++

c7cbf11

This allows to write code for FCI and non-FCI using templates.

Add iter_pnts function

f9cad29

Directly iterate over the points

Add example on how to replace RangeIterator with YBoundary

e6f004b

Move documentation to sphinx

3a0f214

Fix iter_pnts

4e21aac

Add more documentation on YBoundary

3a850b1

github-actions bot reviewed Apr 15, 2026

View reviewed changes

Comment thread tests/integrated/test-fci-boundary/get_par_bndry.cxx Outdated

dschwoerer added 2 commits April 15, 2026 09:57

More clang-tidy fixes

30298ca

Move getYBoundary to a free friend function

db5d3aa

github-actions bot reviewed Apr 15, 2026

View reviewed changes

dschwoerer added 10 commits April 15, 2026 10:27

Declare function inline

c67ae87

Avoid warning for unused push

9558cca

move reference value

8d7fbd1

Fix documentation

6729983

Thanks to Mike Kryjak for noticing.

Set enum size for performance

38e3edb

Also declare size in forward definition

6051461

Change BOUT_ENUM_CLASS to 8bit integer

d170774

Ensure int8_t is defined

20d61c4

Ensure int8_t is defined

399d219

int8_t is defined in cstdint not cstddef

eaac193

github-actions bot reviewed Apr 15, 2026

View reviewed changes

dschwoerer added 2 commits April 15, 2026 14:05

Header is used by macro, add comment for IWYU

372fbe7

Do not move trivially-copyable variable

2dfa950

ZedThree reviewed Apr 15, 2026

View reviewed changes

dschwoerer added 11 commits April 16, 2026 12:35

Do not use custom pow implementation.

1ff62a4

Use bout namespace

afe1909

Use func for functions, and f for field

6d18997

Use point rather then pnt

ad0f14e

fixup std pow

f0e6f28

Only check for FCI if CHECK > 0

92569dd

getBoundariesPar can be const

42836ef

Make dir a function and store a reference to mesh, rather then a pointer

ceb557a

use std::copysign

72cd8a4

Update test

b38fdf3

Make members private, provide public accessors where needed

898d191

dschwoerer marked this pull request as draft April 16, 2026 13:28

	namespace parallel_stencil {
	namespace bout::parallel_stencil {

Conversation

dschwoerer commented Mar 5, 2026

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!

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!

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!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

github-actions bot Apr 15, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Apr 15, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Apr 15, 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!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

ZedThree left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!