From ae1f285b10d073d9109e8b8b6ede9f6aba0f2418 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 27 Mar 2024 13:29:58 +0100 Subject: [PATCH 01/62] Move stencils to separte header Makes it easier to reuse for other code --- include/bout/parallel_boundary_region.hxx | 39 +---------------------- include/bout/sys/parallel_stencils.hxx | 39 +++++++++++++++++++++++ 2 files changed, 40 insertions(+), 38 deletions(-) create mode 100644 include/bout/sys/parallel_stencils.hxx diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 4d5278d00f..44e5cc562e 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -5,6 +5,7 @@ #include "bout/bout_types.hxx" #include +#include "bout/sys/parallel_stencils.hxx" #include #include @@ -14,44 +15,6 @@ * */ -namespace parallel_stencil { -// generated by src/mesh/parallel_boundary_stencil.cxx.py -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; -} -inline BoutReal dirichlet_o1(BoutReal UNUSED(spacing0), BoutReal value0) { - return value0; -} -inline BoutReal dirichlet_o2(BoutReal spacing0, BoutReal value0, BoutReal spacing1, - BoutReal value1) { - return (spacing0 * value1 - spacing1 * value0) / (spacing0 - spacing1); -} -inline BoutReal neumann_o2(BoutReal UNUSED(spacing0), BoutReal value0, BoutReal spacing1, - BoutReal value1) { - return -spacing1 * value0 + value1; -} -inline BoutReal dirichlet_o3(BoutReal spacing0, BoutReal value0, BoutReal spacing1, - BoutReal value1, BoutReal spacing2, BoutReal value2) { - return (pow(spacing0, 2) * spacing1 * value2 - pow(spacing0, 2) * spacing2 * value1 - - spacing0 * pow(spacing1, 2) * value2 + spacing0 * pow(spacing2, 2) * value1 - + pow(spacing1, 2) * spacing2 * value0 - spacing1 * pow(spacing2, 2) * value0) - / ((spacing0 - spacing1) * (spacing0 - spacing2) * (spacing1 - spacing2)); -} -inline BoutReal neumann_o3(BoutReal spacing0, BoutReal value0, BoutReal spacing1, - BoutReal value1, BoutReal spacing2, BoutReal value2) { - return (2 * spacing0 * spacing1 * value2 - 2 * spacing0 * spacing2 * value1 - + pow(spacing1, 2) * spacing2 * value0 - pow(spacing1, 2) * value2 - - spacing1 * pow(spacing2, 2) * value0 + pow(spacing2, 2) * value1) - / ((spacing1 - spacing2) * (2 * spacing0 - spacing1 - spacing2)); -} -} // namespace parallel_stencil - class BoundaryRegionPar : public BoundaryRegionBase { struct RealPoint { diff --git a/include/bout/sys/parallel_stencils.hxx b/include/bout/sys/parallel_stencils.hxx new file mode 100644 index 0000000000..34a51c5285 --- /dev/null +++ b/include/bout/sys/parallel_stencils.hxx @@ -0,0 +1,39 @@ +#pragma once + +namespace parallel_stencil { +// generated by src/mesh/parallel_boundary_stencil.cxx.py +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; +} +inline BoutReal dirichlet_o1(BoutReal UNUSED(spacing0), BoutReal value0) { + return value0; +} +inline BoutReal dirichlet_o2(BoutReal spacing0, BoutReal value0, BoutReal spacing1, + BoutReal value1) { + return (spacing0 * value1 - spacing1 * value0) / (spacing0 - spacing1); +} +inline BoutReal neumann_o2(BoutReal UNUSED(spacing0), BoutReal value0, BoutReal spacing1, + BoutReal value1) { + return -spacing1 * value0 + value1; +} +inline BoutReal dirichlet_o3(BoutReal spacing0, BoutReal value0, BoutReal spacing1, + BoutReal value1, BoutReal spacing2, BoutReal value2) { + return (pow(spacing0, 2) * spacing1 * value2 - pow(spacing0, 2) * spacing2 * value1 + - spacing0 * pow(spacing1, 2) * value2 + spacing0 * pow(spacing2, 2) * value1 + + pow(spacing1, 2) * spacing2 * value0 - spacing1 * pow(spacing2, 2) * value0) + / ((spacing0 - spacing1) * (spacing0 - spacing2) * (spacing1 - spacing2)); +} +inline BoutReal neumann_o3(BoutReal spacing0, BoutReal value0, BoutReal spacing1, + BoutReal value1, BoutReal spacing2, BoutReal value2) { + return (2 * spacing0 * spacing1 * value2 - 2 * spacing0 * spacing2 * value1 + + pow(spacing1, 2) * spacing2 * value0 - pow(spacing1, 2) * value2 + - spacing1 * pow(spacing2, 2) * value0 + pow(spacing2, 2) * value1) + / ((spacing1 - spacing2) * (2 * spacing0 - spacing1 - spacing2)); +} +} // namespace parallel_stencil From b9b42fff9cc68a1524fe801b70109df6a009be80 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 27 Mar 2024 13:40:07 +0100 Subject: [PATCH 02/62] Add basic boundary region iterator Mimicks the parallel case, to write FCI independent code --- include/bout/boundary_region.hxx | 93 ++++++++++++++++++++++++++++++++ src/mesh/boundary_region.cxx | 3 ++ 2 files changed, 96 insertions(+) diff --git a/include/bout/boundary_region.hxx b/include/bout/boundary_region.hxx index 58de12045e..6e7a939d3c 100644 --- a/include/bout/boundary_region.hxx +++ b/include/bout/boundary_region.hxx @@ -4,6 +4,9 @@ class BoundaryRegion; #ifndef BOUT_BNDRY_REGION_H #define BOUT_BNDRY_REGION_H +#include "bout/mesh.hxx" +#include "bout/region.hxx" +#include "bout/sys/parallel_stencils.hxx" #include #include @@ -62,6 +65,7 @@ public: isDone() = 0; ///< Returns true if outside domain. Can use this with nested nextX, nextY }; +class BoundaryRegionIter; /// Describes a region of the boundary, and a means of iterating over it class BoundaryRegion : public BoundaryRegionBase { public: @@ -80,6 +84,95 @@ public: virtual void next1d() = 0; ///< Loop over the innermost elements virtual void nextX() = 0; ///< Just loop over X virtual void nextY() = 0; ///< Just loop over Y + + BoundaryRegionIter begin(); + BoundaryRegionIter end(); +}; + +class BoundaryRegionIter { +public: + BoundaryRegionIter(BoundaryRegion* rgn, bool is_end) + : rgn(rgn), is_end(is_end), dir(rgn->bx + rgn->by) { + //static_assert(std::is_base_of, "BoundaryRegionIter only works on BoundaryRegion"); + + // Ensure only one is non-zero + ASSERT3(rgn->bx * rgn->by == 0); + if (!is_end) { + rgn->first(); + } + } + bool operator!=(const BoundaryRegionIter& rhs) { + if (is_end) { + if (rhs.is_end || rhs.rgn->isDone()) { + return false; + } else { + return true; + } + } + if (rhs.is_end) { + return !rgn->isDone(); + } + return ind() != rhs.ind(); + } + + Ind3D ind() const { return xyz2ind(rgn->x - rgn->bx, rgn->y - rgn->by, z); } + BoundaryRegionIter& operator++() { + ASSERT3(z < nz()); + z++; + if (z == nz()) { + z = 0; + rgn->next(); + } + return *this; + } + BoundaryRegionIter& operator*() { return *this; } + + void dirichlet_o2(Field3D& f, BoutReal value) const { + ynext(f) = parallel_stencil::dirichlet_o2(1, f[ind()], 0.5, value); + } + + BoutReal extrapolate_grad_o2(const Field3D& f) const { return f[ind()] - yprev(f); } + + BoutReal extrapolate_sheath_o2(const Field3D& f) const { + return (f[ind()] * 3 - yprev(f)) * 0.5; + } + + BoutReal extrapolate_next_o2(const Field3D& f) const { return 2 * f[ind()] - yprev(f); } + + BoutReal + extrapolate_next_o2(const std::function& f) const { + return 2 * f(0, ind()) - f(0, ind().yp(-rgn->by).xp(-rgn->bx)); + } + + BoutReal interpolate_sheath(const Field3D& f) const { + return (f[ind()] + ynext(f)) * 0.5; + } + + BoutReal& ynext(Field3D& f) const { return f[ind().yp(rgn->by).xp(rgn->bx)]; } + const BoutReal& ynext(const Field3D& f) const { + return f[ind().yp(rgn->by).xp(rgn->bx)]; + } + BoutReal& yprev(Field3D& f) const { return f[ind().yp(-rgn->by).xp(-rgn->bx)]; } + const BoutReal& yprev(const Field3D& f) const { + return f[ind().yp(-rgn->by).xp(-rgn->bx)]; + } + +private: + BoundaryRegion* rgn; + const bool is_end; + int z{0}; + +public: + const int dir; + +private: + int nx() const { return rgn->localmesh->LocalNx; } + int ny() const { return rgn->localmesh->LocalNy; } + int nz() const { return rgn->localmesh->LocalNz; } + + Ind3D xyz2ind(int x, int y, int z) const { + return Ind3D{(x * ny() + y) * nz() + z, ny(), nz()}; + } }; class BoundaryRegionXIn : public BoundaryRegion { diff --git a/src/mesh/boundary_region.cxx b/src/mesh/boundary_region.cxx index 700ef8a91f..ef4aa13a66 100644 --- a/src/mesh/boundary_region.cxx +++ b/src/mesh/boundary_region.cxx @@ -202,3 +202,6 @@ void BoundaryRegionYUp::nextY() { bool BoundaryRegionYUp::isDone() { return (x > xe) || (y >= localmesh->LocalNy); // Return true if gone out of the boundary } + +BoundaryRegionIter BoundaryRegion::begin() { return BoundaryRegionIter(this, false); } +BoundaryRegionIter BoundaryRegion::end() { return BoundaryRegionIter(this, true); } From 9bdfef5f553b1a783b9f4295be05ccbe8cba247a Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 9 Apr 2024 14:40:17 +0200 Subject: [PATCH 03/62] Provide a new boundary iterator based on RangeIterator The boundary region does not match what is done for the parallel case, thus porting it to a iterator does not work. --- include/bout/boundary_iterator.hxx | 117 +++++++++++++++++++++++++++++ include/bout/boundary_region.hxx | 90 ---------------------- 2 files changed, 117 insertions(+), 90 deletions(-) create mode 100644 include/bout/boundary_iterator.hxx diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx new file mode 100644 index 0000000000..93f02c004d --- /dev/null +++ b/include/bout/boundary_iterator.hxx @@ -0,0 +1,117 @@ +#pragma once + +#include "bout/mesh.hxx" +#include "bout/sys/parallel_stencils.hxx" +#include "bout/sys/range.hxx" + +class BoundaryRegionIter { +public: + BoundaryRegionIter(int x, int y, int bx, int by, Mesh* mesh) + : dir(bx + by), x(x), y(y), bx(bx), by(by), localmesh(mesh) { + ASSERT3(bx * by == 0); + } + bool operator!=(const BoundaryRegionIter& rhs) { return ind() != rhs.ind(); } + + Ind3D ind() const { return xyz2ind(x, y, z); } + BoundaryRegionIter& operator++() { + ASSERT3(z < nz()); + z++; + if (z == nz()) { + z = 0; + _next(); + } + return *this; + } + virtual void _next() = 0; + BoundaryRegionIter& operator*() { return *this; } + + void dirichlet_o2(Field3D& f, BoutReal value) const { + ynext(f) = parallel_stencil::dirichlet_o2(1, f[ind()], 0.5, value); + } + + BoutReal extrapolate_grad_o2(const Field3D& f) const { return f[ind()] - yprev(f); } + + BoutReal extrapolate_sheath_o2(const Field3D& f) const { + return (f[ind()] * 3 - yprev(f)) * 0.5; + } + + BoutReal extrapolate_next_o2(const Field3D& f) const { return 2 * f[ind()] - yprev(f); } + + BoutReal + extrapolate_next_o2(const std::function& f) const { + return 2 * f(0, ind()) - f(0, ind().yp(-by).xp(-bx)); + } + + BoutReal interpolate_sheath(const Field3D& f) const { + return (f[ind()] + ynext(f)) * 0.5; + } + + 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)]; } + + const int dir; + +protected: + int z{0}; + int x; + int y; + const int bx; + const int by; + +private: + Mesh* localmesh; + int nx() const { return localmesh->LocalNx; } + int ny() const { return localmesh->LocalNy; } + int nz() const { return localmesh->LocalNz; } + + Ind3D xyz2ind(int x, int y, int z) const { + return Ind3D{(x * ny() + y) * nz() + z, ny(), nz()}; + } +}; + +class BoundaryRegionIterY : public BoundaryRegionIter { +public: + BoundaryRegionIterY(RangeIterator r, int y, int dir, bool is_end, Mesh* mesh) + : BoundaryRegionIter(r.ind, y, 0, dir, mesh), r(r), is_end(is_end) {} + + bool operator!=(const BoundaryRegionIterY& rhs) { + ASSERT2(y == rhs.y); + if (is_end) { + if (rhs.is_end) { + return false; + } + return !rhs.r.isDone(); + } + if (rhs.is_end) { + return !r.isDone(); + } + return x != rhs.x; + } + + virtual void _next() override { + ++r; + x = r.ind; + } + +private: + RangeIterator r; + bool is_end; +}; + +class NewBoundaryRegionY { +public: + NewBoundaryRegionY(Mesh* mesh, bool lower, RangeIterator r) + : mesh(mesh), lower(lower), r(std::move(r)) {} + BoundaryRegionIterY begin(bool begin = true) { + return BoundaryRegionIterY(r, lower ? mesh->ystart : mesh->yend, lower ? -1 : +1, + !begin, mesh); + } + BoundaryRegionIterY end() { return begin(false); } + +private: + Mesh* mesh; + bool lower; + RangeIterator r; +}; diff --git a/include/bout/boundary_region.hxx b/include/bout/boundary_region.hxx index 6e7a939d3c..22956d1d4a 100644 --- a/include/bout/boundary_region.hxx +++ b/include/bout/boundary_region.hxx @@ -65,7 +65,6 @@ public: isDone() = 0; ///< Returns true if outside domain. Can use this with nested nextX, nextY }; -class BoundaryRegionIter; /// Describes a region of the boundary, and a means of iterating over it class BoundaryRegion : public BoundaryRegionBase { public: @@ -84,95 +83,6 @@ public: virtual void next1d() = 0; ///< Loop over the innermost elements virtual void nextX() = 0; ///< Just loop over X virtual void nextY() = 0; ///< Just loop over Y - - BoundaryRegionIter begin(); - BoundaryRegionIter end(); -}; - -class BoundaryRegionIter { -public: - BoundaryRegionIter(BoundaryRegion* rgn, bool is_end) - : rgn(rgn), is_end(is_end), dir(rgn->bx + rgn->by) { - //static_assert(std::is_base_of, "BoundaryRegionIter only works on BoundaryRegion"); - - // Ensure only one is non-zero - ASSERT3(rgn->bx * rgn->by == 0); - if (!is_end) { - rgn->first(); - } - } - bool operator!=(const BoundaryRegionIter& rhs) { - if (is_end) { - if (rhs.is_end || rhs.rgn->isDone()) { - return false; - } else { - return true; - } - } - if (rhs.is_end) { - return !rgn->isDone(); - } - return ind() != rhs.ind(); - } - - Ind3D ind() const { return xyz2ind(rgn->x - rgn->bx, rgn->y - rgn->by, z); } - BoundaryRegionIter& operator++() { - ASSERT3(z < nz()); - z++; - if (z == nz()) { - z = 0; - rgn->next(); - } - return *this; - } - BoundaryRegionIter& operator*() { return *this; } - - void dirichlet_o2(Field3D& f, BoutReal value) const { - ynext(f) = parallel_stencil::dirichlet_o2(1, f[ind()], 0.5, value); - } - - BoutReal extrapolate_grad_o2(const Field3D& f) const { return f[ind()] - yprev(f); } - - BoutReal extrapolate_sheath_o2(const Field3D& f) const { - return (f[ind()] * 3 - yprev(f)) * 0.5; - } - - BoutReal extrapolate_next_o2(const Field3D& f) const { return 2 * f[ind()] - yprev(f); } - - BoutReal - extrapolate_next_o2(const std::function& f) const { - return 2 * f(0, ind()) - f(0, ind().yp(-rgn->by).xp(-rgn->bx)); - } - - BoutReal interpolate_sheath(const Field3D& f) const { - return (f[ind()] + ynext(f)) * 0.5; - } - - BoutReal& ynext(Field3D& f) const { return f[ind().yp(rgn->by).xp(rgn->bx)]; } - const BoutReal& ynext(const Field3D& f) const { - return f[ind().yp(rgn->by).xp(rgn->bx)]; - } - BoutReal& yprev(Field3D& f) const { return f[ind().yp(-rgn->by).xp(-rgn->bx)]; } - const BoutReal& yprev(const Field3D& f) const { - return f[ind().yp(-rgn->by).xp(-rgn->bx)]; - } - -private: - BoundaryRegion* rgn; - const bool is_end; - int z{0}; - -public: - const int dir; - -private: - int nx() const { return rgn->localmesh->LocalNx; } - int ny() const { return rgn->localmesh->LocalNy; } - int nz() const { return rgn->localmesh->LocalNz; } - - Ind3D xyz2ind(int x, int y, int z) const { - return Ind3D{(x * ny() + y) * nz() + z, ny(), nz()}; - } }; class BoundaryRegionXIn : public BoundaryRegion { From c51a85f94a6590189749f3421b08b337a560de52 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 2 Jul 2024 14:41:42 +0200 Subject: [PATCH 04/62] Only check hasBndry*Y if they would be included 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. --- src/invert/laplace/invert_laplace.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/invert/laplace/invert_laplace.cxx b/src/invert/laplace/invert_laplace.cxx index e9182042de..9af6fd20f9 100644 --- a/src/invert/laplace/invert_laplace.cxx +++ b/src/invert/laplace/invert_laplace.cxx @@ -218,10 +218,10 @@ Field3D Laplacian::solve(const Field3D& b, const Field3D& x0) { // Setting the start and end range of the y-slices int ys = localmesh->ystart, ye = localmesh->yend; - if (localmesh->hasBndryLowerY() && include_yguards) { + if (include_yguards && localmesh->hasBndryLowerY()) { ys = 0; // Mesh contains a lower boundary } - if (localmesh->hasBndryUpperY() && include_yguards) { + if (include_yguards && localmesh->hasBndryUpperY()) { ye = localmesh->LocalNy - 1; // Contains upper boundary } From 06d581bcac128616d696dccb4fad9dc9d208ede3 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 5 Jul 2024 11:42:18 +0200 Subject: [PATCH 05/62] ensure we dont mix non-fci BCs with fci --- src/mesh/impls/bout/boutmesh.cxx | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index d0bdcf33d8..99ec63ff25 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -2905,6 +2905,9 @@ void BoutMesh::addBoundaryRegions() { } RangeIterator BoutMesh::iterateBndryLowerInnerY() const { + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } int xs = 0; int xe = LocalNx - 1; @@ -2940,6 +2943,9 @@ RangeIterator BoutMesh::iterateBndryLowerInnerY() const { } RangeIterator BoutMesh::iterateBndryLowerOuterY() const { + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } int xs = 0; int xe = LocalNx - 1; @@ -2974,6 +2980,10 @@ RangeIterator BoutMesh::iterateBndryLowerOuterY() const { } RangeIterator BoutMesh::iterateBndryLowerY() const { + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } + int xs = 0; int xe = LocalNx - 1; if ((DDATA_INDEST >= 0) && (DDATA_XSPLIT > xstart)) { @@ -3003,6 +3013,10 @@ RangeIterator BoutMesh::iterateBndryLowerY() const { } RangeIterator BoutMesh::iterateBndryUpperInnerY() const { + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } + int xs = 0; int xe = LocalNx - 1; @@ -3037,6 +3051,10 @@ RangeIterator BoutMesh::iterateBndryUpperInnerY() const { } RangeIterator BoutMesh::iterateBndryUpperOuterY() const { + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } + int xs = 0; int xe = LocalNx - 1; @@ -3071,6 +3089,10 @@ RangeIterator BoutMesh::iterateBndryUpperOuterY() const { } RangeIterator BoutMesh::iterateBndryUpperY() const { + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } + int xs = 0; int xe = LocalNx - 1; if ((UDATA_INDEST >= 0) && (UDATA_XSPLIT > xstart)) { From 7f1d26fbaea4edd33af2a8d2e158a8fb2d6eef4d Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 9 Aug 2024 13:45:23 +0200 Subject: [PATCH 06/62] rename to include order --- include/bout/parallel_boundary_region.hxx | 60 +++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 44e5cc562e..ad66c992a6 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -98,6 +98,66 @@ public: return f[ind()] * (1 + length()) - f.ynext(-dir)[ind().yp(-dir)] * length(); } + inline BoutReal + extrapolate_sheath_o1(const std::function& f) const { + return f(0, ind()); + } + inline BoutReal + extrapolate_sheath_o2(const std::function& f) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return extrapolate_sheath_o1(f); + } + return f(0, ind()) * (1 + length()) - f(-dir, ind().yp(-dir)) * length(); + } + + inline BoutReal interpolate_sheath_o1(const Field3D& f) const { + return f[ind()] * (1 - length()) + ynext(f) * length(); + } + + inline BoutReal extrapolate_next_o1(const Field3D& f) const { return f[ind()]; } + inline BoutReal extrapolate_next_o2(const Field3D& f) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return extrapolate_next_o1(f); + } + return f[ind()] * 2 - f.ynext(-dir)[ind().yp(-dir)]; + } + + inline BoutReal + extrapolate_next_o1(const std::function& f) const { + return f(0, ind()); + } + inline BoutReal + extrapolate_next_o2(const std::function& f) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return extrapolate_sheath_o1(f); + } + return f(0, ind()) * 2 - f(-dir, ind().yp(-dir)); + } + + // extrapolate the gradient into the boundary + inline BoutReal extrapolate_grad_o1(const Field3D& f) const { return 0; } + inline BoutReal extrapolate_grad_o2(const Field3D& f) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return extrapolate_grad_o1(f); + } + return f[ind()] - f.ynext(-dir)[ind().yp(-dir)]; + } + + BoundaryRegionParIterBase& operator*() { return *this; } + + BoundaryRegionParIterBase& operator++() { + ++bndry_position; + return *this; + } + + bool operator!=(const BoundaryRegionParIterBase& rhs) { + return bndry_position != rhs.bndry_position; + } + // dirichlet boundary code void dirichlet_o1(Field3D& f, BoutReal value) const { f.ynext(dir)[ind().yp(dir)] = value; From 84f8c8301b8fac38934740aefa8f0dc0a185baec Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 9 Aug 2024 13:50:34 +0200 Subject: [PATCH 07/62] Fix: remove broken code BoundaryRegionIter has been delted --- src/mesh/boundary_region.cxx | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/mesh/boundary_region.cxx b/src/mesh/boundary_region.cxx index ef4aa13a66..700ef8a91f 100644 --- a/src/mesh/boundary_region.cxx +++ b/src/mesh/boundary_region.cxx @@ -202,6 +202,3 @@ void BoundaryRegionYUp::nextY() { bool BoundaryRegionYUp::isDone() { return (x > xe) || (y >= localmesh->LocalNy); // Return true if gone out of the boundary } - -BoundaryRegionIter BoundaryRegion::begin() { return BoundaryRegionIter(this, false); } -BoundaryRegionIter BoundaryRegion::end() { return BoundaryRegionIter(this, true); } From ea05bae1f4d1f3192e8e88ddf409b456744f3a5e Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 27 Sep 2024 11:06:33 +0200 Subject: [PATCH 08/62] Convert test to new iterator scheeme --- tests/integrated/test-fci-boundary/get_par_bndry.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/integrated/test-fci-boundary/get_par_bndry.cxx b/tests/integrated/test-fci-boundary/get_par_bndry.cxx index 8e3cfac2f7..52d07f067e 100644 --- a/tests/integrated/test-fci-boundary/get_par_bndry.cxx +++ b/tests/integrated/test-fci-boundary/get_par_bndry.cxx @@ -24,11 +24,11 @@ int main(int argc, char** argv) { const auto boundary = static_cast(i); const auto boundary_name = toString(boundary); mesh->communicate(fields[i]); - for (const auto& bndry_par : mesh->getBoundariesPar(boundary)) { - output.write("{:s} region\n", boundary_name); - for (bndry_par->first(); !bndry_par->isDone(); bndry_par->next()) { - fields[i][bndry_par->ind()] += 1; - output.write("{:s} increment\n", boundary_name); + for (auto& bndry_par : mesh->getBoundariesPar(static_cast(i))) { + output.write("{:s} region\n", toString(static_cast(i))); + for (const auto& pnt : *bndry_par) { + fields[i][pnt.ind()] += 1; + output.write("{:s} increment\n", toString(static_cast(i))); } } output.write("{:s} done\n", boundary_name); From 17f8efa868b0195e0c81d35e99a2a4eab028f969 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 1 Oct 2024 15:27:21 +0200 Subject: [PATCH 09/62] Add Field2D version for 2D metrics --- include/bout/parallel_boundary_region.hxx | 107 +++++++++++++++++++++- 1 file changed, 105 insertions(+), 2 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index ad66c992a6..9eedc7847d 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -216,8 +216,111 @@ private: // BoutReal get(const Field3D& f, int off) const BoutReal& ynext(const Field3D& f) const { return f.ynext(dir)[ind().yp(dir)]; } BoutReal& ynext(Field3D& f) const { return f.ynext(dir)[ind().yp(dir)]; } - const BoutReal& yprev(const Field3D& f) const { return f.ynext(-dir)[ind().yp(-dir)]; } - BoutReal& yprev(Field3D& f) const { return f.ynext(-dir)[ind().yp(-dir)]; } + + const BoutReal& yprev(const Field3D& f) const { + ASSERT3(valid() > 0); + return f.ynext(-dir)[ind().yp(-dir)]; + } + BoutReal& yprev(Field3D& f) const { + ASSERT3(valid() > 0); + return f.ynext(-dir)[ind().yp(-dir)]; + } + +#if BOUT_USE_METRIC_3D == 0 + const BoutReal& ynext(const Field2D& f) const { return f.ynext(dir)[ind().yp(dir)]; } + BoutReal& ynext(Field2D& f) const { return f.ynext(dir)[ind().yp(dir)]; } + + const BoutReal& yprev(const Field2D& f) const { + ASSERT3(valid() > 0); + return f.ynext(-dir)[ind().yp(-dir)]; + } + BoutReal& yprev(Field2D& f) const { + ASSERT3(valid() > 0); + return f.ynext(-dir)[ind().yp(-dir)]; + } +#endif + +private: + const IndicesVec& bndry_points; + IndicesIter bndry_position; + + constexpr static BoutReal small_value = 1e-2; + +public: + const int dir; + Mesh* localmesh; +}; +} // namespace parallel_boundary_region +} // namespace bout +using BoundaryRegionParIter = bout::parallel_boundary_region::BoundaryRegionParIterBase< + bout::parallel_boundary_region::IndicesVec, + bout::parallel_boundary_region::IndicesIter>; +using BoundaryRegionParIterConst = + bout::parallel_boundary_region::BoundaryRegionParIterBase< + const bout::parallel_boundary_region::IndicesVec, + bout::parallel_boundary_region::IndicesIterConst>; + +class BoundaryRegionPar : public BoundaryRegionBase { +public: + BoundaryRegionPar(const std::string& name, int dir, Mesh* passmesh) + : BoundaryRegionBase(name, passmesh), dir(dir) { + ASSERT0(std::abs(dir) == 1); + BoundaryRegionBase::isParallel = true; + } + BoundaryRegionPar(const std::string& name, BndryLoc loc, int dir, Mesh* passmesh) + : BoundaryRegionBase(name, loc, passmesh), dir(dir) { + BoundaryRegionBase::isParallel = true; + ASSERT0(std::abs(dir) == 1); + } + + /// Add a point to the boundary + void add_point(Ind3D ind, BoutReal x, BoutReal y, BoutReal z, BoutReal length, + char valid) { + bndry_points.push_back({ind, {x, y, z}, length, valid}); + } + void add_point(int ix, int iy, int iz, BoutReal x, BoutReal y, BoutReal z, + BoutReal length, char valid) { + bndry_points.push_back({xyz2ind(ix, iy, iz, localmesh), {x, y, z}, length, valid}); + } + + // final, so they can be inlined + void first() final { bndry_position = std::begin(bndry_points); } + void next() final { ++bndry_position; } + bool isDone() final { return (bndry_position == std::end(bndry_points)); } + + bool contains(const BoundaryRegionPar& bndry) const { + return std::binary_search(std::begin(bndry_points), std::end(bndry_points), + *bndry.bndry_position, + [](const bout::parallel_boundary_region::Indices& i1, + const bout::parallel_boundary_region::Indices& i2) { + return i1.index < i2.index; + }); + } + + // setter + void setValid(char val) { bndry_position->valid = val; } + + // BoundaryRegionParIterConst begin() const { + // return BoundaryRegionParIterConst(bndry_points, bndry_points.begin(), dir); + // } + // BoundaryRegionParIterConst end() const { + // return BoundaryRegionParIterConst(bndry_points, bndry_points.begin(), dir); + // } + BoundaryRegionParIter begin() { + return BoundaryRegionParIter(bndry_points, bndry_points.begin(), dir, localmesh); + } + BoundaryRegionParIter end() { + return BoundaryRegionParIter(bndry_points, bndry_points.end(), dir, localmesh); + } + + const int dir; + +private: + /// Vector of points in the boundary + bout::parallel_boundary_region::IndicesVec bndry_points; + /// Current position in the boundary points + bout::parallel_boundary_region::IndicesIter bndry_position; + static Ind3D xyz2ind(int x, int y, int z, Mesh* mesh) { const int ny = mesh->LocalNy; const int nz = mesh->LocalNz; From 0db3052e66f96f19047edb72b46af073b79b79c8 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 1 Oct 2024 15:52:40 +0200 Subject: [PATCH 10/62] Add setYPrevIfValid --- include/bout/parallel_boundary_region.hxx | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 9eedc7847d..f370295481 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -225,6 +225,11 @@ private: ASSERT3(valid() > 0); return f.ynext(-dir)[ind().yp(-dir)]; } + void setYPrevIfValid(Field3D& f, BoutReal val) const { + if (valid() > 0) { + yprev(f) = val; + } + } #if BOUT_USE_METRIC_3D == 0 const BoutReal& ynext(const Field2D& f) const { return f.ynext(dir)[ind().yp(dir)]; } From d2f2cc80383ce87f3a89031c619c6046e64d5a7d Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 13 Jan 2025 09:53:15 +0100 Subject: [PATCH 11/62] Add offset to parallel boundary region This allows to extend the boundary code to place the boundary further away from the boundary. --- include/bout/parallel_boundary_region.hxx | 17 +++++++++++++---- src/mesh/parallel/fci.cxx | 15 +++++++++------ src/mesh/parallel/shiftedmetricinterp.cxx | 8 ++++---- 3 files changed, 26 insertions(+), 14 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index f370295481..ff7ae14685 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -34,6 +34,8 @@ class BoundaryRegionPar : public BoundaryRegionBase { // BoutReal angle; // How many points we can go in the opposite direction signed char valid; + signed char offset; + unsigned char abs_offset; }; using IndicesVec = std::vector; @@ -78,6 +80,8 @@ public: BoutReal s_z() const { return bndry_position->intersection.s_z; } BoutReal length() const { return bndry_position->length; } signed char valid() const { return bndry_position->valid; } + signed char offset() const { return bndry_position->offset; } + unsigned char abs_offset() const { return bndry_position->abs_offset; } // setter void setValid(signed char val) { bndry_position->valid = val; } @@ -280,12 +284,17 @@ public: /// Add a point to the boundary void add_point(Ind3D ind, BoutReal x, BoutReal y, BoutReal z, BoutReal length, - char valid) { - bndry_points.push_back({ind, {x, y, z}, length, valid}); + char valid, signed char offset) { + bndry_points.push_back({ind, + {x, y, z}, + length, + valid, + offset, + static_cast(std::abs(offset))}); } void add_point(int ix, int iy, int iz, BoutReal x, BoutReal y, BoutReal z, - BoutReal length, char valid) { - bndry_points.push_back({xyz2ind(ix, iy, iz, localmesh), {x, y, z}, length, valid}); + BoutReal length, char valid, signed char offset) { + add_point(xyz2ind(ix, iy, iz, localmesh), x, y, z, length, valid, offset); } // final, so they can be inlined diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 39654519af..3394a3e2af 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -182,7 +182,9 @@ FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy, BoutMask to_remove(map_mesh); const int xend = map_mesh->xstart - + ((map_mesh->xend - map_mesh->xstart + 1) * map_mesh->getNXPE()) - 1; + + (map_mesh->xend - map_mesh->xstart + 1) * map_mesh->getNXPE() - 1; + // Default to the maximum number of points + const int defValid{map_mesh->ystart - 1 + std::abs(offset)}; // Serial loop because call to BoundaryRegionPar::addPoint // (probably?) can't be done in parallel BOUT_FOR_SERIAL(i, xt_prime.getRegion("RGN_NOBNDRY")) { @@ -252,11 +254,12 @@ FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy, // need at least 2 points in the domain. ASSERT2(map_mesh->xend - map_mesh->xstart >= 2); auto boundary = (xt_prime[i] < map_mesh->xstart) ? inner_boundary : outer_boundary; - boundary->add_point(x, y, z, x + dx, y + (0.5 * offset_), - z + dz, // Intersection point in local index space - 0.5, // Distance to intersection - 1 // Default to that there is a point in the other direction - ); + if (!boundary->contains(x, y, z)) { + boundary->add_point(x, y, z, x + dx, y + offset - sgn(offset) * 0.5, + z + dz, // Intersection point in local index space + std::abs(offset) - 0.5, // Distance to intersection + defValid, offset); + } } region_no_boundary = region_no_boundary.mask(to_remove); diff --git a/src/mesh/parallel/shiftedmetricinterp.cxx b/src/mesh/parallel/shiftedmetricinterp.cxx index c71618ab19..22a8f0748a 100644 --- a/src/mesh/parallel/shiftedmetricinterp.cxx +++ b/src/mesh/parallel/shiftedmetricinterp.cxx @@ -146,7 +146,7 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in, 0.25 * (1 // dy/2 + dy(it.ind, mesh.yend + 1) / dy(it.ind, mesh.yend)), // length - yvalid); + yvalid, 1); } } auto backward_boundary_xin = std::make_shared( @@ -162,7 +162,7 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in, 0.25 * (1 // dy/2 + dy(it.ind, mesh.ystart - 1) / dy(it.ind, mesh.ystart)), - yvalid); + yvalid, -1); } } // Create regions for parallel boundary conditions @@ -179,7 +179,7 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in, 0.25 * (1 // dy/2 + dy(it.ind, mesh.yend + 1) / dy(it.ind, mesh.yend)), - yvalid); + yvalid, 1); } } auto backward_boundary_xout = std::make_shared( @@ -195,7 +195,7 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in, 0.25 * (dy(it.ind, mesh.ystart - 1) / dy(it.ind, mesh.ystart) // dy/2 + 1), - yvalid); + yvalid, -1); } } From 1f48a59e5b39235864169de4b0cc468de21fcdb9 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 13 Jan 2025 09:54:15 +0100 Subject: [PATCH 12/62] Add function to check wehther point is in the boundary --- include/bout/parallel_boundary_region.hxx | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index ff7ae14685..2d2e6829a8 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -311,6 +311,15 @@ public: }); } + bool contains(const int ix, const int iy, const int iz) const { + const auto i2 = xyz2ind(ix, iy, iz, localmesh); + for (auto i1 : bndry_points) { + if (i1.index == i2) { + return true; + } + } + return false; + } // setter void setValid(char val) { bndry_position->valid = val; } From 4af5a4d93888f180678b0bd3c5d4c23b48dddc4d Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 13 Jan 2025 09:55:59 +0100 Subject: [PATCH 13/62] Add setValid to BoundaryRegionParIterBase --- include/bout/parallel_boundary_region.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 2d2e6829a8..9af24770f7 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -84,7 +84,7 @@ public: unsigned char abs_offset() const { return bndry_position->abs_offset; } // setter - void setValid(signed char val) { bndry_position->valid = val; } + void setValid(signed char valid) { bndry_position->valid = valid; } bool contains(const BoundaryRegionPar& bndry) const { return std::binary_search( From e5fbcaf08ed7016e708622946bf5a98185f8097f Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 13 Jan 2025 09:57:32 +0100 Subject: [PATCH 14/62] Reimplement ynext for the boundary region 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 --- include/bout/parallel_boundary_region.hxx | 56 +++++++++++++++++------ 1 file changed, 41 insertions(+), 15 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 9af24770f7..b03a0e2f01 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -212,23 +212,49 @@ public: parallel_stencil::neumann_o3(1 - length(), value, 1, f[ind()], 2, yprev(f)); } - const int dir; - -private: - constexpr static BoutReal small_value = 1e-2; - - // BoutReal get(const Field3D& f, int off) - const BoutReal& ynext(const Field3D& f) const { return f.ynext(dir)[ind().yp(dir)]; } - BoutReal& ynext(Field3D& f) const { return f.ynext(dir)[ind().yp(dir)]; } - - const BoutReal& yprev(const Field3D& f) const { - ASSERT3(valid() > 0); - return f.ynext(-dir)[ind().yp(-dir)]; + template + BoutReal& getAt(Field3D& f, int off) const { + if constexpr (check) { + ASSERT3(valid() > -off - 2); + } + auto _off = offset() + off * dir; + return f.ynext(_off)[ind().yp(_off)]; } - BoutReal& yprev(Field3D& f) const { - ASSERT3(valid() > 0); - return f.ynext(-dir)[ind().yp(-dir)]; + template + const BoutReal& getAt(const Field3D& f, int off) const { + if constexpr (check) { + ASSERT3(valid() > -off - 2); + } + auto _off = offset() + off * dir; + return f.ynext(_off)[ind().yp(_off)]; + } + + const BoutReal& ynext(const Field3D& f) const { return getAt(f, 0); } + BoutReal& ynext(Field3D& f) const { return getAt(f, 0); } + const BoutReal& ythis(const Field3D& f) const { return getAt(f, -1); } + BoutReal& ythis(Field3D& f) const { return getAt(f, -1); } + const BoutReal& yprev(const Field3D& f) const { return getAt(f, -2); } + BoutReal& yprev(Field3D& f) const { return getAt(f, -2); } + + template + BoutReal getAt(const std::function& f, + int off) const { + if constexpr (check) { + ASSERT3(valid() > -off - 2); + } + auto _off = offset() + off * dir; + return f(_off, ind().yp(_off)); } + BoutReal ynext(const std::function& f) const { + return getAt(f, 0); + } + BoutReal ythis(const std::function& f) const { + return getAt(f, -1); + } + BoutReal yprev(const std::function& f) const { + return getAt(f, -2); + } + void setYPrevIfValid(Field3D& f, BoutReal val) const { if (valid() > 0) { yprev(f) = val; From ef7b03df582ed1ba825f8ac0dbeedbd680e4c840 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 13 Jan 2025 09:58:41 +0100 Subject: [PATCH 15/62] Fix extrapolaton / interpolation --- include/bout/parallel_boundary_region.hxx | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index b03a0e2f01..0ced58d65e 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -93,18 +93,18 @@ public: } // extrapolate a given point to the boundary - BoutReal extrapolate_o1(const Field3D& f) const { return f[ind()]; } - BoutReal extrapolate_o2(const Field3D& f) const { + BoutReal extrapolate_sheath_o1(const Field3D& f) const { return ythis(f); } + BoutReal extrapolate_sheath_o2(const Field3D& f) const { ASSERT3(valid() >= 0); if (valid() < 1) { return extrapolate_o1(f); } - return f[ind()] * (1 + length()) - f.ynext(-dir)[ind().yp(-dir)] * length(); + return ythis(f) * (1 + length()) - yprev(f) * length(); } inline BoutReal extrapolate_sheath_o1(const std::function& f) const { - return f(0, ind()); + return ythis(f); } inline BoutReal extrapolate_sheath_o2(const std::function& f) const { @@ -112,25 +112,25 @@ public: if (valid() < 1) { return extrapolate_sheath_o1(f); } - return f(0, ind()) * (1 + length()) - f(-dir, ind().yp(-dir)) * length(); + return ythis(f) * (1 + length()) - yprev(f) * length(); } inline BoutReal interpolate_sheath_o1(const Field3D& f) const { - return f[ind()] * (1 - length()) + ynext(f) * length(); + return ythis(f) * (1 - length()) + ynext(f) * length(); } - inline BoutReal extrapolate_next_o1(const Field3D& f) const { return f[ind()]; } + inline BoutReal extrapolate_next_o1(const Field3D& f) const { return ythis(f); } inline BoutReal extrapolate_next_o2(const Field3D& f) const { ASSERT3(valid() >= 0); if (valid() < 1) { return extrapolate_next_o1(f); } - return f[ind()] * 2 - f.ynext(-dir)[ind().yp(-dir)]; + return ythis(f) * 2 - yprev(f); } inline BoutReal extrapolate_next_o1(const std::function& f) const { - return f(0, ind()); + return ythis(f); } inline BoutReal extrapolate_next_o2(const std::function& f) const { @@ -138,7 +138,7 @@ public: if (valid() < 1) { return extrapolate_sheath_o1(f); } - return f(0, ind()) * 2 - f(-dir, ind().yp(-dir)); + return ythis(f) * 2 - yprev(f); } // extrapolate the gradient into the boundary @@ -148,7 +148,7 @@ public: if (valid() < 1) { return extrapolate_grad_o1(f); } - return f[ind()] - f.ynext(-dir)[ind().yp(-dir)]; + return ythis(f) - ynext(f); } BoundaryRegionParIterBase& operator*() { return *this; } @@ -346,6 +346,7 @@ public: } return false; } + // setter void setValid(char val) { bndry_position->valid = val; } From 6dba751e085d55bdba46ff84b70115a1f2a0034d Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 13 Jan 2025 11:12:38 +0100 Subject: [PATCH 16/62] Fix parallel boundary to interpolate into the boundary Previously only the first boundary point was set --- include/bout/parallel_boundary_region.hxx | 32 ++++++++++++++++------- 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 0ced58d65e..d3082f4d45 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -162,17 +162,20 @@ public: return bndry_position != rhs.bndry_position; } +#define ITER() for (int i = 0; i < localmesh->ystart - abs_offset(); ++i) // dirichlet boundary code void dirichlet_o1(Field3D& f, BoutReal value) const { - f.ynext(dir)[ind().yp(dir)] = value; + ITER() { getAt(f, i) = value; } } void dirichlet_o2(Field3D& f, BoutReal value) const { if (length() < small_value) { return dirichlet_o1(f, value); } - ynext(f) = parallel_stencil::dirichlet_o2(1, f[ind()], 1 - length(), value); - // ynext(f) = f[ind()] * (1 + 1/length()) + value / length(); + ITER() { + getAt(f, i) = + parallel_stencil::dirichlet_o2(i + 1, ythis(f), i + 1 - length(), value); + } } void dirichlet_o3(Field3D& f, BoutReal value) const { @@ -181,17 +184,24 @@ public: return dirichlet_o2(f, value); } if (length() < small_value) { - ynext(f) = parallel_stencil::dirichlet_o2(2, yprev(f), 1 - length(), value); + ITER() { + getAt(f, i) = + parallel_stencil::dirichlet_o2(i + 2, yprev(f), i + 1 - length(), value); + } } else { - ynext(f) = - parallel_stencil::dirichlet_o3(2, yprev(f), 1, f[ind()], 1 - length(), value); + ITER() { + getAt(f, i) = parallel_stencil::dirichlet_o3(i + 2, yprev(f), i + 1, ythis(f), + i + 1 - length(), value); + } } } // NB: value needs to be scaled by dy // neumann_o1 is actually o2 if we would use an appropriate one-sided stencil. // But in general we do not, and thus for normal C2 stencils, this is 1st order. - void neumann_o1(Field3D& f, BoutReal value) const { ynext(f) = f[ind()] + value; } + void neumann_o1(Field3D& f, BoutReal value) const { + ITER() { getAt(f, i) = ythis(f) + value; } + } // NB: value needs to be scaled by dy void neumann_o2(Field3D& f, BoutReal value) const { @@ -199,7 +209,7 @@ public: if (valid() < 1) { return neumann_o1(f, value); } - ynext(f) = yprev(f) + 2 * value; + ITER() { getAt(f, i) = yprev(f) + 2 * value; } } // NB: value needs to be scaled by dy @@ -208,8 +218,10 @@ public: if (valid() < 1) { return neumann_o1(f, value); } - ynext(f) = - parallel_stencil::neumann_o3(1 - length(), value, 1, f[ind()], 2, yprev(f)); + ITER() { + getAt(f, i) = parallel_stencil::neumann_o3(i + 1 - length(), value, i + 1, ythis(f), + 2, yprev(f)); + } } template From f105c3f725e6678597f444271c8957f8a71ddbeb Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 13 Jan 2025 11:13:14 +0100 Subject: [PATCH 17/62] Ensure data is sorted --- include/bout/parallel_boundary_region.hxx | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index d3082f4d45..5d6150d94c 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -323,6 +323,9 @@ public: /// Add a point to the boundary void add_point(Ind3D ind, BoutReal x, BoutReal y, BoutReal z, BoutReal length, char valid, signed char offset) { + if (!bndry_points.empty() && bndry_points.back().index > ind) { + is_sorted = false; + } bndry_points.push_back({ind, {x, y, z}, length, @@ -341,6 +344,7 @@ public: bool isDone() final { return (bndry_position == std::end(bndry_points)); } bool contains(const BoundaryRegionPar& bndry) const { + ASSERT2(is_sorted); return std::binary_search(std::begin(bndry_points), std::end(bndry_points), *bndry.bndry_position, [](const bout::parallel_boundary_region::Indices& i1, @@ -388,6 +392,7 @@ private: const int nz = mesh->LocalNz; return Ind3D{(x * ny + y) * nz + z, ny, nz}; } + bool is_sorted{true}; }; #endif // BOUT_PAR_BNDRY_H From f7c1acaf2fcf3b0f488da73351aea98d7d801d65 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 13 Jan 2025 11:13:59 +0100 Subject: [PATCH 18/62] Calculate valid for the general case --- src/mesh/parallel/fci.cxx | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 3394a3e2af..d3fcc8d04c 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -357,17 +357,20 @@ FCITransform::FCITransform(Mesh& mesh, const Coordinates::FieldMetric& dy, bool field_line_maps.emplace_back(mesh, dy, options, -offset, backward_boundary_xin, backward_boundary_xout, zperiodic); } - ASSERT0(mesh.ystart == 1); const std::array bndries = {forward_boundary_xin, forward_boundary_xout, backward_boundary_xin, backward_boundary_xout}; - for (const auto& bndry : bndries) { + for (auto& bndry : bndries) { for (const auto& bndry2 : bndries) { if (bndry->dir == bndry2->dir) { continue; } - for (bndry->first(); !bndry->isDone(); bndry->next()) { - if (bndry2->contains(*bndry)) { - bndry->setValid(0); + for (auto pnt : *bndry) { + for (auto pnt2 : *bndry2) { +#warning this could likely be done faster + if (pnt.ind() == pnt2.ind()) { + pnt.setValid( + static_cast(std::abs((pnt2.offset() - pnt.offset())) - 2)); + } } } } From 55bdbcfd0eedca426be35fbdb173bc85611ea062 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 20 Jan 2025 16:29:05 +0100 Subject: [PATCH 19/62] fix parallel neumann BC --- include/bout/parallel_boundary_region.hxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 5d6150d94c..936a57e9d9 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -200,7 +200,7 @@ public: // neumann_o1 is actually o2 if we would use an appropriate one-sided stencil. // But in general we do not, and thus for normal C2 stencils, this is 1st order. void neumann_o1(Field3D& f, BoutReal value) const { - ITER() { getAt(f, i) = ythis(f) + value; } + ITER() { getAt(f, i) = ythis(f) + value * (i + 1); } } // NB: value needs to be scaled by dy @@ -209,14 +209,14 @@ public: if (valid() < 1) { return neumann_o1(f, value); } - ITER() { getAt(f, i) = yprev(f) + 2 * value; } + ITER() { getAt(f, i) = yprev(f) + (2 + i) * value; } } // NB: value needs to be scaled by dy void neumann_o3(Field3D& f, BoutReal value) const { ASSERT3(valid() >= 0); if (valid() < 1) { - return neumann_o1(f, value); + return neumann_o2(f, value); } ITER() { getAt(f, i) = parallel_stencil::neumann_o3(i + 1 - length(), value, i + 1, ythis(f), From 6f40b4e6f6905d3fb61512fc69e7c354b9e3d2db Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 20 Jan 2025 16:30:54 +0100 Subject: [PATCH 20/62] Add limitFree Taken from hermes-3, adopted for higher order --- include/bout/parallel_boundary_region.hxx | 27 +++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 936a57e9d9..701298977a 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -15,6 +15,22 @@ * */ +inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc) { + if (fm < fc) { + return 1; // Neumann rather than increasing into boundary + } + if (fm < 1e-10) { + return 1; // Low / no density condition + } + BoutReal fp = fc / fm; +#if CHECKLEVEL >= 2 + if (!std::isfinite(fp)) { + throw BoutException("SheathBoundaryParallel limitFree: {}, {} -> {}", fm, fc, fp); + } +#endif + return fp; +} + class BoundaryRegionPar : public BoundaryRegionBase { struct RealPoint { @@ -224,6 +240,17 @@ public: } } + // extrapolate into the boundary using only monotonic decreasing values. + // f needs to be positive + void limitFree(Field3D& f) const { + const auto fac = valid() > 0 ? limitFreeScale(yprev(f), ythis(f)) : 1; + auto val = ythis(f); + ITER() { + val *= fac; + getAt(f, i) = val; + } + } + template BoutReal& getAt(Field3D& f, int off) const { if constexpr (check) { From 36c1f94f4f33983df959a7c30f52bab1fb376511 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 20 Jan 2025 16:31:18 +0100 Subject: [PATCH 21/62] add setAll to parallel BC --- include/bout/parallel_boundary_region.hxx | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 701298977a..1ae49360f2 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -251,6 +251,12 @@ public: } } + void setAll(Field3D& f, const BoutReal val) const { + for (int i = -localmesh->ystart; i <= localmesh->ystart; ++i) { + f.ynext(i)[ind().yp(i)] = val; + } + } + template BoutReal& getAt(Field3D& f, int off) const { if constexpr (check) { From c7cbf117f9db1c3c7264f5716b8d205e3d781e7e Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 31 Jan 2025 10:07:22 +0100 Subject: [PATCH 22/62] Move yboundary iterator from hermes to BOUT++ This allows to write code for FCI and non-FCI using templates. --- include/bout/boundary_iterator.hxx | 51 ++++++++++++++++++++++- include/bout/yboundary_regions.hxx | 66 ++++++++++++++++++++++++++++++ 2 files changed, 116 insertions(+), 1 deletion(-) create mode 100644 include/bout/yboundary_regions.hxx diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx index 93f02c004d..601f6a8a28 100644 --- a/include/bout/boundary_iterator.hxx +++ b/include/bout/boundary_iterator.hxx @@ -1,6 +1,7 @@ #pragma once #include "bout/mesh.hxx" +#include "bout/parallel_boundary_region.hxx" #include "bout/sys/parallel_stencils.hxx" #include "bout/sys/range.hxx" @@ -42,14 +43,62 @@ public: return 2 * f(0, ind()) - f(0, ind().yp(-by).xp(-bx)); } - BoutReal interpolate_sheath(const Field3D& f) const { + BoutReal interpolate_sheath_o1(const Field3D& f) const { return (f[ind()] + ynext(f)) * 0.5; } + BoutReal + extrapolate_sheath_o2(const std::function& f) const { + return 0.5 * (3 * f(0, ind()) - f(0, ind().yp(-by).xp(-bx))); + } + + void limitFree(Field3D& f) const { + const BoutReal fac = + bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f)); + BoutReal val = ythis(f); + for (int i = 1; i <= localmesh->ystart; ++i) { + val *= fac; + f[ind().yp(by * i).xp(bx * i)] = val; + } + } + + void neumann_o1(Field3D& f, BoutReal grad) const { + BoutReal val = ythis(f); + for (int i = 1; i <= localmesh->ystart; ++i) { + val += grad; + f[ind().yp(by * i).xp(bx * i)] = val; + } + } + + void neumann_o2(Field3D& f, BoutReal grad) const { + BoutReal val = yprev(f) + grad; + for (int i = 1; i <= localmesh->ystart; ++i) { + val += grad; + f[ind().yp(by * i).xp(bx * i)] = val; + } + } 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()]; } + const BoutReal& ythis(const Field3D& f) const { return f[ind()]; } + + void setYPrevIfValid(Field3D& f, BoutReal val) const { yprev(f) = val; } + void setAll(Field3D& f, const BoutReal val) const { + for (int i = -localmesh->ystart; i <= localmesh->ystart; ++i) { + f[ind().yp(by * i).xp(bx * i)] = val; + } + } + + int abs_offset() const { return 1; } + +#if BOUT_USE_METRIC_3D == 0 + BoutReal& ynext(Field2D& f) const { return f[ind().yp(by).xp(bx)]; } + const BoutReal& ynext(const Field2D& f) const { return f[ind().yp(by).xp(bx)]; } + BoutReal& yprev(Field2D& f) const { return f[ind().yp(-by).xp(-bx)]; } + const BoutReal& yprev(const Field2D& f) const { return f[ind().yp(-by).xp(-bx)]; } +#endif const int dir; diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx new file mode 100644 index 0000000000..e0e93e17f9 --- /dev/null +++ b/include/bout/yboundary_regions.hxx @@ -0,0 +1,66 @@ +#pragma once + +#include "./boundary_iterator.hxx" +#include "bout/parallel_boundary_region.hxx" + +class YBoundary { +public: + template + void iter_regions(const T& f) { + ASSERT1(is_init); + for (auto& region : boundary_regions) { + f(*region); + } + for (auto& region : boundary_regions_par) { + f(*region); + } + } + + template + void iter(const F& f) { + return iter_regions(f); + } + + void init(Options& options, Mesh* mesh = nullptr) { + if (mesh == nullptr) { + mesh = bout::globals::mesh; + } + + bool lower_y = options["lower_y"].doc("Boundary on lower y?").withDefault(true); + bool upper_y = options["upper_y"].doc("Boundary on upper y?").withDefault(true); + bool outer_x = options["outer_x"].doc("Boundary on outer x?").withDefault(true); + bool inner_x = + options["inner_x"].doc("Boundary on inner x?").withDefault(false); + + if (mesh->isFci()) { + if (outer_x) { + for (auto& bndry : mesh->getBoundariesPar(BoundaryParType::xout)) { + boundary_regions_par.push_back(bndry); + } + } + if (inner_x) { + for (auto& bndry : mesh->getBoundariesPar(BoundaryParType::xin)) { + boundary_regions_par.push_back(bndry); + } + } + } else { + if (lower_y) { + boundary_regions.push_back( + std::make_shared(mesh, true, mesh->iterateBndryLowerY())); + } + if (upper_y) { + boundary_regions.push_back(std::make_shared( + mesh, false, mesh->iterateBndryUpperY())); + } + } + is_init = true; + } + +private: + std::vector> boundary_regions_par; + std::vector> boundary_regions; + + bool is_init{false}; +}; + +extern YBoundary yboundary; From f9cad29209ce213da7bd78b2406de29aec829a69 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 3 Mar 2025 10:01:55 +0100 Subject: [PATCH 23/62] Add iter_pnts function Directly iterate over the points --- include/bout/yboundary_regions.hxx | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index e0e93e17f9..b72f8d32c0 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -5,8 +5,8 @@ class YBoundary { public: - template - void iter_regions(const T& f) { + template + void iter_regions(const F& f) { ASSERT1(is_init); for (auto& region : boundary_regions) { f(*region); @@ -15,6 +15,14 @@ public: f(*region); } } + template + void iter_pnts(const F& f) { + iter_regions([&](auto& region) { + for (auto& pnt : region) { + f(pnt); + } + } + } template void iter(const F& f) { From e6f004b33634f50dc9680bf7c100bcef6237b8fb Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 3 Mar 2025 10:02:46 +0100 Subject: [PATCH 24/62] Add example on how to replace RangeIterator with YBoundary --- include/bout/example_yboundary_regions.cxx | 65 ++++++++++++++++++++++ include/bout/yboundary_regions.hxx | 12 ++++ 2 files changed, 77 insertions(+) create mode 100644 include/bout/example_yboundary_regions.cxx diff --git a/include/bout/example_yboundary_regions.cxx b/include/bout/example_yboundary_regions.cxx new file mode 100644 index 0000000000..0e8d9b07dc --- /dev/null +++ b/include/bout/example_yboundary_regions.cxx @@ -0,0 +1,65 @@ +#include + +class yboundary_example_legacy { +public: + yboundary_example_legacy(Options* opt, const Field3D& N, const Field3D& V) + : N(N), V(V) { + Options& options = *opt; + lower_y = options["lower_y"].doc("Boundary on lower y?").withDefault(lower_y); + upper_y = options["upper_y"].doc("Boundary on upper y?").withDefault(upper_y); + } + + void rhs() { + BoutReal totalFlux = 0; + if (lower_y) { + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Calculate flux through surface [normalised m^-2 s^-1], + // should be positive since V < 0.0 + BoutReal flux = + -0.5 * (N(r.ind, mesh->ystart, jz) + N(r.ind, mesh->ystart - 1, jz)) * 0.5 + * (V(r.ind, mesh->ystart, jz) + V(r.ind, mesh->ystart - 1, jz)); + totalFlux += flux; + } + } + } + if (upper_y) { + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Calculate flux through surface [normalised m^-2 s^-1], + // should be positive since V < 0.0 + BoutReal flux = -0.5 * (N(r.ind, mesh->yend, jz) + N(r.ind, mesh->yend + 1, jz)) + * 0.5 + * (V(r.ind, mesh->yend, jz) + V(r.ind, mesh->yend + 1, jz)); + totalFlux += flux; + } + } + } + } + +private: + bool lower_y{true}; + bool upper_y{true}; + const Field3D& N; + const Field3D& V; +} + +class yboundary_example { +public: + yboundary_example(Options* opt, const Field3D& N, const Field3D& V) : N(N), V(V) { + // Set what kind of yboundaries you want to include + yboundary.init(opt); + } + + void rhs() { + BoutReal totalFlux = 0; + yboundary.iter_pnts([&](auto& pnt) { + BoutReal flux = pnt.interpolate_sheath_o1(N) * pnt.interpolate_sheath_o1(V); + }); + } + +private: + YBoundary ybounday; + const Field3D& N; + const Field3D& V; +}; diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index b72f8d32c0..833e988eb9 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -2,6 +2,18 @@ #include "./boundary_iterator.hxx" #include "bout/parallel_boundary_region.hxx" +/*! + * This class allows to simplify iterating over y-boundaries. + * + * It makes it easier to write code for FieldAligned boundaries, but if a bit + * care is taken the code also works with FluxCoordinateIndependent code. + * + * An example how to replace old code is given here: + * + * \example example_yboundary_regions.hxx + * This is an example how to use the YBoundary class to replace RangeIterator + * boundaries. + */ class YBoundary { public: From 3a0f2146636c05d57fd85b8bfdb1da883bedad54 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 3 Mar 2025 11:25:50 +0100 Subject: [PATCH 25/62] Move documentation to sphinx --- include/bout/example_yboundary_regions.cxx | 65 ---------------- include/bout/yboundary_regions.hxx | 20 ++--- manual/sphinx/user_docs/boundary_options.rst | 79 ++++++++++++++++++++ 3 files changed, 87 insertions(+), 77 deletions(-) delete mode 100644 include/bout/example_yboundary_regions.cxx diff --git a/include/bout/example_yboundary_regions.cxx b/include/bout/example_yboundary_regions.cxx deleted file mode 100644 index 0e8d9b07dc..0000000000 --- a/include/bout/example_yboundary_regions.cxx +++ /dev/null @@ -1,65 +0,0 @@ -#include - -class yboundary_example_legacy { -public: - yboundary_example_legacy(Options* opt, const Field3D& N, const Field3D& V) - : N(N), V(V) { - Options& options = *opt; - lower_y = options["lower_y"].doc("Boundary on lower y?").withDefault(lower_y); - upper_y = options["upper_y"].doc("Boundary on upper y?").withDefault(upper_y); - } - - void rhs() { - BoutReal totalFlux = 0; - if (lower_y) { - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - // Calculate flux through surface [normalised m^-2 s^-1], - // should be positive since V < 0.0 - BoutReal flux = - -0.5 * (N(r.ind, mesh->ystart, jz) + N(r.ind, mesh->ystart - 1, jz)) * 0.5 - * (V(r.ind, mesh->ystart, jz) + V(r.ind, mesh->ystart - 1, jz)); - totalFlux += flux; - } - } - } - if (upper_y) { - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - // Calculate flux through surface [normalised m^-2 s^-1], - // should be positive since V < 0.0 - BoutReal flux = -0.5 * (N(r.ind, mesh->yend, jz) + N(r.ind, mesh->yend + 1, jz)) - * 0.5 - * (V(r.ind, mesh->yend, jz) + V(r.ind, mesh->yend + 1, jz)); - totalFlux += flux; - } - } - } - } - -private: - bool lower_y{true}; - bool upper_y{true}; - const Field3D& N; - const Field3D& V; -} - -class yboundary_example { -public: - yboundary_example(Options* opt, const Field3D& N, const Field3D& V) : N(N), V(V) { - // Set what kind of yboundaries you want to include - yboundary.init(opt); - } - - void rhs() { - BoutReal totalFlux = 0; - yboundary.iter_pnts([&](auto& pnt) { - BoutReal flux = pnt.interpolate_sheath_o1(N) * pnt.interpolate_sheath_o1(V); - }); - } - -private: - YBoundary ybounday; - const Field3D& N; - const Field3D& V; -}; diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index 833e988eb9..6d12ab3541 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -2,18 +2,14 @@ #include "./boundary_iterator.hxx" #include "bout/parallel_boundary_region.hxx" -/*! - * This class allows to simplify iterating over y-boundaries. - * - * It makes it easier to write code for FieldAligned boundaries, but if a bit - * care is taken the code also works with FluxCoordinateIndependent code. - * - * An example how to replace old code is given here: - * - * \example example_yboundary_regions.hxx - * This is an example how to use the YBoundary class to replace RangeIterator - * boundaries. - */ + +/// This class allows to simplify iterating over y-boundaries. +/// +/// It makes it easier to write code for FieldAligned boundaries, but if a bit +/// care is taken the code also works with FluxCoordinateIndependent code. +/// +/// An example how to replace old code is given here: +/// ../../manual/sphinx/user_docs/boundary_options.rst class YBoundary { public: diff --git a/manual/sphinx/user_docs/boundary_options.rst b/manual/sphinx/user_docs/boundary_options.rst index 8493049516..579411b669 100644 --- a/manual/sphinx/user_docs/boundary_options.rst +++ b/manual/sphinx/user_docs/boundary_options.rst @@ -436,6 +436,85 @@ the upper Y boundary of a 2D variable ``var``:: The `BoundaryRegion` class is defined in ``include/boundary_region.hxx`` +Y-Boundaries +------------ + +The sheath boundaries are often implemented in the physics model. +Previously of they where implemented using a `RangeIterator`:: + + class yboundary_example_legacy { + public: + yboundary_example_legacy(Options* opt, const Field3D& N, const Field3D& V) + : N(N), V(V) { + Options& options = *opt; + lower_y = options["lower_y"].doc("Boundary on lower y?").withDefault(lower_y); + upper_y = options["upper_y"].doc("Boundary on upper y?").withDefault(upper_y); + } + + void rhs() { + BoutReal totalFlux = 0; + if (lower_y) { + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Calculate flux through surface [normalised m^-2 s^-1], + // should be positive since V < 0.0 + BoutReal flux = + -0.5 * (N(r.ind, mesh->ystart, jz) + N(r.ind, mesh->ystart - 1, jz)) * 0.5 + * (V(r.ind, mesh->ystart, jz) + V(r.ind, mesh->ystart - 1, jz)); + totalFlux += flux; + } + } + } + if (upper_y) { + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Calculate flux through surface [normalised m^-2 s^-1], + // should be positive since V < 0.0 + BoutReal flux = -0.5 * (N(r.ind, mesh->yend, jz) + N(r.ind, mesh->yend + 1, jz)) + * 0.5 + * (V(r.ind, mesh->yend, jz) + V(r.ind, mesh->yend + 1, jz)); + totalFlux += flux; + } + } + } + } + + private: + bool lower_y{true}; + bool upper_y{true}; + const Field3D& N; + const Field3D& V; + } + + +This can be replaced using the `YBoundary` class, which not only simplifies the +code, but also allows to have the same code working with non-field-aligned +geometries, as flux coordinate independent (FCI) method:: + + #include + + class yboundary_example { + public: + yboundary_example(Options* opt, const Field3D& N, const Field3D& V) : N(N), V(V) { + // Set what kind of yboundaries you want to include + yboundary.init(opt); + } + + void rhs() { + BoutReal totalFlux = 0; + yboundary.iter_pnts([&](auto& pnt) { + BoutReal flux = pnt.interpolate_sheath_o1(N) * pnt.interpolate_sheath_o1(V); + }); + } + + private: + YBoundary ybounday; + const Field3D& N; + const Field3D& V; + }; + + + Boundary regions ---------------- From 4e21aac84aa1c084944045bfe941755660febcfe Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 4 Mar 2025 15:51:40 +0100 Subject: [PATCH 26/62] Fix iter_pnts --- include/bout/yboundary_regions.hxx | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index 6d12ab3541..1d434e2420 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -29,7 +29,7 @@ public: for (auto& pnt : region) { f(pnt); } - } + }); } template @@ -78,5 +78,3 @@ private: bool is_init{false}; }; - -extern YBoundary yboundary; From 3a850b11e3559607ba16d5cb10cce882177d19b5 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 4 Mar 2025 15:52:00 +0100 Subject: [PATCH 27/62] Add more documentation on YBoundary --- manual/sphinx/user_docs/boundary_options.rst | 45 ++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/manual/sphinx/user_docs/boundary_options.rst b/manual/sphinx/user_docs/boundary_options.rst index 579411b669..18535c91b4 100644 --- a/manual/sphinx/user_docs/boundary_options.rst +++ b/manual/sphinx/user_docs/boundary_options.rst @@ -515,6 +515,51 @@ geometries, as flux coordinate independent (FCI) method:: +There are several member functions of ``pnt``. ``pnt`` is of type +`BoundaryRegionParIterBase` and `BoundaryRegionIter`, and both should provide +the same interface. If they don't that is a bug, as the above code is a +template, that gets instantiated for both types, and thus requires both +classes to provide the same interface, one for FCI-like boundaries and one for +field aligned boundaries. + +Here is a short summary of some members of ``pnt``, where ``f`` is a : + +.. list-table:: Members for boundary operation + :widths: 15 70 + :header-rows: 1 + + * - Function + - Description + * - ``pnt.ythis(f)`` + - Returns the value at the last point in the domain + * - ``pnt.ynext(f)`` + - Returns the value at the first point in the domain + * - ``pnt.yprev(f)`` + - Returns the value at the second to last point in the domain, if it is + valid. NB: this point may not be valid. + * - ``pnt.interpolate_sheath_o1(f)`` + - Returns the value at the boundary, assuming the bounday value has been set + * - ``pnt.extrapolate_sheath_o1(f)`` + - Returns the value at the boundary, extrapolating from the bulk, first order + * - ``pnt.extrapolate_sheath_o2(f)`` + - Returns the value at the boundary, extrapolating from the bulk, second order + * - ``pnt.extrapolate_next_o{1,2}(f)`` + - Extrapolate into the boundary from the bulk, first or second order + * - ``pnt.extrapolate_grad_o{1,2}(f)`` + - Extrapolate the gradient into the boundary, first or second order + * - ``pnt.dirichlet_o{1,2,3}(f, v)`` + - Apply dirichlet boundary conditions with value ``v`` and given order + * - ``pnt.neumann_o{1,2,3}(f, v)`` + - Applies a gradient of ``v / dy`` boundary condition. + * - ``pnt.limitFree(f)`` + - Extrapolate into the boundary using only monotonic decreasing values. + ``f`` needs to be positive. + * - ``pnt.dir`` + - The direction of the boundary. + + + + Boundary regions ---------------- From e0f80cf51473648753148b0efbf7e1a1e9edbf37 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 7 Mar 2025 10:39:28 +0100 Subject: [PATCH 28/62] Ensure the field has parallel slices --- include/bout/parallel_boundary_region.hxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 1ae49360f2..27d8c34ccf 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -259,6 +259,7 @@ public: template BoutReal& getAt(Field3D& f, int off) const { + ASSERT4(f.hasParallelSlices()); if constexpr (check) { ASSERT3(valid() > -off - 2); } @@ -267,6 +268,7 @@ public: } template const BoutReal& getAt(const Field3D& f, int off) const { + ASSERT4(f.hasParallelSlices()); if constexpr (check) { ASSERT3(valid() > -off - 2); } From 101f33df885c535c5dc21e54cf64967ea767afc8 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 7 Mar 2025 11:30:36 +0100 Subject: [PATCH 29/62] Lower check level --- include/bout/parallel_boundary_region.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 27d8c34ccf..862670b529 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -259,7 +259,7 @@ public: template BoutReal& getAt(Field3D& f, int off) const { - ASSERT4(f.hasParallelSlices()); + ASSERT3(f.hasParallelSlices()); if constexpr (check) { ASSERT3(valid() > -off - 2); } @@ -268,7 +268,7 @@ public: } template const BoutReal& getAt(const Field3D& f, int off) const { - ASSERT4(f.hasParallelSlices()); + ASSERT3(f.hasParallelSlices()); if constexpr (check) { ASSERT3(valid() > -off - 2); } From 4ea1657fea06c13f37e7ac8714182d566fde0308 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 21 Aug 2025 10:10:16 +0200 Subject: [PATCH 30/62] Prefer std::vector::emplace_back to work around cuda compiler issue --- include/bout/parallel_boundary_region.hxx | 66 +++++++++++++---------- 1 file changed, 37 insertions(+), 29 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 862670b529..3ddb8ed8fe 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -15,6 +15,38 @@ * */ +namespace bout { +namespace parallel_boundary_region { + +struct RealPoint { + BoutReal s_x; + BoutReal s_y; + BoutReal s_z; +}; + +struct Indices { + // Indices of the boundary point + Ind3D index; + // Intersection with boundary in index space + RealPoint intersection; + // Distance to intersection + BoutReal length; + // Angle between field line and boundary + // BoutReal angle; + // How many points we can go in the opposite direction + signed char valid; + signed char offset; + unsigned char abs_offset; + Indices(Ind3D index, RealPoint&& intersection, BoutReal length, signed char valid, + signed char offset, unsigned char abs_offset) + : index(index), intersection(std::move(intersection)), length(length), valid(valid), + offset(offset), abs_offset(abs_offset) {}; +}; + +using IndicesVec = std::vector; +using IndicesIter = IndicesVec::iterator; +using IndicesIterConst = IndicesVec::const_iterator; + inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc) { if (fm < fc) { return 1; // Neumann rather than increasing into boundary @@ -33,27 +65,6 @@ inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc) { class BoundaryRegionPar : public BoundaryRegionBase { - struct RealPoint { - BoutReal s_x; - BoutReal s_y; - BoutReal s_z; - }; - - struct Indices { - // Indices of the boundary point - Ind3D index; - // Intersection with boundary in index space - RealPoint intersection; - // Distance to intersection - BoutReal length; - // Angle between field line and boundary - // BoutReal angle; - // How many points we can go in the opposite direction - signed char valid; - signed char offset; - unsigned char abs_offset; - }; - using IndicesVec = std::vector; using IndicesIter = IndicesVec::iterator; @@ -77,11 +88,11 @@ public: /// Add a point to the boundary void add_point(Ind3D ind, BoutReal x, BoutReal y, BoutReal z, BoutReal length, signed char valid) { - bndry_points.push_back({ind, {x, y, z}, length, valid}); + bndry_points.emplace_back({ind, {x, y, z}, length, valid}); } void add_point(int ix, int iy, int iz, BoutReal x, BoutReal y, BoutReal z, BoutReal length, signed char valid) { - bndry_points.push_back({xyz2ind(ix, iy, iz, localmesh), {x, y, z}, length, valid}); + bndry_points.emplace_back({xyz2ind(ix, iy, iz, localmesh), {x, y, z}, length, valid}); } // final, so they can be inlined @@ -361,12 +372,9 @@ public: if (!bndry_points.empty() && bndry_points.back().index > ind) { is_sorted = false; } - bndry_points.push_back({ind, - {x, y, z}, - length, - valid, - offset, - static_cast(std::abs(offset))}); + bndry_points.emplace_back(ind, bout::parallel_boundary_region::RealPoint{x, y, z}, + length, valid, offset, + static_cast(std::abs(offset))); } void add_point(int ix, int iy, int iz, BoutReal x, BoutReal y, BoutReal z, BoutReal length, char valid, signed char offset) { From c78dd6ec3fd7fc6cec5be9f96cfa18c1792998a5 Mon Sep 17 00:00:00 2001 From: dschwoerer <5637662+dschwoerer@users.noreply.github.com> Date: Thu, 21 Aug 2025 08:16:32 +0000 Subject: [PATCH 31/62] Apply clang-format changes From 6fadc066e5aacc022f2e0a0961e4c87fd8e74e75 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 21 Aug 2025 11:37:58 +0200 Subject: [PATCH 32/62] remove unneeded std::move --- include/bout/parallel_boundary_region.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 3ddb8ed8fe..76a5cbe5d1 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -39,7 +39,7 @@ struct Indices { unsigned char abs_offset; Indices(Ind3D index, RealPoint&& intersection, BoutReal length, signed char valid, signed char offset, unsigned char abs_offset) - : index(index), intersection(std::move(intersection)), length(length), valid(valid), + : index(index), intersection(intersection), length(length), valid(valid), offset(offset), abs_offset(abs_offset) {}; }; From 212982033706fd6b5a6e5c182942de64e4cbbd5b Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 1 Sep 2025 12:21:49 +0200 Subject: [PATCH 33/62] Add limit_at_least function to boundary_iterator --- include/bout/boundary_iterator.hxx | 6 ++++++ include/bout/parallel_boundary_region.hxx | 8 ++++++++ 2 files changed, 14 insertions(+) diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx index 601f6a8a28..118036256e 100644 --- a/include/bout/boundary_iterator.hxx +++ b/include/bout/boundary_iterator.hxx @@ -77,6 +77,12 @@ public: } } + void limit_at_least(Field3D& f, BoutReal value) const { + if (ynext(f) < value) { + ynext(f) = value; + } + } + 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)]; } diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 76a5cbe5d1..6ce078435c 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -223,6 +223,14 @@ public: } } + void limit_at_least(Field3D& f, BoutReal value) const { + ITER() { + if (getAt(f, i) < value) { + getAt(f, i) = value; + } + } + } + // NB: value needs to be scaled by dy // neumann_o1 is actually o2 if we would use an appropriate one-sided stencil. // But in general we do not, and thus for normal C2 stencils, this is 1st order. From 7fc08a34de5cfb65ef8849328297dcbb5eaba96d Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 21 Oct 2025 11:11:33 +0200 Subject: [PATCH 34/62] Interpolate sheath is at least o2 --- include/bout/parallel_boundary_region.hxx | 6 ++++- src/field/field3d.cxx | 27 ++++++++++++++++++++++- 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 6ce078435c..f74bf35bd8 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -142,7 +142,11 @@ public: return ythis(f) * (1 + length()) - yprev(f) * length(); } - inline BoutReal interpolate_sheath_o1(const Field3D& f) const { + inline BoutReal interpolate_sheath_o2(const Field3D& f) const { + return ythis(f) * (1 - length()) + ynext(f) * length(); + } + inline BoutReal + interpolate_sheath_o2(const std::function& f) const { return ythis(f) * (1 - length()) + ynext(f) * length(); } diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index fe9fdc9263..146b36da42 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -445,7 +445,32 @@ void Field3D::setBoundaryTo(const Field3D& f3d) { allocate(); // Make sure data allocated - /// Loop over boundary regions + if (isFci()) { + ASSERT1(f3d.hasParallelSlices()); + if (copyParallelSlices) { + splitParallelSlices(); + for (int i = 0; i < fieldmesh->ystart; ++i) { + yup(i) = f3d.yup(i); + ydown(i) = f3d.ydown(i); + } + } else { + // Set yup/ydown using midpoint values from f3d + ASSERT1(hasParallelSlices()); + + for (auto& region : fieldmesh->getBoundariesPar()) { + for (const auto& pnt : *region) { + // Interpolate midpoint value in f3d + const BoutReal val = pnt.interpolate_sheath_o2(f3d); + // Set the same boundary value in this field + pnt.dirichlet_o1(*this, val); + } + } + } + } + + // Non-FCI. + // Transform to field-aligned coordinates? + // Loop over boundary regions for (const auto& reg : fieldmesh->getBoundaries()) { /// Loop within each region for (reg->first(); !reg->isDone(); reg->next()) { From e135545f2f8d9fcbe997e21fc4892216c98f046e Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 21 Oct 2025 14:00:31 +0200 Subject: [PATCH 35/62] Interpolate sheath is at least o2 --- include/bout/boundary_iterator.hxx | 2 +- manual/sphinx/user_docs/boundary_options.rst | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx index 118036256e..3504e2f967 100644 --- a/include/bout/boundary_iterator.hxx +++ b/include/bout/boundary_iterator.hxx @@ -43,7 +43,7 @@ public: return 2 * f(0, ind()) - f(0, ind().yp(-by).xp(-bx)); } - BoutReal interpolate_sheath_o1(const Field3D& f) const { + BoutReal interpolate_sheath_o2(const Field3D& f) const { return (f[ind()] + ynext(f)) * 0.5; } BoutReal diff --git a/manual/sphinx/user_docs/boundary_options.rst b/manual/sphinx/user_docs/boundary_options.rst index 18535c91b4..9f9df6736c 100644 --- a/manual/sphinx/user_docs/boundary_options.rst +++ b/manual/sphinx/user_docs/boundary_options.rst @@ -503,7 +503,7 @@ geometries, as flux coordinate independent (FCI) method:: void rhs() { BoutReal totalFlux = 0; yboundary.iter_pnts([&](auto& pnt) { - BoutReal flux = pnt.interpolate_sheath_o1(N) * pnt.interpolate_sheath_o1(V); + BoutReal flux = pnt.interpolate_sheath_o2(N) * pnt.interpolate_sheath_o2(V); }); } @@ -537,7 +537,7 @@ Here is a short summary of some members of ``pnt``, where ``f`` is a : * - ``pnt.yprev(f)`` - Returns the value at the second to last point in the domain, if it is valid. NB: this point may not be valid. - * - ``pnt.interpolate_sheath_o1(f)`` + * - ``pnt.interpolate_sheath_o2(f)`` - Returns the value at the boundary, assuming the bounday value has been set * - ``pnt.extrapolate_sheath_o1(f)`` - Returns the value at the boundary, extrapolating from the bulk, first order From 73e418302ec46c74fad42da0dca1f0807e747126 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 21 Oct 2025 16:16:39 +0200 Subject: [PATCH 36/62] Add limitFreeScale and SheathLimitMode --- include/bout/parallel_boundary_region.hxx | 67 ++++++++++++++--------- 1 file changed, 42 insertions(+), 25 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index f74bf35bd8..370175dd7e 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -15,6 +15,8 @@ * */ +BOUT_ENUM_CLASS(SheathLimitMode, limit_free, exponential_free, linear_free); + namespace bout { namespace parallel_boundary_region { @@ -47,6 +49,42 @@ using IndicesVec = std::vector; using IndicesIter = IndicesVec::iterator; using IndicesIterConst = IndicesVec::const_iterator; +/// Limited free gradient of log of a quantity +/// This ensures that the guard cell values remain positive +/// while also ensuring that the quantity never increases +/// +/// fm fc | fp +/// ^ boundary +/// +/// exp( 2*log(fc) - log(fm) ) +inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc, SheathLimitMode mode) { + if ((fm < fc) && (mode == SheathLimitMode::limit_free)) { + return fc; // Neumann rather than increasing into boundary + } + if (fm < 1e-10) { + return fc; // Low / no density condition + } + + BoutReal fp = 0; + switch (mode) { + case SheathLimitMode::limit_free: + case SheathLimitMode::exponential_free: + fp = SQ(fc) / fm; // Exponential + break; + case SheathLimitMode::linear_free: + fp = 2.0 * fc - fm; // Linear + break; + } + +#if CHECKLEVEL >= 2 + if (!std::isfinite(fp)) { + throw BoutException("SheathBoundary limitFree: {}, {} -> {}", fm, fc, fp); + } +#endif + + return fp; +} + inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc) { if (fm < fc) { return 1; // Neumann rather than increasing into boundary @@ -74,31 +112,10 @@ class BoundaryRegionPar : public BoundaryRegionBase { IndicesIter bndry_position; public: - BoundaryRegionPar(const std::string& name, int dir, Mesh* passmesh) - : BoundaryRegionBase(name, passmesh), dir(dir) { - ASSERT0(std::abs(dir) == 1); - BoundaryRegionBase::isParallel = true; - } - BoundaryRegionPar(const std::string& name, BndryLoc loc, int dir, Mesh* passmesh) - : BoundaryRegionBase(name, loc, passmesh), dir(dir) { - BoundaryRegionBase::isParallel = true; - ASSERT0(std::abs(dir) == 1); - } - - /// Add a point to the boundary - void add_point(Ind3D ind, BoutReal x, BoutReal y, BoutReal z, BoutReal length, - signed char valid) { - bndry_points.emplace_back({ind, {x, y, z}, length, valid}); - } - void add_point(int ix, int iy, int iz, BoutReal x, BoutReal y, BoutReal z, - BoutReal length, signed char valid) { - bndry_points.emplace_back({xyz2ind(ix, iy, iz, localmesh), {x, y, z}, length, valid}); - } - - // final, so they can be inlined - void first() final { bndry_position = begin(bndry_points); } - void next() final { ++bndry_position; } - bool isDone() final { return (bndry_position == end(bndry_points)); } + BoundaryRegionParIterBase(IndicesVec& bndry_points, IndicesIter bndry_position, int dir, + Mesh* localmesh) + : bndry_points(bndry_points), bndry_position(bndry_position), dir(dir), + localmesh(localmesh) {}; // getter Ind3D ind() const { return bndry_position->index; } From 4c105b833f7282c43184f31422207c8128b4aad7 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 21 Oct 2025 16:17:23 +0200 Subject: [PATCH 37/62] Add extrapolate_sheath_free --- include/bout/boundary_iterator.hxx | 8 ++++++++ include/bout/parallel_boundary_region.hxx | 8 ++++++++ 2 files changed, 16 insertions(+) diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx index 3504e2f967..2f05f0c82e 100644 --- a/include/bout/boundary_iterator.hxx +++ b/include/bout/boundary_iterator.hxx @@ -51,6 +51,14 @@ public: return 0.5 * (3 * f(0, ind()) - f(0, ind().yp(-by).xp(-bx))); } + BoutReal extrapolate_sheath_free(const Field3D& f, SheathLimitMode mode) const { + const BoutReal fac = + bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f), mode); + BoutReal val = ythis(f); + BoutReal next = mode == SheathLimitMode::linear_free ? val + fac : val * fac; + return 0.5 * (val + next); + } + void limitFree(Field3D& f) const { const BoutReal fac = bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f)); diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 370175dd7e..8b10eccca6 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -291,6 +291,14 @@ public: } } + BoutReal extrapolate_sheath_free(const Field3D& f, SheathLimitMode mode) const { + const auto fac = valid() > 0 ? limitFreeScale(yprev(f), ythis(f), mode) + : (mode == SheathLimitMode::linear_free ? 0 : 1); + auto val = ythis(f); + BoutReal next = mode == SheathLimitMode::linear_free ? val + fac : val * fac; + return val * length() + next * (1 - length()); + } + void setAll(Field3D& f, const BoutReal val) const { for (int i = -localmesh->ystart; i <= localmesh->ystart; ++i) { f.ynext(i)[ind().yp(i)] = val; From 41f9b2844d3de90b924486d0946c5e89e56ff8f1 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 21 Oct 2025 16:17:34 +0200 Subject: [PATCH 38/62] Add set_free --- include/bout/boundary_iterator.hxx | 17 +++++++++++++++++ include/bout/parallel_boundary_region.hxx | 17 +++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx index 2f05f0c82e..7c162ab0bd 100644 --- a/include/bout/boundary_iterator.hxx +++ b/include/bout/boundary_iterator.hxx @@ -59,6 +59,23 @@ public: return 0.5 * (val + next); } + void set_free(Field3D& f, SheathLimitMode mode) const { + const BoutReal fac = + bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f), mode); + BoutReal val = ythis(f); + if (mode == SheathLimitMode::linear_free) { + for (int i = 1; i <= localmesh->ystart; ++i) { + val += fac; + f[ind().yp(by * i).xp(bx * i)] = val; + } + } else { + for (int i = 1; i <= localmesh->ystart; ++i) { + val *= fac; + f[ind().yp(by * i).xp(bx * i)] = val; + } + } + } + void limitFree(Field3D& f) const { const BoutReal fac = bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f)); diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 8b10eccca6..ccefc277a4 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -299,6 +299,23 @@ public: return val * length() + next * (1 - length()); } + void set_free(Field3D& f, SheathLimitMode mode) const { + const auto fac = valid() > 0 ? limitFreeScale(yprev(f), ythis(f), mode) + : (mode == SheathLimitMode::linear_free ? 0 : 1); + auto val = ythis(f); + if (mode == SheathLimitMode::linear_free) { + ITER() { + val += fac; + getAt(f, i) = val; + } + } else { + ITER() { + val *= fac; + getAt(f, i) = val; + } + } + } + void setAll(Field3D& f, const BoutReal val) const { for (int i = -localmesh->ystart; i <= localmesh->ystart; ++i) { f.ynext(i)[ind().yp(i)] = val; From e1364b3c81a8d642ed0caedaec3d3ef00298ddbd Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 22 Oct 2025 11:11:34 +0200 Subject: [PATCH 39/62] add missing interpolate_sheath_o2(func) to FA --- include/bout/boundary_iterator.hxx | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx index 7c162ab0bd..9e2ff04588 100644 --- a/include/bout/boundary_iterator.hxx +++ b/include/bout/boundary_iterator.hxx @@ -46,6 +46,12 @@ public: BoutReal interpolate_sheath_o2(const Field3D& f) const { return (f[ind()] + ynext(f)) * 0.5; } + + BoutReal + interpolate_sheath_o2(const std::function& f) const { + return (f(0, ind()) + f(0, ind().yp(-by).xp(-bx))) * 0.5; + } + BoutReal extrapolate_sheath_o2(const std::function& f) const { return 0.5 * (3 * f(0, ind()) - f(0, ind().yp(-by).xp(-bx))); From 2edebb3a19d29b40d8e985557013a92957d1a917 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 22 Oct 2025 11:12:06 +0200 Subject: [PATCH 40/62] add is_lower() to check for which direction the boundary is --- include/bout/boundary_iterator.hxx | 5 +++++ include/bout/parallel_boundary_region.hxx | 2 ++ 2 files changed, 7 insertions(+) diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx index 9e2ff04588..728c2e1cb0 100644 --- a/include/bout/boundary_iterator.hxx +++ b/include/bout/boundary_iterator.hxx @@ -92,6 +92,11 @@ public: } } + bool is_lower() const { + ASSERT2(bx == 0); + return by == -1; + } + void neumann_o1(Field3D& f, BoutReal grad) const { BoutReal val = ythis(f); for (int i = 1; i <= localmesh->ystart; ++i) { diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index ccefc277a4..c93c5b680d 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -252,6 +252,8 @@ public: } } + bool is_lower() const { return dir == -1; } + // NB: value needs to be scaled by dy // neumann_o1 is actually o2 if we would use an appropriate one-sided stencil. // But in general we do not, and thus for normal C2 stencils, this is 1st order. From 382120525f105cf3595bea311a519cb99a049b28 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 21 Jan 2026 13:04:44 +0100 Subject: [PATCH 41/62] Only call iterateBndryLowerY() for non-FCI --- src/mesh/impls/bout/boutmesh.cxx | 33 ++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 99ec63ff25..2f88e86cb4 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -620,23 +620,28 @@ int BoutMesh::load() { // Add boundary regions addBoundaryRegions(); - // Set cached values - { - int mybndry = static_cast(!(iterateBndryLowerY().isDone())); - int allbndry = 0; - mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(yend)); - has_boundary_lower_y = static_cast(allbndry); - } - { - int mybndry = static_cast(!(iterateBndryUpperY().isDone())); - int allbndry = 0; - mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(ystart)); - has_boundary_upper_y = static_cast(allbndry); - } - // Initialize default coordinates getCoordinates(); + // Set cached values + if (isFci()) { + has_boundary_lower_y = false; + has_boundary_upper_y = false; + } else { + { + int mybndry = static_cast(!(iterateBndryLowerY().isDone())); + int allbndry = 0; + mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(yend)); + has_boundary_lower_y = static_cast(allbndry); + } + { + int mybndry = static_cast(!(iterateBndryUpperY().isDone())); + int allbndry = 0; + mpi->MPI_Allreduce(&mybndry, &allbndry, 1, MPI_INT, MPI_BOR, getXcomm(ystart)); + has_boundary_upper_y = static_cast(allbndry); + } + } + output_info.write(_("\tdone\n")); return 0; From 38293b9467c74c59911c9baf8574c409cd3dfe83 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 4 Mar 2026 10:57:15 +0100 Subject: [PATCH 42/62] Break include cycle --- include/bout/boundary_region.hxx | 1 - 1 file changed, 1 deletion(-) diff --git a/include/bout/boundary_region.hxx b/include/bout/boundary_region.hxx index 22956d1d4a..b73c2e5a02 100644 --- a/include/bout/boundary_region.hxx +++ b/include/bout/boundary_region.hxx @@ -4,7 +4,6 @@ class BoundaryRegion; #ifndef BOUT_BNDRY_REGION_H #define BOUT_BNDRY_REGION_H -#include "bout/mesh.hxx" #include "bout/region.hxx" #include "bout/sys/parallel_stencils.hxx" #include From 6ae2351403dfb782dabec87f34446ebe45deb3ef Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 10 Mar 2026 08:12:08 +0100 Subject: [PATCH 43/62] Fix Field3D::setBoundaryTo for FCI methods Without this fix, boundary conditions set on yup/down fields are not applied when a boundary is copied from one field to another. --- src/field/field3d.cxx | 27 ++++++++++----------------- 1 file changed, 10 insertions(+), 17 deletions(-) diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 146b36da42..419fb05a2b 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -446,26 +446,19 @@ void Field3D::setBoundaryTo(const Field3D& f3d) { allocate(); // Make sure data allocated if (isFci()) { + // Set yup/ydown using midpoint values from f3d ASSERT1(f3d.hasParallelSlices()); - if (copyParallelSlices) { - splitParallelSlices(); - for (int i = 0; i < fieldmesh->ystart; ++i) { - yup(i) = f3d.yup(i); - ydown(i) = f3d.ydown(i); - } - } else { - // Set yup/ydown using midpoint values from f3d - ASSERT1(hasParallelSlices()); - - for (auto& region : fieldmesh->getBoundariesPar()) { - for (const auto& pnt : *region) { - // Interpolate midpoint value in f3d - const BoutReal val = pnt.interpolate_sheath_o2(f3d); - // Set the same boundary value in this field - pnt.dirichlet_o1(*this, val); - } + ASSERT1(hasParallelSlices()); + + for (auto& region : fieldmesh->getBoundariesPar()) { + for (const auto& pnt : *region) { + // Interpolate midpoint value in f3d + const BoutReal val = pnt.interpolate_sheath_o2(f3d); + // Set the same boundary value in this field + pnt.dirichlet_o2(*this, val); } } + return; } // Non-FCI. From 576c609d779ddaa6733d24bb1558d9c3e28ad897 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 25 Nov 2024 16:44:29 +0100 Subject: [PATCH 44/62] Copy BCs in x-direction also for FCI --- src/field/field3d.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 419fb05a2b..9c59d32bd0 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -458,13 +458,15 @@ void Field3D::setBoundaryTo(const Field3D& f3d) { pnt.dirichlet_o2(*this, val); } } - return; } // Non-FCI. // Transform to field-aligned coordinates? // Loop over boundary regions for (const auto& reg : fieldmesh->getBoundaries()) { + if (isFci() && reg->by != 0) { + continue; + } /// Loop within each region for (reg->first(); !reg->isDone(); reg->next()) { for (int z = 0; z < nz; z++) { From 1b99b34fb771022d7512c44fe4dfd533b809612e Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 25 Nov 2024 16:48:06 +0100 Subject: [PATCH 45/62] Use consistently first order interpolation --- src/field/field3d.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 9c59d32bd0..9bda64dd32 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -455,7 +455,7 @@ void Field3D::setBoundaryTo(const Field3D& f3d) { // Interpolate midpoint value in f3d const BoutReal val = pnt.interpolate_sheath_o2(f3d); // Set the same boundary value in this field - pnt.dirichlet_o2(*this, val); + pnt.dirichlet_o1(*this, val); } } } From ff190e8f19fa8fa99544f08281abaa0fade5b2a9 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 10 Mar 2026 08:12:54 +0100 Subject: [PATCH 46/62] Add option to copy parallel slices at setBoundaryTo --- include/bout/field3d.hxx | 3 ++- src/field/field3d.cxx | 28 ++++++++++++++++++---------- 2 files changed, 20 insertions(+), 11 deletions(-) diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index c6e1d1191d..57f8c559bd 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -496,7 +496,8 @@ public: /// This uses 2nd order central differences to set the value /// on the boundary to the value on the boundary in field \p f3d. /// Note: does not just copy values in boundary region. - void setBoundaryTo(const Field3D& f3d); + void setBoundaryTo(const Field3D& f3d) { setBoundaryTo(f3d, true); } + void setBoundaryTo(const Field3D& f3d, bool copyParallelSlices); using FieldData::applyParallelBoundary; void applyParallelBoundary() override; diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 9bda64dd32..6171cc1061 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -439,23 +439,31 @@ void Field3D::applyTDerivBoundary() { } } -void Field3D::setBoundaryTo(const Field3D& f3d) { +void Field3D::setBoundaryTo(const Field3D& f3d, bool copyParallelSlices) { checkData(f3d); allocate(); // Make sure data allocated if (isFci()) { - // Set yup/ydown using midpoint values from f3d ASSERT1(f3d.hasParallelSlices()); - ASSERT1(hasParallelSlices()); - - for (auto& region : fieldmesh->getBoundariesPar()) { - for (const auto& pnt : *region) { - // Interpolate midpoint value in f3d - const BoutReal val = pnt.interpolate_sheath_o2(f3d); - // Set the same boundary value in this field - pnt.dirichlet_o1(*this, val); + if (copyParallelSlices) { + splitParallelSlices(); + for (int i = 0; i < fieldmesh->ystart; ++i) { + yup(i) = f3d.yup(i); + ydown(i) = f3d.ydown(i); + } + } else { + // Set yup/ydown using midpoint values from f3d + ASSERT1(hasParallelSlices()); + + for (auto& region : fieldmesh->getBoundariesPar()) { + for (const auto& pnt : *region) { + // Interpolate midpoint value in f3d + const BoutReal val = pnt.interpolate_sheath_o2(f3d); + // Set the same boundary value in this field + pnt.dirichlet_o1(*this, val); + } } } } From 6aac8345434cb9301b982e3fa573a59e8738b5ac Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 4 Mar 2026 11:43:22 +0100 Subject: [PATCH 47/62] Prefer [[maybe_unused]] --- include/bout/sys/parallel_stencils.hxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/bout/sys/parallel_stencils.hxx b/include/bout/sys/parallel_stencils.hxx index 34a51c5285..1d55931ff3 100644 --- a/include/bout/sys/parallel_stencils.hxx +++ b/include/bout/sys/parallel_stencils.hxx @@ -11,15 +11,15 @@ inline BoutReal pow(BoutReal val, int exp) { ASSERT3(exp == 3); return val * val * val; } -inline BoutReal dirichlet_o1(BoutReal UNUSED(spacing0), BoutReal value0) { +inline BoutReal dirichlet_o1([[maybe_unused]] BoutReal spacing0, BoutReal value0) { return value0; } inline BoutReal dirichlet_o2(BoutReal spacing0, BoutReal value0, BoutReal spacing1, BoutReal value1) { return (spacing0 * value1 - spacing1 * value0) / (spacing0 - spacing1); } -inline BoutReal neumann_o2(BoutReal UNUSED(spacing0), BoutReal value0, BoutReal spacing1, - BoutReal value1) { +inline BoutReal neumann_o2([[maybe_unused]] BoutReal spacing0, BoutReal value0, + BoutReal spacing1, BoutReal value1) { return -spacing1 * value0 + value1; } inline BoutReal dirichlet_o3(BoutReal spacing0, BoutReal value0, BoutReal spacing1, From d5fb4c4f1899061fae855068042d28d1097659c9 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 4 Mar 2026 12:42:44 +0100 Subject: [PATCH 48/62] Fixup parallel_boundary --- include/bout/parallel_boundary_op.hxx | 41 +++++++++++++---------- include/bout/parallel_boundary_region.hxx | 23 +++---------- src/mesh/parallel/fci.cxx | 7 ++++ src/mesh/parallel_boundary_op.cxx | 9 +++-- 4 files changed, 40 insertions(+), 40 deletions(-) diff --git a/include/bout/parallel_boundary_op.hxx b/include/bout/parallel_boundary_op.hxx index d8620e892b..9e551ebc17 100644 --- a/include/bout/parallel_boundary_op.hxx +++ b/include/bout/parallel_boundary_op.hxx @@ -49,7 +49,7 @@ protected: enum class ValueType { GEN, FIELD, REAL }; const ValueType value_type{ValueType::REAL}; - BoutReal getValue(const BoundaryRegionPar& bndry, BoutReal t); + BoutReal getValue(const BoundaryRegionParIter& bndry, BoutReal t); }; template @@ -95,12 +95,13 @@ public: auto dy = f.getCoordinates()->dy; - for (bndry->first(); !bndry->isDone(); bndry->next()) { - BoutReal value = getValue(*bndry, t); + for (auto pnt : *bndry) { + //for (bndry->first(); !bndry->isDone(); bndry->next()) { + BoutReal value = getValue(pnt, t); if (isNeumann) { - value *= dy[bndry->ind()]; + value *= dy[pnt.ind()]; } - static_cast(this)->apply_stencil(f, bndry, value); + static_cast(this)->apply_stencil(f, pnt, value); } } }; @@ -111,24 +112,27 @@ public: class BoundaryOpPar_dirichlet_o1 : public BoundaryOpParTemp { public: using BoundaryOpParTemp::BoundaryOpParTemp; - static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) { - bndry->dirichlet_o1(f, value); + static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt, + BoutReal value) { + pnt.dirichlet_o1(f, value); } }; class BoundaryOpPar_dirichlet_o2 : public BoundaryOpParTemp { public: using BoundaryOpParTemp::BoundaryOpParTemp; - static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) { - bndry->dirichlet_o2(f, value); + static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt, + BoutReal value) { + pnt.dirichlet_o2(f, value); } }; class BoundaryOpPar_dirichlet_o3 : public BoundaryOpParTemp { public: using BoundaryOpParTemp::BoundaryOpParTemp; - static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) { - bndry->dirichlet_o3(f, value); + static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt, + BoutReal value) { + pnt.dirichlet_o3(f, value); } }; @@ -136,8 +140,9 @@ class BoundaryOpPar_neumann_o1 : public BoundaryOpParTemp { public: using BoundaryOpParTemp::BoundaryOpParTemp; - static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) { - bndry->neumann_o1(f, value); + static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt, + BoutReal value) { + pnt.neumann_o1(f, value); } }; @@ -145,8 +150,9 @@ class BoundaryOpPar_neumann_o2 : public BoundaryOpParTemp { public: using BoundaryOpParTemp::BoundaryOpParTemp; - static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) { - bndry->neumann_o2(f, value); + static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt, + BoutReal value) { + pnt.neumann_o2(f, value); } }; @@ -154,8 +160,9 @@ class BoundaryOpPar_neumann_o3 : public BoundaryOpParTemp { public: using BoundaryOpParTemp::BoundaryOpParTemp; - static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) { - bndry->neumann_o3(f, value); + static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt, + BoutReal value) { + pnt.neumann_o3(f, value); } }; diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index c93c5b680d..3caa98767e 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -3,6 +3,7 @@ #include "bout/boundary_region.hxx" #include "bout/bout_types.hxx" +#include #include #include "bout/sys/parallel_stencils.hxx" @@ -101,15 +102,8 @@ inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc) { return fp; } -class BoundaryRegionPar : public BoundaryRegionBase { - - using IndicesVec = std::vector; - using IndicesIter = IndicesVec::iterator; - - /// Vector of points in the boundary - IndicesVec bndry_points; - /// Current position in the boundary points - IndicesIter bndry_position; +template +class BoundaryRegionParIterBase { public: BoundaryRegionParIterBase(IndicesVec& bndry_points, IndicesIter bndry_position, int dir, @@ -130,22 +124,15 @@ public: // setter void setValid(signed char valid) { bndry_position->valid = valid; } - bool contains(const BoundaryRegionPar& bndry) const { - return std::binary_search( - begin(bndry_points), end(bndry_points), *bndry.bndry_position, - [](const Indices& i1, const Indices& i2) { return i1.index < i2.index; }); - } - // extrapolate a given point to the boundary BoutReal extrapolate_sheath_o1(const Field3D& f) const { return ythis(f); } BoutReal extrapolate_sheath_o2(const Field3D& f) const { ASSERT3(valid() >= 0); if (valid() < 1) { - return extrapolate_o1(f); + return extrapolate_sheath_o1(f); } return ythis(f) * (1 + length()) - yprev(f) * length(); } - inline BoutReal extrapolate_sheath_o1(const std::function& f) const { return ythis(f); @@ -184,7 +171,7 @@ public: extrapolate_next_o2(const std::function& f) const { ASSERT3(valid() >= 0); if (valid() < 1) { - return extrapolate_sheath_o1(f); + return extrapolate_next_o1(f); } return ythis(f) * 2 - yprev(f); } diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index d3fcc8d04c..28b8db9326 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -62,6 +62,13 @@ #include #include +namespace { +template +int sgn(T val) { + return (T(0) < val) - (val < T(0)); +} +} // namespace + using namespace std::string_view_literals; FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy, diff --git a/src/mesh/parallel_boundary_op.cxx b/src/mesh/parallel_boundary_op.cxx index ebd9852791..644331221e 100644 --- a/src/mesh/parallel_boundary_op.cxx +++ b/src/mesh/parallel_boundary_op.cxx @@ -1,21 +1,20 @@ #include "bout/parallel_boundary_op.hxx" +#include "bout/bout_types.hxx" #include "bout/constants.hxx" #include "bout/field_factory.hxx" #include "bout/globals.hxx" #include "bout/mesh.hxx" #include "bout/output.hxx" +#include "bout/parallel_boundary_region.hxx" -BoutReal BoundaryOpPar::getValue(const BoundaryRegionPar& bndry, BoutReal t) { - BoutReal value; - +BoutReal BoundaryOpPar::getValue(const BoundaryRegionParIter& bndry, BoutReal t) { switch (value_type) { case ValueType::GEN: return gen_values->generate(bout::generator::Context( bndry.s_x(), bndry.s_y(), bndry.s_z(), CELL_CENTRE, bndry.localmesh, t)); case ValueType::FIELD: // FIXME: Interpolate to s_x, s_y, s_z... - value = (*field_values)[bndry.ind()]; - return value; + return (*field_values)[bndry.ind()]; case ValueType::REAL: return real_value; default: From 4b503c5cf0af017935ef6188638b3d8c4e30b81f Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 6 Mar 2026 14:47:34 +0100 Subject: [PATCH 49/62] Ensure init is only called once --- include/bout/yboundary_regions.hxx | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index 1d434e2420..deef2257ba 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -38,6 +38,12 @@ public: } void init(Options& options, Mesh* mesh = nullptr) { + if (is_init) { + if (optptr == &options) { + return; + } + throw BoutException("YBoundary is already initialised!"); + } if (mesh == nullptr) { mesh = bout::globals::mesh; } @@ -70,6 +76,7 @@ public: } } is_init = true; + optptr = &options; } private: @@ -77,4 +84,5 @@ private: std::vector> boundary_regions; bool is_init{false}; + Options* optptr; }; From e9a0944c65b043db4dfbc2f8211d6f09847cccd9 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 9 Mar 2026 11:58:36 +0100 Subject: [PATCH 50/62] Add headers from clang-tidy --- include/bout/boundary_iterator.hxx | 8 ++++++++ include/bout/parallel_boundary_region.hxx | 6 ++++++ include/bout/yboundary_regions.hxx | 6 ++++++ 3 files changed, 20 insertions(+) diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx index 728c2e1cb0..e3ecf8234b 100644 --- a/include/bout/boundary_iterator.hxx +++ b/include/bout/boundary_iterator.hxx @@ -1,9 +1,17 @@ #pragma once +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/build_defines.hxx" +#include "bout/field2d.hxx" +#include "bout/field3d.hxx" #include "bout/mesh.hxx" #include "bout/parallel_boundary_region.hxx" #include "bout/sys/parallel_stencils.hxx" #include "bout/sys/range.hxx" +#include +#include +#include class BoundaryRegionIter { public: diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 3caa98767e..ca1437782f 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -3,9 +3,15 @@ #include "bout/boundary_region.hxx" #include "bout/bout_types.hxx" +#include +#include +#include #include #include +#include "bout/build_defines.hxx" +#include "bout/field2d.hxx" +#include "bout/region.hxx" #include "bout/sys/parallel_stencils.hxx" #include #include diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index deef2257ba..74a68377c9 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -1,7 +1,13 @@ #pragma once #include "./boundary_iterator.hxx" +#include "bout/assert.hxx" +#include "bout/boutexception.hxx" +#include "bout/globals.hxx" +#include "bout/options.hxx" #include "bout/parallel_boundary_region.hxx" +#include +#include /// This class allows to simplify iterating over y-boundaries. /// From c59626a0498b0cfd1be46f84f0fa38b4855a43cc Mon Sep 17 00:00:00 2001 From: dschwoerer <5637662+dschwoerer@users.noreply.github.com> Date: Tue, 10 Mar 2026 07:25:19 +0000 Subject: [PATCH 51/62] [bot] Apply format changes --- include/bout/parallel_boundary_region.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index ca1437782f..dbefeb60e6 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -49,7 +49,7 @@ struct Indices { Indices(Ind3D index, RealPoint&& intersection, BoutReal length, signed char valid, signed char offset, unsigned char abs_offset) : index(index), intersection(intersection), length(length), valid(valid), - offset(offset), abs_offset(abs_offset) {}; + offset(offset), abs_offset(abs_offset){}; }; using IndicesVec = std::vector; @@ -115,7 +115,7 @@ public: BoundaryRegionParIterBase(IndicesVec& bndry_points, IndicesIter bndry_position, int dir, Mesh* localmesh) : bndry_points(bndry_points), bndry_position(bndry_position), dir(dir), - localmesh(localmesh) {}; + localmesh(localmesh){}; // getter Ind3D ind() const { return bndry_position->index; } From 2c4179d9510854f45408d57ba5e9ba230e2e368d Mon Sep 17 00:00:00 2001 From: dschwoerer <5637662+dschwoerer@users.noreply.github.com> Date: Tue, 10 Mar 2026 07:30:38 +0000 Subject: [PATCH 52/62] [bot] Add last format changes commit to ignore file --- .git-blame-ignore-revs | 1 + 1 file changed, 1 insertion(+) diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index 84be8e384a..44192ad222 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -14,3 +14,4 @@ c8f38049359170a34c915e209276238ea66b9a1e 8d5cb39e03c2644715a50684f8cd0318b4e82676 ec69e8838be2dde140a915e50506f8e1ce3cb534 f2bc0488a298f136294c523bc5ab4086d090059b +c59626a0498b0cfd1be46f84f0fa38b4855a43cc From 717864631ccf1417c41acbfde9a347e6888764a6 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 10 Mar 2026 12:35:32 +0100 Subject: [PATCH 53/62] Apply clang-tidy fixes --- include/bout/boundary_iterator.hxx | 32 ++++++++++++----------- include/bout/boundary_region.hxx | 2 -- include/bout/parallel_boundary_region.hxx | 26 ++++++++++-------- include/bout/sys/parallel_stencils.hxx | 2 +- include/bout/yboundary_regions.hxx | 11 +++++--- src/mesh/parallel/fci.cxx | 6 ++--- 6 files changed, 43 insertions(+), 36 deletions(-) diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx index e3ecf8234b..7d8e7ed027 100644 --- a/include/bout/boundary_iterator.hxx +++ b/include/bout/boundary_iterator.hxx @@ -7,19 +7,20 @@ #include "bout/field3d.hxx" #include "bout/mesh.hxx" #include "bout/parallel_boundary_region.hxx" +#include "bout/region.hxx" #include "bout/sys/parallel_stencils.hxx" #include "bout/sys/range.hxx" -#include #include #include class BoundaryRegionIter { public: + virtual ~BoundaryRegionIter() = default; BoundaryRegionIter(int x, int y, int bx, int by, Mesh* mesh) : dir(bx + by), x(x), y(y), bx(bx), by(by), localmesh(mesh) { ASSERT3(bx * by == 0); } - bool operator!=(const BoundaryRegionIter& rhs) { return ind() != rhs.ind(); } + bool operator!=(const BoundaryRegionIter& rhs) const { return ind() != rhs.ind(); } Ind3D ind() const { return xyz2ind(x, y, z); } BoundaryRegionIter& operator++() { @@ -44,11 +45,13 @@ public: return (f[ind()] * 3 - yprev(f)) * 0.5; } - BoutReal extrapolate_next_o2(const Field3D& f) const { return 2 * f[ind()] - yprev(f); } + BoutReal extrapolate_next_o2(const Field3D& f) const { + return (2 * f[ind()]) - yprev(f); + } BoutReal extrapolate_next_o2(const std::function& f) const { - return 2 * f(0, ind()) - f(0, ind().yp(-by).xp(-bx)); + return (2 * f(0, ind())) - f(0, ind().yp(-by).xp(-bx)); } BoutReal interpolate_sheath_o2(const Field3D& f) const { @@ -68,8 +71,8 @@ public: BoutReal extrapolate_sheath_free(const Field3D& f, SheathLimitMode mode) const { const BoutReal fac = bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f), mode); - BoutReal val = ythis(f); - BoutReal next = mode == SheathLimitMode::linear_free ? val + fac : val * fac; + const BoutReal val = ythis(f); + const BoutReal next = mode == SheathLimitMode::linear_free ? val + fac : val * fac; return 0.5 * (val + next); } @@ -122,9 +125,7 @@ public: } void limit_at_least(Field3D& f, BoutReal value) const { - if (ynext(f) < value) { - ynext(f) = value; - } + ynext(f) = std::max(ynext(f), value); } BoutReal& ynext(Field3D& f) const { return f[ind().yp(by).xp(bx)]; } @@ -141,7 +142,7 @@ public: } } - int abs_offset() const { return 1; } + static int abs_offset() { return 1; } #if BOUT_USE_METRIC_3D == 0 BoutReal& ynext(Field2D& f) const { return f[ind().yp(by).xp(bx)]; } @@ -166,13 +167,14 @@ private: int nz() const { return localmesh->LocalNz; } Ind3D xyz2ind(int x, int y, int z) const { - return Ind3D{(x * ny() + y) * nz() + z, ny(), nz()}; + return Ind3D{((x * ny() + y) * nz()) + z, ny(), nz()}; } }; class BoundaryRegionIterY : public BoundaryRegionIter { public: - BoundaryRegionIterY(RangeIterator r, int y, int dir, bool is_end, Mesh* mesh) + virtual ~BoundaryRegionIterY() = default; + BoundaryRegionIterY(const RangeIterator& r, int y, int dir, bool is_end, Mesh* mesh) : BoundaryRegionIter(r.ind, y, 0, dir, mesh), r(r), is_end(is_end) {} bool operator!=(const BoundaryRegionIterY& rhs) { @@ -189,7 +191,7 @@ public: return x != rhs.x; } - virtual void _next() override { + void _next() override { ++r; x = r.ind; } @@ -201,8 +203,8 @@ private: class NewBoundaryRegionY { public: - NewBoundaryRegionY(Mesh* mesh, bool lower, RangeIterator r) - : mesh(mesh), lower(lower), r(std::move(r)) {} + NewBoundaryRegionY(Mesh* mesh, bool lower, const RangeIterator& r) + : mesh(mesh), lower(lower), r(r) {} BoundaryRegionIterY begin(bool begin = true) { return BoundaryRegionIterY(r, lower ? mesh->ystart : mesh->yend, lower ? -1 : +1, !begin, mesh); diff --git a/include/bout/boundary_region.hxx b/include/bout/boundary_region.hxx index b73c2e5a02..58de12045e 100644 --- a/include/bout/boundary_region.hxx +++ b/include/bout/boundary_region.hxx @@ -4,8 +4,6 @@ class BoundaryRegion; #ifndef BOUT_BNDRY_REGION_H #define BOUT_BNDRY_REGION_H -#include "bout/region.hxx" -#include "bout/sys/parallel_stencils.hxx" #include #include diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index dbefeb60e6..1d06bc65cb 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -1,18 +1,22 @@ #ifndef BOUT_PAR_BNDRY_H #define BOUT_PAR_BNDRY_H +#include "bout/assert.hxx" #include "bout/boundary_region.hxx" +#include "bout/bout_enum_class.hxx" #include "bout/bout_types.hxx" #include #include #include #include +#include #include #include "bout/build_defines.hxx" #include "bout/field2d.hxx" #include "bout/region.hxx" #include "bout/sys/parallel_stencils.hxx" +#include "bout/utils.hxx" #include #include @@ -79,7 +83,7 @@ inline BoutReal limitFreeScale(BoutReal fm, BoutReal fc, SheathLimitMode mode) { fp = SQ(fc) / fm; // Exponential break; case SheathLimitMode::linear_free: - fp = 2.0 * fc - fm; // Linear + fp = (2.0 * fc) - fm; // Linear break; } @@ -139,11 +143,11 @@ public: } return ythis(f) * (1 + length()) - yprev(f) * length(); } - inline BoutReal + BoutReal extrapolate_sheath_o1(const std::function& f) const { return ythis(f); } - inline BoutReal + BoutReal extrapolate_sheath_o2(const std::function& f) const { ASSERT3(valid() >= 0); if (valid() < 1) { @@ -152,16 +156,16 @@ public: return ythis(f) * (1 + length()) - yprev(f) * length(); } - inline BoutReal interpolate_sheath_o2(const Field3D& f) const { + BoutReal interpolate_sheath_o2(const Field3D& f) const { return ythis(f) * (1 - length()) + ynext(f) * length(); } - inline BoutReal + BoutReal interpolate_sheath_o2(const std::function& f) const { return ythis(f) * (1 - length()) + ynext(f) * length(); } - inline BoutReal extrapolate_next_o1(const Field3D& f) const { return ythis(f); } - inline BoutReal extrapolate_next_o2(const Field3D& f) const { + BoutReal extrapolate_next_o1(const Field3D& f) const { return ythis(f); } + BoutReal extrapolate_next_o2(const Field3D& f) const { ASSERT3(valid() >= 0); if (valid() < 1) { return extrapolate_next_o1(f); @@ -169,11 +173,11 @@ public: return ythis(f) * 2 - yprev(f); } - inline BoutReal + BoutReal extrapolate_next_o1(const std::function& f) const { return ythis(f); } - inline BoutReal + BoutReal extrapolate_next_o2(const std::function& f) const { ASSERT3(valid() >= 0); if (valid() < 1) { @@ -183,8 +187,8 @@ public: } // extrapolate the gradient into the boundary - inline BoutReal extrapolate_grad_o1(const Field3D& f) const { return 0; } - inline BoutReal extrapolate_grad_o2(const Field3D& f) const { + BoutReal extrapolate_grad_o1([[maybe_unused]] const Field3D& f) const { return 0; } + BoutReal extrapolate_grad_o2(const Field3D& f) const { ASSERT3(valid() >= 0); if (valid() < 1) { return extrapolate_grad_o1(f); diff --git a/include/bout/sys/parallel_stencils.hxx b/include/bout/sys/parallel_stencils.hxx index 1d55931ff3..ed6daa5558 100644 --- a/include/bout/sys/parallel_stencils.hxx +++ b/include/bout/sys/parallel_stencils.hxx @@ -20,7 +20,7 @@ inline BoutReal dirichlet_o2(BoutReal spacing0, BoutReal value0, BoutReal spacin } inline BoutReal neumann_o2([[maybe_unused]] BoutReal spacing0, BoutReal value0, BoutReal spacing1, BoutReal value1) { - return -spacing1 * value0 + value1; + return (-spacing1 * value0) + value1; } inline BoutReal dirichlet_o3(BoutReal spacing0, BoutReal value0, BoutReal spacing1, BoutReal value1, BoutReal spacing2, BoutReal value2) { diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index 74a68377c9..e5a9086397 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -54,10 +54,13 @@ public: mesh = bout::globals::mesh; } - bool lower_y = options["lower_y"].doc("Boundary on lower y?").withDefault(true); - bool upper_y = options["upper_y"].doc("Boundary on upper y?").withDefault(true); - bool outer_x = options["outer_x"].doc("Boundary on outer x?").withDefault(true); - bool inner_x = + const bool lower_y = + options["lower_y"].doc("Boundary on lower y?").withDefault(true); + const bool upper_y = + options["upper_y"].doc("Boundary on upper y?").withDefault(true); + const bool outer_x = + options["outer_x"].doc("Boundary on outer x?").withDefault(true); + const bool inner_x = options["inner_x"].doc("Boundary on inner x?").withDefault(false); if (mesh->isFci()) { diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index 28b8db9326..a14522b3d9 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -189,7 +189,7 @@ FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy, BoutMask to_remove(map_mesh); const int xend = map_mesh->xstart - + (map_mesh->xend - map_mesh->xstart + 1) * map_mesh->getNXPE() - 1; + + ((map_mesh->xend - map_mesh->xstart + 1) * map_mesh->getNXPE()) - 1; // Default to the maximum number of points const int defValid{map_mesh->ystart - 1 + std::abs(offset)}; // Serial loop because call to BoundaryRegionPar::addPoint @@ -262,7 +262,7 @@ FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy, ASSERT2(map_mesh->xend - map_mesh->xstart >= 2); auto boundary = (xt_prime[i] < map_mesh->xstart) ? inner_boundary : outer_boundary; if (!boundary->contains(x, y, z)) { - boundary->add_point(x, y, z, x + dx, y + offset - sgn(offset) * 0.5, + boundary->add_point(x, y, z, x + dx, y + offset - (sgn(offset) * 0.5), z + dz, // Intersection point in local index space std::abs(offset) - 0.5, // Distance to intersection defValid, offset); @@ -366,7 +366,7 @@ FCITransform::FCITransform(Mesh& mesh, const Coordinates::FieldMetric& dy, bool } const std::array bndries = {forward_boundary_xin, forward_boundary_xout, backward_boundary_xin, backward_boundary_xout}; - for (auto& bndry : bndries) { + for (const auto& bndry : bndries) { for (const auto& bndry2 : bndries) { if (bndry->dir == bndry2->dir) { continue; From 7975518600c1930a1c98dfd92ffec91e34b474a6 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 10 Mar 2026 12:38:49 +0100 Subject: [PATCH 54/62] Add missing header --- include/bout/sys/parallel_stencils.hxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/bout/sys/parallel_stencils.hxx b/include/bout/sys/parallel_stencils.hxx index ed6daa5558..2cbec9c8b0 100644 --- a/include/bout/sys/parallel_stencils.hxx +++ b/include/bout/sys/parallel_stencils.hxx @@ -1,4 +1,6 @@ #pragma once +#include +#include namespace parallel_stencil { // generated by src/mesh/parallel_boundary_stencil.cxx.py From 4ee3d9d2a728f32f247bc314c7fb71d133f222a1 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 11 Mar 2026 10:47:00 +0100 Subject: [PATCH 55/62] Add contains to YBoundary --- include/bout/yboundary_regions.hxx | 33 ++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index e5a9086397..50df198539 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -86,12 +86,45 @@ public: } is_init = true; optptr = &options; + // Cache boundary regions + _contains.emplace_back(mesh, false); + _contains.emplace_back(mesh, false); + iter_pnts([&](const auto& pnt) { + if (pnt.dir == 1) { + _contains[1][pnt.ind()] = true; + } else if (pnt.dir == -1) { + _contains[0][pnt.ind()] = true; + } + }); + } + + bool contains_ylow(Ind3D ind) const { return _contains[0][ind] } + bool contains_yhigh(Ind3D ind) const { return _contains[1][ind] } + template + bool contains(Ind3D ind) { + if constexpr (dir == 1) { + return _contains[1][ind]; + } else if constexpr (dir == -1) { + return _contains[0][ind]; + } else { + static_assert(false); // Only for +1 and -1 + } + } + bool contains(int dir, Ind3D ind) const { + if (dir == 1) { + return contains<+1>(ind); + } else if (dir == -1) { + return contains<-1>(ind); + } else { + throw BoutException("only dir == 1 and dir == -1 are implemented, not {}", dir); + } } private: std::vector> boundary_regions_par; std::vector> boundary_regions; + std::vector _contains; bool is_init{false}; Options* optptr; }; From 44fd33161c9ab7d2e24c82d45c7c133914610a5a Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 11 Mar 2026 11:05:14 +0100 Subject: [PATCH 56/62] Do not use else after return --- include/bout/yboundary_regions.hxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index 50df198539..96d0410dad 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -104,20 +104,20 @@ public: bool contains(Ind3D ind) { if constexpr (dir == 1) { return _contains[1][ind]; - } else if constexpr (dir == -1) { + } + if constexpr (dir == -1) { return _contains[0][ind]; - } else { - static_assert(false); // Only for +1 and -1 } + static_assert(false); // Only for +1 and -1 } bool contains(int dir, Ind3D ind) const { if (dir == 1) { return contains<+1>(ind); - } else if (dir == -1) { + } + if (dir == -1) { return contains<-1>(ind); - } else { - throw BoutException("only dir == 1 and dir == -1 are implemented, not {}", dir); } + throw BoutException("only dir == 1 and dir == -1 are implemented, not {}", dir); } private: From e94bb27003733544cd0a518d698d4825aaecad80 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 11 Mar 2026 11:26:24 +0100 Subject: [PATCH 57/62] Ensure YBoundary can compile --- include/bout/yboundary_regions.hxx | 9 +++++---- tests/integrated/test-fci-boundary/get_par_bndry.cxx | 6 ++++++ 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index 96d0410dad..68b10a1b52 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -6,6 +6,7 @@ #include "bout/globals.hxx" #include "bout/options.hxx" #include "bout/parallel_boundary_region.hxx" +#include "bout/mask.hxx" #include #include @@ -98,17 +99,17 @@ public: }); } - bool contains_ylow(Ind3D ind) const { return _contains[0][ind] } - bool contains_yhigh(Ind3D ind) const { return _contains[1][ind] } + bool contains_ylow(Ind3D ind) const { return _contains[0][ind]; } + bool contains_yhigh(Ind3D ind) const { return _contains[1][ind]; } template - bool contains(Ind3D ind) { + bool contains(Ind3D ind) const { + static_assert(dir == 1 || dir == -1); if constexpr (dir == 1) { return _contains[1][ind]; } if constexpr (dir == -1) { return _contains[0][ind]; } - static_assert(false); // Only for +1 and -1 } bool contains(int dir, Ind3D ind) const { if (dir == 1) { diff --git a/tests/integrated/test-fci-boundary/get_par_bndry.cxx b/tests/integrated/test-fci-boundary/get_par_bndry.cxx index 52d07f067e..7b2fd53c60 100644 --- a/tests/integrated/test-fci-boundary/get_par_bndry.cxx +++ b/tests/integrated/test-fci-boundary/get_par_bndry.cxx @@ -7,6 +7,8 @@ #include "bout/output.hxx" #include "bout/parallel_boundary_region.hxx" +#include "bout/yboundary_regions.hxx" + #include #include @@ -36,6 +38,10 @@ int main(int argc, char** argv) { dump[fmt::format("field_{:s}", boundary_name)] = fields[i]; } + Options dummy; + YBoundary ybndry; + ybndry.init(dummy, mesh); + bout::writeDefaultOutputFile(dump); BoutFinalise(); From f26bbfbe34e2b2ef8c80069c51ebb8faaf48a983 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 11 Mar 2026 12:24:35 +0100 Subject: [PATCH 58/62] Add test for YBoundary --- include/bout/yboundary_regions.hxx | 2 +- .../test-fci-boundary/get_par_bndry.cxx | 24 ++++++++++-- tests/integrated/test-fci-boundary/runtest | 39 ++++++++++++------- 3 files changed, 47 insertions(+), 18 deletions(-) diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index 68b10a1b52..8450b7c4a8 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -4,9 +4,9 @@ #include "bout/assert.hxx" #include "bout/boutexception.hxx" #include "bout/globals.hxx" +#include "bout/mask.hxx" #include "bout/options.hxx" #include "bout/parallel_boundary_region.hxx" -#include "bout/mask.hxx" #include #include diff --git a/tests/integrated/test-fci-boundary/get_par_bndry.cxx b/tests/integrated/test-fci-boundary/get_par_bndry.cxx index 7b2fd53c60..15530a340b 100644 --- a/tests/integrated/test-fci-boundary/get_par_bndry.cxx +++ b/tests/integrated/test-fci-boundary/get_par_bndry.cxx @@ -37,10 +37,28 @@ int main(int argc, char** argv) { dump[fmt::format("field_{:s}", boundary_name)] = fields[i]; } + for (const auto& name : {"forward_xt_prime", "backward_xt_prime"}) { + Field3D tmp; + mesh->get(tmp, name); + dump[name] = tmp; + } + + { + Options dummy; + YBoundary ybndry; + ybndry.init(dummy, mesh); + + std::vector fields(mesh->ystart * 2 + 1, Field3D{0.0}); + for (auto& field : fields) { + field.allocate(); + } + ybndry.iter_pnts( + [&](const auto& pnt) { fields[pnt.dir + mesh->ystart][pnt.ind()] += 1; }); - Options dummy; - YBoundary ybndry; - ybndry.init(dummy, mesh); + for (int i = -mesh->ystart; i <= mesh->ystart; ++i) { + dump[fmt::format("ybndry_{}", i)] = fields[i + mesh->ystart]; + } + } bout::writeDefaultOutputFile(dump); diff --git a/tests/integrated/test-fci-boundary/runtest b/tests/integrated/test-fci-boundary/runtest index e749055185..a54b13bd01 100755 --- a/tests/integrated/test-fci-boundary/runtest +++ b/tests/integrated/test-fci-boundary/runtest @@ -4,9 +4,20 @@ from boututils.run_wrapper import launch_safe from boututils.datafile import DataFile -from boutdata.collect import collect +from boutdata.collect import collect as _collect import numpy as np +from numpy.testing import assert_allclose + + +def collect(*args): + return _collect( + *args, + info=False, + path=directory, + xguards=False, + yguards=False, + ) nprocs = [1] @@ -29,27 +40,27 @@ regions = { } regions = {k: v.astype(int) for k, v in regions.items()} + for x in "xout", "xin": regions[x] = regions[f"{x}_fwd"] + regions[f"{x}_bwd"] for x in "fwd", "bwd": regions[x] = regions[f"xin_{x}"] + regions[f"xout_{x}"] regions["all"] = regions["xin"] + regions["xout"] +bndrys = { + "ybndry_-1": regions["xout_bwd"], + "ybndry_0": regions["xout_fwd"] * 0, + "ybndry_1": regions["xout_fwd"], +} + for nproc in nprocs: cmd = "./get_par_bndry" _, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) for k, v in regions.items(): - data = collect( - f"field_{k}", - info=False, - path=directory, - xguards=False, - yguards=False, - ) - assert np.allclose(data, v), ( - f"{k} does not match", - np.sum(data), - np.sum(v), - np.max(data), - ) + data = collect(f"field_{k}") + assert_allclose(data, v) + for i in range(-1, 2): + name = f"ybndry_{i}" + data = collect(name) + assert_allclose(bndrys[name], data) From 52cb85abf1df70206979a127f64a0881f7a9ec88 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 11 Mar 2026 12:34:07 +0100 Subject: [PATCH 59/62] Apply clang-tidy suggestions --- include/bout/boundary_iterator.hxx | 4 ++-- include/bout/yboundary_regions.hxx | 1 + tests/integrated/test-fci-boundary/get_par_bndry.cxx | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx index 7d8e7ed027..765f83e694 100644 --- a/include/bout/boundary_iterator.hxx +++ b/include/bout/boundary_iterator.hxx @@ -10,8 +10,8 @@ #include "bout/region.hxx" #include "bout/sys/parallel_stencils.hxx" #include "bout/sys/range.hxx" +#include #include -#include class BoundaryRegionIter { public: @@ -173,7 +173,7 @@ private: class BoundaryRegionIterY : public BoundaryRegionIter { public: - virtual ~BoundaryRegionIterY() = default; + ~BoundaryRegionIterY() override = default; BoundaryRegionIterY(const RangeIterator& r, int y, int dir, bool is_end, Mesh* mesh) : BoundaryRegionIter(r.ind, y, 0, dir, mesh), r(r), is_end(is_end) {} diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx index 8450b7c4a8..6254ef7862 100644 --- a/include/bout/yboundary_regions.hxx +++ b/include/bout/yboundary_regions.hxx @@ -7,6 +7,7 @@ #include "bout/mask.hxx" #include "bout/options.hxx" #include "bout/parallel_boundary_region.hxx" +#include "bout/region.hxx" #include #include diff --git a/tests/integrated/test-fci-boundary/get_par_bndry.cxx b/tests/integrated/test-fci-boundary/get_par_bndry.cxx index 15530a340b..0871577b86 100644 --- a/tests/integrated/test-fci-boundary/get_par_bndry.cxx +++ b/tests/integrated/test-fci-boundary/get_par_bndry.cxx @@ -48,7 +48,7 @@ int main(int argc, char** argv) { YBoundary ybndry; ybndry.init(dummy, mesh); - std::vector fields(mesh->ystart * 2 + 1, Field3D{0.0}); + std::vector fields((mesh->ystart * 2) + 1, Field3D{0.0}); for (auto& field : fields) { field.allocate(); } From 6d72a7e990c58c322183459955119870b1ec9fcd Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 12 Mar 2026 11:30:43 +0100 Subject: [PATCH 60/62] Use std::ranges::any_of --- include/bout/parallel_boundary_region.hxx | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 1d06bc65cb..fa8e8efbba 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -451,12 +451,7 @@ public: bool contains(const int ix, const int iy, const int iz) const { const auto i2 = xyz2ind(ix, iy, iz, localmesh); - for (auto i1 : bndry_points) { - if (i1.index == i2) { - return true; - } - } - return false; + return std::ranges::any_of(bndry_points, [](auto i1) { return i2.index == i2; }); } // setter From 5968f1387372f688d8647040330553f3888283f0 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 12 Mar 2026 12:14:30 +0100 Subject: [PATCH 61/62] Fixes for prek --- include/bout/field3d.hxx | 4 ++-- include/bout/parallel_boundary_region.hxx | 4 ++-- manual/sphinx/user_docs/boundary_options.rst | 12 ++++++------ src/field/field3d.cxx | 6 +++--- src/mesh/impls/bout/boutmesh.cxx | 6 +++--- src/mesh/parallel/shiftedmetricinterp.cxx | 2 +- tests/integrated/test-fci-boundary/runtest | 1 - 7 files changed, 17 insertions(+), 18 deletions(-) diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index f607445705..09d9c93f1b 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -2,7 +2,7 @@ * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu * * Contact: Ben Dudson, bd512@york.ac.uk - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -361,7 +361,7 @@ public: * Direct access to the underlying data array * * If CHECK > 2 then bounds checking is performed - * + * * If CHECK <= 2 then no checks are performed, to * allow inlining and optimisation of inner loops */ diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index fa8e8efbba..05338c8479 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -53,7 +53,7 @@ struct Indices { Indices(Ind3D index, RealPoint&& intersection, BoutReal length, signed char valid, signed char offset, unsigned char abs_offset) : index(index), intersection(intersection), length(length), valid(valid), - offset(offset), abs_offset(abs_offset){}; + offset(offset), abs_offset(abs_offset) {}; }; using IndicesVec = std::vector; @@ -119,7 +119,7 @@ public: BoundaryRegionParIterBase(IndicesVec& bndry_points, IndicesIter bndry_position, int dir, Mesh* localmesh) : bndry_points(bndry_points), bndry_position(bndry_position), dir(dir), - localmesh(localmesh){}; + localmesh(localmesh) {}; // getter Ind3D ind() const { return bndry_position->index; } diff --git a/manual/sphinx/user_docs/boundary_options.rst b/manual/sphinx/user_docs/boundary_options.rst index 9f9df6736c..f04d09cabc 100644 --- a/manual/sphinx/user_docs/boundary_options.rst +++ b/manual/sphinx/user_docs/boundary_options.rst @@ -450,7 +450,7 @@ Previously of they where implemented using a `RangeIterator`:: lower_y = options["lower_y"].doc("Boundary on lower y?").withDefault(lower_y); upper_y = options["upper_y"].doc("Boundary on upper y?").withDefault(upper_y); } - + void rhs() { BoutReal totalFlux = 0; if (lower_y) { @@ -478,14 +478,14 @@ Previously of they where implemented using a `RangeIterator`:: } } } - + private: bool lower_y{true}; bool upper_y{true}; const Field3D& N; const Field3D& V; } - + This can be replaced using the `YBoundary` class, which not only simplifies the code, but also allows to have the same code working with non-field-aligned @@ -499,20 +499,20 @@ geometries, as flux coordinate independent (FCI) method:: // Set what kind of yboundaries you want to include yboundary.init(opt); } - + void rhs() { BoutReal totalFlux = 0; yboundary.iter_pnts([&](auto& pnt) { BoutReal flux = pnt.interpolate_sheath_o2(N) * pnt.interpolate_sheath_o2(V); }); } - + private: YBoundary ybounday; const Field3D& N; const Field3D& V; }; - + There are several member functions of ``pnt``. ``pnt`` is of type diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 1f04d3e099..7f5e4d4b03 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -7,7 +7,7 @@ * Copyright 2010 - 2025 BOUT++ developers * * Contact: Ben Dudson, dudson2@llnl.gov - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -230,7 +230,7 @@ const Region& Field3D::getRegion2D(const std::string& region_name) const } /*************************************************************** - * OPERATORS + * OPERATORS ***************************************************************/ /////////////////// ASSIGNMENT //////////////////// @@ -734,7 +734,7 @@ Field3D lowPass(const Field3D& var, int zmax, bool keep_zonal, const std::string return result; } -/* +/* * Use FFT to shift by an angle in the Z direction */ void shiftZ(Field3D& var, int jx, int jy, double zangle) { diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index ae7690aff7..a19e416beb 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -2223,9 +2223,9 @@ void BoutMesh::topology() { } for (int i = 0; i < limiter_count; ++i) { - int const yind = limiter_yinds[i]; - int const xstart = limiter_xstarts[i]; - int const xend = limiter_xends[i]; + const int yind = limiter_yinds[i]; + const int xstart = limiter_xstarts[i]; + const int xend = limiter_xends[i]; output_info.write("Adding a limiter between y={} and {}. X indices {} to {}\n", yind, yind + 1, xstart, xend); add_target(yind, xstart, xend); diff --git a/src/mesh/parallel/shiftedmetricinterp.cxx b/src/mesh/parallel/shiftedmetricinterp.cxx index 22a8f0748a..4b470da2ad 100644 --- a/src/mesh/parallel/shiftedmetricinterp.cxx +++ b/src/mesh/parallel/shiftedmetricinterp.cxx @@ -215,7 +215,7 @@ void ShiftedMetricInterp::checkInputGrid() { "Should be 'orthogonal'."); } } // else: coordinate_system variable not found in grid input, indicates older input - // file so must rely on the user having ensured the type is correct + // file so must rely on the user having ensured the type is correct } /*! diff --git a/tests/integrated/test-fci-boundary/runtest b/tests/integrated/test-fci-boundary/runtest index a54b13bd01..74a8abf741 100755 --- a/tests/integrated/test-fci-boundary/runtest +++ b/tests/integrated/test-fci-boundary/runtest @@ -6,7 +6,6 @@ from boututils.run_wrapper import launch_safe from boututils.datafile import DataFile from boutdata.collect import collect as _collect -import numpy as np from numpy.testing import assert_allclose From ade82b73d1eb4ea7e364094fbff7bddb0e354db8 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 12 Mar 2026 12:47:36 +0100 Subject: [PATCH 62/62] Fix usage of any_of --- include/bout/parallel_boundary_region.hxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 05338c8479..d42a2b2366 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -451,7 +451,8 @@ public: bool contains(const int ix, const int iy, const int iz) const { const auto i2 = xyz2ind(ix, iy, iz, localmesh); - return std::ranges::any_of(bndry_points, [](auto i1) { return i2.index == i2; }); + return std::ranges::any_of(bndry_points.begin(), bndry_points.end(), + [&i2](auto i1) { return i1.index == i2; }); } // setter