Merge pull request #2892 from boutproject/boutpp-fixes

Fixes for the python interface

examples/CMakeLists.txt

@@ -11,6 +11,7 @@ add_subdirectory(backtrace)
 add_subdirectory(blob2d)
 add_subdirectory(blob2d-outerloop)
 add_subdirectory(blob2d-laplacexz)
+add_subdirectory(boutpp)
 add_subdirectory(boundary-conditions/advection)
 add_subdirectory(conducting-wall-mode)
 add_subdirectory(conduction)

examples/blob2d/blob2d.cxx

@@ -25,7 +25,6 @@ class Blob2D : public PhysicsModel {
   BoutReal rho_s;   ///< Bohm gyro radius
   BoutReal Omega_i; ///< Ion cyclotron frequency
   BoutReal c_s;     ///< Bohm sound speed
-  BoutReal n0;      ///< Reference density
 
   // Constants to calculate the parameters
   BoutReal Te0; ///< Isothermal temperature [eV]
@@ -61,7 +60,6 @@ class Blob2D : public PhysicsModel {
     m_i = options["m_i"].withDefault(2 * 1.667e-27);
     m_e = options["m_e"].withDefault(9.11e-31);
 
-    n0 = options["n0"].doc("Background density in cubic m").withDefault(1e19);
     D_vort = options["D_vort"].doc("Viscous diffusion coefficient").withDefault(0.0);
     D_n = options["D_n"].doc("Density diffusion coefficient").withDefault(0.0);
 

examples/blob2d/delta_0.25/BOUT.inp

@@ -87,8 +87,6 @@ flags = 49152  # set_rhs i.e. identity matrix in boundaries
 
 Te0 = 5    # Electron Temperature (eV)
 
-n0 = 2e+18  # Background plasma density (m^-3)
-
 compressible = false  # Compressibility?
 
 boussinesq = true  # Boussinesq approximation (no perturbed n in vorticity)

examples/blob2d/delta_1/BOUT.inp

@@ -87,8 +87,6 @@ flags = 49152  # set_rhs i.e. identity matrix in boundaries
 
 Te0 = 5    # Electron Temperature (eV)
 
-n0 = 2e+18  # Background plasma density (m^-3)
-
 compressible = false  # Compressibility?
 
 boussinesq = true  # Boussinesq approximation (no perturbed n in vorticity)

examples/blob2d/delta_10/BOUT.inp

@@ -87,8 +87,6 @@ flags = 49152  # set_rhs i.e. identity matrix in boundaries
 
 Te0 = 5    # Electron Temperature (eV)
 
-n0 = 2e+18  # Background plasma density (m^-3)
-
 compressible = false  # Compressibility?
 
 boussinesq = true  # Boussinesq approximation (no perturbed n in vorticity)

examples/boutpp/CMakeLists.txt

@@ -0,0 +1,10 @@
+cmake_minimum_required(VERSION 3.13)
+
+if (NOT TARGET bout++::bout++)
+  find_package(bout++ REQUIRED)
+endif()
+
+bout_copy_file(runexample)
+bout_copy_file(blob2d.py)
+bout_copy_file(simulation.py)
+bout_copy_file(data/BOUT.inp)

examples/boutpp/blob2d.py

@@ -24,7 +24,7 @@ def init(self, restart):
 
         self.phiSolver = bc.Laplacian()
 
-        options = bc.Options("model")
+        options = bc.Options.root("model")
         # Temperature in eV
         Te0 = options.get("Te0", 30)
         e = options.get("e", 1.602e-19)
@@ -70,12 +70,20 @@ def init(self, restart):
 
         # /************ Create a solver for potential ********/
 
+        opts_boussinesq = bc.Options.root("phiBoussinesq")
+        opts_non_boussinesq = bc.Options.root("phiSolver")
+
         if self.boussinesq:
             # BOUT.inp section "phiBoussinesq"
-            self.phiSolver = bc.Laplacian(bc.Options("phiBoussinesq"))
+            opts_used = opts_boussinesq
+            opts_unused = opts_non_boussinesq
         else:
             # BOUT.inp section "phiSolver"
-            self.phiSolver = bc.Laplacian(bc.Options("phiSolver"))
+            opts_used = opts_non_boussinesq
+            opts_unused = opts_boussinesq
+
+        self.phiSolver = bc.Laplacian(opts_used)
+        opts_unused.setConditionallyUsed()
 
         # Starting guess for first solve (if iterative)
         self.phi = bc.create3D("0")
@@ -165,8 +173,8 @@ def ensure_blob():
 
 # settings used by the core code
 
-NOUT = 50      # number of time-steps
-TIMESTEP = 50  # time between outputs [1/wci]
+nout = 50      # number of time-steps
+timestep = 50  # time between outputs [1/wci]
 
 
 MXG = 2      # Number of X guard cells
@@ -198,17 +206,17 @@ def ensure_blob():
 
 [mesh:ddz]
 
-first = FFT
-second = FFT
+first = C2
+second = C2
 upwind = W3
 
 ###################################################
 # Time-integration solver
 
 [solver]
 
-ATOL = 1.0e-10  # absolute tolerance
-RTOL = 1.0e-5   # relative tolerance
+atol = 1e-10  # absolute tolerance
+rtol = 1e-05   # relative tolerance
 mxstep = 10000  # Maximum internal steps per output
 
 ###################################################
@@ -221,22 +229,20 @@ def ensure_blob():
 
 fourth_order = true  # 4th order or 2nd order
 
-flags = 0  # inversion flags for phi
-             # 0  = Zero value
-             # 10 = Zero gradient AC inner & outer
-             # 15 = Zero gradient AC and DC
-             # 768 = Zero laplace inner & outer
+# 0  = Zero value
+# 10 = Zero gradient AC inner & outer
+# 15 = Zero gradient AC and DC
+# 768 = Zero laplace inner & outer
 
 [phiSolver:precon]  # Preconditioner (if pctype=user)
-filter     = 0.     # Must not filter solution
-flags      = 49152  # set_rhs i.e. identity matrix in boundaries
+filter = 0.0     # Must not filter solution
+flags = 49152  # set_rhs i.e. identity matrix in boundaries
 
 ###################################################
 # Electrostatic potential solver (Boussinesq)
 
 [phiBoussinesq]
 # By default type is tri (serial) or spt (parallel)
-flags = 0
 
 ##################################################
 # general settings for the model
@@ -245,14 +251,12 @@ def ensure_blob():
 
 Te0 = 5    # Electron Temperature (eV)
 
-n0 = 2e18  # Background plasma density (m^-3)
-
 compressible = false  # Compressibility?
 
 boussinesq = true  # Boussinesq approximation (no perturbed n in vorticity)
 
-D_vort = 1e-6  # Viscosity
-D_n = 1e-6    # Diffusion
+D_vort = 1e-06  # Viscosity
+D_n = 1e-06    # Diffusion
 
 R_c = 1.5  # Radius of curvature (m)
 
@@ -261,7 +265,7 @@ def ensure_blob():
 # These can be overridden for individual variables in
 # a section of that name.
 
-[All]
+[all]
 scale = 0.0 # default size of initial perturbations
 
 bndry_all = neumann # Zero-gradient on all boundaries
@@ -278,9 +282,8 @@ def ensure_blob():
 
 
 if __name__ == "__main__":
-    if "--create" in sys.argv:
-        sys.argv.remove("--create")
-        ensure_blob()
+    ensure_blob()
+
     bc.init("-d blob".split(" ") + sys.argv[1:])
 
     # Create an instance

examples/boutpp/data/BOUT.inp

@@ -0,0 +1,9 @@
+nout=10
+timestep=10
+
+[mesh]
+nx=160
+ny=1
+nz=n/n
+
+MYG=0

tools/pylib/_boutpp_build/bout_options.pxd

@@ -43,6 +43,7 @@ cdef extern from "bout/options.hxx":
         void get(string, double&, double)
         void get(string, bool&, bool)
         void cleanCache()
+        void setConditionallyUsed()
 
 
 cdef extern from "bout/optionsreader.hxx":

tools/pylib/_boutpp_build/boutcpp.pxd.jinja

@@ -148,10 +148,10 @@ cdef extern from "bout/physicsmodel.hxx":
 ctypedef void (*Method)(void *param, void *user_data)
 cdef extern from "helper.h":
     cppclass PythonModel(PhysicsModel):
-        int rhs(double t)
+        int rhs(double t) except +raise_bout_py_error
         void pyinit()
         void free()
-        void solve()
+        void solve() except +raise_bout_py_error
         Solver * getSolver()
         void set_rhs_func(PythonModelCallback*)
         void set_init_func(PythonModelCallback*)

tools/pylib/_boutpp_build/boutpp.pyx.jinja

@@ -583,9 +583,9 @@ cdef class {{ field.field_type }}:
 
 {% endfor %}
     def __dealloc__(self):
-        self.__boutpp_dealloc()
+        self._boutpp_dealloc()
 
-    def __boutpp_dealloc(self):
+    def _boutpp_dealloc(self):
         if self.isSelfOwned and self.cobj != NULL:
             del self.cobj
             self.cobj = NULL
@@ -645,9 +645,9 @@ cdef class {{ vec }}:
 
 
     def __dealloc__(self):
-        self.__boutpp_dealloc()
+        self._boutpp_dealloc()
 
-    def __boutpp_dealloc(self):
+    def _boutpp_dealloc(self):
         if self.isSelfOwned and self.cobj != NULL:
             del self.cobj
             self.cobj=NULL
@@ -742,9 +742,9 @@ cdef class Mesh:
         return msh
 
     def __dealloc__(self):
-        self.__boutpp_dealloc()
+        self._boutpp_dealloc()
 
-    def __boutpp_dealloc(self):
+    def _boutpp_dealloc(self):
         if self.cobj and self.isSelfOwned:
             del self.cobj
             self.cobj = NULL
@@ -850,9 +850,9 @@ cdef class Coordinates:
 {% endfor %}
 
     def __dealloc__(self):
-        self.__boutpp_dealloc()
+        self._boutpp_dealloc()
 
-    def __boutpp_dealloc(self):
+    def _boutpp_dealloc(self):
         if self.cobj and self.isSelfOwned:
             del self.cobj
             self.cobj = NULL
@@ -931,9 +931,9 @@ cdef class FieldFactory:
         checkInit()
         cobj=< c.FieldFactory*>0
     def __dealloc__(self):
-        self.__boutpp_dealloc()
+        self._boutpp_dealloc()
 
-    def __boutpp_dealloc(self):
+    def _boutpp_dealloc(self):
         if self.cobj != NULL:
             del self.cobj
             self.cobj = NULL
@@ -965,9 +965,9 @@ cdef class PythonModelCallback:
         self.cobj = new c.PythonModelCallback(callback, <void*>method)
 
     def __dealloc__(self):
-        self.__boutpp_dealloc()
+        self._boutpp_dealloc()
 
-    def __boutpp_dealloc(self):
+    def _boutpp_dealloc(self):
         if self.cobj:
             del self.cobj
             self.cobj = NULL
@@ -1037,12 +1037,12 @@ cdef class PhysicsModelBase(object):
         self.cmodel.set_init_func(self.callbackinit)
 
     def __dealloc__(self):
-        if hasattr(self, "__boutpp_dealloc"):
-            self.__boutpp_dealloc()
+        if hasattr(self, "_boutpp_dealloc"):
+            self._boutpp_dealloc()
         else:
-            PhysicsModelBase.__boutpp_dealloc(self)
+            PhysicsModelBase._boutpp_dealloc(self)
 
-    def __boutpp_dealloc(self):
+    def _boutpp_dealloc(self):
         if self.cmodel != <c.PythonModel *> 0:
             self.cmodel.free()
             del self.cmodel
@@ -1123,8 +1123,8 @@ class PhysicsModel(PhysicsModelBase):
     def __dealloc__(self):
         super(PhysicsModel,self).__dealloc__()
 
-    def __boutpp_dealloc(self):
-        super(PhysicsModel,self).__boutpp_dealloc()
+    def _boutpp_dealloc(self):
+        super(PhysicsModel,self)._boutpp_dealloc()
 
 cdef extern from "bout/bout.hxx":
     int BoutInitialise(int&, char **&) except +raise_bout_py_error
@@ -1204,13 +1204,14 @@ def finalise():
                   PythonModelCallback)
     for obj in objects:
         if isinstance(obj, ourClasses):
-            if hasattr(obj, "__boutpp_dealloc"):
-                obj.__boutpp_dealloc()
+            if hasattr(obj, "_boutpp_dealloc"):
+                obj._boutpp_dealloc()
             else:
                 for ourClass in ourClasses:
                     if isinstance(obj, ourClass):
-                         ourClass.__boutpp_dealloc(obj)
-                         break
+                        if hasattr(ourClass, "_boutpp_dealloc"):
+                            ourClass._boutpp_dealloc(obj)
+                            break
     del objects
     # Actually finalise
     if wasInit:
@@ -1715,10 +1716,19 @@ cdef class Options:
             opt.get(key, ret_str, default_)
             return ret_str.decode()
 
+    def setConditionallyUsed(self):
+        """Set the attribute "conditionally used" to be true for \p options
+        and all its children/sections, causing `Options::getUnused` to
+        assume those options have been used. This is useful to ignore
+        options when checking for typos etc.
+        """
+        cdef c.Options* opt = self.cobj
+        opt.setConditionallyUsed()
+
     def __dealloc__(self):
-        self.__boutpp_dealloc()
+        self._boutpp_dealloc()
 
-    def __boutpp_dealloc(self):
+    def _boutpp_dealloc(self):
         if self.isSelfOwned and self.cobj != NULL:
             del self.cobj
             self.cobj = NULL