diff --git a/timflow/steady/aquifer.py b/timflow/steady/aquifer.py
index 89713d2..9bcdec7 100644
--- a/timflow/steady/aquifer.py
+++ b/timflow/steady/aquifer.py
@@ -9,8 +9,6 @@
 
 """
 
-import inspect  # Used for storing the input
-
 import numpy as np
 import pandas as pd
 
@@ -132,9 +130,6 @@ def add_element(self, e):
     def isinside(self, x, y):
         raise Exception("Must overload AquiferData.isinside()")
 
-    def storeinput(self, frame):
-        self.inputargs, _, _, self.inputvalues = inspect.getargvalues(frame)
-
     def findlayer(self, z):
         """Returns layer-number, layer-type and model-layer-number."""
         if z > self.z[0]:
diff --git a/timflow/steady/circareasink.py b/timflow/steady/circareasink.py
index 0cf7c25..4b92e85 100644
--- a/timflow/steady/circareasink.py
+++ b/timflow/steady/circareasink.py
@@ -7,8 +7,6 @@
     CircAreaSink(ml, xc=0, yc=0, R=50, N=0.001, layer=0)
 """
 
-import inspect  # Used for storing the input
-
 import numpy as np
 from scipy.special import i0, i1, k0, k1
 
@@ -43,7 +41,6 @@ class CircAreaSink(Element):
     def __init__(
         self, model, xc=0, yc=0, R=1, N=0.001, layer=0, name="CircAreasink", label=None
     ):
-        self.storeinput(inspect.currentframe())
         Element.__init__(
             self, model, nparam=1, nunknowns=0, layers=layer, name=name, label=label
         )
diff --git a/timflow/steady/constant.py b/timflow/steady/constant.py
index 89d662b..dde0bab 100644
--- a/timflow/steady/constant.py
+++ b/timflow/steady/constant.py
@@ -8,8 +8,6 @@
     Constant(ml, xr=0, yr=0, hr=10.0, layer=0)
 """
 
-import inspect  # Used for storing the input
-
 import numpy as np
 
 from timflow.steady.element import Element
@@ -96,7 +94,6 @@ class Constant(ConstantBase, PotentialEquation):
     """
 
     def __init__(self, model, xr=0, yr=0, hr=0.0, layer=0, label=None):
-        self.storeinput(inspect.currentframe())
         ConstantBase.__init__(
             self, model, xr=xr, yr=yr, hr=hr, layer=layer, name="Constant", label=label
         )
diff --git a/timflow/steady/element.py b/timflow/steady/element.py
index bd93a5c..8095f98 100644
--- a/timflow/steady/element.py
+++ b/timflow/steady/element.py
@@ -9,8 +9,6 @@ def initialize(self):
             pass
 """
 
-import inspect  # Used for storing the input
-
 import numpy as np
 
 __all__ = ["Element"]
@@ -140,9 +138,6 @@ def headinside(self):
     def setparams(self, sol):
         raise Exception("Must overload Element.setparams()")
 
-    def storeinput(self, frame):
-        self.inputargs, _, _, self.inputvalues = inspect.getargvalues(frame)
-
     # def stoptrace(self, xyz, layer, ltype, step, direction):
     #    return False, 0
 
@@ -161,20 +156,3 @@ def qztop(self, x, y, aq):
 
     def plot(self, **kwargs):
         pass
-
-    def write(self):
-        rv = self.name + "(" + self.model.modelname + ",\n"
-        for key in self.inputargs[2:]:  # The first two are ignored
-            if isinstance(self.inputvalues[key], np.ndarray):
-                rv += (
-                    key
-                    + " = "
-                    + np.array2string(self.inputvalues[key], separator=",")
-                    + ",\n"
-                )
-            elif isinstance(self.inputvalues[key], str):
-                rv += key + " = '" + self.inputvalues[key] + "',\n"
-            else:
-                rv += key + " = " + str(self.inputvalues[key]) + ",\n"
-        rv += ")\n"
-        return rv
diff --git a/timflow/steady/inhomogeneity.py b/timflow/steady/inhomogeneity.py
index b5caa55..9f94158 100644
--- a/timflow/steady/inhomogeneity.py
+++ b/timflow/steady/inhomogeneity.py
@@ -13,7 +13,6 @@
     )
 """
 
-import inspect  # Used for storing the input
 from warnings import warn
 
 import numpy as np
