New Y-Boundaries for unified sheath code

.git-blame-ignore-revs

@@ -16,6 +16,7 @@ d1cfb8abd6aa5c76e6c1a4d7ab20929c65f8afc2
 8d5cb39e03c2644715a50684f8cd0318b4e82676
 ec69e8838be2dde140a915e50506f8e1ce3cb534
 f2bc0488a298f136294c523bc5ab4086d090059b
+c59626a0498b0cfd1be46f84f0fa38b4855a43cc
 1b4707187a3a85126338303dc766280b8fb2dc56
 b2e2f3575e68f771be2a4341af5fd14e2870469e
 64e4285ec38569f66d31e589a4610cefd16d8076

include/bout/boundary_iterator.hxx

@@ -0,0 +1,218 @@
+#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 <algorithm>
+#include <functional>
+
+class BoundaryRegionIter {
+public:
+  virtual ~BoundaryRegionIter() = default;
+  BoundaryRegionIter(int x, int y, int bx, int by, Mesh* mesh)
+      : dir(bx + by), x(x), y(y), bx(bx), by(by), localmesh(mesh) {
+    ASSERT3(bx * by == 0);
+  }
+  bool operator!=(const BoundaryRegionIter& rhs) 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;
+  }
+  virtual void _next() = 0;
+  BoundaryRegionIter& operator*() { return *this; }
+
+  void dirichlet_o2(Field3D& f, BoutReal value) const {
+    ynext(f) = parallel_stencil::dirichlet_o2(1, f[ind()], 0.5, value);
+  }
+
+  BoutReal extrapolate_grad_o2(const Field3D& f) const { return f[ind()] - yprev(f); }
+
+  BoutReal extrapolate_sheath_o2(const Field3D& f) const {
+    return (f[ind()] * 3 - yprev(f)) * 0.5;
+  }
+
+  BoutReal extrapolate_next_o2(const Field3D& f) const {
+    return (2 * f[ind()]) - yprev(f);
+  }
+
+  BoutReal
+  extrapolate_next_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& 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<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return (f(0, ind()) + f(0, ind().yp(-by).xp(-bx))) * 0.5;
+  }
+
+  BoutReal
+  extrapolate_sheath_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& 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; }
+
+#if BOUT_USE_METRIC_3D == 0
+  BoutReal& ynext(Field2D& f) const { return f[ind().yp(by).xp(bx)]; }
+  const BoutReal& ynext(const Field2D& f) const { return f[ind().yp(by).xp(bx)]; }
+  BoutReal& yprev(Field2D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+  const BoutReal& yprev(const Field2D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+#endif
+
+  const int dir;
+
+protected:
+  int z{0};
+  int x;
+  int y;
+  const int bx;
+  const int by;
+
+private:
+  Mesh* localmesh;
+  int nx() const { return localmesh->LocalNx; }
+  int ny() const { return localmesh->LocalNy; }
+  int nz() const { return localmesh->LocalNz; }
+
+  Ind3D xyz2ind(int x, int y, int z) const {
+    return Ind3D{((x * ny() + y) * nz()) + z, ny(), nz()};
+  }
+};
+
+class BoundaryRegionIterY : public BoundaryRegionIter {
+public:
+  ~BoundaryRegionIterY() override = default;
+  BoundaryRegionIterY(const RangeIterator& r, int y, int dir, bool is_end, Mesh* mesh)
+      : BoundaryRegionIter(r.ind, y, 0, dir, mesh), r(r), is_end(is_end) {}
+
+  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;
+    x = r.ind;
+  }
+
+private:
+  RangeIterator r;
+  bool is_end;
+};
+
+class NewBoundaryRegionY {
+public:
+  NewBoundaryRegionY(Mesh* mesh, bool lower, const RangeIterator& r)
+      : mesh(mesh), lower(lower), r(r) {}
+  BoundaryRegionIterY begin(bool begin = true) {
+    return BoundaryRegionIterY(r, lower ? mesh->ystart : mesh->yend, lower ? -1 : +1,
+                               !begin, mesh);
+  }
+  BoundaryRegionIterY end() { return begin(false); }
+
+private:
+  Mesh* mesh;
+  bool lower;
+  RangeIterator r;
+};

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;

include/bout/parallel_boundary_op.hxx

@@ -49,7 +49,7 @@ protected:
   enum class ValueType { GEN, FIELD, REAL };
   const ValueType value_type{ValueType::REAL};
 
-  BoutReal getValue(const BoundaryRegionPar& bndry, BoutReal t);
+  BoutReal getValue(const BoundaryRegionParIter& bndry, BoutReal t);
 };
 
 template <class T, bool isNeumann = false>
@@ -95,12 +95,13 @@ public:
 
     auto dy = f.getCoordinates()->dy;
 
-    for (bndry->first(); !bndry->isDone(); bndry->next()) {
-      BoutReal value = getValue(*bndry, t);
+    for (auto pnt : *bndry) {
+      //for (bndry->first(); !bndry->isDone(); bndry->next()) {
+      BoutReal value = getValue(pnt, t);
       if (isNeumann) {
-        value *= dy[bndry->ind()];
+        value *= dy[pnt.ind()];
       }
-      static_cast<T*>(this)->apply_stencil(f, bndry, value);
+      static_cast<T*>(this)->apply_stencil(f, pnt, value);
     }
   }
 };
@@ -111,51 +112,57 @@ public:
 class BoundaryOpPar_dirichlet_o1 : public BoundaryOpParTemp<BoundaryOpPar_dirichlet_o1> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->dirichlet_o1(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.dirichlet_o1(f, value);
   }
 };
 
 class BoundaryOpPar_dirichlet_o2 : public BoundaryOpParTemp<BoundaryOpPar_dirichlet_o2> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->dirichlet_o2(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.dirichlet_o2(f, value);
   }
 };
 
 class BoundaryOpPar_dirichlet_o3 : public BoundaryOpParTemp<BoundaryOpPar_dirichlet_o3> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->dirichlet_o3(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.dirichlet_o3(f, value);
   }
 };
 
 class BoundaryOpPar_neumann_o1
     : public BoundaryOpParTemp<BoundaryOpPar_neumann_o1, true> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->neumann_o1(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.neumann_o1(f, value);
   }
 };
 
 class BoundaryOpPar_neumann_o2
     : public BoundaryOpParTemp<BoundaryOpPar_neumann_o2, true> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->neumann_o2(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.neumann_o2(f, value);
   }
 };
 
 class BoundaryOpPar_neumann_o3
     : public BoundaryOpParTemp<BoundaryOpPar_neumann_o3, true> {
 public:
   using BoundaryOpParTemp::BoundaryOpParTemp;
-  static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
-    bndry->neumann_o3(f, value);
+  static void apply_stencil(Field3D& f, const BoundaryRegionParIter& pnt,
+                            BoutReal value) {
+    pnt.neumann_o3(f, value);
   }
 };