Fix #51

tests/test_notebooks.py

@@ -10,7 +10,8 @@
 testdir = tempfile.mkdtemp()
 
 def get_notebooks():
-    return [f for f in os.listdir(nbdir) if f.endswith('.ipynb')]
+    skip = ["benchmarking_besselaes.ipynb"]
+    return [f for f in os.listdir(nbdir) if f.endswith('.ipynb') and f not in skip]
 
 def get_jupyter_kernel():
     try:

timml/__init__.py

@@ -1,28 +1,38 @@
-'''
+"""
 Copyright (C), 2015, Mark Bakker.
 Mark Bakker, Delft University of Technology
 mark dot bakker at tudelft dot nl
 
 TimML is a computer program for the simulation of steady-state
 multiaquifer flow with analytic elements and consists of a
 library of Python scripts and FORTRAN extensions.
-'''
-#from __future__ import division, print_function, absolute_import
+"""
+# from __future__ import division, print_function, absolute_import
 
-#--version number
-__name__='timml'
-__author__='Mark Bakker'
+# --version number
+__name__ = "timml"
+__author__ = "Mark Bakker"
 from .version import __version__
 
 # Import all classes and functions
 from .model import ModelMaq, Model3D, Model
 from .well import WellBase, Well, HeadWell
 from .constant import Constant, ConstantStar
-from .linesink import LineSinkBase, HeadLineSinkZero, HeadLineSink, \
-                      LineSinkDitch, HeadLineSinkString, LineSinkDitchString, \
-                      HeadLineSinkContainer
-from .linedoublet import ImpLineDoublet, ImpLineDoubletString, \
-                         LeakyLineDoublet, LeakyLineDoubletString
+from .linesink import (
+    LineSinkBase,
+    HeadLineSinkZero,
+    HeadLineSink,
+    LineSinkDitch,
+    HeadLineSinkString,
+    LineSinkDitchString,
+    HeadLineSinkContainer,
+)
+from .linedoublet import (
+    ImpLineDoublet,
+    ImpLineDoubletString,
+    LeakyLineDoublet,
+    LeakyLineDoubletString,
+)
 from .circareasink import CircAreaSink
 from .inhomogeneity import PolygonInhomMaq, PolygonInhom3D, BuildingPit
 from .inhomogeneity1d import StripInhomMaq, StripInhom3D
@@ -31,5 +41,9 @@
 from .linesink1d import LineSink1D, HeadLineSink1D
 from .linedoublet1d import ImpLineDoublet1D, LeakyLineDoublet1D
 from .stripareasink import StripAreaSink
+from . import bessel
 
 __all__ = [s for s in dir() if not s.startswith("_")]
+
+# default bessel module is numba
+bessel.set_bessel_method(method="numba")

timml/bessel.py

@@ -0,0 +1,21 @@
+from importlib import import_module
+from warnings import warn
+
+
+def set_bessel_method(method="numba"):
+    global bessel
+    if method == "fortran":
+        try:
+            besselaesnew = import_module("timml.src.besselaesnew")
+            bessel = besselaesnew.besselaesnew
+            bessel.initialize()
+        except ImportError:
+            warn("Cannot import compiled fortran bessel module! Defaulting to numba!")
+            bessel = import_module("timml.besselaesnumba.besselaesnumba")
+    elif method == "numba":
+        bessel = import_module("timml.besselaesnumba.besselaesnumba")
+    else:
+        raise ValueError("method must be one of ['fortran', 'numba']")
+
+
+bessel = None  # is set in timml.__init__ or modified by set_bessel_method()

timml/linedoublet.py

@@ -2,16 +2,9 @@
 import matplotlib.pyplot as plt
 import inspect  # Used for storing the input
 from .element import Element
-from .besselaesnumba import besselaesnumba
-besselaesnumba.initialize()
-try:
-    from .src import besselaesnew
-    besselaesnew.besselaesnew.initialize()
-    #print('succes on f2py')
-except:
-    pass
 from .controlpoints import controlpoints
 from .equation import DisvecEquation, LeakyWallEquation
+from . import bessel
 
 __all__ = ['ImpLineDoublet', 'ImpLineDoubletString', 'LeakyLineDoublet', 
            'LeakyLineDoubletString']
@@ -35,10 +28,6 @@ def __init__(self, model, x1=-1, y1=0, x2=1, y2=0, delp=0.0, res=0.0,
         if addtomodel: self.model.add_element(self)
         self.aq = aq
         self.zcinout = zcinout
-        if self.model.f2py:
-            self.bessel = besselaesnew.besselaesnew
-        else:
-            self.bessel = besselaesnumba
 
     def __repr__(self):
         return self.name + ' from ' + str((self.x1, self.y1)) + ' to ' + str(
@@ -90,7 +79,7 @@ def potinf(self, x, y, aq=None):
             potrv = rv.reshape((self.order + 1, self.nlayers,
                                 aq.naq))  # clever way of using a reshaped rv here
             pot = np.zeros((self.order + 1, aq.naq))
-            pot[:, :] = self.bessel.potbesldv(float(x), float(y), self.z1, self.z2, aq.lab,
+            pot[:, :] = bessel.bessel.potbesldv(float(x), float(y), self.z1, self.z2, aq.lab,
                                               self.order, aq.ilap, aq.naq)
             potrv[:] = self.aq.coef[self.layers] * pot[:, np.newaxis, :]
         return rv
@@ -111,7 +100,7 @@ def disvecinf(self, x, y, aq=None):
             qxqyrv = rv.reshape((2, self.order + 1, self.nlayers,
                                  aq.naq))  # clever way of using a reshaped rv here
             qxqy = np.zeros((2 * (self.order + 1), aq.naq))
-            qxqy[:, :] = self.bessel.disbesldv(float(x), float(y), self.z1, self.z2, aq.lab,
+            qxqy[:, :] = bessel.bessel.disbesldv(float(x), float(y), self.z1, self.z2, aq.lab,
                                                self.order, aq.ilap, aq.naq)
             qxqyrv[0, :] = self.aq.coef[self.layers] * qxqy[:self.order + 1,
                                                        np.newaxis, :]

timml/linesink.py

@@ -3,20 +3,14 @@
 import inspect  # Used for storing the input
 from .element import Element
 from .equation import HeadEquation, PotentialEquation
-from .besselaesnumba import besselaesnumba
-besselaesnumba.initialize()
-try:
-    from .src import besselaesnew
-    besselaesnew.besselaesnew.initialize()
-    #print('succes on f2py')
-except:
-    pass
-
 from .controlpoints import controlpoints, strengthinf_controlpoints
+from . import bessel
+
 
 __all__ = ['LineSinkBase', 'HeadLineSinkZero', 'HeadLineSink', 'LineSinkDitch',
            'HeadLineSinkString', 'LineSinkDitchString']
 
+
 class LineSinkChangeTrace:
     def changetrace(self, xyzt1, xyzt2, aq, layer, ltype, modellayer, direction, hstepmax, verbose=False):
         changed = False
@@ -142,10 +136,6 @@ def __init__(self, model, x1=-1, y1=0, x2=1, y2=0, Qls=100.0, \
         self.wh = wh
         self.addtomodel = addtomodel
         if self.addtomodel: self.model.add_element(self)
-        if self.model.f2py:
-            self.bessel = besselaesnew.besselaesnew
-        else:
-            self.bessel = besselaesnumba
 
     def __repr__(self):
         return self.name + ' from ' + str((self.x1, self.y1)) + ' to ' + str(
@@ -179,7 +169,7 @@ def potinf(self, x, y, aq=None):
         rv = np.zeros((self.nparam, aq.naq))
         if aq == self.aq:
             pot = np.zeros(aq.naq)
-            pot[:] = self.bessel.potbeslsho(float(x), float(y), self.z1,
+            pot[:] = bessel.bessel.potbeslsho(float(x), float(y), self.z1,
                 self.z2, aq.lab, 0, aq.ilap, aq.naq)
             rv[:] = self.aq.coef[self.layers] * pot
         return rv
@@ -190,7 +180,7 @@ def disvecinf(self, x, y, aq=None):
         rv = np.zeros((2, self.nparam, aq.naq))
         if aq == self.aq:
             qxqy = np.zeros((2, aq.naq))
-            qxqy[:, :] = self.bessel.disbeslsho(float(x), float(y), self.z1,
+            qxqy[:, :] = bessel.bessel.disbeslsho(float(x), float(y), self.z1,
                 self.z2, aq.lab, 0, aq.ilap, aq.naq)
             rv[0] = self.aq.coef[self.layers] * qxqy[0]
             rv[1] = self.aq.coef[self.layers] * qxqy[1]
@@ -242,10 +232,6 @@ def __init__(self, model, x1=-1, y1=0, x2=1, y2=0, \
         if addtomodel: self.model.add_element(self)
         self.aq = aq
         self.zcinout = zcinout
-        if self.model.f2py:
-            self.bessel = besselaesnew.besselaesnew
-        else:
-            self.bessel = besselaesnumba
 
     def __repr__(self):
         return self.name + ' from ' + str((self.x1, self.y1)) + ' to ' + str(
@@ -299,7 +285,7 @@ def potinf(self, x, y, aq=None):
             # clever way of using a reshaped rv here
             potrv = rv.reshape((self.order + 1, self.nlayers, aq.naq))
             pot = np.zeros((self.order + 1, aq.naq))
-            pot[:, :] = self.bessel.potbeslsv(float(x), float(y), self.z1,
+            pot[:, :] = bessel.bessel.potbeslsv(float(x), float(y), self.z1,
                 self.z2, aq.lab, self.order, aq.ilap, aq.naq)
             potrv[:] = self.aq.coef[self.layers] * pot[:, np.newaxis, :]
         return rv
@@ -319,7 +305,7 @@ def disvecinf(self, x, y, aq=None):
         if aq == self.aq:
             qxqyrv = rv.reshape((2, self.order + 1, self.nlayers, aq.naq))
             qxqy = np.zeros((2 * (self.order + 1), aq.naq))
-            qxqy[:, :] = self.bessel.disbeslsv(float(x), float(y), self.z1,
+            qxqy[:, :] = bessel.bessel.disbeslsv(float(x), float(y), self.z1,
                 self.z2, aq.lab, self.order, aq.ilap, aq.naq)
             qxqyrv[0, :] = self.aq.coef[self.layers] * qxqy[:self.order + 1,
                                                        np.newaxis, :]

timml/model.py

@@ -41,21 +41,13 @@ class Model(PlotTim):
 
     """
 