@@ -218,7 +217,6 @@ def __init__(
             assert topboundary[:4] == "conf", (
                 "Error: infiltration can only be added if topboundary='conf'"
             )
-        self.storeinput(inspect.currentframe())
         (
             kaq,
             c,
@@ -315,7 +313,6 @@ def __init__(
             assert topboundary[:4] == "conf", (
                 "Error: infiltration can only be added if topboundary='conf'"
             )
-        self.storeinput(inspect.currentframe())
         kaq, kzoverkh, c, npor, ltype = param_3d(
             kaq, z, kzoverkh, npor, topboundary, topres
         )
diff --git a/timflow/steady/inhomogeneity1d.py b/timflow/steady/inhomogeneity1d.py
index 74665f3..0131ed8 100644
--- a/timflow/steady/inhomogeneity1d.py
+++ b/timflow/steady/inhomogeneity1d.py
@@ -10,7 +10,6 @@
 
 """
 
-import inspect  # user for storing the input
 import warnings
 
 import matplotlib.pyplot as plt
@@ -354,7 +353,6 @@ def __init__(
             c = []
         if z is None:
             z = [1, 0]
-        self.storeinput(inspect.currentframe())
         (
             kaq,
             c,
@@ -433,7 +431,6 @@ def __init__(
     ):
         if z is None:
             z = [1, 0]
-        self.storeinput(inspect.currentframe())
         (
             kaq,
             kzoverkh,
diff --git a/timflow/steady/linedoublet.py b/timflow/steady/linedoublet.py
index b9997c1..29ad827 100644
--- a/timflow/steady/linedoublet.py
+++ b/timflow/steady/linedoublet.py
@@ -7,7 +7,6 @@
     ImpLineDoublet(ml, x1=-10, y1=0, x2=10, y2=0, layers=0)
 """
 
-import inspect  # Used for storing the input
 import warnings
 
 import matplotlib.pyplot as plt
@@ -228,7 +227,6 @@ def __init__(
         label=None,
         addtomodel=True,
     ):
-        self.storeinput(inspect.currentframe())
         LineDoubletHoBase.__init__(
             self,
             model,
@@ -302,7 +300,6 @@ def __init__(
         label=None,
         addtomodel=True,
     ):
-        self.storeinput(inspect.currentframe())
         LineDoubletHoBase.__init__(
             self,
             model,
@@ -454,7 +451,6 @@ class ImpermeableWallString(LineDoubletStringBase, DisvecEquation):
     def __init__(self, model, xy=None, layers=0, order=0, label=None):
         if xy is None:
             xy = [(-1, 0), (1, 0)]
-        self.storeinput(inspect.currentframe())
         LineDoubletStringBase.__init__(
             self,
             model,
@@ -507,7 +503,6 @@ class LeakyWallString(LineDoubletStringBase, LeakyWallEquation):
     def __init__(self, model, xy=None, res=np.inf, layers=0, order=0, label=None):
         if xy is None:
             xy = [(-1, 0), (1, 0)]
-        self.storeinput(inspect.currentframe())
         LineDoubletStringBase.__init__(
             self,
             model,
diff --git a/timflow/steady/linedoublet1d.py b/timflow/steady/linedoublet1d.py
index e08f8a5..3fcd5e6 100644
--- a/timflow/steady/linedoublet1d.py
+++ b/timflow/steady/linedoublet1d.py
@@ -8,7 +8,6 @@
 
 """
 
-import inspect  # Used for storing the input
 import warnings
 
 import matplotlib.pyplot as plt
@@ -127,7 +126,6 @@ class ImpermeableWall1D(LineDoublet1D, DisvecEquation):
     """Create 1D impermeable wall."""
 
     def __init__(self, model, xld=0, layers=0, label=None):
-        self.storeinput(inspect.currentframe())
         LineDoublet1D.__init__(
             self,
             model,
@@ -173,7 +171,6 @@ class LeakyWall1D(LineDoublet1D, LeakyWallEquation):
     tiny = 1e-6
 
     def __init__(self, model, xld=0, res=np.inf, layers=0, label=None):
-        self.storeinput(inspect.currentframe())
         LineDoublet1D.__init__(
             self,
             model,
diff --git a/timflow/steady/linesink.py b/timflow/steady/linesink.py
index e540fb1..145accd 100644
--- a/timflow/steady/linesink.py
+++ b/timflow/steady/linesink.py
@@ -7,7 +7,6 @@
     River(ml, x1=-10, y1=0, x2=10, y2=0, hls=5, layers=0)
 """
 
-import inspect  # Used for storing the input
 import warnings
 
 import matplotlib.pyplot as plt
@@ -286,7 +285,6 @@ def __init__(
         label=None,
         addtomodel=True,
     ):
-        self.storeinput(inspect.currentframe())
         LineSinkBase.__init__(
             self,
             model,
@@ -566,7 +564,6 @@ def __init__(
         addtomodel=True,
     ):
         """Initialize a steady head-specified line-sink."""
-        self.storeinput(inspect.currentframe())
         LineSinkHoBase.__init__(
             self,
             model,
@@ -672,7 +669,6 @@ def __init__(
         addtomodel=True,
     ):
         """Initialize a steady discharge-specified line-sink."""
-        self.storeinput(inspect.currentframe())
         River.__init__(
             self,
             model,
@@ -933,7 +929,6 @@ def __init__(
         """Initialize a steady string of head-specified line-sinks."""
         if xy is None:
             xy = [(-1, 0), (1, 0)]
-        self.storeinput(inspect.currentframe())
         LineSinkStringBase2.__init__(
             self,
             model,
@@ -1086,7 +1081,6 @@ def __init__(
         """Initialize a steady string of discharge-specified line-sinks."""
         if xy is None:
             xy = [(-1, 0), (1, 0)]
-        self.storeinput(inspect.currentframe())
         RiverString.__init__(
             self,
             model,
diff --git a/timflow/steady/linesink1d.py b/timflow/steady/linesink1d.py
index c5ca4d4..bfff0c1 100644
--- a/timflow/steady/linesink1d.py
+++ b/timflow/steady/linesink1d.py
@@ -7,7 +7,6 @@
     River1D(ml, xls=0, hls=1, layers=0)
 """
 
-import inspect  # Used for storing the input
 import warnings
 
 import matplotlib.pyplot as plt
@@ -171,7 +170,6 @@ class LineSink1D(LineSink1DBase, MscreenWellEquation):
 
     def __init__(self, model, xls=0, sigls=1, layers=0, label=None):
         """Initialize a steady 1D line-sink with specified discharge per length."""
-        self.storeinput(inspect.currentframe())
         LineSink1DBase.__init__(
             self,
             model,
@@ -228,7 +226,6 @@ class River1D(LineSink1DBase, HeadEquation):
 
     def __init__(self, model, xls=0, hls=1, res=0, wh=1, layers=0, label=None):
         """Initialize a steady 1D line-sink with a specified head."""
-        self.storeinput(inspect.currentframe())
         LineSink1DBase.__init__(
             self,
             model,
diff --git a/timflow/steady/model.py b/timflow/steady/model.py
index 79ae990..d013f89 100644
--- a/timflow/steady/model.py
+++ b/timflow/steady/model.py
@@ -10,7 +10,6 @@
     ml.solve()
 """
 
-import inspect  # Used for storing the input
 import multiprocessing as mp
 import warnings
 
@@ -27,6 +26,22 @@
 __all__ = ["Model", "ModelMaq", "Model3D", "ModelXsection"]
 
 
+def _compute_head_mp(args):
+    """Helper function for parallel computation of head_array."""
+    model, xi, yi, layers, i = args
+    return i, model.head(xi, yi, layers=layers)
+
+
+def _compute_velocity_mp(args):
+    """Helper function for parallel computation of velocity_array."""
+    model, xi, yi, zi, i = args
+    try:
+        vv = model.velocomp(xi, yi, zi)
+    except (ZeroDivisionError, ValueError):
+        vv = np.full((3,), np.nan)
+    return i, vv
+
+
 class Model:
     """Create a model consisting of an arbitrary sequence of aquifers and leaky layers.
 
@@ -60,7 +75,6 @@ def __init__(self, kaq, z, c, npor, ltype, model3d=False):
         self.elementlist = []
         self.elementdict = {}  # only elements that have a label
         self.aq = Aquifer(self, kaq, c, z, npor, ltype, model3d=model3d)
-        self.modelname = "ml"  # Used for writing out input
         self.name = "Model"
         self.model_type = "steady"  # Model type for plotting and other purposes
 
@@ -84,9 +98,6 @@ def remove_element(self, e):
             self.elementdict.pop(e.label)
         self.elementlist.remove(e)
 
-    def storeinput(self, frame):
-        self.inputargs, _, _, self.inputvalues = inspect.getargvalues(frame)
-
     def potential(self, x, y, aq=None):
         if aq is None:
             aq = self.aq.find_aquifer_data(x, y)
@@ -259,6 +270,64 @@ def head(self, x, y, layers=None, aq=None):
         else:
             return rv[layers]
 
+    def head_array(self, x, y, layers=None, show_progress=False, parallel=False):
+        """Head for array of points.
+
+        Parameters
+        ----------
+        x : 1D array or list
+            x values of points
+        y : 1D array or list
+            y values of points
+        layers : integer, list or array, optional
+            layers for which grid is returned
+        show_progress : bool
+            show computation progress, by printing dots per row or with tqdm progressbar
+            when parallel is True. Default is False.
+        parallel : bool or int, optional
+            if `True`, computes head_array in parallel using multiprocessing,
+            by default `False`. If an integer is provided, it is interpreted as the
+            number of processes to use.
+
+        Returns
+        -------
+        h : array
+            heads array with size (nlayers, npoints)
+        """
+        parallel, process_map, tqdm = check_tqdm_parallel(parallel)
+        x = np.atleast_1d(x)
+        y = np.atleast_1d(y)
+        npts = len(x)
+        assert npts == len(y), "x and y must have the same length"
+        if layers is None:
+            nlayers = self.aq.find_aquifer_data(x[0], y[0]).naq
+        else:
+            nlayers = len(np.atleast_1d(layers))
+        h = np.empty((nlayers, npts))
+        if not parallel:
+            for i in (
+                tqdm(range(npts), disable=not show_progress) if tqdm else range(npts)
+            ):
+                h[:, i] = self.head(x[i], y[i], layers)
+        else:
+            nproc = mp.cpu_count() if parallel is True else int(parallel)
+            chunksize = max(1, npts // (4 * nproc)) if nproc > 0 else 1
+            tasks = [(self, x[i], y[i], layers, i) for i in range(npts)]
+            results = process_map(
+                _compute_head_mp,
+                tasks,
+                total=npts,
+                desc="head array",
+                disable=not show_progress,
+                tqdm_class=tqdm,
+                max_workers=nproc,
+                chunksize=chunksize,
+            )
+
+            for i, result in results:
+                h[:, i] = result
+        return h
+
     def headgrid(
         self, xg, yg, layers=None, printrow=False, show_progress=False, parallel=False
     ):
@@ -275,9 +344,10 @@ def headgrid(
         show_progress : bool
             show computation progress, by printing dots per row or with tqdm progressbar
             when parallel is True. Default is False.
-        parallel : bool, optional
-            if `True`, computes headgrid in parallel using multi threading,
-            by default `False`
+        parallel : bool or int, optional
+            if `True`, computes headgrid in parallel using multiprocessing,
+            by default `False`. If an integer is provided, it is interpreted as the
+            number of processes to use.
         printrow : bool, optional
 
             .. deprecated:: 0.2.0
@@ -289,6 +359,7 @@ def headgrid(
 
         See Also
         --------
+        :func:`~timflow.steady.Model.head_array`
         :func:`~timflow.steady.Model.headgrid2`
         """
         if printrow:
@@ -298,43 +369,15 @@ def headgrid(
                 stacklevel=2,
             )
             show_progress = printrow
-
-        parallel, thread_map, tqdm = check_tqdm_parallel(parallel)
-
-        xg = np.atleast_1d(xg)
-        yg = np.atleast_1d(yg)
-        nx, ny = len(xg), len(yg)
-        if layers is None:
-            Nlayers = self.aq.find_aquifer_data(xg[0], yg[0]).naq
-        else:
-            Nlayers = len(np.atleast_1d(layers))
-        h = np.empty((Nlayers, ny, nx))
-        if not parallel:
-            for j in range(ny):
-                if show_progress:
-                    print(".", end="", flush=True)
-                for i in range(nx):
-                    h[:, j, i] = self.head(xg[i], yg[j], layers)
-            if show_progress:
-                print("", flush=True)
-        else:
-
-            def compute(ij):
-                i, j = ij
-                return i, j, self.head(xg[i], yg[j], layers)
-
-            results = thread_map(
-                compute,
-                [(i, j) for j in range(ny) for i in range(nx)],
-                total=nx * ny,
-                desc="headgrid",
-                disable=not show_progress,
-                tqdm_class=tqdm,
-            )
-            for i, j, result in results:
-                h[:, j, i] = result
-
-        return h
+        x, y = np.meshgrid(xg, yg)
+        h = self.head_array(
+            x.ravel(),
+            y.ravel(),
+            layers=layers,
+            show_progress=show_progress,
+            parallel=parallel,
+        )
+        return h.reshape((h.shape[0], len(yg), len(xg)))
 
     def headgrid2(
         self,
@@ -362,6 +405,10 @@ def headgrid2(
         show_progress : bool
             show computation progress, by printing dots per row or with tqdm progressbar
             when parallel is True. Default is False.
+        parallel : bool or int, optional
+            if `True`, computes headgrid in parallel using multiprocessing,
+            by default `False`. If an integer is provided, it is interpreted as the
+            number of processes to use.
         printrow : boolean, optional
 
             .. deprecated:: 0.2.0
@@ -373,6 +420,7 @@ def headgrid2(
 
         See Also
         --------
+        :func:`~timflow.steady.Model.head_array`
         :func:`~timflow.steady.Model.headgrid`
         """
         xg, yg = np.linspace(x1, x2, nx), np.linspace(y1, y2, ny)
@@ -436,7 +484,7 @@ def disvecgrid(
             show computation progress, by printing dots per row or with tqdm progressbar
             when parallel is True. Default is True.
         parallel : bool, optional
-            if `True`, computes discharge vector grid in parallel using multi threading,
+            if `True`, computes discharge vector grid in parallel using multiprocessing,
             by default `False`
 
         Returns
@@ -570,9 +618,69 @@ def velocity(self, x, y, z):
         """
         return self.velocomp(x, y, z)
 
-    def velocity_grid(self, xg, yg, zg, show_progress=True, parallel=False):
+    def velocity_array(self, x, y, z, show_progress=True, parallel=False):
         """Compute velocity grid.
 
+        Parameters
+        ----------
+        x : 1d-array
+            x values
+        y : 1d-array
+            y values
+        z : 1d-array
+            z values
+        show_progress : bool
+            show computation progress with tqdm progressbar if tqdm is installed.
+            Default is True.
+        parallel : bool or int, optional
+            if `True`, computes velocity grid in parallel using multi processing,
+            by default `False`. If an integer is provided, it is interpreted
+            as the number of processes to use.
+
+        Returns
+        -------
+        velocity : array
+            velocity vector (vx, vy, vz) at each point in grid,
+            size (3, len(x))
+        """
+        parallel, process_map, tqdm = check_tqdm_parallel(parallel)
+        x = np.atleast_1d(x)
+        y = np.atleast_1d(y)
+        z = np.atleast_1d(z)
+        npts = len(x)
+        assert npts == len(y) == len(z), "x, y and z must have the same length"
+        v = np.empty((3, npts))
+        if not parallel:
+            for i in (
+                tqdm(range(npts), disable=not show_progress) if tqdm else range(npts)
+            ):
+                try:
+                    vv = self.velocomp(x[i], y[i], z[i])
+                except (ZeroDivisionError, ValueError):
+                    vv = np.full((3,), np.nan)
+                v[:, i] = vv
+        else:
+            nproc = mp.cpu_count() if parallel is True else int(parallel)
+            chunksize = max(1, npts // (4 * nproc)) if nproc > 0 else 1
+            tasks = [(self, x[i], y[i], z[i], i) for i in range(npts)]
+            results = process_map(
+                _compute_velocity_mp,
+                tasks,
+                total=npts,
+                desc="velocity array",
+                disable=not show_progress,
+                tqdm_class=tqdm,
+                max_workers=nproc,
+                chunksize=chunksize,
+            )
+            for i, result in results:
+                v[:, i] = result
+
+        return v
+
+    def velocity_grid(self, xg, yg, zg, show_progress=True, parallel=False):
+        """Compute velocities for an array of points.
+
         Parameters
         ----------
         xg : 1d-array
@@ -581,57 +689,31 @@ def velocity_grid(self, xg, yg, zg, show_progress=True, parallel=False):
             y values of grid
         zg : 1d-array
             z values of grid
-        parallel : bool, optional
-            if `True`, computes velocity grid in parallel using multi threading,
-            by default `False`
+        show_progress : bool, optional
+            if `True`, shows progress bar when computing velocity grid, by default `True`
+        parallel : bool or int, optional
+            if `True`, computes velocity grid in parallel using multiprocessing,
+            by default `False`. If an integer is provided, it specifies the number of
+            processes to use.
 
         Returns
         -------
         velocity : array
-            velocity vector (vz, vy, vx) at each point in grid,
-            size (3, len(z), len(y), len(x))
+            velocity vector (vx, vy, vz) at each point in grid,
+            size (3, len(zg), len(yg), len(xg))
         """
-        parallel, thread_map, tqdm = check_tqdm_parallel(parallel)
-
-        def compute(kij):
-            k, i, j = kij
-            try:
-                vv = self.velocomp(xg[j], yg[i], zg[k])
-            except ZeroDivisionError:
-                vv = np.full((3,), np.nan)
-            return k, i, j, vv
-
         xg = np.atleast_1d(xg)
         yg = np.atleast_1d(yg)
         zg = np.atleast_1d(zg)
-        nz, ny, nx = len(zg), len(yg), len(xg)
-        v = np.empty((3, nz, ny, nx))
-        if not parallel:
-            for k in range(nz):
-                if show_progress:
-                    print(".", end="", flush=True)
-                for i in range(ny):
-                    for j in range(nx):
-                        try:
-                            vv = self.velocomp(xg[j], yg[i], zg[k])
-                        except ZeroDivisionError:
-                            vv = np.full((3,), np.nan)
-                        v[:, k, i, j] = vv
-            if show_progress:
-                print("", flush=True)
-        else:
-            results = thread_map(
-                compute,
-                [(k, i, j) for k in range(nz) for i in range(ny) for j in range(nx)],
-                total=nz * nx * ny,
-                desc="velocity grid",
-                disable=not show_progress,
-                tqdm_class=tqdm,
-            )
-            for k, i, j, result in results:
-                v[:, k, i, j] = result
-
-        return v
+        x, y, z = np.meshgrid(xg, yg, zg)
+        v = self.velocity_array(
+            x.ravel(),
+            y.ravel(),
+            z.ravel(),
+            show_progress=show_progress,
+            parallel=parallel,
+        )
+        return v.reshape((3, len(yg), len(xg), len(zg))).transpose((0, 3, 1, 2))
 
     def velocomp(self, x, y, z, aq=None, layer_ltype=None):
         if aq is None:
@@ -792,32 +874,6 @@ def solve_mp(self, nproc=4, printmat=0, sendback=0, silent=False):
             return sol
         return
 
-    def write(self):
-        rv = self.modelname + " = " + self.name + "(\n"
-        for key in self.inputargs[1:]:  # The first argument (self) is ignored
-            if isinstance(self.inputvalues[key], np.ndarray):
-                rv += (
-                    key
-                    + " = "
-                    + np.array2string(self.inputvalues[key], separator=",")
-                    + ",\n"
-                )
-            elif isinstance(self.inputvalues[key], str):
-                rv += key + " = '" + self.inputvalues[key] + "',\n"
-            else:
-                rv += key + " = " + str(self.inputvalues[key]) + ",\n"
-        rv += ")\n"
-        return rv
-
-    def writemodel(self, fname):
-        self.initialize()  # So that the model can be written without solving first
-        f = open(fname, "w")
-        f.write("from timflow.steady import *\n")
-        f.write(self.write())
-        for e in self.elementlist:
-            f.write(e.write())
-        f.close()
-
     def aquifer_summary(self):
         """Return DataFrame with summary of aquifer(s) parameters in model.
 
@@ -925,7 +981,6 @@ def __init__(self, kaq=1, z=None, c=None, npor=0.3, topboundary="conf", hstar=No
             c = []
         if z is None:
             z = [1, 0]
-        self.storeinput(inspect.currentframe())
         kaq, c, npor, ltype = param_maq(kaq, z, c, npor, topboundary)
         super().__init__(kaq=kaq, z=z, c=c, npor=npor, ltype=ltype)
         self.name = "ModelMaq"
@@ -1004,7 +1059,6 @@ def __init__(
         """
         if z is None:
             z = [1, 0]
-        self.storeinput(inspect.currentframe())
         kaq, kzoverkh, c, npor, ltype = param_3d(
             kaq, z, kzoverkh, npor, topboundary, topres
         )
@@ -1041,7 +1095,6 @@ def __init__(self, naq=1):
         self.elementlist = []
         self.elementdict = {}  # only elements that have a label
         self.aq = SimpleAquifer(naq)
-        self.modelname = "ml"  # Used for writing out input
 
         self.plots = PlotSteady(self)
         self.name = "ModelXsection"
diff --git a/timflow/steady/plots.py b/timflow/steady/plots.py
index 3db4baa..4f9dcd8 100644
--- a/timflow/steady/plots.py
+++ b/timflow/steady/plots.py
@@ -98,9 +98,10 @@ def contour(
             if list of strings: use strings as names in legend
         return_contours : bool, optional
             if True, return list of contour sets for each contoured layer
-        parallel : bool, optional
-            if True, compute headgrid in parallel using multiple threads,
-            default is False
+        parallel : bool or int, optional
+            if True, compute head grid in parallel using multiprocessing,
+            default is False. If int is provided, it is interpreted as the number of
+            processes to use.
         **kwargs
             additional keyword arguments passed to ax.contour()
 
@@ -783,7 +784,7 @@ def quiver_xy(
         figsize : tuple of 2 values
             size of figure
         parallel : bool
-            if True, compute velocity grid in parallel using multiple threads,
+            if True, compute velocity grid in parallel using multiprocessing,
             default is False
         **kwargs
             additional keyword arguments passed to ax.quiver()
@@ -819,7 +820,7 @@ def quiver_z(
         figsize : tuple of 2 values
             size of figure
         parallel : bool
-            if True, compute velocity grid in parallel using multiple threads,
+            if True, compute velocity grid in parallel using multiprocessing,
             default is False
         **kwargs
             additional keyword arguments passed to ax.quiver()
diff --git a/timflow/steady/uflow.py b/timflow/steady/uflow.py
index 83b36d5..1242426 100644
--- a/timflow/steady/uflow.py
+++ b/timflow/steady/uflow.py
@@ -7,8 +7,6 @@
     Uflow(ml, slope=1e-4, angle=0.0)
 """
 
-import inspect  # Used for storing the input
-
 import numpy as np
 
 from timflow.steady.element import Element
@@ -42,7 +40,6 @@ def __init__(self, model, slope, angle, label=None):
             "timflow error: Uflow can only be added to "
             "model with background confined aquifer"
         )
-        self.storeinput(inspect.currentframe())
         Element.__init__(
             self, model, nparam=2, nunknowns=0, layers=0, name="Uflow", label=label
         )
diff --git a/timflow/steady/well.py b/timflow/steady/well.py
index ea128ac..a5cd27d 100644
--- a/timflow/steady/well.py
+++ b/timflow/steady/well.py
@@ -8,7 +8,6 @@
     Well(ml, xw=100, yw=200, Qw=1000, layers=[0])
 """
 
-import inspect  # Used for storing the input
 import warnings
 
 import matplotlib.pyplot as plt
@@ -405,7 +404,6 @@ def __init__(
         yc=None,
     ):
         """Initialize a steady well with a specified discharge."""
-        self.storeinput(inspect.currentframe())
         WellBase.__init__(
             self,
             model,
@@ -502,7 +500,6 @@ def __init__(
         addtomodel=True,
     ):
         """Initialize a steady well with a specified head."""
-        self.storeinput(inspect.currentframe())
         WellBase.__init__(
             self,
             model,
@@ -584,7 +581,6 @@ def __init__(
         addtomodel=True,
     ):
         """Initialize a steady well with a target head at a control point."""
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             xw,
@@ -687,7 +683,6 @@ def __init__(
         xc=None,
         yc=None,
     ):
-        self.storeinput(inspect.currentframe())
         WellBase.__init__(
             self,
             model,
diff --git a/timflow/transient/circareasink.py b/timflow/transient/circareasink.py
index 83797c2..2b55039 100644
--- a/timflow/transient/circareasink.py
+++ b/timflow/transient/circareasink.py
@@ -8,8 +8,6 @@
     CircAreaSink(ml, xc=0, yc=0, R=50, tsandN=[(0, 0.001)], layers=0)
 """
 
-import inspect  # Used for storing the input
-
 import matplotlib.pyplot as plt
 import numpy as np
 from scipy.special import iv, kv
@@ -40,7 +38,6 @@ class CircAreaSink(Element):
     def __init__(
         self, model, xc=0, yc=0, R=0.1, tsandN=[(0, 1)], name="CircAreaSink", label=None
     ):
-        self.storeinput(inspect.currentframe())
         Element.__init__(
             self,
             model,
diff --git a/timflow/transient/element.py b/timflow/transient/element.py
index 255559c..f9f33ba 100644
--- a/timflow/transient/element.py
+++ b/timflow/transient/element.py
@@ -9,7 +9,6 @@ def initialize(self):
             pass
 """
 
-import inspect  # Used for storing the input
 from abc import ABC, abstractmethod
 
 import numpy as np
@@ -324,26 +323,6 @@ def headinside(self, t):
         print("This function not implemented for this element")
         return
 
-    def storeinput(self, frame):
-        self.inputargs, _, _, self.inputvalues = inspect.getargvalues(frame)
-
-    def write(self):
-        rv = self.name + "(" + self.model.modelname + ",\n"
-        for key in self.inputargs[2:]:  # The first two are ignored
-            if isinstance(self.inputvalues[key], np.ndarray):
-                rv += (
-                    key
-                    + " = "
-                    + np.array2string(self.inputvalues[key], separator=",")
-                    + ",\n"
-                )
-            elif isinstance(self.inputvalues[key], str):
-                rv += key + " = '" + self.inputvalues[key] + "',\n"
-            else:
-                rv += key + " = " + str(self.inputvalues[key]) + ",\n"
-        rv += ")\n"
-        return rv
-
     def run_after_solve(self):  # noqa: B027
         """Function to run after a solution is completed.
 
diff --git a/timflow/transient/linedoublet.py b/timflow/transient/linedoublet.py
index 29c2234..414fff5 100644
--- a/timflow/transient/linedoublet.py
+++ b/timflow/transient/linedoublet.py
@@ -7,7 +7,6 @@
     LeakyLineDoublet(ml, x1=-10, y1=0, x2=10, y2=0, layers=0)
 """
 
-import inspect  # Used for storing the input
 import warnings
 
 import matplotlib.pyplot as plt
@@ -243,7 +242,6 @@ def __init__(
         label=None,
         addtomodel=True,
     ):
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             x1=x1,
@@ -299,7 +297,6 @@ class LeakyWallString(Element, LeakyWallEquation):
     def __init__(
         self, model, xy=[(-1, 0), (1, 0)], res="imp", order=0, layers=0, label=None
     ):
-        self.storeinput(inspect.currentframe())
         Element.__init__(
             self,
             model,
diff --git a/timflow/transient/linesink.py b/timflow/transient/linesink.py
index fb74428..98bdfc5 100644
--- a/timflow/transient/linesink.py
+++ b/timflow/transient/linesink.py
@@ -8,7 +8,6 @@
     River(ml, x1=-10, y1=0, x2=10, y2=0, tsandbc=[(0, 5)], layers=0)
 """
 
-import inspect  # Used for storing the input
 import warnings
 
 import matplotlib.pyplot as plt
@@ -223,7 +222,6 @@ def __init__(
         addtomodel=True,
     ):
         """Initialize a transient line-sink with a specified discharge."""
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             x1=x1,
@@ -245,7 +243,6 @@ def __init__(
 #'''HeadLineSink that remains zero and constant through time'''
 # def __init__(self, model, x1=-1, y1=0, x2=1, y2=0, res=0.0, wh='H', \
 #             layers=0, label=None, addtomodel=True):
-#    self.storeinput(inspect.currentframe())
 #    LineSinkBase.__init__(self, model, x1=x1, y1=y1, x2=x2, y2=y2, \
 #                          tsandbc=[(0, 0)], res=res, wh=wh, layers=layers,\
 #                          type='z', name='ZeroHeadLineSink', label=label, \
@@ -325,7 +322,6 @@ def __init__(
         addtomodel=True,
     ):
         """Initialize a transient head-specified line-sink."""
-        self.storeinput(inspect.currentframe())
         if tsandh == "fixed":
             tsandh = [(0, 0)]
             etype = "z"
@@ -635,7 +631,6 @@ def __init__(
     ):
         """Initialize a transient multi-screen line-sink."""
         # assert len(layers) > 1, "number of layers must be at least 2"
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             x1=x1,
@@ -725,7 +720,6 @@ def __init__(
         label=None,
     ):
         """Initialize a transient ditch string with specified discharge."""
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             tsandbc=tsandQ,
@@ -820,7 +814,6 @@ def __init__(
         Astorage=None,
         label=None,
     ):
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             tsandbc=tsandQ,
@@ -1071,7 +1064,6 @@ def __init__(
         addtomodel=True,
     ):
         """Initialize a transient higher-order discharge line-sink."""
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             x1=x1,
@@ -1110,7 +1102,6 @@ def __init__(
         addtomodel=True,
     ):
         """Initialize a transient higher-order head-specified line-sink."""
-        self.storeinput(inspect.currentframe())
         if tsandh == "fixed":
             tsandh = [(0, 0)]
             etype = "z"
diff --git a/timflow/transient/model.py b/timflow/transient/model.py
index 9829ce1..acb806c 100644
--- a/timflow/transient/model.py
+++ b/timflow/transient/model.py
@@ -10,7 +10,7 @@
     ml.solve()
 """
 
-import inspect  # Used for storing the input
+import multiprocessing as mp
 from warnings import warn
 
 import numpy as np
@@ -27,6 +27,22 @@
 from timflow.version import check_tqdm_parallel
 
 
+def _compute_head_mp(args):
+    """Helper function for parallel computation of head_array."""
+    model, xi, yi, t, layers, i = args
+    return i, model.head(xi, yi, t, layers)
+
+
+def _compute_velocity_mp(args):
+    """Helper function for parallel computation of velocity_array."""
+    model, xi, yi, zi, t, i = args
+    try:
+        vv = model.velocomp(xi, yi, zi, t)
+    except (ZeroDivisionError, ValueError):
+        vv = np.full((3,), np.nan)
+    return i, vv
+
+
 class Model:
     def __init__(
         self,
@@ -80,8 +96,7 @@ def __init__(
             model3d,
         )
         self.compute_laplace_parameters()
-        self.name = "TimModel"
-        self.modelname = "ml"  # Used for writing out input
+        self.name = "Model"
         self.model_type = "transient"  # Model type for plotting and other purposes
         self.steady = steady
         if self.steady is not None:
@@ -446,6 +461,69 @@ def velocomp(self, x, y, z, t, aq=None, layer_ltype=None):
 
         return velo
 
+    def velocity_array(self, x, y, z, t, show_progress=True, parallel=False):
+        """Compute velocity for an array of points.
+
+        Parameters
+        ----------
+        x : 1d-array
+            x values
+        y : 1d-array
+            y values
+        z : 1d-array
+            z values
+        t : float
+            time at which velocity computed
+        show_progress : bool, optional
+            if `True`, shows progress bar when computing velocity grid, by default `True`
+        parallel : bool or int, optional
+            if `True`, computes velocity grid in parallel using multiprocessing,
+            by default `False`. If an integer is provided, it specifies the number of
+            processes to use.
+
+        Returns
+        -------
+        velocity : array
+            velocity vector (vx, vy, vz) at each point in grid,
+            size (3, len(x))
+        """
+        parallel, thread_map, tqdm = check_tqdm_parallel(parallel)
+
+        x = np.atleast_1d(x)
+        y = np.atleast_1d(y)
+        z = np.atleast_1d(z)
+        npts = len(x)
+        v = np.empty((3, npts))
+        if not parallel:
+            for i in (
+                tqdm(range(npts), desc="velocity array", disable=not show_progress)
+                if tqdm
+                else range(npts)
+            ):
+                try:
+                    vv = self.velocomp(x[i], y[i], z[i], t)
+                except ZeroDivisionError:
+                    vv = np.full((3,), np.nan)
+                v[:, i] = vv
+        else:
+            nproc = mp.cpu_count() if parallel is True else int(parallel)
+            chunksize = max(1, npts // (4 * nproc)) if nproc > 0 else 1
+            tasks = [(self, x[i], y[i], z[i], t, i) for i in range(npts)]
+            results = thread_map(
+                _compute_velocity_mp,
+                tasks,
+                total=npts,
+                desc="velocity array",
+                disable=not show_progress,
+                tqdm_class=tqdm,
+                max_workers=nproc,
+                chunksize=chunksize,
+            )
+            for i, result in results:
+                v[:, i] = result
+
+        return v
+
     def velocity_grid(self, xg, yg, zg, t, show_progress=True, parallel=False):
         """Compute velocity grid.
 
@@ -461,57 +539,30 @@ def velocity_grid(self, xg, yg, zg, t, show_progress=True, parallel=False):
             time for which grid is returned
         show_progress : bool, optional
             if `True`, shows progress bar when computing velocity grid, by default `True`
-        parallel : bool, optional
-            if `True`, computes velocity grid in parallel using multi threading,
-            by default `False`
+        parallel : bool or int, optional
+            if `True`, computes velocity grid in parallel using multiprocessing,
+            by default `False`. If an integer is provided, it specifies the number of
+            processes to use.
 
         Returns
         -------
         velocity : array
-            velocity vector (vz, vy, vx) at each point in grid,
+            velocity vector (vx, vy, vz) at each point in grid,
             size (3, len(z), len(y), len(x))
         """
-        parallel, thread_map, tqdm = check_tqdm_parallel(parallel)
-
-        def compute(kij):
-            k, i, j = kij
-            try:
-                vv = self.velocomp(xg[j], yg[i], zg[k], t)
-            except ZeroDivisionError:
-                vv = np.full((3,), np.nan)
-            return k, i, j, vv
-
         xg = np.atleast_1d(xg)
         yg = np.atleast_1d(yg)
         zg = np.atleast_1d(zg)
-        nz, ny, nx = len(zg), len(yg), len(xg)
-        v = np.empty((3, nz, ny, nx))
-        if not parallel:
-            for k in range(nz):
-                if show_progress:
-                    print(".", end="", flush=True)
-                for i in range(ny):
-                    for j in range(nx):
-                        try:
-                            vv = self.velocomp(xg[j], yg[i], zg[k], t)
-                        except ZeroDivisionError:
-                            vv = np.full((3,), np.nan)
-                        v[:, k, i, j] = vv
-            if show_progress:
-                print("", flush=True)
-        else:
-            results = thread_map(
-                compute,
-                [(k, i, j) for k in range(nz) for i in range(ny) for j in range(nx)],
-                total=nz * nx * ny,
-                desc="velocity grid",
-                disable=not show_progress,
-                tqdm_class=tqdm,
-            )
-            for k, i, j, result in results:
-                v[:, k, i, j] = result
-
-        return v
+        x, y, z = np.meshgrid(xg, yg, zg)
+        v = self.velocity_array(
+            x.ravel(),
+            y.ravel(),
+            z.ravel(),
+            t,
+            show_progress=show_progress,
+            parallel=parallel,
+        )
+        return v.reshape((3, len(yg), len(xg), len(zg))).transpose((0, 3, 1, 2))
 
     def velo_one(self, x, y, z, t, aq=None, layer_ltype=[0, 0]):
         # implemented for one layer and one time
@@ -591,8 +642,78 @@ def disvecalongline(self, x, y, t, layers=None):
             qx[:, :, i], qy[:, :, i] = self.disvec(xg[i], yg[i], t, layers)
         return qx, qy
 
+    def head_array(self, x, y, t, layers=None, show_progress=False, parallel=False):
+        """Head for array of points.
+
+        Parameters
+        ----------
+        x : 1D array or list
+            x values of points
+        y : 1D array or list
+            y values of points
+        t : float or 1D array or list
+            times for which grid is returned
+        layers : integer, list or array, optional
+            layers for which grid is returned
+        show_progress : bool
+            show computation progress, by printing dots per row or with tqdm progressbar
+            when parallel is True. Default is False.
+        parallel : bool or int, optional
+            if `True`, computes head_array in parallel using multiprocessing,
+            by default `False`. If an integer is provided, it specifies the number of
+            processes to use.
+
+        Returns
+        -------
+        h : array size `nlayers, ntimes, npoints`
+        """
+        parallel, process_map, tqdm = check_tqdm_parallel(parallel)
+        x = np.atleast_1d(x)
+        y = np.atleast_1d(y)
+        t = np.atleast_1d(t)
+        npts = len(x)
+        assert npts == len(y), "x and y must have the same length"
+        ntimes = len(t)
+        if layers is None:
+            nlayers = self.aq.find_aquifer_data(x[0], y[0]).naq
+        else:
+            nlayers = len(np.atleast_1d(layers))
+        h = np.empty((nlayers, ntimes, npts))
+        if not parallel:
+            for i in (
+                tqdm(range(npts), disable=not show_progress, desc="head array")
+                if tqdm
+                else range(npts)
+            ):
+                h[:, :, i] = self.head(x[i], y[i], t, layers)
+        else:
+            nproc = mp.cpu_count() if parallel is True else int(parallel)
+            chunksize = max(1, npts // (4 * nproc)) if nproc > 0 else 1
+            tasks = [(self, x[i], y[i], t, layers, i) for i in range(npts)]
+            results = process_map(
+                _compute_head_mp,
+                tasks,
+                total=npts,
+                desc="head array",
+                disable=not show_progress,
+                tqdm_class=tqdm,
+                max_workers=nproc,
+                chunksize=chunksize,
+            )
+
+            for i, result in results:
+                h[:, :, i] = result
+        return h
+
     def headgrid(
-        self, xg, yg, t, layers=None, printrow=False, show_progress=False, parallel=False
+        self,
+        xg,
+        yg,
+        t,
+        layers=None,
+        printrow=False,
+        show_progress=False,
+        parallel=False,
     ):
         """Grid of heads.
 
@@ -609,9 +730,10 @@ def headgrid(
         show_progress : bool
             show computation progress, by printing dots per row or with tqdm progressbar
             when parallel is True. Default is False.
-        parallel : bool, optional
+        parallel : bool or int, optional
             if `True`, computes headgrid in parallel using multithreading,
-            by default `False`
+            by default `False`. If an integer is provided, it specifies the number of
+            processes to use.
         printrow : bool, optional
 
             .. deprecated:: 0.2.0
@@ -633,44 +755,18 @@ def headgrid(
             )
             show_progress = printrow
 
-        parallel, thread_map, tqdm = check_tqdm_parallel(parallel)
-
-        xg = np.atleast_1d(xg)
-        yg = np.atleast_1d(yg)
-        t = np.atleast_1d(t)
         nx = len(xg)
         ny = len(yg)
-        ntimes = len(t)
-        if layers is None:
-            nlayers = self.aq.find_aquifer_data(xg[0], yg[0]).naq
-        else:
-            nlayers = len(np.atleast_1d(layers))
-        t = np.atleast_1d(t)
-        h = np.empty((nlayers, ntimes, ny, nx))
-        if not parallel:
-            for j in range(ny):
-                if show_progress:
-                    print(".", end="", flush=True)
-                for i in range(nx):
-                    h[:, :, j, i] = self.head(xg[i], yg[j], t, layers)
-        else:
-
-            def compute(ij):
-                i, j = ij
-                return i, j, self.head(xg[i], yg[j], t, layers)
-
-            results = thread_map(
-                compute,
-                [(i, j) for j in range(ny) for i in range(nx)],
-                total=nx * ny,
-                desc="headgrid",
-                disable=not show_progress,
-                tqdm_class=tqdm,
-            )
-
-            for i, j, result in results:
-                h[:, :, j, i] = result
-        return h
+        x, y = np.meshgrid(xg, yg)
+        h = self.head_array(
+            x.ravel(),
+            y.ravel(),
+            t,
+            layers=layers,
+            show_progress=show_progress,
+            parallel=parallel,
+        )
+        return h.reshape((h.shape[0], h.shape[1], ny, nx))
 
     def headgrid2(
         self,
@@ -701,9 +797,10 @@ def headgrid2(
         show_progress : bool
             show computation progress, by printing dots per row or with tqdm progressbar
             when parallel is True. Default is False.
-        parallel : bool, optional
-            if `True`, computes headgrid in parallel using multi threading,
-            by default `False`
+        parallel : bool or int, optional
+            if `True`, computes headgrid in parallel using multiprocessing,
+            by default `False`. If an integer is provided, it specifies the number of
+            processes to use.
         printrow : boolean, optional
 
             .. deprecated:: 0.2.0
@@ -754,7 +851,7 @@ def disvecgrid(
             show computation progress, by printing dots per row or with tqdm progressbar
             when parallel is True. Default is True.
         parallel : bool, optional
-            if `True`, computes discharge vector grid in parallel using multi threading,
+            if `True`, computes discharge vector grid in parallel using multiprocessing,
             by default `False`
 
         Returns
@@ -873,35 +970,6 @@ def solve(self, printmat=0, sendback=0, silent=False):
             return sol
         return
 
-    def storeinput(self, frame):
-        self.inputargs, _, _, self.inputvalues = inspect.getargvalues(frame)
-
-    def write(self):
-        rv = self.modelname + " = " + self.name + "(\n"
-        for key in self.inputargs[1:]:  # The first argument (self) is ignored
-            if isinstance(self.inputvalues[key], np.ndarray):
-                rv += (
-                    key
-                    + " = "
-                    + np.array2string(self.inputvalues[key], separator=",")
-                    + ",\n"
-                )
-            elif isinstance(self.inputvalues[key], str):
-                rv += key + " = '" + self.inputvalues[key] + "',\n"
-            else:
-                rv += key + " = " + str(self.inputvalues[key]) + ",\n"
-        rv += ")\n"
-        return rv
-
-    def writemodel(self, fname):
-        self.initialize()  # So that model can be written without solving first
-        f = open(fname, "w")
-        f.write("from timflow.transient import *\n")
-        f.write(self.write())
-        for e in self.elementlist:
-            f.write(e.write())
-        f.close()
-
     def aquifer_summary(self):
         """Return DataFrame with summary of aquifer(s) parameters in model.
 
@@ -995,7 +1063,6 @@ def __init__(
         M=10,
         steady=None,
     ):
-        self.storeinput(inspect.currentframe())
         if phreatictop is None:
             phreatictop = False
             if topboundary[:3] == "phr":
@@ -1111,7 +1178,6 @@ def __init__(
         steady=None,
     ):
         """Z must have the length of the number of layers + 1."""
-        self.storeinput(inspect.currentframe())
         if phreatictop is None:
             phreatictop = False
             if topboundary[:3] == "phr":
@@ -1211,8 +1277,7 @@ def __init__(
         self.M = M
         self.aq = SimpleAquifer(naq)
         self.compute_laplace_parameters()
-        self.name = "TimModel"
-        self.modelname = "ml"  # Used for writing out input
+        self.name = "ModelXsection"
         self.steady = steady
         if self.steady is not None:
             self.steady.solve()
diff --git a/timflow/transient/plots.py b/timflow/transient/plots.py
index 370fffb..9483ec0 100644
--- a/timflow/transient/plots.py
+++ b/timflow/transient/plots.py
@@ -229,9 +229,10 @@ def contour(
             if list of strings: use strings as names in legend
         return_contours : bool, optional
             if True, return list of contour sets for each contoured layer
-        parallel : bool, optional
-            if True, compute headgrid in parallel using multiple threads,
-            default is False
+        parallel : bool or int, optional
+            if True, compute headgrid in parallel using multiprocessing,
+            default is False. If int is provided, it is interpreted as the number of
+            processes to use.
         show_progress : bool, optional
             if True, show progress bar when computing headgrid in parallel,
             default is False.
@@ -373,7 +374,7 @@ def quiver_xy(
         figsize : tuple of 2 values
             size of figure
         parallel : bool
-            if True, compute velocity grid in parallel using multiple threads,
+            if True, compute velocity grid in parallel using multiprocessing,
             default is False
         **kwargs
             additional keyword arguments passed to ax.quiver()
@@ -411,7 +412,7 @@ def quiver_z(
         figsize : tuple of 2 values
             size of figure
         parallel : bool
-            if True, compute velocity grid in parallel using multiple threads,
+            if True, compute velocity grid in parallel using multiprocessing,
             default is False
         **kwargs
             additional keyword arguments passed to ax.quiver()
diff --git a/timflow/transient/well.py b/timflow/transient/well.py
index 45c201f..8449bb4 100644
--- a/timflow/transient/well.py
+++ b/timflow/transient/well.py
@@ -8,8 +8,6 @@
     Well(ml, xw=100, yw=200, tsandbc=[(0, 1000)], layers=[0])
 """
 
-import inspect  # Used for storing the input
-
 import matplotlib.pyplot as plt
 import numpy as np
 
@@ -258,7 +256,6 @@ class DischargeWell(WellBase):
     def __init__(
         self, model, xw=0, yw=0, tsandQ=[(0, 1)], rw=0.1, res=0, layers=0, label=None
     ):
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             xw,
@@ -327,7 +324,6 @@ def __init__(
         label=None,
     ):
         """Initialize a transient well with a specified discharge."""
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             xw,
@@ -406,7 +402,6 @@ def __init__(
         self, model, xw=0, yw=0, rw=0.1, tsandh=[(0, 1)], res=0, layers=0, label=None
     ):
         """Initialize a transient well with a specified head."""
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             xw,
@@ -443,7 +438,6 @@ def __init__(
         label=None,
         fp=None,
     ):
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             xw,
@@ -698,7 +692,6 @@ def __init__(
         label=None,
     ):
         """Initialize a transient string of wells with a specified discharge."""
-        self.storeinput(inspect.currentframe())
         super().__init__(
             model,
             xy,
diff --git a/timflow/version.py b/timflow/version.py
index deabe24..8a2deb8 100644
--- a/timflow/version.py
+++ b/timflow/version.py
@@ -45,8 +45,8 @@ def check_tqdm_parallel(parallel):
     -------
     parallel : bool
         Whether parallel processing was requested and can be used.
-    thread_map : function or None
-        The thread_map function from tqdm if parallel processing is available, else None.
+    process_map : function or None
+        The process_map function from tqdm if parallel processing is available, else None.
     tqdm : class or None
         The tqdm class from tqdm if parallel processing is available, else None.
     """
@@ -54,7 +54,7 @@ def check_tqdm_parallel(parallel):
         return parallel, None, None
     try:
         from tqdm import tqdm
-        from tqdm.contrib.concurrent import thread_map
+        from tqdm.contrib.concurrent import process_map
     except ImportError:
         warnings.warn(
             "Parallel requires 'tqdm'. Install 'timflow[parallel]' or 'tqdm' to"
@@ -63,6 +63,6 @@ def check_tqdm_parallel(parallel):
             stacklevel=2,
         )
         parallel = False
-        thread_map = None
+        process_map = None
         tqdm = None
-    return parallel, thread_map, tqdm
+    return parallel, process_map, tqdm