diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs index 43a5292e9d..945f769e70 100644 --- a/.git-blame-ignore-revs +++ b/.git-blame-ignore-revs @@ -16,6 +16,7 @@ d1cfb8abd6aa5c76e6c1a4d7ab20929c65f8afc2 8d5cb39e03c2644715a50684f8cd0318b4e82676 ec69e8838be2dde140a915e50506f8e1ce3cb534 f2bc0488a298f136294c523bc5ab4086d090059b +c59626a0498b0cfd1be46f84f0fa38b4855a43cc 1b4707187a3a85126338303dc766280b8fb2dc56 b2e2f3575e68f771be2a4341af5fd14e2870469e 64e4285ec38569f66d31e589a4610cefd16d8076 diff --git a/include/bout/boundary_iterator.hxx b/include/bout/boundary_iterator.hxx new file mode 100644 index 0000000000..66c66d45b6 --- /dev/null +++ b/include/bout/boundary_iterator.hxx @@ -0,0 +1,231 @@ +#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/region.hxx" +#include "bout/sys/parallel_stencils.hxx" +#include "bout/sys/range.hxx" +#include +#include + +class BoundaryRegionIter { +public: + virtual ~BoundaryRegionIter() = default; + BoundaryRegionIter(int x, int y, int bx, int by, const Mesh& mesh) + : _dir(bx + by), _x(x), _y(y), _bx(bx), _by(by), localmesh(&mesh) { + ASSERT3(_bx * _by == 0); + } + bool operator!=(const BoundaryRegionIter& rhs) const { 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; + } + +protected: + virtual void _next() = 0; + +public: + BoundaryRegionIter& operator*() { return *this; } + + void dirichlet_o2(Field3D& f, BoutReal value) const { + ynext(f) = bout::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_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))); + } + + BoutReal extrapolate_sheath_free(const Field3D& f, SheathLimitMode mode) const { + const BoutReal fac = + bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f), mode); + const BoutReal val = ythis(f); + const BoutReal next = mode == SheathLimitMode::linear_free ? val + fac : val * fac; + 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)); + BoutReal val = ythis(f); + for (int i = 1; i <= localmesh->ystart; ++i) { + val *= fac; + f[ind().yp(_by * i).xp(_bx * i)] = val; + } + } + + 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) { + 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; + } + } + + void limit_at_least(Field3D& f, BoutReal value) const { + ynext(f) = std::max(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)]; } + 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; + } + } + + static int abs_offset() { return 1; } + + 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)]; } + + int dir() const { return _dir; } + +protected: + int x() const { return _x; } + int y() const { return _y; } + int z() const { return _z; } + +protected: + void setx(int x) { _x = x; } + void sety(int y) { _y = y; } + +private: + int _dir; + + int _z{0}; + int _x; + int _y; + int _bx; + int _by; + + const 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() override = default; + BoundaryRegionIterY(const RangeIterator& r, int y, int dir, bool is_end, + const 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(); + } + + void _next() override { + ++r; + setx(r.ind); + } + +private: + RangeIterator r; + bool is_end; +}; + +class NewBoundaryRegionY { +public: + NewBoundaryRegionY(const 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); + } + BoundaryRegionIterY end() { return begin(false); } + +private: + const Mesh* mesh; + bool lower; + RangeIterator r; +}; diff --git a/include/bout/boundary_region.hxx b/include/bout/boundary_region.hxx index 58de12045e..2a89cc1a50 100644 --- a/include/bout/boundary_region.hxx +++ b/include/bout/boundary_region.hxx @@ -4,6 +4,7 @@ class BoundaryRegion; #ifndef BOUT_BNDRY_REGION_H #define BOUT_BNDRY_REGION_H +#include #include #include @@ -15,7 +16,7 @@ extern Mesh* mesh; ///< Global mesh } // namespace bout /// Location of boundary -enum class BndryLoc { +enum class BndryLoc : std::int8_t { xin, xout, ydown, @@ -36,6 +37,9 @@ constexpr BndryLoc BNDRY_PAR_BKWD_XIN = BndryLoc::par_bkwd_xin; constexpr BndryLoc BNDRY_PAR_FWD_XOUT = BndryLoc::par_fwd_xout; constexpr BndryLoc BNDRY_PAR_BKWD_XOUT = BndryLoc::par_bkwd_xout; +/// Physical type of y boundary +enum class YBndryType : std::int8_t { sheath, not_sheath, all }; + class BoundaryRegionBase { public: BoundaryRegionBase() = delete; diff --git a/include/bout/boundary_region_iter.hxx b/include/bout/boundary_region_iter.hxx new file mode 100644 index 0000000000..ee0132c527 --- /dev/null +++ b/include/bout/boundary_region_iter.hxx @@ -0,0 +1,547 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace bout { +namespace boundary { + +/// Physical type of y boundary +enum class BndryType : std::int8_t { + sheath, + not_sheath_par, + core, + sol, + sol_perp, + sol_par, + all, + num +}; + +/// Types of free boundary condition +/// ================== =================================================== +/// Name Description +/// ================== =================================================== +/// ``limited`` use exponential if decreasing, otherwise Neumanm +/// ``exponential`` use exponential extrapolation +/// ``linear`` use linear extrapolation +/// ================== =================================================== +enum class BoundaryFreeExtrapolation : std::int8_t { limited, exponential, linear }; +//BOUT_ENUM_CLASS(BoundaryFreeExtrapolation, limited, exponential, linear); + +template +class BoundaryRegionIterBase { /// get the index at the last point in domain +public: + Ind3D ind() const { return static_cast(this)->_ind(); } + /// get the length from the point in the domain to the boundary in index + /// space. It is in the range [0, 1] + BoutReal length() const { return static_cast(this)->_length(); } + /// Lower bound of how many points are between the first point in the domain + /// and the boundary in the other direction. + signed char valid() const { return static_cast(this)->_valid(); } + /// Get the width of the boundary at the current point + int boundary_width() const { return static_cast(this)->_boundary_width(); } + /// Is this the lower boundary? + bool is_lower() const { return static_cast(this)->_is_lower(); } + + /* + * FIELD3D ACCESSORS + */ + + /// get the value at a given offset `off` of a field `f`. + /// off = -1 is the second point in the boundary + /// off = 0 is the first point in the boundary + /// off = 1 is the last point in the domain + /// off = 2 is the second to last point in the domain + template + BoutReal& getAt(Field3D& f, int off) const { + return static_cast(this)->template _getAt(f, off); + } + /// get the value at a given offset `off` of a field `f`. + template + BoutReal& getAt(const Field3D& f, int off) const { + return static_cast(this)->template _getAt(f, off); + } + + /// Get the first point in the boundary + const BoutReal& next(const Field3D& f) const { + return static_cast(this)->_getAt(f, 0); + } + /// Get the first point in the boundary + BoutReal& next(Field3D& f) const { + return static_cast(this)->_getAt(f, 0); + } + /// Get the last point in the domain + const BoutReal& current(const Field3D& f) const { + return static_cast(this)->_getAt(f, 1); + } + /// Get the last point in the domain + BoutReal& current(Field3D& f) const { + return static_cast(this)->_getAt(f, 1); + } + /// Get the second to last point in the domain - this may not be valid and thus throw + const BoutReal& prev(const Field3D& f) const { + return static_cast(this)->_getAt(f, 2); + } + /// Get the offset from the last point in the domain + /// For FA this is always ±1, for FCI this can be up to ±MYG, excluding 0 + int offset() const { return static_cast(this)->_offset(); } + /* + * FIELD2D ACCESSORS + */ + + /// get the value at a given offset `off` of a field `f`. + /// off = -1 is the second point in the boundary + /// off = 0 is the first point in the boundary + /// off = 1 is the last point in the domain + /// off = 2 is the second to last point in the domain + template + BoutReal& getAt(Field2D& f, int off) const { + return static_cast(this)->template _getAt(f, off); + } + /// get the value at a given offset `off` of a field `f`. + template + BoutReal& getAt(const Field2D& f, int off) const { + return static_cast(this)->template _getAt(f, off); + } + + /// Get the first point in the boundary + const BoutReal& next(const Field2D& f) const { + return static_cast(this)->_getAt(f, 0); + } + /// Get the first point in the boundary + BoutReal& next(Field2D& f) const { + return static_cast(this)->_getAt(f, 0); + } + /// Get the last point in the domain + const BoutReal& current(const Field2D& f) const { + return static_cast(this)->_getAt(f, 1); + } + /// Get the last point in the domain + BoutReal& current(Field2D& f) const { + return static_cast(this)->_getAt(f, 1); + } + /// Get the second to last point in the domain - this may not be valid and thus throw + const BoutReal& prev(const Field2D& f) const { + return static_cast(this)->_getAt(f, 2); + } + + /* + * FUNCTIONS ACCESSORS + */ + + /// get the value at a given offset `off` of a field `f`. + /// off = -1 is the second point in the boundary + /// off = 0 is the first point in the boundary + /// off = 1 is the last point in the domain + /// off = 2 is the second to last point in the domain + template + BoutReal& getAt(const std::function& func, + int off) const { + return static_cast(this)->template _getAt(func, off); + } + /// Get the first point in the boundary + const BoutReal& + next(const std::function& func) const { + return static_cast(this)->_getAt(func, 0); + } + /// Get the last point in the domain + const BoutReal& + current(const std::function& func) const { + return static_cast(this)->_getAt(func, 1); + } + /// Get the second to last point in the domain - this may not be valid and thus throw + const BoutReal& + prev(const std::function& func) const { + return static_cast(this)->_getAt(func, 2); + } + + /* + * INTERPOLATION and EXTRAPOLATION + */ + + // extrapolate a given field to the boundary + BoutReal extrapolate_boundary_o1(const Field3D& f) const { return current(f); } + // extrapolate a given field to the boundary + BoutReal extrapolate_boundary_o2(const Field3D& f) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return extrapolate_boundary_o1(f); + } + return current(f) * (1 + length()) - prev(f) * length(); + } + /// Extrapolate a given function to the boundary + BoutReal extrapolate_bounday_o1( + const std::function& func) const { + return current(func); + } + /// Extrapolate a given function to the boundary + BoutReal extrapolate_boundary_o2( + const std::function& func) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return extrapolate_boundary_o1(func); + } + return current(func) * (1 + length()) - prev(func) * length(); + } + + /// Interpolate a field to the boundary, using the boundary values + BoutReal interpolate_boundary_o2(const Field3D& f) const { + return current(f) * (1 - length()) + next(f) * length(); + } + /// Interpolate a field to the boundary, using the boundary values + BoutReal interpolate_boundary_o2( + const std::function& func) const { + return current(func) * (1 - length()) + next(func) * length(); + } + /// Extrapolate to the first boundary value freely + BoutReal extrapolate_next_o1(const Field3D& f) const { return current(f); } + /// Extrapolate to the first boundary value freely + BoutReal extrapolate_next_o2(const Field3D& f) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return extrapolate_next_o1(f); + } + return current(f) * 2 - prev(f); + } + + /// Extrapolate to the first boundary value freely + BoutReal + extrapolate_next_o1(const std::function& func) const { + return current(func); + } + /// Extrapolate to the first boundary value freely + BoutReal + extrapolate_next_o2(const std::function& func) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return extrapolate_next_o1(func); + } + return current(func) * 2 - prev(func); + } + + /// extrapolate the gradient into the boundary + BoutReal extrapolate_grad_o1([[maybe_unused]] const Field3D& f) const { return 0; } + /// extrapolate the gradient into the boundary + BoutReal extrapolate_grad_o2(const Field3D& f) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return extrapolate_grad_o1(f); + } + return current(f) - next(f); + } + + BoutReal extrapolate_boundary_free(const Field3D& f, + BoundaryFreeExtrapolation mode) const { + const auto fac = valid() > 0 ? limitFreeScale(prev(f), current(f), mode) + : (mode == BoundaryFreeExtrapolation::linear ? 0 : 1); + auto val = current(f); + BoutReal next = mode == BoundaryFreeExtrapolation::linear ? val + fac : val * fac; + return val * length() + next * (1 - length()); + } + + /* + * APPLY BOUNDARY CONDITIONS + */ + + /// Apply a dirichlet boundary condition + void dirichlet_o1(Field3D& f, BoutReal value) const { + for (int i = 0; i < boundary_width(); ++i) { + getAt(f, i) = value; + } + } + + /// Apply a dirichlet boundary condition + void dirichlet_o2(Field3D& f, BoutReal value) const { + if (length() < small_value) { + return dirichlet_o1(f, value); + } + for (int i = 0; i < boundary_width(); ++i) { + getAt(f, i) = + parallel_stencil::dirichlet_o2(i + 1, current(f), i + 1 - length(), value); + } + } + + /// Apply a dirichlet boundary condition + void dirichlet_o3(Field3D& f, BoutReal value) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return dirichlet_o2(f, value); + } + if (length() < small_value) { + for (int i = 0; i < boundary_width(); ++i) { + getAt(f, i) = + parallel_stencil::dirichlet_o2(i + 2, prev(f), i + 1 - length(), value); + } + } else { + for (int i = 0; i < boundary_width(); ++i) { + getAt(f, i) = parallel_stencil::dirichlet_o3(i + 2, prev(f), i + 1, current(f), + i + 1 - length(), value); + } + } + } + + /// Ensure the value in the boundary is at least `value` + void limit_at_least(Field3D& f, BoutReal value) const { + for (int i = 0; i < boundary_width(); ++i) { + if (getAt(f, i) < value) { + getAt(f, i) = value; + } + } + } + + /// Apply neumann boundary condition, where `value` is the gradient in index space + + // neumann_o1 would give second order convergence, given 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 { + for (int i = 0; i < boundary_width(); ++i) { + getAt(f, i) = current(f) + value * (i + 1); + } + } + + /// Apply neumann boundary condition, where `value` is the gradient in index space + void neumann_o2(Field3D& f, BoutReal value) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return neumann_o1(f, value); + } + for (int i = 0; i < boundary_width(); ++i) { + getAt(f, i) = prev(f) + (2 + i) * value; + } + } + + /// Apply neumann boundary condition, where `value` is the gradient in index space + void neumann_o3(Field3D& f, BoutReal value) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return neumann_o2(f, value); + } + for (int i = 0; i < boundary_width(); ++i) { + getAt(f, i) = parallel_stencil::neumann_o3(i + 1 - length(), value, i + 1, + current(f), 2, prev(f)); + } + } + + void set_free(Field3D& f, BoundaryFreeExtrapolation mode) const { + int fac; + if (valid() > 0) { + fac = limitFreeScale(prev(f), current(f), mode); + } else { + fac = mode == BoundaryFreeExtrapolation::linear ? 0 : 1; + } + auto val = current(f); + if (mode == BoundaryFreeExtrapolation::linear) { + for (int i = 0; i < boundary_width(); ++i) { + val += fac; + getAt(f, i) = val; + } + } else { + for (int i = 0; i < boundary_width(); ++i) { + val *= fac; + getAt(f, i) = val; + } + } + } + void setSmallValue(BoutReal val) { + ASSERT2(val > 0); + ASSERT2(val < 0.5); + small_value = val; + } + +private: + BoutReal small_value = 1e-4; +}; + +namespace { +/// 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, BoundaryFreeExtrapolation mode) { + if ((fm < fc) && (mode == BoundaryFreeExtrapolation::limited)) { + return fc; // Neumann rather than increasing into boundary + } + if (fm < 1e-10) { + return fc; // Low / no density condition + } + + BoutReal fp = 0; + switch (mode) { + case BoundaryFreeExtrapolation::limited: + case BoundaryFreeExtrapolation::exponential: + fp = SQ(fc) / fm; // Exponential + break; + case BoundaryFreeExtrapolation::linear: + fp = (2.0 * fc) - fm; // Linear + break; + } + +#if CHECKLEVEL >= 2 + if (!std::isfinite(fp)) { + throw BoutException("SheathBoundary limitFree {}: {}, {} -> {}", + static_cast(mode), fm, fc, fp); + } +#endif + + return fp; +} +} // namespace + +class BoundaryRegionFCI : public BoundaryRegionBase { +public: + BoundaryRegionFCI(const std::string& name, [[maybe_unused]] const BndryLoc& loc, + int dir, Mesh* mesh) + : BoundaryRegionBase(name, mesh), _dir(dir), localmesh(mesh) {}; + /// 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.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) { + add_point(xyz2ind(ix, iy, iz), x, y, z, length, valid, offset); + } + bool contains(int ix, int iy, int iz) { + const auto ind = xyz2ind(ix, iy, iz); + ensureSorted(); + const auto found = + std::lower_bound(std::begin(bndry_points), std::end(bndry_points), ind); + return found != std::end(bndry_points) and found->index == ind; + } + int dir() { return _dir; } + // legacy interface + void first() override { throw BoutException("Legacy interface is not suppored"); } + void next() override { throw BoutException("Legacy interface is not suppored"); } + bool isDone() override { throw BoutException("Legacy interface is not suppored"); } + +private: + friend class BoundaryRegionIterFCI; + int _dir; + // Vector of points in the boundary + bout::parallel_boundary_region::IndicesVec bndry_points; + bool is_sorted{true}; + Mesh* localmesh; + Ind3D xyz2ind(int x, int y, int z) const { + const int ny = localmesh->LocalNy; + const int nz = localmesh->LocalNz; + return Ind3D{(x * ny + y) * nz + z, ny, nz}; + } + void ensureSorted() { + if (is_sorted) { + return; + } + std::sort(std::begin(bndry_points), std::end(bndry_points)); + } +}; + +class BoundaryRegionIterFCI : public BoundaryRegionIterBase { +private: + // TODO(dave) make non-const? + const BoundaryRegionFCI* region; + size_t pos{0}; + +public: + BoundaryRegionIterFCI() = delete; + BoundaryRegionIterFCI(const BoundaryRegionFCI* reg, bool isstart) + : region(reg), pos(isstart ? 0 : reg->bndry_points.size()) {} + void setValid(char valid) { + const_cast(region)->bndry_points[pos].valid = valid; + }; + BoutReal s_x() const { return region->bndry_points[pos].intersection.s_x; }; + BoutReal s_y() const { return region->bndry_points[pos].intersection.s_y; }; + BoutReal s_z() const { return region->bndry_points[pos].intersection.s_z; }; + Mesh* localmesh() const { return region->localmesh; }; + int dir() const { return region->_dir; } + template + BoutReal& _getAt(Field3D& f, int off) const { + ASSERT3(f.hasParallelSlices()); + if constexpr (check) { + ASSERT3(_valid() > -off - 2); + } + auto _off = _offset() + off * region->_dir; + return f.ynext(_off)[ind().yp(_off)]; + } + template + const BoutReal& _getAt(const Field3D& f, int off) const { + ASSERT3(f.hasParallelSlices()); + if constexpr (check) { + ASSERT3(_valid() > -off - 2); + } + auto _off = _offset() + off * region->_dir; + return f.ynext(_off)[ind().yp(_off)]; + } + template + BoutReal& _getAt(Field2D& f, int off) const { + ASSERT3(f.hasParallelSlices()); + if constexpr (check) { + ASSERT3(_valid() > -off - 2); + } + auto _off = _offset() + off * region->_dir; + return f.ynext(_off)[ind().yp(_off)]; + } + template + const BoutReal& _getAt(const Field2D& f, int off) const { + ASSERT3(f.hasParallelSlices()); + if constexpr (check) { + ASSERT3(_valid() > -off - 2); + } + auto _off = _offset() + off * region->_dir; + return f.ynext(_off)[ind().yp(_off)]; + } + template + BoutReal getAt(const std::function& f, + int off) const { + if constexpr (check) { + ASSERT3(valid() > -off - 2); + } + auto _off = _offset() + off * region->_dir; + return f(_off, ind().yp(_off)); + } + signed char _offset() const { return region->bndry_points[pos].offset; } + signed char _valid() const { return region->bndry_points[pos].valid; } + Ind3D _ind() const { return region->bndry_points[pos].index; } + signed char _boundary_width() const { + return region->localmesh->ystart - region->bndry_points[pos].abs_offset; + } + const BoutReal& _length() const { return region->bndry_points[pos].length; } + bool operator!=(BoundaryRegionIterFCI lhs) { + ASSERT3(region == lhs.region); + return pos != lhs.pos; + } + BoundaryRegionIterFCI& operator++() { + ++pos; + return *this; + } + // No-op for compatibility + BoundaryRegionIterFCI& operator*() { return *this; } +}; +} // namespace boundary +} // namespace bout + +inline bout::boundary::BoundaryRegionIterFCI +begin(const bout::boundary::BoundaryRegionFCI reg) { + return bout::boundary::BoundaryRegionIterFCI(®, true); +} +inline bout::boundary::BoundaryRegionIterFCI +end(const bout::boundary::BoundaryRegionFCI reg) { + return bout::boundary::BoundaryRegionIterFCI(®, false); +} diff --git a/include/bout/bout_enum_class.hxx b/include/bout/bout_enum_class.hxx index 6f632747ed..334fa657da 100644 --- a/include/bout/bout_enum_class.hxx +++ b/include/bout/bout_enum_class.hxx @@ -2,7 +2,7 @@ * Copyright 2019 B.D.Dudson, J.T.Omotani, P.Hill * * Contact Ben Dudson, bd512@york.ac.uk - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -26,8 +26,9 @@ #include "bout/macro_for_each.hxx" #include "bout/options.hxx" // IWYU pragma: keep -#include // IWYU pragma: keep -#include // IWYU pragma: keep +#include // IWYU pragma: keep +#include // IWYU pragma: keep +#include // IWYU pragma: keep /// Create some macro magic similar to bout/macro_for_each.hxx, but allowing for the enum /// class name to be passed through to each _call @@ -66,7 +67,7 @@ /// Create an enum class with toString and FromString functions, and an /// Options::as overload to read the enum #define BOUT_ENUM_CLASS(enumname, ...) \ - enum class enumname { __VA_ARGS__ }; \ + enum class enumname : std::int8_t { __VA_ARGS__ }; \ \ inline std::string toString(enumname e) { \ \ diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index e7ead42ee5..be7de0f5b6 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -3,16 +3,16 @@ * * ChangeLog * ========= - * + * * 2014-11-10 Ben Dudson * * Created by separating metric from Mesh * - * + * ************************************************************************** * Copyright 2014-2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -33,13 +33,18 @@ #ifndef BOUT_COORDINATES_H #define BOUT_COORDINATES_H +#include "bout/field_data.hxx" #include #include #include #include #include +#include +#include + class Mesh; +class YBoundary; /*! * Represents a coordinate system, and associated operators @@ -221,9 +226,13 @@ public: // solver Field2D Laplace_perpXY(const Field2D& A, const Field2D& f); + friend std::shared_ptr getYBoundary(Coordinates* coords, YBndryType type); + private: + std::shared_ptr makeYBoundary(YBndryType type) const; int nz; // Size of mesh in Z. This is mesh->ngz-1 Mesh* localmesh; + Options* localoptions; CELL_LOC location; /// Handles calculation of yup and ydown @@ -249,6 +258,8 @@ private: void checkCovariant(); // check that contravariant tensors are positive (if expected) and finite (always) void checkContravariant(); + + mutable std::array, 3> ybndrys; }; /* @@ -256,10 +267,10 @@ private: class TokamakCoordinates : public Coordinates { public: TokamakCoordinates(Mesh *mesh) : Coordinates(mesh) { - + } private: - + }; */ diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 785b139fbb..81cba6d809 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -498,7 +498,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/include/bout/field_data.hxx b/include/bout/field_data.hxx index 185dcabf2d..4b8519d54b 100644 --- a/include/bout/field_data.hxx +++ b/include/bout/field_data.hxx @@ -6,7 +6,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 @@ -33,6 +33,7 @@ class FieldData; #include "bout/bout_types.hxx" #include "bout/unused.hxx" +#include #include #include #include @@ -40,12 +41,15 @@ class FieldData; class BoundaryOp; class BoundaryOpPar; +namespace bout::boundary { +class BoundaryRegionFCI; +} class Coordinates; class Mesh; #include "bout/boundary_region.hxx" class BoundaryRegionPar; -enum class BndryLoc; +enum class BndryLoc : std::int8_t; #include "bout/sys/expressionparser.hxx" @@ -77,7 +81,7 @@ public: /// Number of BoutReals in one element virtual int elementSize() const { return 1; } - virtual void doneComms(){}; // Notifies that communications done + virtual void doneComms() {}; // Notifies that communications done // Boundary conditions void setBoundary(const std::string& name); ///< Set the boundary conditions @@ -86,13 +90,13 @@ public: copyBoundary(const FieldData& f); ///< Copy the boundary conditions from another field virtual void applyBoundary(bool UNUSED(init) = false) {} - virtual void applyTDerivBoundary(){}; + virtual void applyTDerivBoundary() {}; - virtual void applyParallelBoundary(){}; - virtual void applyParallelBoundary(BoutReal UNUSED(t)){}; - virtual void applyParallelBoundary(const std::string& UNUSED(condition)){}; + virtual void applyParallelBoundary() {}; + virtual void applyParallelBoundary(BoutReal UNUSED(t)) {}; + virtual void applyParallelBoundary(const std::string& UNUSED(condition)) {}; virtual void applyParallelBoundary(const std::string& UNUSED(region), - const std::string& UNUSED(condition)){}; + const std::string& UNUSED(condition)) {}; // JMAD void addBndryFunction(FuncPtr userfunc, BndryLoc location); void addBndryGenerator(FieldGeneratorPtr gen, BndryLoc location); diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index a1ed6a9011..f982bcef1f 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -3,7 +3,7 @@ * * Interface for mesh classes. Contains standard variables and useful * routines. - * + * * Changelog * ========= * @@ -11,7 +11,7 @@ * * Removing coordinate system into separate * Coordinates class * * Adding index derivative functions from derivs.cxx - * + * * 2010-06 Ben Dudson, Sean Farley * * Initial version, adapted from GridData class * * Incorporates code from topology.cpp and Communicator @@ -20,7 +20,7 @@ * Copyright 2010-2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -67,7 +67,9 @@ class Mesh; #include class BoundaryRegion; -class BoundaryRegionPar; +namespace bout::boundary { +class BoundaryRegionFCI; +} class GridDataSource; class MeshFactory : public Factory { @@ -500,12 +502,13 @@ public: /// mesh->getBoundariesPar(Mesh::BoundaryParType::all) /// get only xout: /// mesh->getBoundariesPar(Mesh::BoundaryParType::xout) - virtual std::vector> - getBoundariesPar(BoundaryParType type = BoundaryParType::all) = 0; + virtual std::vector> + getBoundariesPar(BoundaryParType type = BoundaryParType::all) const = 0; /// Add a parallel(Y) boundary to this processor - virtual void addBoundaryPar(std::shared_ptr UNUSED(bndry), - BoundaryParType UNUSED(type)) {} + virtual void + addBoundaryPar(std::shared_ptr UNUSED(bndry), + BoundaryParType UNUSED(type)) {} virtual BoutReal GlobalX(int jx) const = 0; ///< Continuous X index between 0 and 1 virtual BoutReal GlobalY(int jy) const = 0; ///< Continuous Y index (0 -> 1) diff --git a/include/bout/msg_stack.hxx b/include/bout/msg_stack.hxx index ddc0ccad98..aa19939058 100644 --- a/include/bout/msg_stack.hxx +++ b/include/bout/msg_stack.hxx @@ -77,7 +77,8 @@ public: /// Dummy functions which should be optimised out int push(const std::string&) { return 0; } template - int push(fmt::format_string format, const Args&... args) { + int push([[maybe_unused]] fmt::format_string format, + [[maybe_unused]] const Args&... args) { return 0; } diff --git a/include/bout/parallel_boundary_op.hxx b/include/bout/parallel_boundary_op.hxx index d8620e892b..b0f74b9f9f 100644 --- a/include/bout/parallel_boundary_op.hxx +++ b/include/bout/parallel_boundary_op.hxx @@ -2,6 +2,7 @@ #define BOUT_PAR_BNDRY_OP_H #include "bout/boundary_op.hxx" +#include "bout/boundary_region_iter.hxx" #include "bout/bout_types.hxx" #include "bout/field_factory.hxx" #include "bout/parallel_boundary_region.hxx" @@ -16,28 +17,46 @@ class BoundaryOpPar : public BoundaryOpBase { public: BoundaryOpPar() = default; - BoundaryOpPar(BoundaryRegionPar* region, std::shared_ptr value) + BoundaryOpPar(bout::boundary::BoundaryRegionFCI* region, + std::shared_ptr value) : bndry(region), gen_values(std::move(value)), value_type(ValueType::GEN) {} - BoundaryOpPar(BoundaryRegionPar* region, Field3D* value) + BoundaryOpPar(bout::boundary::BoundaryRegionFCI* region, Field3D* value) : bndry(region), field_values(value), value_type(ValueType::FIELD) {} - BoundaryOpPar(BoundaryRegionPar* region, BoutReal value) + BoundaryOpPar(bout::boundary::BoundaryRegionFCI* region, BoutReal value) : bndry(region), real_value(value), value_type(ValueType::REAL) {} - BoundaryOpPar(BoundaryRegionPar* region) + BoundaryOpPar(bout::boundary::BoundaryRegionFCI* region) : bndry(region), real_value(0.), value_type(ValueType::REAL) {} + BoundaryOpPar(BoundaryOpPar* region, std::shared_ptr value) + : bndry(region->bndry), gen_values(std::move(value)), value_type(ValueType::GEN) {} + BoundaryOpPar(BoundaryOpPar* region, Field3D* value) + : bndry(region->bndry), field_values(value), value_type(ValueType::FIELD) {} + BoundaryOpPar(BoundaryOpPar* region, BoutReal value) + : bndry(region->bndry), real_value(value), value_type(ValueType::REAL) {} + BoundaryOpPar(BoundaryOpPar* region) + : bndry(region->bndry), real_value(0.), value_type(ValueType::REAL) {} ~BoundaryOpPar() override = default; // Note: All methods must implement clone, except for modifiers (see below) - virtual BoundaryOpPar* clone(BoundaryRegionPar* region, + virtual BoundaryOpPar* clone(bout::boundary::BoundaryRegionFCI* region, const std::list& args) = 0; - virtual BoundaryOpPar* clone(BoundaryRegionPar* region, Field3D* f) = 0; + virtual BoundaryOpPar* clone(BoundaryOpPar* region, + const std::list& args) = 0; + virtual BoundaryOpPar* clone(bout::boundary::BoundaryRegionFCI* region, Field3D* f) = 0; + virtual BoundaryOpPar* clone(BoundaryOpPar* region, Field3D* f) = 0; + virtual BoundaryOpPar* + clone(bout::boundary::BoundaryRegionFCI* region, const std::list& args, + const std::map& UNUSED(keywords)) { + // If not implemented, call two-argument version + return clone(region, args); + } virtual BoundaryOpPar* - clone(BoundaryRegionPar* region, const std::list& args, + clone(BoundaryOpPar* region, const std::list& args, const std::map& UNUSED(keywords)) { // If not implemented, call two-argument version return clone(region, args); } - BoundaryRegionPar* bndry{nullptr}; + bout::boundary::BoundaryRegionFCI* bndry{nullptr}; protected: /// Possible ways to get boundary values @@ -49,7 +68,7 @@ protected: enum class ValueType { GEN, FIELD, REAL }; const ValueType value_type{ValueType::REAL}; - BoutReal getValue(const BoundaryRegionPar& bndry, BoutReal t); + BoutReal getValue(const bout::boundary::BoundaryRegionIterFCI& bndry, BoutReal t); }; template @@ -59,8 +78,23 @@ public: using BoundaryOpPar::clone; - // Note: All methods must implement clone, except for modifiers (see below) - BoundaryOpPar* clone(BoundaryRegionPar* region, + BoundaryOpPar* clone(BoundaryOpPar* region, + const std::list& args) override { + if (!args.empty()) { + try { + real_value = stringToReal(args.front()); + return new T(region, real_value); + } catch (const BoutException&) { + std::shared_ptr newgen = nullptr; + // First argument should be an expression + newgen = FieldFactory::get()->parse(args.front()); + return new T(region, newgen); + } + } + + return new T(region); + } + BoundaryOpPar* clone(bout::boundary::BoundaryRegionFCI* region, const std::list& args) override { if (!args.empty()) { try { @@ -77,7 +111,10 @@ public: return new T(region); } - BoundaryOpPar* clone(BoundaryRegionPar* region, Field3D* f) override { + BoundaryOpPar* clone(bout::boundary::BoundaryRegionFCI* region, Field3D* f) override { + return new T(region, f); + } + BoundaryOpPar* clone(BoundaryOpPar* region, Field3D* f) override { return new T(region, f); } @@ -91,16 +128,17 @@ public: void apply(Field3D& f) override { return apply(f, 0); } void apply(Field3D& f, BoutReal t) override { - f.ynext(bndry->dir).allocate(); // Ensure unique before modifying + f.ynext(bndry->dir()).allocate(); // Ensure unique before modifying 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 +149,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 bout::boundary::BoundaryRegionIterFCI& 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 bout::boundary::BoundaryRegionIterFCI& 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 bout::boundary::BoundaryRegionIterFCI& pnt, + BoutReal value) { + pnt.dirichlet_o3(f, value); } }; @@ -136,8 +177,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 bout::boundary::BoundaryRegionIterFCI& pnt, + BoutReal value) { + pnt.neumann_o1(f, value); } }; @@ -145,8 +187,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 bout::boundary::BoundaryRegionIterFCI& pnt, + BoutReal value) { + pnt.neumann_o2(f, value); } }; @@ -154,8 +197,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 bout::boundary::BoundaryRegionIterFCI& 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 4d5278d00f..b078d66ef6 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -1,10 +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 @@ -14,99 +26,105 @@ * */ -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 +BOUT_ENUM_CLASS(SheathLimitMode, limit_free, exponential_free, linear_free); -class BoundaryRegionPar : public BoundaryRegionBase { +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; - }; - - using IndicesVec = std::vector; - using IndicesIter = IndicesVec::iterator; +struct RealPoint { + BoutReal s_x; + BoutReal s_y; + BoutReal s_z; +}; - /// Vector of points in the boundary - IndicesVec bndry_points; - /// Current position in the boundary points - IndicesIter bndry_position; +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(intersection), length(length), valid(valid), + offset(offset), abs_offset(abs_offset) {}; +}; -public: - BoundaryRegionPar(const std::string& name, int dir, Mesh* passmesh) - : BoundaryRegionBase(name, passmesh), dir(dir) { - ASSERT0(std::abs(dir) == 1); - BoundaryRegionBase::isParallel = true; +inline bool operator<(const Indices& lhs, const Indices& rhs) { + return lhs.index < rhs.index; +} +inline bool operator<(const Indices& lhs, const Ind3D& rhs) { return lhs.index < rhs; } + +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 } - 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); + if (fm < 1e-10) { + return fc; // Low / no density condition } - /// 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}); + 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; } - 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}); + +#if CHECKLEVEL >= 2 + if (!std::isfinite(fp)) { + throw BoutException("SheathBoundary limitFree: {}, {} -> {}", fm, fc, fp); } +#endif - // 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)); } + return fp; +} + +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; +} + +template +class BoundaryRegionParIterBase { + +public: + 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; } @@ -115,37 +133,99 @@ 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; } + 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_sheath_o1(f); + } + return ythis(f) * (1 + length()) - yprev(f) * length(); + } + BoutReal + extrapolate_sheath_o1(const std::function& f) const { + return ythis(f); + } + BoutReal + extrapolate_sheath_o2(const std::function& f) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return extrapolate_sheath_o1(f); + } + return ythis(f) * (1 + length()) - yprev(f) * length(); } - // 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 interpolate_sheath_o2(const Field3D& f) const { + return ythis(f) * (1 - length()) + ynext(f) * length(); + } + BoutReal + interpolate_sheath_o2(const std::function& f) const { + return ythis(f) * (1 - length()) + ynext(f) * length(); + } + + 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_o1(f); + return extrapolate_next_o1(f); } - return f[ind()] * (1 + length()) - f.ynext(-dir)[ind().yp(-dir)] * length(); + return ythis(f) * 2 - yprev(f); } + BoutReal + extrapolate_next_o1(const std::function& f) const { + return ythis(f); + } + BoutReal + extrapolate_next_o2(const std::function& f) const { + ASSERT3(valid() >= 0); + if (valid() < 1) { + return extrapolate_next_o1(f); + } + return ythis(f) * 2 - yprev(f); + } + + // extrapolate the gradient into the boundary + 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); + } + return ythis(f) - ynext(f); + } + + BoundaryRegionParIterBase& operator*() { return *this; } + + BoundaryRegionParIterBase& operator++() { + ++bndry_position; + return *this; + } + + bool operator!=(const BoundaryRegionParIterBase& rhs) { + 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 { @@ -154,17 +234,34 @@ 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); + } } } + void limit_at_least(Field3D& f, BoutReal value) const { + ITER() { + if (getAt(f, i) < value) { + getAt(f, i) = value; + } + } + } + + 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. - 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 * (i + 1); } + } // NB: value needs to be scaled by dy void neumann_o2(Field3D& f, BoutReal value) const { @@ -172,34 +269,229 @@ public: if (valid() < 1) { return neumann_o1(f, value); } - ynext(f) = 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), + 2, yprev(f)); + } + } + + // 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; + } + } + + 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 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; + } + } + + template + BoutReal& getAt(Field3D& f, int off) const { + ASSERT3(f.hasParallelSlices()); + if constexpr (check) { + ASSERT3(valid() > -off - 2); + } + auto _off = offset() + off * _dir; + return f.ynext(_off)[ind().yp(_off)]; + } + template + const BoutReal& getAt(const Field3D& f, int off) const { + ASSERT3(f.hasParallelSlices()); + if constexpr (check) { + ASSERT3(valid() > -off - 2); } - ynext(f) = - parallel_stencil::neumann_o3(1 - length(), value, 1, f[ind()], 2, yprev(f)); + auto _off = offset() + off * _dir; + return f.ynext(_off)[ind().yp(_off)]; } - const int dir; + 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; + } + } + +#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; + int _dir; + +public: + int dir() const { return _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, signed char offset) { + if (!bndry_points.empty() && bndry_points.back().index > ind) { + is_sorted = false; + } + 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) { + add_point(xyz2ind(ix, iy, iz, localmesh), x, y, z, length, valid, offset); + } + + // 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 { + 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, + const bout::parallel_boundary_region::Indices& i2) { + return i1.index < i2.index; + }); + } + + 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.begin(), bndry_points.end(), + [&i2](auto i1) { return i1.index == i2; }); + } + + // 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); + } + + int dir() const { return _dir; } + +private: + int _dir; + /// 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; - // 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)]; } static Ind3D xyz2ind(int x, int y, int z, Mesh* mesh) { const int ny = mesh->LocalNy; const int nz = mesh->LocalNz; return Ind3D{(x * ny + y) * nz + z, ny, nz}; } + bool is_sorted{true}; }; #endif // BOUT_PAR_BNDRY_H diff --git a/include/bout/sys/parallel_stencils.hxx b/include/bout/sys/parallel_stencils.hxx new file mode 100644 index 0000000000..651ebc3ca1 --- /dev/null +++ b/include/bout/sys/parallel_stencils.hxx @@ -0,0 +1,34 @@ +#pragma once +#include +#include + +namespace bout::parallel_stencil { +// generated by src/mesh/parallel_boundary_stencil.cxx.py +using std::pow; + +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([[maybe_unused]] BoutReal 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 bout::parallel_stencil diff --git a/include/bout/sys/uuid.h b/include/bout/sys/uuid.h index d3d954e2fe..973103ef70 100644 --- a/include/bout/sys/uuid.h +++ b/include/bout/sys/uuid.h @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -75,7 +76,7 @@ namespace uuids { namespace detail { template -constexpr inline unsigned char hex2char(TChar const ch) { +constexpr inline unsigned char hex2char(const TChar ch) { if (ch >= static_cast('0') && ch <= static_cast('9')) return ch - static_cast('0'); if (ch >= static_cast('a') && ch <= static_cast('f')) @@ -86,14 +87,14 @@ constexpr inline unsigned char hex2char(TChar const ch) { } template -constexpr inline bool is_hex(TChar const ch) { +constexpr inline bool is_hex(const TChar ch) { return (ch >= static_cast('0') && ch <= static_cast('9')) || (ch >= static_cast('a') && ch <= static_cast('f')) || (ch >= static_cast('A') && ch <= static_cast('F')); } template -constexpr inline unsigned char hexpair2char(TChar const a, TChar const b) { +constexpr inline unsigned char hexpair2char(const TChar a, const TChar b) { return (hex2char(a) << 4) | hex2char(b); } @@ -104,7 +105,7 @@ class sha1 { static constexpr unsigned int block_bytes = 64; - inline static uint32_t left_rotate(uint32_t value, size_t const count) { + inline static uint32_t left_rotate(uint32_t value, const size_t count) { return (value << count) ^ (value >> (32 - count)); } @@ -129,7 +130,7 @@ class sha1 { } } - void process_block(void const* const start, void const* const end) { + void process_block(const void* const start, const void* const end) { const uint8_t* begin = static_cast(start); const uint8_t* finish = static_cast(end); while (begin != finish) { @@ -138,13 +139,13 @@ class sha1 { } } - void process_bytes(void const* const data, size_t const len) { + void process_bytes(const void* const data, const size_t len) { const uint8_t* block = static_cast(data); process_block(block, block + len); } - uint32_t const* get_digest(digest32_t digest) { - size_t const bitCount = this->m_byteCount * 8; + const uint32_t* get_digest(digest32_t digest) { + const size_t bitCount = this->m_byteCount * 8; process_byte(0x80); if (this->m_blockByteIndex > 56) { while (m_blockByteIndex != 0) { @@ -165,13 +166,13 @@ class sha1 { process_byte(static_cast((bitCount >> 24) & 0xFF)); process_byte(static_cast((bitCount >> 16) & 0xFF)); process_byte(static_cast((bitCount >> 8) & 0xFF)); - process_byte(static_cast((bitCount)&0xFF)); + process_byte(static_cast((bitCount) & 0xFF)); memcpy(digest, m_digest, 5 * sizeof(uint32_t)); return digest; } - uint8_t const* get_digest_bytes(digest8_t digest) { + const uint8_t* get_digest_bytes(digest8_t digest) { digest32_t d32; get_digest(d32); size_t di = 0; @@ -343,13 +344,13 @@ class uuid { public: using value_type = uint8_t; - constexpr uuid() noexcept : data({}){}; + constexpr uuid() noexcept : data({}) {}; uuid(value_type (&arr)[16]) noexcept { std::copy(std::cbegin(arr), std::cend(arr), std::begin(data)); } - uuid(std::array const& arr) noexcept { + uuid(const std::array& arr) noexcept { std::copy(std::cbegin(arr), std::cend(arr), std::begin(data)); } @@ -394,12 +395,12 @@ class uuid { void swap(uuid& other) noexcept { data.swap(other.data); } - inline char const* as_bytes() const { - return reinterpret_cast(data.data()); + inline const char* as_bytes() const { + return reinterpret_cast(data.data()); } template - static bool is_valid_uuid(CharT const* str) noexcept { + static bool is_valid_uuid(const CharT* str) noexcept { bool firstDigit = true; int hasBraces = 0; size_t index = 0; @@ -454,12 +455,12 @@ class uuid { template , class Allocator = std::allocator> static bool - is_valid_uuid(std::basic_string const& str) noexcept { + is_valid_uuid(const std::basic_string& str) noexcept { return is_valid_uuid(str.c_str()); } template - static uuid from_string(CharT const* str) noexcept { + static uuid from_string(const CharT* str) noexcept { CharT digit = 0; bool firstDigit = true; int hasBraces = 0; @@ -511,40 +512,40 @@ class uuid { template , class Allocator = std::allocator> static uuid - from_string(std::basic_string const& str) noexcept { + from_string(const std::basic_string& str) noexcept { return from_string(str.c_str()); } private: std::array data{{0}}; - friend bool operator==(uuid const& lhs, uuid const& rhs) noexcept; - friend bool operator<(uuid const& lhs, uuid const& rhs) noexcept; + friend bool operator==(const uuid& lhs, const uuid& rhs) noexcept; + friend bool operator<(const uuid& lhs, const uuid& rhs) noexcept; template friend std::basic_ostream& operator<<(std::basic_ostream& s, - uuid const& id); + const uuid& id); }; // -------------------------------------------------------------------------------------------------------------------------- // operators and non-member functions // -------------------------------------------------------------------------------------------------------------------------- -inline bool operator==(uuid const& lhs, uuid const& rhs) noexcept { +inline bool operator==(const uuid& lhs, const uuid& rhs) noexcept { return lhs.data == rhs.data; } -inline bool operator!=(uuid const& lhs, uuid const& rhs) noexcept { +inline bool operator!=(const uuid& lhs, const uuid& rhs) noexcept { return !(lhs == rhs); } -inline bool operator<(uuid const& lhs, uuid const& rhs) noexcept { +inline bool operator<(const uuid& lhs, const uuid& rhs) noexcept { return lhs.data < rhs.data; } template std::basic_ostream& operator<<(std::basic_ostream& s, - uuid const& id) { + const uuid& id) { // save current flags std::ios_base::fmtflags f(s.flags()); @@ -567,7 +568,7 @@ std::basic_ostream& operator<<(std::basic_ostream& s template , class Allocator = std::allocator> -inline std::basic_string to_string(uuid const& id) { +inline std::basic_string to_string(const uuid& id) { std::basic_stringstream sstr; sstr << id; return sstr.str(); @@ -691,11 +692,11 @@ using uuid_random_generator = basic_uuid_random_generator; class uuid_name_generator { public: - explicit uuid_name_generator(uuid const& namespace_uuid) noexcept + explicit uuid_name_generator(const uuid& namespace_uuid) noexcept : nsuuid(namespace_uuid) {} template - uuid operator()(CharT const* name) { + uuid operator()(const CharT* name) { size_t size = 0; if (std::is_same::value) size = strlen(name); @@ -709,7 +710,7 @@ class uuid_name_generator { template , class Allocator = std::allocator> - uuid operator()(std::basic_string const& name) { + uuid operator()(const std::basic_string& name) { reset(); process_characters(name.data(), name.size()); return make_uuid(); @@ -726,7 +727,7 @@ class uuid_name_generator { template ::value>> - void process_characters(char_type const* const characters, size_t const count) { + void process_characters(const char_type* const characters, const size_t count) { for (size_t i = 0; i < count; i++) { uint32_t c = characters[i]; hasher.process_byte(static_cast((c >> 0) & 0xFF)); @@ -736,7 +737,7 @@ class uuid_name_generator { } } - void process_characters(const char* const characters, size_t const count) { + void process_characters(const char* const characters, const size_t count) { hasher.process_bytes(characters, count); } @@ -846,7 +847,7 @@ struct hash { using argument_type = uuids::uuid; using result_type = std::size_t; - result_type operator()(argument_type const& uuid) const { + result_type operator()(const argument_type& uuid) const { std::hash hasher; return static_cast(hasher(uuids::to_string(uuid))); } diff --git a/include/bout/yboundary_regions.hxx b/include/bout/yboundary_regions.hxx new file mode 100644 index 0000000000..b1aec6bc2e --- /dev/null +++ b/include/bout/yboundary_regions.hxx @@ -0,0 +1,154 @@ +#pragma once + +#include + +#include "./boundary_iterator.hxx" +#include "bout/assert.hxx" +#include "bout/boutexception.hxx" +#include "bout/field_data.hxx" +#include "bout/globals.hxx" +#include "bout/mask.hxx" +#include "bout/options.hxx" +#include "bout/parallel_boundary_region.hxx" +#include "bout/region.hxx" +#include +#include + +/// 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: + template + void iter_regions(const Func& func) { + for (auto& region : boundary_regions) { + func(*region); + } + for (auto& region : boundary_regions_par) { + func(*region); + } + } + template + void iter_points(const Func& func) { + iter_regions([&](auto& region) { + for (auto& point : region) { + func(point); + } + }); + } + + template + void iter(const Func& func) { + return iter_points(func); + } + + YBoundary(YBndryType type, Options* options_ptr, const Mesh& mesh) { + bool lower_y = true; + bool upper_y = true; + bool outer_x = true; + bool inner_x = false; + if (options_ptr != nullptr) { + auto& options = *options_ptr; + if (!mesh.isFci()) { + 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); + } else { + outer_x = + options["outer_x"].doc("Boundary on outer x?").withDefault(outer_x); + inner_x = + options["inner_x"].doc("Boundary on inner x?").withDefault(inner_x); + } + } + switch (type) { + case YBndryType::sheath: + break; + case YBndryType::not_sheath: + lower_y = !lower_y; + upper_y = !upper_y; + outer_x = !outer_x; + inner_x = !inner_x; + break; + case YBndryType::all: + lower_y = true; + upper_y = true; + outer_x = true; + inner_x = true; + } + + 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())); + } + } + // Cache boundary regions + _contains.emplace_back(&mesh, false); + _contains.emplace_back(&mesh, false); + iter_points([&](const auto& point) { + if (point.dir() == 1) { + _contains[1][point.ind()] = true; + } else if (point.dir() == -1) { + _contains[0][point.ind()] = true; + } + }); + } + + bool contains_low(Ind3D ind) const { return _contains[0][ind]; } + bool contains_high(Ind3D ind) const { return _contains[1][ind]; } + template + 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]; + } + } + bool contains(int dir, Ind3D ind) const { + if (dir == 1) { + return contains<+1>(ind); + } + if (dir == -1) { + return contains<-1>(ind); + } + 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; +}; + +inline std::shared_ptr getYBoundary(Coordinates* coords, + YBndryType type = YBndryType::sheath) { + auto itype = static_cast(type); + if (coords->ybndrys[itype] == nullptr) { + coords->ybndrys[itype] = coords->makeYBoundary(type); + } + return coords->ybndrys[itype]; +} diff --git a/manual/sphinx/user_docs/boundary_options.rst b/manual/sphinx/user_docs/boundary_options.rst index 8493049516..1aba0321d2 100644 --- a/manual/sphinx/user_docs/boundary_options.rst +++ b/manual/sphinx/user_docs/boundary_options.rst @@ -436,6 +436,128 @@ 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) {} + + void rhs() { + BoutReal totalFlux = 0; + mesh->getCoordinates()->getYBoundary()->iter_points([&](auto& point) { + BoutReal flux = point.interpolate_sheath_o2(N) * point.interpolate_sheath_o2(V); + totalFlux += flux; + }); + } + + private: + const Field3D& N; + const Field3D& V; + }; + + + +There are several member functions of ``point``. ``point`` 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 ``point``, where ``f`` is a : + +.. list-table:: Members for boundary operation + :widths: 15 70 + :header-rows: 1 + + * - Function + - Description + * - ``point.ythis(f)`` + - Returns the value at the last point in the domain + * - ``point.ynext(f)`` + - Returns the value at the first point in the boundary, i.e. one beyond the domain. + * - ``point.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. + * - ``point.interpolate_sheath_o2(f)`` + - Returns the value at the boundary, assuming the bounday value has been set + * - ``point.extrapolate_sheath_o1(f)`` + - Returns the value at the boundary, extrapolating from the bulk, first order + * - ``point.extrapolate_sheath_o2(f)`` + - Returns the value at the boundary, extrapolating from the bulk, second order + * - ``point.extrapolate_next_o{1,2}(f)`` + - Extrapolate into the boundary from the bulk, first or second order + * - ``point.extrapolate_grad_o{1,2}(f)`` + - Extrapolate the gradient into the boundary, first or second order + * - ``point.dirichlet_o{1,2,3}(f, v)`` + - Apply dirichlet boundary conditions with value ``v`` and given order + * - ``point.neumann_o{1,2,3}(f, v)`` + - Applies a gradient of ``v / dy`` boundary condition. + * - ``point.limitFree(f)`` + - Extrapolate into the boundary using only monotonic decreasing values. + ``f`` needs to be positive. + * - ``point.dir()`` + - The direction of the boundary. + + + + Boundary regions ---------------- diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 9154f78a4f..c0c01fefc5 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -43,6 +43,7 @@ #include #include #include +#include #include #include #include @@ -445,14 +446,42 @@ void Field3D::applyTDerivBoundary() { } } -void Field3D::setBoundaryTo(const Field3D& f3d) { +void Field3D::setBoundaryTo(const Field3D& f3d, bool copyParallelSlices) { checkData(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_boundary_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()) { + if (isFci() && reg->by != 0) { + continue; + } /// Loop within each region for (reg->first(); !reg->isDone(); reg->next()) { for (int z = 0; z < nz; z++) { diff --git a/src/field/field_data.cxx b/src/field/field_data.cxx index b925cb1426..154380dcb6 100644 --- a/src/field/field_data.cxx +++ b/src/field/field_data.cxx @@ -3,6 +3,7 @@ #include "bout/parallel_boundary_region.hxx" #include "bout/unused.hxx" #include +#include #include #include #include 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 } diff --git a/src/mesh/boundary_factory.cxx b/src/mesh/boundary_factory.cxx index fd5bbdc70d..a9c278a4d6 100644 --- a/src/mesh/boundary_factory.cxx +++ b/src/mesh/boundary_factory.cxx @@ -1,6 +1,7 @@ #include "bout/parallel_boundary_op.hxx" #include "bout/parallel_boundary_region.hxx" #include +#include #include #include #include @@ -105,7 +106,8 @@ BoundaryOpBase* BoundaryFactory::create(const string& name, BoundaryRegionBase* // Clone the boundary operation, passing the region to operate over, // an empty args list and empty keyword map list args; - return pop->clone(dynamic_cast(region), args, {}); + return pop->clone(dynamic_cast(region), args, + {}); } else { // Perpendicular boundary BoundaryOp* op = findBoundaryOp(trim(name)); @@ -197,7 +199,8 @@ BoundaryOpBase* BoundaryFactory::create(const string& name, BoundaryRegionBase* BoundaryOpPar* pop = findBoundaryOpPar(trim(func)); if (pop != nullptr) { // An operation with arguments - return pop->clone(dynamic_cast(region), arglist, keywords); + return pop->clone(dynamic_cast(region), arglist, + keywords); } } else { // Perpendicular boundary diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index d2634d3f88..634db26f3a 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -4,20 +4,22 @@ * given the contravariant metric tensor terms **************************************************************************/ +#include "bout/field_data.hxx" #include #include #include #include -#include -#include -#include - #include #include +#include #include +#include #include +#include +#include +#include -#include +#include #include "invert3x3.hxx" #include "parallel/fci.hxx" @@ -369,7 +371,7 @@ Coordinates::Coordinates(Mesh* mesh, FieldMetric dx, FieldMetric dy, FieldMetric g_22(std::move(g_22)), g_33(std::move(g_33)), g_12(std::move(g_12)), g_13(std::move(g_13)), g_23(std::move(g_23)), ShiftTorsion(std::move(ShiftTorsion)), IntShiftTorsion(std::move(IntShiftTorsion)), nz(mesh->LocalNz), localmesh(mesh), - location(CELL_CENTRE) {} + localoptions(nullptr), location(CELL_CENTRE) {} Coordinates::Coordinates(Mesh* mesh, Options* options) : dx(1., mesh), dy(1., mesh), dz(1., mesh), d1_dx(mesh), d1_dy(mesh), d1_dz(mesh), @@ -381,7 +383,8 @@ Coordinates::Coordinates(Mesh* mesh, Options* options) G1_13(mesh), G1_23(mesh), G2_11(mesh), G2_22(mesh), G2_33(mesh), G2_12(mesh), G2_13(mesh), G2_23(mesh), G3_11(mesh), G3_22(mesh), G3_33(mesh), G3_12(mesh), G3_13(mesh), G3_23(mesh), G1(mesh), G2(mesh), G3(mesh), ShiftTorsion(mesh), - IntShiftTorsion(mesh), localmesh(mesh), location(CELL_CENTRE) { + IntShiftTorsion(mesh), localmesh(mesh), localoptions(options), + location(CELL_CENTRE) { if (options == nullptr) { options = Options::getRoot()->getSection("mesh"); @@ -604,7 +607,7 @@ Coordinates::Coordinates(Mesh* mesh, Options* options, const CELL_LOC loc, G1_13(mesh), G1_23(mesh), G2_11(mesh), G2_22(mesh), G2_33(mesh), G2_12(mesh), G2_13(mesh), G2_23(mesh), G3_11(mesh), G3_22(mesh), G3_33(mesh), G3_12(mesh), G3_13(mesh), G3_23(mesh), G1(mesh), G2(mesh), G3(mesh), ShiftTorsion(mesh), - IntShiftTorsion(mesh), localmesh(mesh), location(loc) { + IntShiftTorsion(mesh), localmesh(mesh), localoptions(options), location(loc) { std::string suffix = getLocationSuffix(location); @@ -2007,3 +2010,7 @@ void Coordinates::checkContravariant() { } } } + +std::shared_ptr Coordinates::makeYBoundary(YBndryType type) const { + return std::make_shared(type, localoptions, *localmesh); +} diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index ea1b6a41b2..e713ca4f6c 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -28,6 +28,7 @@ #include #include +#include #include #include #include @@ -643,23 +644,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; @@ -2218,9 +2224,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); @@ -2928,6 +2934,11 @@ void BoutMesh::addBoundaryRegions() { } RangeIterator BoutMesh::iterateBndryLowerInnerY() const { +#if CHECK > 0 + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } +#endif int xs = 0; int xe = LocalNx - 1; @@ -2963,6 +2974,11 @@ RangeIterator BoutMesh::iterateBndryLowerInnerY() const { } RangeIterator BoutMesh::iterateBndryLowerOuterY() const { +#if CHECK > 0 + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } +#endif int xs = 0; int xe = LocalNx - 1; @@ -2997,6 +3013,12 @@ RangeIterator BoutMesh::iterateBndryLowerOuterY() const { } RangeIterator BoutMesh::iterateBndryLowerY() const { +#if CHECK > 0 + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } +#endif + int xs = 0; int xe = LocalNx - 1; if ((DDATA_INDEST >= 0) && (DDATA_XSPLIT > xstart)) { @@ -3026,6 +3048,12 @@ RangeIterator BoutMesh::iterateBndryLowerY() const { } RangeIterator BoutMesh::iterateBndryUpperInnerY() const { +#if CHECK > 0 + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } +#endif + int xs = 0; int xe = LocalNx - 1; @@ -3060,6 +3088,12 @@ RangeIterator BoutMesh::iterateBndryUpperInnerY() const { } RangeIterator BoutMesh::iterateBndryUpperOuterY() const { +#if CHECK > 0 + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } +#endif + int xs = 0; int xe = LocalNx - 1; @@ -3094,6 +3128,12 @@ RangeIterator BoutMesh::iterateBndryUpperOuterY() const { } RangeIterator BoutMesh::iterateBndryUpperY() const { +#if CHECK > 0 + if (this->isFci()) { + throw BoutException("FCI should never use this iterator"); + } +#endif + int xs = 0; int xe = LocalNx - 1; if ((UDATA_INDEST >= 0) && (UDATA_XSPLIT > xstart)) { @@ -3124,12 +3164,12 @@ RangeIterator BoutMesh::iterateBndryUpperY() const { std::vector BoutMesh::getBoundaries() { return boundary; } -std::vector> -BoutMesh::getBoundariesPar(BoundaryParType type) { +using bout::boundary::BoundaryRegionFCI; +std::vector> +BoutMesh::getBoundariesPar(BoundaryParType type) const { return par_boundary[static_cast(type)]; } - -void BoutMesh::addBoundaryPar(std::shared_ptr bndry, +void BoutMesh::addBoundaryPar(std::shared_ptr bndry, BoundaryParType type) { output_info << "Adding new parallel boundary: " << bndry->label << endl; switch (type) { diff --git a/src/mesh/impls/bout/boutmesh.hxx b/src/mesh/impls/bout/boutmesh.hxx index b42bf325b5..49d8b6112a 100644 --- a/src/mesh/impls/bout/boutmesh.hxx +++ b/src/mesh/impls/bout/boutmesh.hxx @@ -163,9 +163,9 @@ public: // Boundary regions std::vector getBoundaries() override; - std::vector> - getBoundariesPar(BoundaryParType type) override; - void addBoundaryPar(std::shared_ptr bndry, + std::vector> + getBoundariesPar(BoundaryParType type) const override; + void addBoundaryPar(std::shared_ptr bndry, BoundaryParType type) override; std::set getPossibleBoundaries() const override; @@ -410,7 +410,7 @@ protected: private: std::vector boundary; // Vector of boundary regions - std::array>, + std::array>, static_cast(BoundaryParType::SIZE)> par_boundary; // Vector of parallel boundary regions diff --git a/src/mesh/parallel/fci.cxx b/src/mesh/parallel/fci.cxx index fa965a335a..21fd4ecc32 100644 --- a/src/mesh/parallel/fci.cxx +++ b/src/mesh/parallel/fci.cxx @@ -39,6 +39,7 @@ #include "fci.hxx" #include "bout/assert.hxx" +#include "bout/boundary_region_iter.hxx" #include "bout/bout_types.hxx" #include "bout/boutexception.hxx" #include "bout/field2d.hxx" @@ -63,11 +64,12 @@ #include using namespace std::string_view_literals; +using bout::boundary::BoundaryRegionFCI; FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy, Options& options, int offset, - const std::shared_ptr& inner_boundary, - const std::shared_ptr& outer_boundary, bool zperiodic) + const std::shared_ptr& inner_boundary, + const std::shared_ptr& outer_boundary, bool zperiodic) : map_mesh(&mesh), offset_(offset), region_no_boundary(map_mesh->getRegion("RGN_NOBNDRY")), corner_boundary_mask(map_mesh) { @@ -183,7 +185,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; - // Serial loop because call to BoundaryRegionPar::addPoint + // Default to the maximum number of points + const int defValid{map_mesh->ystart - 1 + std::abs(offset)}; + // Serial loop because call to BoundaryRegionFCI::addPoint // (probably?) can't be done in parallel BOUT_FOR_SERIAL(i, xt_prime.getRegion("RGN_NOBNDRY")) { // z is periodic, so make sure the z-index wraps around @@ -252,11 +256,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 - (std::copysign(0.5, offset)), + 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); @@ -323,13 +328,13 @@ FCITransform::FCITransform(Mesh& mesh, const Coordinates::FieldMetric& dy, bool mesh.get(Z, "Z", 0.0, false); auto forward_boundary_xin = - std::make_shared("FCI_forward", BNDRY_PAR_FWD_XIN, +1, &mesh); + std::make_shared("FCI_forward", BNDRY_PAR_FWD_XIN, +1, &mesh); auto backward_boundary_xin = - std::make_shared("FCI_backward", BNDRY_PAR_BKWD_XIN, -1, &mesh); + std::make_shared("FCI_backward", BNDRY_PAR_BKWD_XIN, -1, &mesh); auto forward_boundary_xout = - std::make_shared("FCI_forward", BNDRY_PAR_FWD_XOUT, +1, &mesh); + std::make_shared("FCI_forward", BNDRY_PAR_FWD_XOUT, +1, &mesh); auto backward_boundary_xout = - std::make_shared("FCI_backward", BNDRY_PAR_BKWD_XOUT, -1, &mesh); + std::make_shared("FCI_backward", BNDRY_PAR_BKWD_XOUT, -1, &mesh); // Add the boundary region to the mesh's vector of parallel boundaries mesh.addBoundaryPar(forward_boundary_xin, BoundaryParType::xin_fwd); @@ -344,17 +349,22 @@ 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 (const auto& bndry2 : bndries) { - if (bndry->dir == bndry2->dir) { + if (bndry->dir() == bndry2->dir()) { continue; } - for (bndry->first(); !bndry->isDone(); bndry->next()) { - if (bndry2->contains(*bndry)) { - bndry->setValid(0); + for (auto point : *bndry) { + for (auto point2 : *bndry2) { + if (point.ind() == point2.ind()) { + // This point has a boundary in both directions. Calculate the + // distance between two points, to check how many non-boundary + // points exist. + point.setValid(static_cast( + std::abs((point2.offset() - point.offset())) - 2)); + } } } } diff --git a/src/mesh/parallel/fci.hxx b/src/mesh/parallel/fci.hxx index 65529a4c4e..f5d0642675 100644 --- a/src/mesh/parallel/fci.hxx +++ b/src/mesh/parallel/fci.hxx @@ -31,6 +31,7 @@ #include "bout/boutexception.hxx" #include "bout/coordinates.hxx" #include "bout/region.hxx" +#include #include #include #include @@ -41,7 +42,6 @@ #include #include -class BoundaryRegionPar; class FieldPerp; class Field2D; class Field3D; @@ -68,8 +68,9 @@ class FCIMap { public: FCIMap() = delete; FCIMap(Mesh& mesh, const Coordinates::FieldMetric& dy, Options& options, int offset, - const std::shared_ptr& inner_boundary, - const std::shared_ptr& outer_boundary, bool zperiodic); + const std::shared_ptr& inner_boundary, + const std::shared_ptr& outer_boundary, + bool zperiodic); /// Direction of map int offset() const { return offset_; } diff --git a/src/mesh/parallel/shiftedmetricinterp.cxx b/src/mesh/parallel/shiftedmetricinterp.cxx index c71618ab19..7171bdd0bd 100644 --- a/src/mesh/parallel/shiftedmetricinterp.cxx +++ b/src/mesh/parallel/shiftedmetricinterp.cxx @@ -31,10 +31,10 @@ #include "shiftedmetricinterp.hxx" +#include "bout/boundary_region_iter.hxx" #include "bout/boutexception.hxx" #include "bout/constants.hxx" #include "bout/field3d.hxx" -#include "bout/parallel_boundary_region.hxx" ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in, Field2D zShift_in, BoutReal zlength_in, @@ -133,7 +133,7 @@ ShiftedMetricInterp::ShiftedMetricInterp(Mesh& mesh, CELL_LOC location_in, // Create regions for parallel boundary conditions Field2D dy; mesh.get(dy, "dy", 1.); - auto forward_boundary_xin = std::make_shared( + auto forward_boundary_xin = std::make_shared( "parallel_forward_xin", BNDRY_PAR_FWD_XIN, +1, &mesh); for (auto it = mesh.iterateBndryUpperY(); not it.isDone(); it.next()) { for (int z = mesh.zstart; z <= mesh.zend; z++) { @@ -146,10 +146,10 @@ 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( + auto backward_boundary_xin = std::make_shared( "parallel_backward_xin", BNDRY_PAR_BKWD_XIN, -1, &mesh); for (auto it = mesh.iterateBndryLowerY(); not it.isDone(); it.next()) { for (int z = mesh.zstart; z <= mesh.zend; z++) { @@ -162,11 +162,11 @@ 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 - auto forward_boundary_xout = std::make_shared( + auto forward_boundary_xout = std::make_shared( "parallel_forward_xout", BNDRY_PAR_FWD_XOUT, +1, &mesh); for (auto it = mesh.iterateBndryUpperY(); not it.isDone(); it.next()) { for (int z = mesh.zstart; z <= mesh.zend; z++) { @@ -179,10 +179,10 @@ 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( + auto backward_boundary_xout = std::make_shared( "parallel_backward_xout", BNDRY_PAR_BKWD_XOUT, -1, &mesh); for (auto it = mesh.iterateBndryLowerY(); not it.isDone(); it.next()) { for (int z = mesh.zstart; z <= mesh.zend; z++) { @@ -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); } } @@ -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/src/mesh/parallel_boundary_op.cxx b/src/mesh/parallel_boundary_op.cxx index ebd9852791..bf0374ddf6 100644 --- a/src/mesh/parallel_boundary_op.cxx +++ b/src/mesh/parallel_boundary_op.cxx @@ -1,21 +1,22 @@ #include "bout/parallel_boundary_op.hxx" +#include "bout/boundary_region_iter.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 bout::boundary::BoundaryRegionIterFCI& 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)); + 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: diff --git a/src/mesh/parallel_boundary_stencil.cxx.py b/src/mesh/parallel_boundary_stencil.cxx.py index d0988ee099..adc163eece 100644 --- a/src/mesh/parallel_boundary_stencil.cxx.py +++ b/src/mesh/parallel_boundary_stencil.cxx.py @@ -26,18 +26,10 @@ def run(cmd): if __name__ == "__main__": with tmpf("w", dir=".", delete=False) as f: - f.write("namespace {\n") + f.write("namespace bout::parallel_stencil {\n") f.write( - """ -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; -} + """using std::pow; + """ ) diff --git a/tests/integrated/test-fci-boundary/get_par_bndry.cxx b/tests/integrated/test-fci-boundary/get_par_bndry.cxx index 8e3cfac2f7..1d5b9a5d63 100644 --- a/tests/integrated/test-fci-boundary/get_par_bndry.cxx +++ b/tests/integrated/test-fci-boundary/get_par_bndry.cxx @@ -6,6 +6,7 @@ #include "bout/options_io.hxx" #include "bout/output.hxx" #include "bout/parallel_boundary_region.hxx" +#include "bout/yboundary_regions.hxx" #include @@ -24,17 +25,38 @@ 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); 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; + } + + { + const Options dummy; + auto ybndry = getYBoundary(mesh->getCoordinates()); + + std::vector fields((mesh->ystart * 2) + 1, Field3D{0.0}); + for (auto& field : fields) { + field.allocate(); + } + ybndry->iter( + [&](const auto& pnt) { fields[pnt.dir() + mesh->ystart][pnt.ind()] += 1; }); + + 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..74a8abf741 100755 --- a/tests/integrated/test-fci-boundary/runtest +++ b/tests/integrated/test-fci-boundary/runtest @@ -4,9 +4,19 @@ 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 +39,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) diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index aa0220609b..6959b5e480 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -1,6 +1,7 @@ #pragma once #include +#include #include #include #include @@ -172,9 +173,9 @@ public: bool hasBndryUpperY() const override { return false; } void addBoundary(BoundaryRegion* region) override { boundaries.push_back(region); } std::vector getBoundaries() override { return boundaries; } - std::vector> - getBoundariesPar(BoundaryParType UNUSED(type)) override { - return std::vector>(); + std::vector> + getBoundariesPar(BoundaryParType UNUSED(type)) const override { + return std::vector>(); } BoutReal GlobalX(int jx) const override { return jx; } BoutReal GlobalY(int jy) const override { return jy; }