-    def __init__(self, kaq, c, z, npor, ltype, f2py=False):
+    def __init__(self, kaq, c, z, npor, ltype):
         # All input variables are numpy arrays
         # That should be checked outside this function
         self.elementlist = []
         self.elementdict = {}  # only elements that have a label
         self.aq = Aquifer(self, kaq, c, z, npor, ltype)
         self.modelname = 'ml'  # Used for writing out input
-        self.f2py = False
-        if f2py:
-            try:
-                from .src import besselaesnew
-                self.f2py = True
-            except:
-                print('FORTRAN extension not found while f2py=True')
-                print('Using Numba instead')
 
     def initialize(self):
         # remove inhomogeneity elements (they are added again)
@@ -462,10 +454,10 @@ class ModelMaq(Model):
     """
 
     def __init__(self, kaq=1, z=[1, 0], c=[], npor=0.3, topboundary='conf',
-                 hstar=None, f2py=False):
+                 hstar=None):
         self.storeinput(inspect.currentframe())
         kaq, c, npor, ltype = param_maq(kaq, z, c, npor, topboundary)
-        Model.__init__(self, kaq, c, z, npor, ltype, f2py)
+        Model.__init__(self, kaq, c, z, npor, ltype)
         self.name = 'ModelMaq'
         if self.aq.ltype[0] == 'l':
             ConstantStar(self, hstar, aq=self.aq)
@@ -513,8 +505,7 @@ class Model3D(Model):
     """
 
     def __init__(self, kaq=1, z=[1, 0], kzoverkh=1, npor=0.3,
-                 topboundary='conf', topres=0, topthick=0, hstar=0,
-                 f2py=False):
+                 topboundary='conf', topres=0, topthick=0, hstar=0):
         '''Model3D
         for semi-confined aquifers, set top equal to 'semi' and provide
         topres: resistance of top
@@ -525,7 +516,7 @@ def __init__(self, kaq=1, z=[1, 0], kzoverkh=1, npor=0.3,
                                        topres)
         if topboundary == 'semi':
             z = np.hstack((z[0] + topthick, z))
-        Model.__init__(self, kaq, c, z, npor, ltype, f2py)
+        Model.__init__(self, kaq, c, z, npor, ltype)
         self.name = 'Model3D'
         if self.aq.ltype[0] == 'l':
             ConstantStar(self, hstar, aq=self.aq)