BoundaryOp refactor

examples/performance/apply-boundary/apply-boundary.cxx

@@ -0,0 +1,132 @@
+#include <bout/physicsmodel.hxx>
+#include <bout/sys/timer.hxx>
+
+class BoundaryOp_timing : public PhysicsModel {
+  int init(bool restarting);
+  int rhs(BoutReal UNUSED(t)) { return 1; }
+  Field3D f_dirichlet, f_neumann, f_dirichlet_o3;
+  Field3D f_dirichlet_val, f_neumann_val, f_dirichlet_o3_val;
+  Field3D f_dirichlet_expr, f_neumann_expr, f_dirichlet_o3_expr;
+};
+
+int BoundaryOp_timing::init(bool UNUSED(restarting)) {
+  SOLVE_FOR3(f_dirichlet, f_neumann, f_dirichlet_o3);
+  SOLVE_FOR3(f_dirichlet_val, f_neumann_val, f_dirichlet_o3_val);
+  SOLVE_FOR3(f_dirichlet_expr, f_neumann_expr, f_dirichlet_o3_expr);
+
+  Options* opt = Options::getRoot();
+  int ntests;
+  OPTION(opt, ntests, 10000);
+
+
+  // test Dirichlet_O3
+
+  {
+    Timer timer("dirichlet_o3");
+    for (int i=0; i<ntests; i++) {
+      f_dirichlet_o3.applyBoundary();
+    }
+  }
+
+
+  // test Dirichlet
+
+  {
+    Timer timer("dirichlet");
+    for (int i=0; i<ntests; i++) {
+      f_dirichlet.applyBoundary();
+    }
+  }
+
+
+  // test Neumann
+
+  {
+    Timer timer("neumann");
+    for (int i=0; i<ntests; i++) {
+      f_neumann.applyBoundary();
+    }
+  }
+
+
+  // test Dirichlet_O3
+
+  {
+    Timer timer("dirichlet_o3_val");
+    for (int i=0; i<ntests; i++) {
+      f_dirichlet_o3_val.applyBoundary();
+    }
+  }
+
+
+  // test Dirichlet
+
+  {
+    Timer timer("dirichlet_val");
+    for (int i=0; i<ntests; i++) {
+      f_dirichlet_val.applyBoundary();
+    }
+  }
+
+
+  // test Neumann
+
+  {
+    Timer timer("neumann_val");
+    for (int i=0; i<ntests; i++) {
+      f_neumann_val.applyBoundary();
+    }
+  }
+
+
+  // test Dirichlet_O3
+
+  {
+    Timer timer("dirichlet_o3_expr");
+    for (int i=0; i<ntests; i++) {
+      f_dirichlet_o3_expr.applyBoundary();
+    }
+  }
+
+
+  // test Dirichlet
+
+  {
+    Timer timer("dirichlet_expr");
+    for (int i=0; i<ntests; i++) {
+      f_dirichlet_expr.applyBoundary();
+    }
+  }
+
+
+  // test Neumann
+
+  {
+    Timer timer("neumann_expr");
+    for (int i=0; i<ntests; i++) {
+      f_neumann_expr.applyBoundary();
+    }
+  }
+
+
+  output<<"dirichlet:\t"<<ntests<<" iterations took "<<Timer::getTime("dirichlet")<<" s"<<endl;
+  output<<"neumann:\t"<<ntests<<" iterations took "<<Timer::getTime("neumann")<<" s"<<endl;
+  output<<"dirichlet_o3:\t"<<ntests<<" iterations took "<<Timer::getTime("dirichlet_o3")<<" s"<<endl;
+  output<<endl;
+
+  output<<"dirichlet_val:\t"<<ntests<<" iterations took "<<Timer::getTime("dirichlet_val")<<" s"<<endl;
+  output<<"neumann_val:\t"<<ntests<<" iterations took "<<Timer::getTime("neumann_val")<<" s"<<endl;
+  output<<"dirichlet_o3_val:\t"<<ntests<<" iterations took "<<Timer::getTime("dirichlet_o3_val")<<" s"<<endl;
+  output<<endl;
+
+  output<<"dirichlet_expr:\t"<<ntests<<" iterations took "<<Timer::getTime("dirichlet_expr")<<" s"<<endl;
+  output<<"neumann_expr:\t"<<ntests<<" iterations took "<<Timer::getTime("neumann_expr")<<" s"<<endl;
+  output<<"dirichlet_o3_expr:\t"<<ntests<<" iterations took "<<Timer::getTime("dirichlet_o3_expr")<<" s"<<endl;
+
+  BoutFinalise();
+  exit(0);
+
+  return 1;
+}
+
+BOUTMAIN(BoundaryOp_timing);

examples/performance/apply-boundary/data/BOUT.inp

