diff --git a/ttim/aquifer.py b/ttim/aquifer.py
index a24101f..b99043e 100644
--- a/ttim/aquifer.py
+++ b/ttim/aquifer.py
@@ -22,21 +22,21 @@ def __init__(
     ):
         """Kzoverkh and model3d only need to be specified when model is model3d."""
         self.model = model
-        self.kaq = np.atleast_1d(kaq).astype("d")
-        self.z = np.atleast_1d(z).astype("d")
+        self.kaq = np.atleast_1d(kaq).astype(float)
+        self.z = np.atleast_1d(z).astype(float)
         self.naq = len(self.kaq)
         self.nlayers = len(self.z) - 1
-        self.Haq = np.atleast_1d(Haq).astype("d")
-        self.Hll = np.atleast_1d(Hll).astype("d")
+        self.Haq = np.atleast_1d(Haq).astype(float)
+        self.Hll = np.atleast_1d(Hll).astype(float)
         self.T = self.kaq * self.Haq
         self.Tcol = self.T.reshape(self.naq, 1)
-        self.c = np.atleast_1d(c).astype("d")
+        self.c = np.atleast_1d(c).astype(float)
         self.c[self.c > 1e100] = 1e100
-        self.Saq = np.atleast_1d(Saq).astype("d")
-        self.Sll = np.atleast_1d(Sll).astype("d")
+        self.Saq = np.atleast_1d(Saq).astype(float)
+        self.Sll = np.atleast_1d(Sll).astype(float)
         self.Sll[self.Sll < 1e-20] = 1e-20  # Cannot be zero
-        self.poraq = np.atleast_1d(poraq).astype("d")
-        self.porll = np.atleast_1d(porll).astype("d")
+        self.poraq = np.atleast_1d(poraq).astype(float)
+        self.porll = np.atleast_1d(porll).astype(float)
         self.ltype = np.atleast_1d(ltype)
         self.zaqtop = self.z[:-1][self.ltype == "a"]
         self.zaqbot = self.z[1:][self.ltype == "a"]
@@ -50,7 +50,7 @@ def __init__(
         self.phreatictop = phreatictop
         self.kzoverkh = kzoverkh
         if self.kzoverkh is not None:
-            self.kzoverkh = np.atleast_1d(self.kzoverkh).astype("d")
+            self.kzoverkh = np.atleast_1d(self.kzoverkh).astype(float)
             if len(self.kzoverkh) == 1:
                 self.kzoverkh = self.kzoverkh * np.ones(self.naq)
         self.model3d = model3d
@@ -96,10 +96,10 @@ def initialize(self):
                 self.kzoverkh[:-1] * self.kaq[:-1]
             ) + 0.5 * self.Haq[1:] / (self.kzoverkh[1:] * self.kaq[1:])
         #
-        self.eigval = np.zeros((self.naq, self.model.npval), "D")
-        self.lab = np.zeros((self.naq, self.model.npval), "D")
-        self.eigvec = np.zeros((self.naq, self.naq, self.model.npval), "D")
-        self.coef = np.zeros((self.naq, self.naq, self.model.npval), "D")
+        self.eigval = np.zeros((self.naq, self.model.npval), dtype=complex)
+        self.lab = np.zeros((self.naq, self.model.npval), dtype=complex)
+        self.eigvec = np.zeros((self.naq, self.naq, self.model.npval), dtype=complex)
+        self.coef = np.zeros((self.naq, self.naq, self.model.npval), dtype=complex)
         b = np.diag(np.ones(self.naq))
         for i in range(self.model.npval):
             w, v = self.compute_lab_eigvec(self.model.p[i])
diff --git a/ttim/aquifer_parameters.py b/ttim/aquifer_parameters.py
index fcbf257..bd66b65 100644
--- a/ttim/aquifer_parameters.py
+++ b/ttim/aquifer_parameters.py
@@ -13,13 +13,13 @@ def param_maq(
     phreatictop=False,
 ):
     # Computes the parameters for a TimModel from input for a maq model
