Cleanup YBoundary

include/bout/boundary_common.hxx

@@ -0,0 +1,57 @@
+#pragma once
+
+#include <bout/bout_types.hxx>
+
+#include <cstdint>
+
+namespace bout {
+namespace boundary {
+/// 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::uint8_t { limited, exponential, linear };
+
+// 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,
+                               bout::boundary::BoundaryFreeExtrapolation mode) {
+  if ((fm < fc) && (mode == bout::boundary::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 bout::boundary::BoundaryFreeExtrapolation::limited:
+  case bout::boundary::BoundaryFreeExtrapolation::exponential:
+    fp = fc * fc / fm; // Exponential
+    break;
+  case bout::boundary::BoundaryFreeExtrapolation::linear:
+    fp = (2.0 * fc) - fm; // Linear
+    break;
+  }
+
+#if CHECKLEVEL >= 2
+  if (!std::isfinite(fp)) {
+    throw BoutException("SheathBoundary limitFree: {}, {} -> {}", fm, fc, fp);
+  }
+#endif
+
+  return fp;
+}
+
+} // namespace boundary
+} // namespace bout

include/bout/boundary_factory.hxx

@@ -13,17 +13,18 @@ class BoundaryRegion;
 class BoundaryModifier;
 
 #include <map>
+#include <memory>
 #include <string>
 
 /// Create BoundaryOp objects on demand
 /*!
  * This implements a simple string parser, used to match boundary condition
  * names like "dirichlet" with a BoundaryOp object.
- * 
+ *
  * Modifiers: Simple modifications of boundary conditions can be performed,
  * for example transforming the coordinate system
- * 
- * This is a singleton, so only one instance can exist. This is 
+ *
+ * This is a singleton, so only one instance can exist. This is
  * enforced by making the constructor private, and having a getInstance() method
  * to return a pointer to the only instance.
  *
@@ -33,7 +34,7 @@ class BoundaryModifier;
  * Boundaries are defined as classes which inherit from BoundaryOp
  * These define a clone() function which creates a new BoundaryOp,
  * given a list of arguments. See boundary_standard.hxx for examples.
- * 
+ *
  *     class MyBoundary : public BoundaryOp {
  *      public:
  *       BoundaryOp* clone(BoundaryRegion *region, const list<string> &args) {
@@ -43,23 +44,23 @@ class BoundaryModifier;
  *       void apply(Field2D &f);
  *       void apply(Field3D &f);
  *     };
- * 
+ *
  * The singleton instance of BoundaryFactory from getInstance():
- * 
+ *
  *     BoundaryFactory* bf = BoundaryFactory::getInstance();
  *
  * New boundary types can be added to the BoundaryFactory
  *
  *     bf->add(new MyBoundary, "myboundary");
- * 
+ *
  *
  * Subsequent calls to create() or createFromOptions() can make use
  * of the boundary type "myboundary".
  *
  * BoundaryOpBase *bndry = bf->create("myboundary()", new BoundaryRegionXOut("xout", 0, 10, localmesh));
- * 
+ *
  * where the region is defined in boundary_region.hxx
- * 
+ *
  */
 class BoundaryFactory {
 public:
@@ -70,13 +71,16 @@ public:
   static void cleanup(); ///< Frees all memory
 
   /// Create a boundary operation object
-  BoundaryOpBase* create(const std::string& name, BoundaryRegionBase* region);
-  BoundaryOpBase* create(const char* name, BoundaryRegionBase* region);
+  BoundaryOpBase* create(const std::string& name,
+                         const std::shared_ptr<BoundaryRegionBase>& region);
+  BoundaryOpBase* create(const char* name,
+                         const std::shared_ptr<BoundaryRegionBase>& region);
 
   /// Create a boundary object using the options file
   BoundaryOpBase* createFromOptions(const std::string& varname,
-                                    BoundaryRegionBase* region);
-  BoundaryOpBase* createFromOptions(const char* varname, BoundaryRegionBase* region);
+                                    const std::shared_ptr<BoundaryRegionBase>& region);
+  BoundaryOpBase* createFromOptions(const char* varname,
+                                    const std::shared_ptr<BoundaryRegionBase>& region);
 
   /*!
    * Add available boundary conditions and modifiers

include/bout/boundary_iterator.hxx

@@ -38,122 +38,6 @@ protected:
 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<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; }
-
-  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; }