@@ -0,0 +1,33 @@
+[mesh]
+nx = 128
+ny = 128
+nz = 128
+
+ixseps1 = -1
+
+[f_dirichlet]
+bndry_all = dirichlet
+
+[f_neumann]
+bndry_all = neumann
+
+[f_dirichlet_o3]
+bndry_all = dirichlet_o3
+
+[f_dirichlet_val]
+bndry_all = dirichlet(42.)
+
+[f_neumann_val]
+bndry_all = neumann(42.)
+
+[f_dirichlet_o3_val]
+bndry_all = dirichlet_o3(42.)
+
+[f_dirichlet_expr]
+bndry_all = dirichlet(x+y+z)
+
+[f_neumann_expr]
+bndry_all = neumann(x+y+z)
+
+[f_dirichlet_o3_expr]
+bndry_all = dirichlet_o3(x+y+z)

examples/performance/apply-boundary/makefile

@@ -0,0 +1,5 @@
+BOUT_TOP = ../../..
+
+SOURCEC = apply-boundary.cxx
+
+include $(BOUT_TOP)/make.config

include/boundary_factory.hxx

@@ -14,13 +14,20 @@ class BoundaryFactory;
 using std::string;
 using std::map;
 
+template<typename T>
+using BoundaryRegionOp = typename std::conditional<std::is_convertible<T*, BoundaryRegionPar*>::value, BoundaryOpPar, BoundaryOp>::type;
+
 /// 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
+ *
+ * Keywords can also be passed, e.g. "dirichlet(3., width=4)" will set a
+ * Dirichlet boundary condition to a constant value 3, and force 4 guard cells
+ * to be used (regardless of the values of xstart/ystart/xend/yend).
  * 
  * This is a singleton, so only one instance can exist. This is 
  * enforced by making the constructor private, and having a getInstance() method
@@ -55,7 +62,7 @@ using std::map;
  * 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));
+ * auto *bndry = bf->create("myboundary()", new BoundaryRegionXOut("xout", 0, 10, localmesh));
  * 
  * where the region is defined in boundary_region.hxx
  * 