-    z = np.atleast_1d(z).astype("d")
-    kaq = np.atleast_1d(kaq).astype("d")
-    Saq = np.atleast_1d(Saq).astype("d")
-    poraq = np.atleast_1d(poraq).astype("d")
-    c = np.atleast_1d(c).astype("d")
-    Sll = np.atleast_1d(Sll).astype("d")
-    porll = np.atleast_1d(porll).astype("d")
+    z = np.atleast_1d(z).astype(float)
+    kaq = np.atleast_1d(kaq).astype(float)
+    Saq = np.atleast_1d(Saq).astype(float)
+    poraq = np.atleast_1d(poraq).astype(float)
+    c = np.atleast_1d(c).astype(float)
+    Sll = np.atleast_1d(Sll).astype(float)
+    porll = np.atleast_1d(porll).astype(float)
     H = z[:-1] - z[1:]
     assert np.all(H >= 0), (
         "Error: Not all layer thicknesses are" + " non-negative" + str(H)
@@ -106,18 +106,18 @@ def param_3d(
     toppor=0.3,
 ):
     # Computes the parameters for a TimModel from input for a 3D model
-    kaq = np.atleast_1d(kaq).astype("d")
-    z = np.atleast_1d(z).astype("d")
+    kaq = np.atleast_1d(kaq).astype(float)
+    z = np.atleast_1d(z).astype(float)
     naq = len(z) - 1
     if len(kaq) == 1:
         kaq = kaq * np.ones(naq)
-    Saq = np.atleast_1d(Saq).astype("d")
+    Saq = np.atleast_1d(Saq).astype(float)
     if len(Saq) == 1:
         Saq = Saq * np.ones(naq)
-    kzoverkh = np.atleast_1d(kzoverkh).astype("d")
+    kzoverkh = np.atleast_1d(kzoverkh).astype(float)
     if len(kzoverkh) == 1:
         kzoverkh = kzoverkh * np.ones(naq)
-    poraq = np.atleast_1d(poraq).astype("d")
+    poraq = np.atleast_1d(poraq).astype(float)
     if len(poraq) == 1:
         poraq = poraq * np.ones(naq)
     Haq = z[:-1] - z[1:]
diff --git a/ttim/aquifernew.py b/ttim/aquifernew.py
deleted file mode 100644
index 641f08e..0000000
--- a/ttim/aquifernew.py
+++ /dev/null
@@ -1,209 +0,0 @@
-# flake8: noqa
-import inspect  # Used for storing the input
-
-import matplotlib.pyplot as plt
-import numpy as np
-
-
-class AquiferData:
-    def __init__(self, model, kaq, c, z, npor, ltype, Saq, Sll, phreatictop):
-        self.model = model
-        self.kaq = np.atleast_1d(kaq).astype("d")
-        self.Naq = len(kaq)
-        self.c = np.atleast_1d(c).astype("d")
-        self.c[self.c > 1e100] = 1e100
-        # Needed for tracing
-        self.z = np.atleast_1d(z).astype("d")
-        self.Hlayer = self.z[:-1] - self.z[1:]  # thickness of all layers
-        self.nlayers = len(self.z) - 1
-        self.npor = np.atleast_1d(npor)
-        self.ltype = np.atleast_1d(ltype)
-        #
-        self.layernumber = np.zeros(self.nlayers, dtype="int")
-        self.layernumber[self.ltype == "a"] = np.arange(self.Naq)
-        self.layernumber[self.ltype == "l"] = np.arange(self.nlayers - self.Naq)
-        if self.ltype[0] == "a":
-            # first leaky layer below first aquifer layer
-            self.layernumber[self.ltype == "l"] += 1
-        self.zaqtop = self.z[:-1][self.ltype == "a"]
-        self.zaqbot = self.z[1:][self.ltype == "a"]
-        self.Haq = self.zaqtop - self.zaqbot
-        self.T = self.kaq * self.Haq
-        self.Tcol = self.T[:, np.newaxis]
-        self.zlltop = self.z[:-1][self.ltype == "l"]
-        self.zllbot = self.z[1:][self.ltype == "l"]
-        if self.ltype[0] == "a":
-            self.zlltop = np.hstack((self.z[0], self.zlltop))
-            self.zllbot = np.hstack((self.z[0], self.zllbot))
-        self.Hll = self.zlltop - self.zllbot
-        self.nporaq = self.npor[self.ltype == "a"]
-        if self.ltype[0] == "a":
-            self.nporll = np.ones(len(self.npor[self.ltype == "l"]) + 1)
-            self.nporll[1:] = self.npor[self.ltype == "l"]
-        else:
-            self.nporll = self.npor[self.ltype == "l"]
-        # Storage coefs
-        self.Saq = np.atleast_1d(Saq).astype("d")
-        self.Sll = np.atleast_1d(Sll).astype("d")
-        self.Sll[self.Sll < 1e-20] = 1e-20  # Cannot be zero
-        # what to do with these?
-        self.phreatictop = phreatictop  # used in calibration
-        self.area = 1e200  # Smaller than default of ml.aq so that inhom is found
-
-    def __repr__(self):
-        return "Inhom T: " + str(self.T)
-
-    def initialize(self):
-        """Eigval[Naq,npval]: Array with eigenvalues lab[Naq,npval]: Array with lambda
-        values lab2[Naq,Nin,npint]: Array with lambda values reorganized per interval
-        eigvec[Naq,Naq,npval]: Array with eigenvector matrices coef[Naq,Naq,npval]:
-
-        Array with coefficients; coef[ilayers,:,np] are the coefficients if the element
-        is in ilayers belonging to Laplace parameter number np.
-        """
-        # Recompute T for when kaq is changed manually
-        self.T = self.kaq * self.Haq
-        self.Tcol = self.T.reshape(self.Naq, 1)
-        self.D = self.T / self.Saq
-        #
-        self.eigval = np.zeros((self.Naq, self.model.Np), "D")
-        self.lab = np.zeros((self.Naq, self.model.Np), "D")
-        self.eigvec = np.zeros((self.Naq, self.Naq, self.model.Np), "D")
-        self.coef = np.zeros((self.Naq, self.Naq, self.model.Np), "D")
-        b = np.diag(np.ones(self.Naq))
-        for i in range(self.model.Np):
-            w, v = self.compute_lab_eigvec(
-                self.model.p[i]
-            )  # Eigenvectors are columns of v
-            ## moved to compute_lab_eigvec routine
-            # index = np.argsort( abs(w) )[::-1]
-            # w = w[index]; v = v[:,index]
-            self.eigval[:, i] = w
-            self.eigvec[:, :, i] = v
-            self.coef[:, :, i] = np.linalg.solve(v, b).T
-        self.lab = 1.0 / np.sqrt(self.eigval)
-        self.lab2 = self.lab.copy()
-        self.lab2.shape = (self.Naq, self.model.Nin, self.model.Npin)
-        self.lababs = np.abs(self.lab2[:, :, 0])  # used to check distances
-
-    def compute_lab_eigvec(self, p, returnA=False, B=None):
-        sqrtpSc = np.sqrt(p * self.Sll * self.c)
-        a, b = np.zeros_like(sqrtpSc), np.zeros_like(sqrtpSc)
-        small = np.abs(sqrtpSc) < 200
-        a[small] = sqrtpSc[small] / np.tanh(sqrtpSc[small])
-        b[small] = sqrtpSc[small] / np.sinh(sqrtpSc[small])
-        a[~small] = sqrtpSc[~small] / (
-            (1.0 - np.exp(-2.0 * sqrtpSc[~small]))
-            / (1.0 + np.exp(-2.0 * sqrtpSc[~small]))
-        )
-        b[~small] = (
-            sqrtpSc[~small]
-            * 2.0
-            * np.exp(-sqrtpSc[~small])
-            / (1.0 - np.exp(-2.0 * sqrtpSc[~small]))
-        )
-        if self.ltype[0] == "l":
-            if abs(sqrtpSc[0]) < 200:
-                dzero = sqrtpSc[0] * np.tanh(sqrtpSc[0])
-            else:
-                dzero = sqrtpSc[0] * cmath_tanh(
-                    sqrtpSc[0]
-                )  # Bug in complex tanh in numpy
-        d0 = p / self.D
-        if B is not None:
-            d0 = d0 * B  # B is vector of load efficiency paramters
-        d0[:-1] += a[1:] / (self.c[1:] * self.T[:-1])
-        d0[1:] += a[1:] / (self.c[1:] * self.T[1:])
-        ## Need to make option 'lea' possible somehow
-        # if self.topboundary[:3] == 'lea':
-        #    d0[0] += dzero / ( self.c[0] * self.T[0] )
-        # elif self.topboundary[:3] == 'sem':
-        #    d0[0] += a[0] / ( self.c[0] * self.T[0] )
-        if self.ltype[0] == "l":
-            d0[0] += a[0] / (self.c[0] * self.T[0])
-
-        dm1 = -b[1:] / (self.c[1:] * self.T[:-1])
-        dp1 = -b[1:] / (self.c[1:] * self.T[1:])
-        A = np.diag(dm1, -1) + np.diag(d0, 0) + np.diag(dp1, 1)
-        if returnA:
-            return A
-        w, v = np.linalg.eig(A)
-        # sorting moved here
-        index = np.argsort(abs(w))[::-1]
-        w = w[index]
-        v = v[:, index]
-        return w, v
-
-    def headToPotential(self, h, layers):
-        return h * self.Tcol[layers]
-
-    def potentialToHead(self, pot, layers):
-        return pot / self.Tcol[layers]
-
-    def isInside(self, x, y):
-        print("Must overload AquiferData.isInside method")
-        return True
-
-    def inWhichLayer(self, z):
-        """Get layer given elevation z.
-
-        Returns -9999 if above top of system, +9999 if below bottom of system, negative
-        for in leaky layer.
-
-        leaky layer -n is on top of aquifer n
-        """
-        if z > self.zt[0]:
-            return -9999
-        for i in range(self.Naquifers - 1):
-            if z >= self.zb[i]:
-                return i
-            if z > self.zt[i + 1]:
-                return -i - 1
-        if z >= self.zb[self.Naquifers - 1]:
-            return self.Naquifers - 1
-        return +9999
-
-    def set_kaq(self, value, layer):
-        self.kaq[layer] = value
-
-    def set_Saq(self, value, layer):
-        self.Saq[layer] = value
-
-    def findlayer(self, z):
-        """Returns layer-number, layer-type and model-layer-number."""
-        if z > self.z[0]:
-            modellayer, ltype = -1, "above"
-            layernumber = None
-        elif z < self.z[-1]:
-            modellayer, ltype = len(self.layernumber), "below"
-            layernumber = None
-        else:
-            modellayer = np.argwhere((z <= self.z[:-1]) & (z >= self.z[1:]))[0, 0]
-            layernumber = self.layernumber[modellayer]
-            ltype = self.ltype[modellayer]
-        return layernumber, ltype, modellayer
-
-
-class Aquifer(AquiferData):
-    def __init__(self, model, kaq, Haq, c, Saq, Sll, topboundary, phreatictop):
-        # AquiferData.__init__(self, model, kaq, Haq, c, Saq, Sll, \
-        #                     topboundary, phreatictop)
-        AquiferData.__init__(self, model, kaq, c, z, npor, ltype, Saq, Sll, phreatictop)
-        self.inhomList = []
-        self.area = 1e300  # Needed to find smallest inhomogeneity
-
-    def __repr__(self):
-        return "Background Aquifer T: " + str(self.T)
-
-    def initialize(self):
-        AquiferData.initialize(self)
-        for inhom in self.inhomList:
-            inhom.initialize()
-
-    def find_aquifer_data(self, x, y):
-        rv = self
-        for aq in self.inhomList:
-            if aq.isInside(x, y):
-                if aq.area < rv.area:
-                    rv = aq
-        return rv
diff --git a/ttim/circareasink.py b/ttim/circareasink.py
index a787584..62a53c3 100644
--- a/ttim/circareasink.py
+++ b/ttim/circareasink.py
@@ -77,10 +77,12 @@ def potinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rv = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
+        rv = np.zeros(
+            (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+        )
         if aq == self.aq:
             r = np.sqrt((x - self.xc) ** 2 + (y - self.yc) ** 2)
-            # pot = np.zeros(self.model.npint, "D")
+            # pot = np.zeros(self.model.npint, dtype=complex)
             if r < self.R:
                 for i in range(self.aq.naq):
                     for j in range(self.model.nint):
@@ -103,10 +105,12 @@ def disvecinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        qx = np.zeros((self.nparam, aq.naq, self.model.npval), "D")
-        qy = np.zeros((self.nparam, aq.naq, self.model.npval), "D")
+        qx = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
+        qy = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
         if aq == self.aq:
-            qr = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
+            qr = np.zeros(
+                (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+            )
             r = np.sqrt((x - self.xc) ** 2 + (y - self.yc) ** 2)
             if r < self.R:
                 for i in range(self.aq.naq):
diff --git a/ttim/circinhom.py b/ttim/circinhom.py
index 71e7ecb..ac73af4 100644
--- a/ttim/circinhom.py
+++ b/ttim/circinhom.py
@@ -92,9 +92,9 @@ def __init__(self, Norder, Nterms):
         self.Nterms = Nterms + 1
         self.krange = np.arange(self.Nterms)
         self.minonek = (-np.ones(self.Nterms)) ** self.krange
-        self.hankeltot = np.ones((self.Norder, 2 * self.Nterms), "d")
-        self.muk = np.ones((self.Norder, self.Nterms), "d")
-        self.nuk = np.ones((self.Norder, self.Nterms), "d")
+        self.hankeltot = np.ones((self.Norder, 2 * self.Nterms), dtype=float)
+        self.muk = np.ones((self.Norder, self.Nterms), dtype=float)
+        self.nuk = np.ones((self.Norder, self.Nterms), dtype=float)
         for n in range(self.Norder):
             mu = 4.0 * n**2
             for k in range(1, self.Nterms):
@@ -114,7 +114,7 @@ def __init__(self, Norder, Nterms):
 
     def ivratio(self, rho, R, lab):
         lab = np.atleast_1d(lab)
-        rv = np.empty((self.Norder, len(lab)), "D")
+        rv = np.empty((self.Norder, len(lab)), dtype=complex)
         for k in range(len(lab)):
             top = np.sum(
                 self.minonek * self.hankelnk / (2.0 * rho / lab[k]) ** self.krange, 1
@@ -127,7 +127,7 @@ def ivratio(self, rho, R, lab):
 
     def kvratio(self, rho, R, lab):
         lab = np.atleast_1d(lab)
-        rv = np.empty((self.Norder, len(lab)), "D")
+        rv = np.empty((self.Norder, len(lab)), dtype=complex)
         for k in range(len(lab)):
             top = np.sum(self.hankelnk / (2.0 * rho / lab[k]) ** self.krange, 1)
             bot = np.sum(self.hankelnk / (2.0 * R / lab[k]) ** self.krange, 1)
@@ -136,7 +136,7 @@ def kvratio(self, rho, R, lab):
 
     def ivratiop(self, rho, R, lab):
         lab = np.atleast_1d(lab)
-        rv = np.empty((self.Norder, len(lab)), "D")
+        rv = np.empty((self.Norder, len(lab)), dtype=complex)
         for k in range(len(lab)):
             top = np.sum(
                 self.muk * self.hankeln2k / (2.0 * rho / lab[k]) ** (2 * self.krange), 1
@@ -154,7 +154,7 @@ def ivratiop(self, rho, R, lab):
 
     def kvratiop(self, rho, R, lab):
         lab = np.atleast_1d(lab)
-        rv = np.empty((self.Norder, len(lab)), "D")
+        rv = np.empty((self.Norder, len(lab)), dtype=complex)
         for k in range(len(lab)):
             top = np.sum(
                 self.muk * self.hankeln2k / (2.0 * rho / lab[k]) ** (2 * self.krange), 1
@@ -227,13 +227,17 @@ def initialize(self):
         #        assert self.R / abs(self.aqout.lab2[i,j,0]) < 900, 'radius too large compared to aqin lab2[i,j,0] '+str((i,j))
         # self.facin = 1.0 / iv(0, self.R / self.aqin.lab2)
         # self.facout = 1.0 / kv(0, self.R / self.aqout.lab2)
-        self.parameters = np.zeros((self.model.Ngvbc, self.Nparam, self.model.Np), "D")
+        self.parameters = np.zeros(
+            (self.model.Ngvbc, self.Nparam, self.model.Np), dtype=complex
+        )
 
     def potinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.findAquiferData(x, y)
-        rv = np.zeros((self.nparam, aq.naq, self.model.nin, self.model.npin), "D")
+        rv = np.zeros(
+            (self.nparam, aq.naq, self.model.nin, self.model.npin), dtype=complex
+        )
         if aq == self.aqin:
             r = np.sqrt((x - self.x0) ** 2 + (y - self.y0) ** 2)
             for i in range(self.aqin.Naq):
@@ -269,10 +273,12 @@ def disinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.findAquiferData(x, y)
-        qx = np.zeros((self.nparam, aq.naq, self.model.np), "D")
-        qy = np.zeros((self.nparam, aq.naq, self.model.np), "D")
+        qx = np.zeros((self.nparam, aq.naq, self.model.np), dtype=complex)
+        qy = np.zeros((self.nparam, aq.naq, self.model.np), dtype=complex)
         if aq == self.aqin:
-            qr = np.zeros((self.nparam, aq.naq, self.model.nin, self.model.npin), "D")
+            qr = np.zeros(
+                (self.nparam, aq.naq, self.model.nin, self.model.npin), dtype=complex
+            )
             r = np.sqrt((x - self.x0) ** 2 + (y - self.y0) ** 2)
             if r < 1e-20:
                 r = 1e-20  # As we divide by that on the return
@@ -296,7 +302,9 @@ def disinf(self, x, y, aq=None):
             qx[:] = qr * (x - self.x0) / r
             qy[:] = qr * (y - self.y0) / r
         if aq == self.aqout:
-            qr = np.zeros((self.Nparam, aq.Naq, self.model.Nin, self.model.Npin), "D")
+            qr = np.zeros(
+                (self.Nparam, aq.Naq, self.model.Nin, self.model.Npin), dtype=complex
+            )
             r = np.sqrt((x - self.x0) ** 2 + (y - self.y0) ** 2)
             for i in range(self.aqout.Naq):
                 for j in range(self.model.Nin):
diff --git a/ttim/element.py b/ttim/element.py
index d23da5d..332eb0e 100644
--- a/ttim/element.py
+++ b/ttim/element.py
@@ -1,11 +1,12 @@
 import inspect  # Used for storing the input
+from abc import ABC, abstractmethod
 
 import numpy as np
 
 from .invlapnumba import invlapcomp
 
 
-class Element:
+class Element(ABC):
     def __init__(
         self,
         model,
@@ -38,7 +39,7 @@ def __init__(
         self.layers = np.atleast_1d(layers)
         self.nlayers = len(self.layers)
         #
-        tsandbc = np.atleast_2d(tsandbc).astype("d")
+        tsandbc = np.atleast_2d(tsandbc).astype(float)
         tsandbc_error = (
             "tsandQ or tsandh need to be 2D lists"
             + " or arrays, like [(0, 1), (2, 5), (8, 0)] "
@@ -68,6 +69,7 @@ def setbc(self):
             self.bc = self.bcin.copy()
         self.ntstart = len(self.tstart)
 
+    @abstractmethod
     def initialize(self):
         """Initialize the element.
 
@@ -78,17 +80,18 @@ def initialize(self):
         are initialized when Model.solve is called The initialization class needs to be
         overloaded by all derived classes
         """
-        pass
 
+    @abstractmethod
     def potinf(self, x, y, aq=None):
         """Returns complex array of size (nparam, naq, npval)."""
-        raise Exception("Must overload Element.potinf()")
 
     def potential(self, x, y, aq=None):
         """Returns complex array of size (ngvbc, naq, npval)."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        return np.sum(self.parameters[:, :, np.newaxis, :] * self.potinf(x, y, aq), 1)
+        return np.sum(
+            self.parameters[:, :, np.newaxis, :] * self.potinf(x, y, aq), axis=1
+        )
 
     def unitpotential(self, x, y, aq=None):
         """Returns complex array of size (naq, npval).
@@ -108,9 +111,9 @@ def unitpotentialone(self, x, y, jtime, aq=None):
             aq = self.model.aq.find_aquifer_data(x, y)
         return np.sum(self.potinfone(x, y, jtime, aq), 0)
 
+    @abstractmethod
     def disvecinf(self, x, y, aq=None):
         """Returns 2 complex arrays of size (nparam, naq, npval)."""
-        raise Exception("Must overload Element.disvecinf()")
 
     def disvec(self, x, y, aq=None):
         """Returns 2 complex arrays of size (ngvbc, naq, npval)."""
@@ -224,7 +227,7 @@ def discharge(self, t, derivative=0):
         """
         # Could potentially be more efficient if s is pre-computed for
         # all elements, but may not be worthwhile to store as it is quick now
-        time = np.atleast_1d(t).astype("d")
+        time = np.atleast_1d(t).astype(float)
         rv = np.zeros((self.nlayers, len(time)))
         if self.type == "g":
             s = self.dischargeinflayers * self.model.p**derivative
@@ -274,7 +277,7 @@ def dischargeold(self, t, derivative=0):
         """
         # Could potentially be more efficient if s is pre-computed for
         # all elements, but may not be worthwhile to store as it is quick now
-        time = np.atleast_1d(t).astype("d")
+        time = np.atleast_1d(t).astype(float)
         if (time[0] < self.model.tmin) or (time[-1] > self.model.tmax):
             print(
                 "Warning, some of the times are smaller than tmin or"
@@ -324,13 +327,13 @@ def write(self):
         rv += ")\n"
         return rv
 
-    def run_after_solve(self):
+    def run_after_solve(self):  # noqa: B027
         """Function to run after a solution is completed.
 
-        for most elements nothing needs to be done, but for strings of elements some
-        arrays may need to be filled
+        For most elements nothing needs to be done, but for strings of elements some
+        arrays may need to be filled.
         """
-        pass
 
+    @abstractmethod
     def plot(self, ax=None):
-        pass
+        """Plot the element."""
diff --git a/ttim/equation.py b/ttim/equation.py
index 8a13fdf..7fc43dd 100644
--- a/ttim/equation.py
+++ b/ttim/equation.py
@@ -14,9 +14,13 @@ def equation(self):
         Well: q_s = Q / (2*pi*r_w*H)
         LineSink: q_s = sigma / H = Q / (L*H)
         """
-        mat = np.empty((self.nunknowns, self.model.neq, self.model.npval), "D")
+        mat = np.empty(
+            (self.nunknowns, self.model.neq, self.model.npval), dtype=complex
+        )
         # rhs needs be initialized zero
-        rhs = np.zeros((self.nunknowns, self.model.ngvbc, self.model.npval), "D")
+        rhs = np.zeros(
+            (self.nunknowns, self.model.ngvbc, self.model.npval), dtype=complex
+        )
         for icp in range(self.ncp):
             istart = icp * self.nlayers
             ieq = 0
@@ -52,8 +56,12 @@ def equation(self):
         Element with total given discharge, uniform but unknown head and
         InternalStorageEquation.
         """
-        mat = np.zeros((self.nunknowns, self.model.neq, self.model.npval), "D")
-        rhs = np.zeros((self.nunknowns, self.model.ngvbc, self.model.npval), "D")
+        mat = np.zeros(
+            (self.nunknowns, self.model.neq, self.model.npval), dtype=complex
+        )
+        rhs = np.zeros(
+            (self.nunknowns, self.model.ngvbc, self.model.npval), dtype=complex
+        )
         ieq = 0
         for e in self.model.elementlist:
             if e.nunknowns > 0:
@@ -110,8 +118,12 @@ def equation(self):
         (really written as constant potential element) Returns matrix part nunknowns,
         neq, npval, complex Returns rhs part nunknowns, nvbc, npval, complex
         """
-        mat = np.empty((self.nunknowns, self.model.neq, self.model.npval), "D")
-        rhs = np.zeros((self.nunknowns, self.model.ngvbc, self.model.npval), "D")
+        mat = np.empty(
+            (self.nunknowns, self.model.neq, self.model.npval), dtype=complex
+        )
+        rhs = np.zeros(
+            (self.nunknowns, self.model.ngvbc, self.model.npval), dtype=complex
+        )
         for icp in range(self.ncp):
             istart = icp * self.nlayers
             ieq = 0
@@ -141,8 +153,12 @@ def equation(self):
         Returns matrix part (nunknowns,neq,npval), complex
         Returns rhs part (nunknowns,nvbc,npval), complex.
         """
-        mat = np.empty((self.nunknowns, self.model.neq, self.model.npval), "D")
-        rhs = np.zeros((self.nunknowns, self.model.ngvbc, self.model.npval), "D")
+        mat = np.empty(
+            (self.nunknowns, self.model.neq, self.model.npval), dtype=complex
+        )
+        rhs = np.zeros(
+            (self.nunknowns, self.model.ngvbc, self.model.npval), dtype=complex
+        )
         for icp in range(self.ncp):
             istart = icp * self.nlayers
             ieq = 0
@@ -194,8 +210,12 @@ def equation(self):
         head_out - c * q_s = h_in
         Set h_i - h_(i + 1) = 0 and Sum Q_i = Q
         """
-        mat = np.zeros((self.nunknowns, self.model.neq, self.model.npval), "D")
-        rhs = np.zeros((self.nunknowns, self.model.ngvbc, self.model.npval), "D")
+        mat = np.zeros(
+            (self.nunknowns, self.model.neq, self.model.npval), dtype=complex
+        )
+        rhs = np.zeros(
+            (self.nunknowns, self.model.ngvbc, self.model.npval), dtype=complex
+        )
         ieq = 0
         for icp in range(self.ncp):
             istart = icp * self.nlayers
@@ -261,8 +281,12 @@ def equation(self):
         In case of storage:
         Sum Q_i - A * p^2 * headin = Q
         """
-        mat = np.zeros((self.nunknowns, self.model.neq, self.model.npval), "D")
-        rhs = np.zeros((self.nunknowns, self.model.ngvbc, self.model.npval), "D")
+        mat = np.zeros(
+            (self.nunknowns, self.model.neq, self.model.npval), dtype=complex
+        )
+        rhs = np.zeros(
+            (self.nunknowns, self.model.ngvbc, self.model.npval), dtype=complex
+        )
         ieq = 0
         for icp in range(self.ncp):
             istart = icp * self.nlayers
@@ -333,8 +357,8 @@ def equation(self):
         mat[-1, istartself : istartself + self.nparam, :] = 1.0
         if self.Astorage is not None:
             # Used to store last equation in case of ditch storage
-            matlast = np.zeros((self.model.neq, self.model.npval), "D")
-            rhslast = np.zeros((self.model.npval), "D")
+            matlast = np.zeros((self.model.neq, self.model.npval), dtype=complex)
+            rhslast = np.zeros((self.model.npval), dtype=complex)
             ieq = 0
             for e in self.model.elementlist:
                 head = (
@@ -376,8 +400,12 @@ def equation(self):
 class InhomEquation:
     def equation(self):
         """Mix-in class that returns matrix rows for inhomogeneity conditions."""
-        mat = np.zeros((self.nunknowns, self.model.neq, self.model.npval), "D")
-        rhs = np.zeros((self.nunknowns, self.model.ngvbc, self.model.npval), "D")
+        mat = np.zeros(
+            (self.nunknowns, self.model.neq, self.model.npval), dtype=complex
+        )
+        rhs = np.zeros(
+            (self.nunknowns, self.model.ngvbc, self.model.npval), dtype=complex
+        )
         for icp in range(self.ncp):
             istart = icp * 2 * self.nlayers
             ieq = 0
diff --git a/ttim/linedoublet.py b/ttim/linedoublet.py
index 27b8b40..cbb9c9b 100644
--- a/ttim/linedoublet.py
+++ b/ttim/linedoublet.py
@@ -30,8 +30,7 @@ def __init__(
         label=None,
         addtomodel=True,
     ):
-        Element.__init__(
-            self,
+        super().__init__(
             model,
             nparam=1,
             nunknowns=0,
@@ -73,7 +72,7 @@ def initialize(self):
         self.sinout = np.sin(self.thetanormOut) * np.ones(self.ncp)
         #
         thetacp = np.arange(np.pi, 0, -np.pi / self.ncp) - 0.5 * np.pi / self.ncp
-        Zcp = np.zeros(self.ncp, "D")
+        Zcp = np.zeros(self.ncp, dtype=complex)
         Zcp.real = np.cos(thetacp)
         # control point just on positive site (this is handy later on)
         Zcp.imag = 1e-6
@@ -110,9 +109,11 @@ def potinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rv = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
+        rv = np.zeros(
+            (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+        )
         if aq == self.aq:
-            pot = np.zeros((self.order + 1, self.model.npint), "D")
+            pot = np.zeros((self.order + 1, self.model.npint), dtype=complex)
             for i in range(self.aq.naq):
                 for j in range(self.model.nint):
                     if besselnumba.isinside(
@@ -141,10 +142,14 @@ def disvecinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rvx = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
-        rvy = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
+        rvx = np.zeros(
+            (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+        )
+        rvy = np.zeros(
+            (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+        )
         if aq == self.aq:
-            qxqy = np.zeros((2 * (self.order + 1), self.model.npint), "D")
+            qxqy = np.zeros((2 * (self.order + 1), self.model.npint), dtype=complex)
             for i in range(self.aq.naq):
                 for j in range(self.model.nint):
                     if besselnumba.isinside(
@@ -230,8 +235,7 @@ def __init__(
         addtomodel=True,
     ):
         self.storeinput(inspect.currentframe())
-        LineDoubletHoBase.__init__(
-            self,
+        super().__init__(
             model,
             x1=x1,
             y1=y1,
@@ -249,9 +253,9 @@ def __init__(
         self.nunknowns = self.nparam
 
     def initialize(self):
-        LineDoubletHoBase.initialize(self)
+        super().initialize()
         self.parameters = np.zeros(
-            (self.model.ngvbc, self.nparam, self.model.npval), "D"
+            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
         )
 
 
@@ -301,7 +305,7 @@ def __init__(
         self.res = res
         self.order = order
         self.ldlist = []
-        xy = np.atleast_2d(xy).astype("d")
+        xy = np.atleast_2d(xy).astype(float)
         self.x, self.y = xy[:, 0], xy[:, 1]
         self.nld = len(self.x) - 1
         for i in range(self.nld):
@@ -340,7 +344,7 @@ def initialize(self):
         self.yldlayout = np.hstack((self.yld[:, 0], self.yld[-1, 1]))
         self.aq = self.model.aq.find_aquifer_data(self.ldlist[0].xc, self.ldlist[0].yc)
         self.parameters = np.zeros(
-            (self.model.ngvbc, self.nparam, self.model.npval), "D"
+            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
         )
         self.setbc()
         # As parameters are only stored for the element not the list,
@@ -361,7 +365,7 @@ def potinf(self, x, y, aq=None):
         """Returns array (nunknowns,nperiods)."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rv = np.zeros((self.nparam, aq.naq, self.model.npval), "D")
+        rv = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
         for i, ld in enumerate(self.ldlist):
             rv[i * ld.nparam : (i + 1) * ld.nparam, :] = ld.potinf(x, y, aq)
         return rv
@@ -370,8 +374,8 @@ def disvecinf(self, x, y, aq=None):
         """Returns array (nunknowns,nperiods)."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rvx = np.zeros((self.nparam, aq.naq, self.model.npval), "D")
-        rvy = np.zeros((self.nparam, aq.naq, self.model.npval), "D")
+        rvx = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
+        rvy = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
         for i, ld in enumerate(self.ldlist):
             qx, qy = ld.disvecinf(x, y, aq)
             rvx[i * ld.nparam : (i + 1) * ld.nparam, :] = qx
diff --git a/ttim/linesink.py b/ttim/linesink.py
index 5fae724..74faf42 100644
--- a/ttim/linesink.py
+++ b/ttim/linesink.py
@@ -35,8 +35,7 @@ def __init__(
         label=None,
         addtomodel=True,
     ):
-        Element.__init__(
-            self,
+        super().__init__(
             model,
             nparam=1,
             nunknowns=0,
@@ -106,9 +105,11 @@ def potinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rv = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
+        rv = np.zeros(
+            (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+        )
         if aq == self.aq:
-            pot = np.zeros(self.model.npint, "D")
+            pot = np.zeros(self.model.npint, dtype=complex)
             for i in range(self.aq.naq):
                 for j in range(self.model.nint):
                     pot[:] = besselnumba.bessellsuniv(
@@ -123,10 +124,14 @@ def disvecinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rvx = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
-        rvy = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
+        rvx = np.zeros(
+            (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+        )
+        rvy = np.zeros(
+            (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+        )
         if aq == self.aq:
-            qxqy = np.zeros((2, self.model.npint), "D")
+            qxqy = np.zeros((2, self.model.npint), dtype=complex)
             for i in range(self.aq.naq):
                 for j in range(self.model.nint):
                     if besselnumba.isinside(
@@ -195,8 +200,7 @@ def __init__(
         addtomodel=True,
     ):
         self.storeinput(inspect.currentframe())
-        LineSinkBase.__init__(
-            self,
+        super().__init__(
             model,
             x1=x1,
             y1=y1,
@@ -302,8 +306,7 @@ def __init__(
             etype = "z"
         else:
             etype = "v"
-        LineSinkBase.__init__(
-            self,
+        super().__init__(
             model,
             x1=x1,
             y1=y1,
@@ -321,9 +324,9 @@ def __init__(
         self.nunknowns = self.nparam
 
     def initialize(self):
-        LineSinkBase.initialize(self)
+        super().initialize()
         self.parameters = np.zeros(
-            (self.model.ngvbc, self.nparam, self.model.npval), "D"
+            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
         )
         # Needed in solving, solve for a unit head
         self.pc = self.aq.T[self.layers]
@@ -339,8 +342,7 @@ def __init__(
         name="LineSinkStringBase",
         label=None,
     ):
-        Element.__init__(
-            self,
+        super().__init__(
             model,
             nparam=1,
             nunknowns=0,
@@ -369,7 +371,7 @@ def initialize(self):
         self.ylslayout = np.hstack((self.yls[:, 0], self.yls[-1, 1]))
         self.aq = self.model.aq.find_aquifer_data(self.lslist[0].xc, self.lslist[0].yc)
         self.parameters = np.zeros(
-            (self.model.ngvbc, self.nparam, self.model.npval), "D"
+            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
         )
         self.setbc()
         # As parameters are only stored for the element not the list
@@ -382,8 +384,12 @@ def initialize(self):
             self.resfacp.extend(ls.resfacp.tolist())  # Needed in solving
         self.resfach = np.array(self.resfach)
         self.resfacp = np.array(self.resfacp)
-        self.dischargeinf = np.zeros((self.nparam, self.aq.naq, self.model.npval), "D")
-        self.dischargeinflayers = np.zeros((self.nparam, self.model.npval), "D")
+        self.dischargeinf = np.zeros(
+            (self.nparam, self.aq.naq, self.model.npval), dtype=complex
+        )
+        self.dischargeinflayers = np.zeros(
+            (self.nparam, self.model.npval), dtype=complex
+        )
         self.xc, self.yc = np.zeros(self.nls), np.zeros(self.nls)
         for i in range(self.nls):
             self.dischargeinf[i * self.nlayers : (i + 1) * self.nlayers, :] = (
@@ -399,7 +405,7 @@ def potinf(self, x, y, aq=None):
         """Returns array (nunknowns, Nperiods)."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rv = np.zeros((self.nparam, aq.naq, self.model.npval), "D")
+        rv = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
         for i in range(self.nls):
             rv[i * self.nlayers : (i + 1) * self.nlayers, :] = self.lslist[i].potinf(
                 x, y, aq
@@ -410,8 +416,8 @@ def disvecinf(self, x, y, aq=None):
         """Returns array (nunknowns,Nperiods)."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rvx = np.zeros((self.nparam, aq.naq, self.model.npval), "D")
-        rvy = np.zeros((self.nparam, aq.naq, self.model.npval), "D")
+        rvx = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
+        rvy = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
         for i in range(self.nls):
             qx, qy = self.lslist[i].disvecinf(x, y, aq)
             rvx[i * self.nlayers : (i + 1) * self.nlayers, :] = qx
@@ -535,8 +541,7 @@ def __init__(
             etype = "z"
         else:
             etype = "v"
-        LineSinkStringBase.__init__(
-            self,
+        super().__init__(
             model,
             tsandbc=tsandh,
             layers=layers,
@@ -544,7 +549,7 @@ def __init__(
             name="HeadLineSinkString",
             label=label,
         )
-        xy = np.atleast_2d(xy).astype("d")
+        xy = np.atleast_2d(xy).astype(float)
         self.x = xy[:, 0]
         self.y = xy[:, 1]
         self.nls = len(self.x) - 1
@@ -571,7 +576,7 @@ def initialize(self):
                     addtomodel=False,
                 )
             )
-        LineSinkStringBase.initialize(self)
+        super().initialize()
         self.pc = np.zeros(self.nls * self.nlayers)
         for i in range(self.nls):
             self.pc[i * self.nlayers : (i + 1) * self.nlayers] = self.lslist[i].pc
@@ -601,8 +606,7 @@ def __init__(
     ):
         # assert len(layers) > 1, "number of layers must be at least 2"
         self.storeinput(inspect.currentframe())
-        LineSinkBase.__init__(
-            self,
+        super().__init__(
             model,
             x1=x1,
             y1=y1,
@@ -625,9 +629,9 @@ def __init__(
             self.vres = self.vres[0] * np.ones(self.nlayers - 1)
 
     def initialize(self):
-        LineSinkBase.initialize(self)
+        super().initialize()
         self.parameters = np.zeros(
-            (self.model.ngvbc, self.nparam, self.model.npval), "D"
+            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
         )
         # Qv = (hn - hn-1) / vresfac[n - 1]
         self.vresfac = self.vres / (self.wv * self.L)
@@ -691,8 +695,7 @@ def __init__(
         label=None,
     ):
         self.storeinput(inspect.currentframe())
-        LineSinkStringBase.__init__(
-            self,
+        super().__init__(
             model,
             tsandbc=tsandQ,
             layers=layers,
@@ -700,7 +703,7 @@ def __init__(
             name="LineSinkDitchString",
             label=label,
         )
-        xy = np.atleast_2d(xy).astype("d")
+        xy = np.atleast_2d(xy).astype(float)
         self.x, self.y = xy[:, 0], xy[:, 1]
         self.nls = len(self.x) - 1
         for i in range(self.nls):
@@ -723,7 +726,7 @@ def __init__(
         self.model.addelement(self)
 
     def initialize(self):
-        LineSinkStringBase.initialize(self)
+        super().initialize()
         # set vresfac to zero, as I don't quite know what it would mean if
         # it is not zero
         self.vresfac = np.zeros_like(self.resfach)
@@ -787,8 +790,7 @@ def __init__(
         label=None,
     ):
         self.storeinput(inspect.currentframe())
-        LineSinkStringBase.__init__(
-            self,
+        super().__init__(
             model,
             tsandbc=tsandQ,
             layers=layers,
@@ -796,7 +798,7 @@ def __init__(
             name="LineSinkDitchString",
             label=label,
         )
-        xy = np.atleast_2d(xy).astype("d")
+        xy = np.atleast_2d(xy).astype(float)
         self.x, self.y = xy[:, 0], xy[:, 1]
         self.nls = len(self.x) - 1
         for i in range(self.nls):
@@ -819,7 +821,7 @@ def __init__(
         self.model.addelement(self)
 
     def initialize(self):
-        LineSinkStringBase.initialize(self)
+        super().initialize()
         # set vresfac to zero, as I don't quite know what it would mean if
         # it is not zero
         self.vresfac = np.zeros_like(self.resfach)
@@ -848,8 +850,7 @@ def __init__(
         label=None,
         addtomodel=True,
     ):
-        Element.__init__(
-            self,
+        super().__init__(
             model,
             nparam=1,
             nunknowns=0,
@@ -886,7 +887,7 @@ def initialize(self):
         self.L = np.abs(self.z1 - self.z2)
         #
         thetacp = np.arange(np.pi, 0, -np.pi / self.ncp) - 0.5 * np.pi / self.ncp
-        Zcp = np.zeros(self.ncp, "D")
+        Zcp = np.zeros(self.ncp, dtype=complex)
         Zcp.real = np.cos(thetacp)
         # control point just on positive site (this is handy later on)
         Zcp.imag = 1e-6
@@ -926,9 +927,11 @@ def potinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rv = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
+        rv = np.zeros(
+            (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+        )
         if aq == self.aq:
-            pot = np.zeros((self.order + 1, self.model.npint), "D")
+            pot = np.zeros((self.order + 1, self.model.npint), dtype=complex)
             for i in range(self.aq.naq):
                 for j in range(self.model.nint):
                     if besselnumba.isinside(
@@ -956,10 +959,14 @@ def disvecinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rvx = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
-        rvy = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
+        rvx = np.zeros(
+            (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+        )
+        rvy = np.zeros(
+            (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+        )
         if aq == self.aq:
-            qxqy = np.zeros((2 * (self.order + 1), self.model.npint), "D")
+            qxqy = np.zeros((2 * (self.order + 1), self.model.npint), dtype=complex)
             for i in range(self.aq.naq):
                 for j in range(self.model.nint):
                     if besselnumba.isinside(
@@ -1031,8 +1038,7 @@ def __init__(
             etype = "z"
         else:
             etype = "v"
-        LineSinkHoBase.__init__(
-            self,
+        super().__init__(
             model,
             x1=x1,
             y1=y1,
@@ -1051,9 +1057,9 @@ def __init__(
         self.nunknowns = self.nparam
 
     def initialize(self):
-        LineSinkHoBase.initialize(self)
+        super().initialize()
         self.parameters = np.zeros(
-            (self.model.ngvbc, self.nparam, self.model.npval), "D"
+            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
         )
         self.pc = np.empty(self.nparam)
         for i, T in enumerate(self.aq.T[self.layers]):
diff --git a/ttim/model.py b/ttim/model.py
index d87c407..6bf1406 100644
--- a/ttim/model.py
+++ b/ttim/model.py
@@ -179,7 +179,7 @@ def potential(self, x, y, t, layers=None, aq=None, derivative=0, returnphi=0):
             layers = range(aq.naq)
         nlayers = len(layers)
         time = np.atleast_1d(t) - self.tstart  # used to be ).copy()
-        pot = np.zeros((self.ngvbc, aq.naq, self.npval), "D")
+        pot = np.zeros((self.ngvbc, aq.naq, self.npval), dtype=complex)
         for i in range(self.ngbc):
             pot[i, :] += self.gbclist[i].unitpotential(x, y, aq)
         for e in self.vzbclist:
@@ -218,7 +218,7 @@ def potentialone(self, x, y, time, layers=None, aq=None, derivative=0, returnphi
         time = np.atleast_1d(time) - self.tstart  # used to be ).copy()
         jtime = np.searchsorted(self.tintervals, time)[0] - 1
         assert 0 <= jtime <= len(self.tintervals), "time not in tintervals"
-        pot = np.zeros((self.ngvbc, aq.naq, self.npint), "D")
+        pot = np.zeros((self.ngvbc, aq.naq, self.npint), dtype=complex)
         for i in range(self.ngbc):
             pot[i, :] += self.gbclist[i].unitpotentialone(x, y, jtime, aq)
         for e in self.vzbclist:
@@ -258,8 +258,8 @@ def disvec(self, x, y, t, layers=None, aq=None, derivative=0):
             layers = range(aq.naq)
         nlayers = len(layers)
         time = np.atleast_1d(t) - self.tstart
-        disx = np.zeros((self.ngvbc, aq.naq, self.npval), "D")
-        disy = np.zeros((self.ngvbc, aq.naq, self.npval), "D")
+        disx = np.zeros((self.ngvbc, aq.naq, self.npval), dtype=complex)
+        disy = np.zeros((self.ngvbc, aq.naq, self.npval), dtype=complex)
         for i in range(self.ngbc):
             qx, qy = self.gbclist[i].unitdisvec(x, y, aq)
             disx[i, :] += qx
@@ -577,8 +577,8 @@ def solve(self, printmat=0, sendback=0, silent=False):
             if silent is False:
                 print("No unknowns. Solution complete")
             return
-        mat = np.empty((self.neq, self.neq, self.npval), "D")
-        rhs = np.empty((self.neq, self.ngvbc, self.npval), "D")
+        mat = np.empty((self.neq, self.neq, self.npval), dtype=complex)
+        rhs = np.empty((self.neq, self.ngvbc, self.npval), dtype=complex)
         ieq = 0
         for e in self.elementlist:
             if e.nunknowns > 0:
diff --git a/ttim/well.py b/ttim/well.py
index 6b083d6..9c54501 100644
--- a/ttim/well.py
+++ b/ttim/well.py
@@ -29,8 +29,7 @@ def __init__(
         name="WellBase",
         label=None,
     ):
-        Element.__init__(
-            self,
+        super().__init__(
             model,
             nparam=1,
             nunknowns=0,
@@ -85,10 +84,12 @@ def potinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rv = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
+        rv = np.zeros(
+            (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+        )
         if aq == self.aq:
             r = np.sqrt((x - self.xw) ** 2 + (y - self.yw) ** 2)
-            pot = np.zeros(self.model.npint, "D")
+            pot = np.zeros(self.model.npint, dtype=complex)
             if r < self.rw:
                 r = self.rw  # If at well, set to at radius
             for i in range(self.aq.naq):
@@ -105,10 +106,10 @@ def potinfone(self, x, y, jtime, aq=None):
         """Can be called with only one x,y value for time interval jtime."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        rv = np.zeros((self.nparam, aq.naq, self.model.npint), "D")
+        rv = np.zeros((self.nparam, aq.naq, self.model.npint), dtype=complex)
         if aq == self.aq:
             r = np.sqrt((x - self.xw) ** 2 + (y - self.yw) ** 2)
-            pot = np.zeros(self.model.npint, "D")
+            pot = np.zeros(self.model.npint, dtype=complex)
             if r < self.rw:
                 r = self.rw  # If at well, set to at radius
             for i in range(self.aq.naq):
@@ -122,12 +123,14 @@ def disvecinf(self, x, y, aq=None):
         """Can be called with only one x,y value."""
         if aq is None:
             aq = self.model.aq.find_aquifer_data(x, y)
-        qx = np.zeros((self.nparam, aq.naq, self.model.npval), "D")
-        qy = np.zeros((self.nparam, aq.naq, self.model.npval), "D")
+        qx = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
+        qy = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
         if aq == self.aq:
-            qr = np.zeros((self.nparam, aq.naq, self.model.nint, self.model.npint), "D")
+            qr = np.zeros(
+                (self.nparam, aq.naq, self.model.nint, self.model.npint), dtype=complex
+            )
             r = np.sqrt((x - self.xw) ** 2 + (y - self.yw) ** 2)
-            # pot = np.zeros(self.model.npint, "D")
+            # pot = np.zeros(self.model.npint, dtype=complex)
             if r < self.rw:
                 r = self.rw  # If at well, set to at radius
             for i in range(self.aq.naq):
@@ -244,8 +247,7 @@ def __init__(
         self, model, xw=0, yw=0, tsandQ=[(0, 1)], rw=0.1, res=0, layers=0, label=None
     ):
         self.storeinput(inspect.currentframe())
-        WellBase.__init__(
-            self,
+        super().__init__(
             model,
             xw,
             yw,
@@ -313,8 +315,7 @@ def __init__(
         label=None,
     ):
         self.storeinput(inspect.currentframe())
-        WellBase.__init__(
-            self,
+        super().__init__(
             model,
             xw,
             yw,
@@ -342,9 +343,9 @@ def __init__(
         self.wbstype = wbstype
 
     def initialize(self):
-        WellBase.initialize(self)
+        super().initialize()
         self.parameters = np.zeros(
-            (self.model.ngvbc, self.nparam, self.model.npval), "D"
+            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
         )
 
     def setflowcoef(self):
@@ -392,8 +393,7 @@ def __init__(
         self, model, xw=0, yw=0, rw=0.1, tsandh=[(0, 1)], res=0, layers=0, label=None
     ):
         self.storeinput(inspect.currentframe())
-        WellBase.__init__(
-            self,
+        super().__init__(
             model,
             xw,
             yw,
@@ -408,9 +408,9 @@ def __init__(
         self.nunknowns = self.nparam
 
     def initialize(self):
-        WellBase.initialize(self)
+        super().initialize()
         self.parameters = np.zeros(
-            (self.model.ngvbc, self.nparam, self.model.npval), "D"
+            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
         )
         # Needed in solving for a unit head
         self.pc = self.aq.T[self.layers]
@@ -430,8 +430,7 @@ def __init__(
         fp=None,
     ):
         self.storeinput(inspect.currentframe())
-        WellBase.__init__(
-            self,
+        super().__init__(
             model,
             xw,
             yw,