@@ -69,12 +76,16 @@ class BoundaryFactory {
   static void cleanup(); ///< Frees all memory
 
   /// Create a boundary operation object
-  BoundaryOpBase* create(const string &name, BoundaryRegionBase *region);
-  BoundaryOpBase* create(const char* name, BoundaryRegionBase *region);
+  template<typename T>
+  BoundaryRegionOp<T>* create(const string &name, T* region);
+  template<typename T>
+  BoundaryRegionOp<T>* create(const char* name, T* region);
 
   /// Create a boundary object using the options file
-  BoundaryOpBase* createFromOptions(const string &varname, BoundaryRegionBase *region);
-  BoundaryOpBase* createFromOptions(const char* varname, BoundaryRegionBase *region);
+  template<typename T>
+  BoundaryRegionOp<T>* createFromOptions(const string &varname, T* region);
+  template<typename T>
+  BoundaryRegionOp<T>* createFromOptions(const char* varname, T* region);
 
   /*!
    * Add available boundary conditions and modifiers
@@ -119,14 +130,18 @@ class BoundaryFactory {
   // map<string, BoundaryModifier*> par_modmap;
 
   // Functions to look up operations and modifiers
-  BoundaryOp* findBoundaryOp(const string &s);
+  // Standard or parallel boundary conditions
+  template<typename T>
+  T* findBoundaryOp(const string &s);
+  // Boundary modifiers
   BoundaryModifier* findBoundaryMod(const string &s);
-  // Parallel boundary conditions
-  BoundaryOpPar* findBoundaryOpPar(const string &s);
-  // To be implemented...
-  // BoundaryModifier* findBoundaryMod(const string &s);
 
 };
 
+template<>
+BoundaryOp* BoundaryFactory::findBoundaryOp<BoundaryOp>(const string &s);
+template<>
+BoundaryOpPar* BoundaryFactory::findBoundaryOp<BoundaryOpPar>(const string &s);
+
 #endif // __BNDRY_FACTORY_H__
 

include/boundary_op.hxx

@@ -8,27 +8,39 @@ class BoundaryModifier;
 #include "boundary_region.hxx"
 #include "field2d.hxx"
 #include "field3d.hxx"
+#include "unused.hxx"
 #include "vector2d.hxx"
 #include "vector3d.hxx"
-#include "unused.hxx"
 
 #include <cmath>
-#include <string>
 #include <list>
 #include <map>
+#include <string>
 using std::string;
 using std::list;
 
-class BoundaryOpBase {
+/// An operation on a boundary
+class BoundaryOp {
 public:
-  BoundaryOpBase() {}
-  virtual ~BoundaryOpBase() {}
+  BoundaryOp(bool apply_ddt = false)
+      : bndry(nullptr), apply_to_ddt(apply_ddt), val(0.), gen(nullptr) {}
+  BoundaryOp(BoundaryRegion *region, bool apply_ddt = false)
+      : bndry(region), apply_to_ddt(apply_ddt), val(0.), gen(nullptr) {}
+  BoundaryOp(BoundaryRegion *region, BoutReal val_in, std::shared_ptr<FieldGenerator> g)
+      : bndry(region), apply_to_ddt(false), val(val_in), gen(std::move(g)) {}
+  virtual ~BoundaryOp() {}
+
+  // Note: All methods must implement clone, except for modifiers (see below)
+  virtual BoundaryOp *clone(BoundaryRegion *UNUSED(region),
+                            const list<string> &UNUSED(args),
+                            const std::map<std::string, std::string> &UNUSED(keywords)) {
+    throw BoutException("BoundaryOp::clone not implemented");
+    return nullptr;
+  }
 
   /// Apply a boundary condition on field f
-  virtual void apply(Field2D &f) = 0;
-  virtual void apply(Field2D &f,BoutReal UNUSED(t)){return apply(f);}//JMAD
-  virtual void apply(Field3D &f) = 0;
-  virtual void apply(Field3D &f,BoutReal UNUSED(t)){return apply(f);}//JMAD
+  virtual void apply(Field2D &f, BoutReal t = 0.) = 0;
+  virtual void apply(Field3D &f, BoutReal t = 0.) = 0;
 
   virtual void apply(Vector2D &f) {
     apply(f.x);
@@ -41,58 +53,54 @@ public:
     apply(f.y);
     apply(f.z);
   }
+
+  /// Apply a boundary condition on ddt(f)
+  virtual void apply_ddt(Field2D &f);
+  virtual void apply_ddt(Field3D &f);
+  virtual void apply_ddt(Vector2D &f) { apply(ddt(f)); }
+  virtual void apply_ddt(Vector3D &f) { apply(ddt(f)); }
+
+  BoundaryRegion *bndry;
+  const bool apply_to_ddt; // True if this boundary condition should be applied on the
+                           // time derivatives, false if it should be applied to the field
+                           // values
+protected:
+  const BoutReal val;                  // constant value for boundary condition
+  std::shared_ptr<FieldGenerator> gen; // Generator
 };
 
 /// An operation on a boundary
-class BoundaryOp : public BoundaryOpBase {
+template <typename Derived, bool needs_delta = false>
+class BoundaryOpWithApply : public BoundaryOp {
 public:
-  BoundaryOp() {
-    bndry = nullptr;
-    apply_to_ddt = false;
-  }
-  BoundaryOp(BoundaryRegion *region) {bndry = region; apply_to_ddt=false;}
-  ~BoundaryOp() override {}
+  using BoundaryOp::BoundaryOp;
 
   // Note: All methods must implement clone, except for modifiers (see below)
-  virtual BoundaryOp* clone(BoundaryRegion *UNUSED(region), const list<string> &UNUSED(args)) {
-    throw BoutException("BoundaryOp::clone not implemented");
-  }
-
-  /// Clone using positional args and keywords
-  /// If not implemented, check if keywords are passed, then call two-argument version
-  virtual BoundaryOp *clone(BoundaryRegion *region, const list<string> &args,
-                            const std::map<std::string, std::string> &keywords) {
-    if (!keywords.empty()) {
-      // Given keywords, but not using
-      throw BoutException("Keywords ignored in boundary : %s", keywords.begin()->first.c_str());
-    }
-
-    return clone(region, args);
+  virtual BoundaryOp *clone(BoundaryRegion *UNUSED(region),
+                            const list<string> &UNUSED(args)) {
+    ASSERT1(false); // this implementation should never get called
+    return nullptr;
   }
 
-  /// Apply a boundary condition on ddt(f)
-  virtual void apply_ddt(Field2D &f) {
-    apply(ddt(f));
-  }
-  virtual void apply_ddt(Field3D &f) {
-    apply(ddt(f));
-  }
-  virtual void apply_ddt(Vector2D &f) {
-    apply(ddt(f));
-  }
-  virtual void apply_ddt(Vector3D &f) {
-    apply(ddt(f));
-  }
+  /// Apply a boundary condition on field f
+  void apply(Field2D &f, BoutReal t = 0.) override { applyTemplate(f, t); }
+  void apply(Field3D &f, BoutReal t = 0.) override { applyTemplate(f, t); }
 
-  BoundaryRegion *bndry;
-  bool apply_to_ddt; // True if this boundary condition should be applied on the time derivatives, false if it should be applied to the field values
+private:
+  template <typename T> void applyTemplate(T &f, BoutReal t);
 };
 
 class BoundaryModifier : public BoundaryOp {
 public:
-  BoundaryModifier() : op(nullptr) {}
-  BoundaryModifier(BoundaryOp *operation) : BoundaryOp(operation->bndry), op(operation) {}
-  virtual BoundaryOp* cloneMod(BoundaryOp *op, const list<string> &args) = 0;
+  BoundaryModifier(bool apply_ddt = false) : BoundaryOp(apply_ddt), op(nullptr) {}
+  BoundaryModifier(BoundaryOp *operation, bool apply_ddt = false)
+      : BoundaryOp(operation->bndry, apply_ddt), op(operation) {}
+  virtual BoundaryOp *cloneMod(BoundaryOp *op, const list<string> &args) = 0;
+  virtual BoundaryOpPar *cloneMod(BoundaryOpPar *UNUSED(op),
+                                  const list<string> &UNUSED(args)) {
+    throw BoutException("BoundaryModifier should not be called on a BoundaryOpPar.");
+  }
+
 protected:
   BoundaryOp *op;
 };