Automate testing by generating a MODFLOW 6 model

[Huite]

Hans and Jan had run into a case where the leaky line doublet seemed to give different answers compared to MWELL and a checking model MODFLOW model. I tried to make a parametrized / reproducible example in flopy... and I quickly found out that @dbrakenhoff had already reported and fixed this, see #71.

Anyway, the script I'd made has most of the logic for converting a TimML ModelMaq instance to a Flopy Modflow6 model. Here's the implementation of a ModflowConverter, it doesn't require anything fancy in terms of dependencies, just numpy, flopy, and matplotlib:

Details

from collections import defaultdict
from typing import Tuple
import tempfile

import timml
import flopy
import numpy as np
from matplotlib.path import Path as MplPath


IntArray = np.ndarray
BoolArray = np.ndarray



def locate_points(points, xgrid, ygrid) -> Tuple[IntArray]:
    x, y = np.atleast_2d(points).T
    i_x = np.minimum(np.searchsorted(xgrid, x), xgrid.size - 1)
    i_y = np.minimum(np.searchsorted(ygrid, y), ygrid.size - 1)
    return i_x, i_y


def points_in_circle(points, center, radius) -> BoolArray:
    dist_squared = ((points - center) ** 2).sum(axis=1)
    return dist_squared < radius**2


def points_in_polygon(points, polygon_xy) -> BoolArray:
    return MplPath(polygon_xy).contains_points(points)


def rasterization_indices(xy, xgrid, ygrid) -> Tuple[IntArray]:
    """
    Vectorized implementation of Bresenham's line algorithm to burn lines into
    a 2D raster.

    Parameters
    ----------
    xy: np.ndarray of shape (n_vertex, 2, 2)
    xgrid: np.ndarray of shape (n_x)
    ygrid: np.ndarray of shape (n_y)

    Returns
    -------
    y_index: np.ndarray of int
    x_index: np.ndarray of int
    """

    def _bresenhamline_nslope(slope) -> np.ndarray:
        scale = np.amax(np.abs(slope), axis=1).reshape(-1, 1)
        zeroslope = (scale == 0).all(1)
        scale[zeroslope] = 1.0
        normalizedslope = slope.astype(float) / scale
        normalizedslope[zeroslope] = 0.0
        return normalizedslope

    def bresenham_lines(start, end) -> IntArray:
        max_iter = np.amax(np.amax(np.abs(end - start), axis=1))
        _, dim = start.shape
        nslope = _bresenhamline_nslope(end - start)

        # steps to iterate on
        stepseq = np.arange(0, max_iter + 1)
        stepmat = np.tile(stepseq, (dim, 1)).T

        # some hacks for broadcasting properly
        bline = start[:, np.newaxis, :] + nslope[:, np.newaxis, :] * stepmat

        # Approximate to nearest int
        return np.rint(bline).astype(int).reshape((-1, dim))

    i_x, i_y = locate_points(xy, xgrid, ygrid)
    start = np.column_stack((i_y[:-1], i_x[:-1]))
    end = np.column_stack((i_y[1:], i_x[1:]))
    return tuple(bresenham_lines(start=start, end=end).T)


def aquifer_z_data(aquifer) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    index = 1 - aquifer.ilap
    top = aquifer.z[index]
    bottom = aquifer.z[index + 1 :]
    thickness = aquifer.Hlayer[index:]
    return top, bottom, thickness


def index3d(i_layer, i_y, i_x):
    return (
        np.repeat(i_layer, i_x.size),
        np.tile(i_y, len(i_layer)),
        np.tile(i_x, len(i_layer)),
    )


class ModflowConverter:
    AQUIFERS = slice(0, None, 2)
    LEAKY_LAYERS = slice(1, None, 2)

    def __init__(self, timml_model, xmin, ymin, xmax, ymax, cellsize):
        """
        Convert a TimML model into a steady-state MODFLOW 6 model.

        The solutions of the analytic elements of TimML are infinite in the x-y
        plane while MODFLOW uses a discrete finite grid. This grid is defined by:

        * xmin
        * ymin
        * xmax
        * ymax
        * cellsize

        Additionally, when no boundaries are specified, MODFLOW will assume a
        no flow boundary on the edges of the grid domain. To get comparable
        results, the TimML model should generally be surrounded by e.g.
        HeadLineSinks.

        This converts the following TimML elements to their MODFLOW equivalent:

        * Aquifer and Polygon inhomogeneities:
          - z is converted into DIS top and bottom.
          - k and c and are converted into NPF k.
          - if topboundary="semi", hstar is converted into a GHB stage, and the
            first c is converted into a GHB conductance.
          - if N is given, it is converted into RCH.
        * AreaSink:
          - is converted to RCH.
        * Well:
          - is converted to WEL.
        * HeadWell
          - if resistance is provided: is converted to GHB.
          - else: is converted to CHD.
        * HeadLineSink:
          - if resistance is provided: is converted to GHB.
          - else: is converted to CHD.
        * LeakyLineDoublet:
          - is burned into the NPF k values, by reducing conductivity.
        * LeakyBuildingPit:
          - is treated like a PolygonInhom.
          - the resistance along the edges are burned into the NPF k values
            like for a LeakyLineDoublet.

        """
        self.timml_model = timml_model
        self.xmin = xmin
        self.xmax = xmax
        self.ymin = ymin
        self.ymax = ymax
        self.cellsize = cellsize
        self.cell_area = cellsize * cellsize
        self.x = np.arange(xmin, xmax, cellsize) + 0.5 * cellsize
        self.y = np.arange(ymin, ymax, cellsize) + 0.5 * cellsize
        self.yy, self.xx = [
            a.ravel() for a in np.meshgrid(self.y, self.x, indexing="ij")
        ]
        self.xy = np.column_stack((self.xx, self.yy))
        self.nx = self.x.size
        self.ny = self.y.size
        # MODFLOW 6 has 3D layers:
        self.nlayer = timml_model.aq.naq * 2 - 1
        self.shape2d = (self.ny, self.nx)
        self.shape = (self.nlayer, self.ny, self.nx)
        self._initialize_top_bottom()
        self._initialize_k()
        self._group_elements()
        # This is for aquifer and inhomogeneities top:
        self.top_recharge = np.zeros(self.shape)
        self.top_cond = np.full(self.shape, np.nan)
        self.top_head = np.full(self.shape, np.nan)
        # We'll store all other boundaries in these arrays:
        self.recharge = np.zeros(self.shape)
        self.bound_q = np.zeros(self.shape)
        self.bound_cond = np.full(self.shape, np.nan)
        self.bound_head = np.full(self.shape, np.nan)
        self.constant_head = np.full(self.shape, np.nan)

    def _group_elements(self):
        self.elements = defaultdict(list)
        for element in self.timml_model.elementlist + self.timml_model.aq.inhomlist:
            self.elements[type(element).__name__].append(element)

    def _initialize_top_bottom(self):
        top, bottom, thickness = aquifer_z_data(self.timml_model.aq)
        self.top = np.full(self.shape2d, top)
        self.bottom = np.broadcast_to(
            bottom[:, np.newaxis, np.newaxis], self.shape
        ).copy()
        self.thickness = np.broadcast_to(
            thickness[:, np.newaxis, np.newaxis], self.shape
        ).copy()

    def _initialize_k(self):
        self.k = np.zeros(self.shape)
        self.k[self.AQUIFERS] = self.timml_model.aq.kaq[:, np.newaxis, np.newaxis]
        if self.nlayer > 1:
            kll = (
                self.thickness[self.LEAKY_LAYERS]
                / self.timml_model.aq.c[1:, np.newaxis, np.newaxis]
            )
            self.k[self.LEAKY_LAYERS, ...] = kll

    def _convert_polygoninhomogeneity(self, inhom):
        polygon_xy = np.column_stack((inhom.x, inhom.y))
        top, bottom, thickness = aquifer_z_data(inhom)
        kll = thickness[self.LEAKY_LAYERS] / inhom.c[1:, np.newaxis, np.newaxis]
        xy_indices = points_in_polygon(self.xy, polygon_xy).reshape(self.shape2d)
        self.top[xy_indices] = top
        self.bottom[:, xy_indices] = bottom
        self.thickness[:, xy_indices] = thickness
        self.k[self.AQUIFERS, xy_indices] = inhom.kaq
        if self.nlayer > 1:
            self.k[self.LEAKY_LAYERS, xy_indices] = kll
        N = getattr(inhom, "N", None)
        if N is None:
            N = 0
        self.top_recharge[0, xy_indices] = N
        hstar = inhom.hstar
        cond = self.cell_area / inhom.c[0]
        if hstar is None:
            hstar = np.nan
            cond = np.nan
        self.top_head[0, xy_indices] = hstar
        self.top_cond[0, xy_indices] = cond
        return

    def convert_polygoninhomogenity(self):
        for inhom in self.elements["PolygonInhomMaq"]:
            self._convert_polygoninhomogeneity(inhom)

    def convert_circareasink(self):
        for circ in self.elements["CircAreaSink"]:
            xy_indices = points_in_circle(self.xy, (circ.xc, circ.yc), circ.R).reshape(
                self.shape2d
            )
            self.recharge[circ.layers * 2, xy_indices] = circ.N
        return

    def _burn_resistance(self, xy, layers, res):
        i_y, i_x = rasterization_indices(xy, self.x, self.y)
        i_layer = layers * 2
        i_layer, i_y, i_x = index3d(i_layer, i_y, i_x)
        k_grid = self.k[i_layer, i_y, i_x]
        effective_k = self.cellsize / (self.cellsize / k_grid + res)
        self.k[i_layer, i_y, i_x] = effective_k
        return

    def _convert_linedoublet(self, doublet):
        xy = [(doublet.x1, doublet.y1), (doublet.x2, doublet.y2)]
        self._burn_resistance(xy, doublet.layers, doublet.res)
        return

    def convert_linedoublet(self):
        for doublet in self.elements["LeakyLineDoublet"]:
            self._convert_linedoublet(doublet)
        return

    def convert_linedoubletstring(self):
        for doubletstring in self.elements["LeakyLineDoubletString"]:
            for doublet in doubletstring.ldlist:
                self._convert_linedoublet(doublet)
        return

    def convert_leakybuildingpit(self):
        for pit in self.elements["LeakyBuildingPit"]:
            self._convert_polygoninhomogeneity(pit)
            xy = np.column_stack((pit.x, pit.y))
            self._burn_resistance(xy, pit.layers, pit.res)
        return

    def _convert_linesink(self, sink):
        xy = [(sink.x1, sink.y1), (sink.x2, sink.y2)]
        i_y, i_x = rasterization_indices(xy, self.x, self.y)
        i_layer = sink.layers * 2
        i_layer, i_y, i_x = index3d(i_layer, i_y, i_x)
        # FUTURE: interpolate head along the length of the linesink
        if sink.res > 0:
            self.bound_cond[i_layer, i_y, i_x] = sink.wh / sink.res
            self.bound_head[i_layer, i_y, i_x] = sink.hls[0]
        else:
            self.constant_head[i_layer, i_y, i_x] = sink.hls[0]
        return

    def convert_linesink(self):
        for sink in self.elements["HeadLineSink"]:
            self._convert_linesink(sink)
        return

    def convert_linesinkstring(self):
        for sinkstring in self.elements["HeadLineSinkString"]:
            for sink in sinkstring.lslist:
                self._convert_linesink(sink)
        return

    def convert_well(self):
        for well in self.elements["Well"]:
            i_x, i_y = locate_points([well.xw, well.yw], self.x, self.y)
            i_layer = well.layers * 2
            self.bound_q[i_layer, i_y, i_x] = -well.Qw
        return

    def convert_headwell(self):
        for well in self.elements["HeadWell"]:
            i_x, i_y = locate_points([well.xw, well.yw], self.x, self.y)
            i_layer = well.layers * 2
            if well.resfac > 0:
                self.bound_head[i_layer, i_y, i_x] = well.hc
                self.bound_cond[i_layer, i_y, i_x] = 1.0 / well.resfac
            else:
                self.constant_head[i_layer, i_y, i_x] = well.hc
        return

    def convert(self):
        self.convert_polygoninhomogenity()
        self.convert_leakybuildingpit()
        self.convert_linesink()
        self.convert_linesinkstring()
        self.convert_linedoublet()
        self.convert_linedoubletstring()
        self.convert_headwell()
        self.convert_well()
        self.convert_circareasink()
        return

    def to_mf6(self):
        def cell_ids(valid):
            return [tuple(i) for i in np.argwhere(valid)]

        constant_head = self.constant_head
        top_recharge = self.top_recharge
        top_head = self.top_head
        top_cond = self.top_cond
        recharge = self.recharge
        bound_q = self.bound_q
        bound_head = self.bound_head
        bound_cond = self.bound_cond
        nlay, nrow, ncol = self.shape

        workspace = tempfile.TemporaryDirectory()
        print(workspace)
        sim = flopy.mf6.MFSimulation(
            # sim_name="TimmlConverted", sim_ws=".", exe_name="mf6")
            sim_name="TimmlConverted",
            sim_ws=workspace.name,
            exe_name="mf6",
        )

        flopy.mf6.ModflowTdis(sim, nper=1)
        flopy.mf6.ModflowIms(
            sim,
            outer_dvclose=1.0e-5,
            outer_maximum=3,  # Linear problem, should be able to solve in one cg-solve...
            inner_maximum=200,
            inner_dvclose=1.0e-6,
        )
        gwf = flopy.mf6.ModflowGwf(sim, modelname="gwf", save_flows=True)
        flopy.mf6.ModflowGwfdis(
            gwf,
            nlay=nlay,
            nrow=nrow,
            ncol=ncol,
            delr=self.cellsize,
            delc=self.cellsize,
            top=self.top,
            botm=self.bottom,
        )
        flopy.mf6.ModflowGwfnpf(
            gwf,
            icelltype=0,  # confined
            k=self.k,
            k33=self.k,
        )
        flopy.mf6.ModflowGwfic(gwf, strt=0.0)

        # Setup recharge for the top
        valid_rch = top_recharge != 0
        if valid_rch.any():
            rch_top_stress = {0: [(i, top_recharge[i]) for i in cell_ids(valid_rch)]}
            flopy.mf6.ModflowGwfrch(gwf, stress_period_data=rch_top_stress)

        # Setup general head boundary for the top
        valid_ghb = ~np.isnan(top_cond)
        if valid_ghb.any():
            ghb_top_stress = {
                0: [(i, top_head[i], top_cond[i]) for i in cell_ids(valid_ghb)]
            }
            flopy.mf6.ModflowGwfghb(gwf, stress_period_data=ghb_top_stress)

        # Setup constant head
        valid_chd = ~np.isnan(constant_head)
        if valid_chd.any():
            chd_stress = {0: [(i, constant_head[i]) for i in cell_ids(valid_chd)]}
            flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_stress)

        # Setup recharge
        valid_rch = recharge != 0
        if valid_rch.any():
            rch_stress = {0: [(i, recharge[i]) for i in cell_ids(valid_rch)]}
            flopy.mf6.ModflowGwfrch(gwf, stress_period_data=rch_stress)

        # Setup general head boundary
        valid_ghb = ~np.isnan(bound_cond)
        if valid_ghb.any():
            ghb_stress = {
                0: [(i, bound_head[i], bound_cond[i]) for i in cell_ids(valid_ghb)]
            }
            flopy.mf6.ModflowGwfghb(gwf, stress_period_data=ghb_stress)

        # Setup wells for fixed discharge boundaries
        valid_wel = bound_q != 0
        if valid_wel.any():
            wel_stress = {0: [(tuple(i), bound_q[i]) for i in cell_ids(valid_wel)]}
            flopy.mf6.ModflowGwfwel(gwf, stress_period_data=wel_stress)

        flopy.mf6.ModflowGwfoc(
            gwf,
            head_filerecord="gwf.hds",
            saverecord=[("HEAD", "ALL")],
        )
        return sim, gwf

<\details>

[Huite]

With this and a running + plotting utility, it becomes quite easy to make automated comparison calculations:

Details

def compare(timml_model, xmin, ymin, xmax, ymax, cellsize):
    timml_model.solve()
    converter = ModflowConverter(
        timml_model,
        xmin,
        ymin,
        xmax,
        ymax,
        cellsize=cellsize,
    )
    converter.convert()
    sim, gwf = converter.to_mf6()
    sim.write_simulation()
    sim.run_simulation()
    mf_head = gwf.output.head().get_alldata()[0]
    tim_head = model.headgrid(xg=converter.x, yg=converter.y)

    _, (ax0, ax1) = plt.subplots(ncols=2, figsize=(12, 5))
    levels = np.linspace(mf_head.min(), mf_head.max(), 11)
    cs = ax0.contour(converter.x, converter.y, mf_head[0], levels=levels)
    ax0.clabel(cs, cs.levels)
    ax0.contour(
        converter.x, converter.y, tim_head[0], levels=levels, linestyles="dashed"
    )
    ax0.set_aspect(1.0)
    ax0.set_ylabel("y (m)")
    ax0.set_xlabel("x (m)")
    ax0.set_title(
        "Isohypses, first layer, plan view.\nSolid: MODFLOW 6, dashed: TimML."
    )
    ax0.legend()
    i_y_mid = converter.ny // 2
    ax1.plot(converter.x, mf_head[0, i_y_mid], label="MODFLOW 6")
    ax1.plot(converter.x, tim_head[0, i_y_mid], linestyle="dashed", label="TimML")
    ax1.set_ylabel("head (m)")
    ax1.set_xlabel("x (m)")
    ax1.set_title(
        f"Head, first layer, along y={converter.y[i_y_mid]}\n"
        "Solid: MODFLOW 6, dashed: TimML."
    )
    return

# Setup some geometry coordinates

xmin = -1000.0
xmax = 1000.0
ymin = -1000.0
ymax = 1000.0
cellsize = 10.0

xy_left = np.array(
    [
        [xmin, ymin - 10 * y_interval],
        [xmin, ymin],
        [xmin, ymax],
        [xmin, ymax + 10 * y_interval],
    ]
)
xy_right = np.array(
    [
        [xmax, ymin - 10 * y_interval],
        [xmax, ymin],
        [xmax, ymax],
        [xmax, ymax + 10 * y_interval],
    ]
)

hls_xy = np.array(
    [
        [xmin, ymin],
        [xmax, ymin],
        [xmax, ymax],
        [xmin, ymax],
        [xmin, ymin],
    ]
)

poly_xy = np.array(
    [
        [-250.0, -250.0],
        [250.0, -250.0],
        [250.0, 250.0],
        [-250.0, 250.0],
        [-250.0, -250.0],
    ]
)

model = timml.ModelMaq(kaq=10.0, z=[10.0, 0.0])
timml.HeadLineSinkString(model, hls_xy, hls=1.0, order=4)
timml.HeadLineSink(model, x1=-250.0, y1=-250.0, x2=250.0, y2=250.0, hls=0.0, order=3)
compare(model, xmin, ymin, xmax, ymax, cellsize)

model = timml.ModelMaq(kaq=10.0, z=[10.0, 0.0])
timml.HeadLineSinkString(model, hls_xy, hls=0.0, order=4)
timml.CircAreaSink(model, xc=0.0, yc=0.0, R=100.0, N=0.001)
compare(model, xmin, ymin, xmax, ymax, cellsize)

model = timml.ModelMaq(kaq=10.0, z=[10.0, 0.0])
timml.HeadLineSinkString(model, hls_xy, hls=1.0, order=4)
timml.LeakyBuildingPit(
    model,
    poly_xy,
    kaq=10.0,
    z=[11.0, 10.0, 0.0],
    hstar=0.0,
    c=1_000.0,
    res=100.0,
    topboundary="semi",
)
compare(model, xmin, ymin, xmax, ymax, cellsize)

[Huite]

I think this could be a useful addition to testing, apart from straightforward analytical solution (such as in: https://github.com/mbakker7/timml/blob/master/notebooks/circular_buildingpit.ipynb)

I could make a PR implementing such testing, which would run through dozens of tests, but some difficulty lies in deciding what the right testing criteria are. E.g. this is the outcome for the leaky line doublet with the incorrect resfac:

The image clearly shows that something is wrong, but what is the right measure? The mean head difference is ~0 due to symmetry. The effect on standard deviation is "diluted".

For comparison, this is with current master:

The image shows that both models now agree, but the standard deviation has only improved by a factor ~4. Taking the heads along y=5, we get an std of 0.058, but it remains 0.013 in current master. They'll never get exactly the same result:

• MODFLOW is discretized: in the above example, the head difference for MODFLOW is spread out because I've burned the resistance of the line doublet into the k field (because horizontal flow barriers are a pain to setup)
• MODFLOW has a finite domain: in the above example, I force flow from left to right in the TimML model using much larger linesinks. In MODFLOW the "northern" and "southern" boundaries are closed.
• TimML shows some head differences when the linestrings are a bit long (see corners in the previous comment).

We can still guard against regressions relatively easily of course: run all test cases with plots, eyeball whether the plots look okay, and use those values for mean differences and standard deviation. If we then get something wrong in a subsequent commit, it should trigger a test failure.

[Huite]

The above obviously requires both flopy and MODFLOW 6, which rather complicate setting up the testing on the CI. Fortunately, we use so little of MODFLOW 6 functionality, that the relevant functionality can be implemented in less than 200 lines:

class HeadBoundary:
    def __init__(self, head, conductance):
        head = head.ravel().copy()
        conductance = conductance.ravel().copy()
        self.head = head
        self.head[np.isnan(head)] = 0.0
        self.conductance = conductance
        self.conductance[np.isnan(conductance)] = 0.0

    def formulate(self, model):
        model.hcof -= self.conductance
        model.b -= self.conductance * self.head
 

class Well:
    def __init__(self, rate):
        rate = rate.ravel().copy()
        self.rate = rate
        self.rate[np.isnan(rate)] = 0.0

    def formulate(self, model):
        model.b -= self.rate


class Recharge:
    def __init__(self, rate):
        rate = rate.ravel().copy()
        self.rate = rate
        self.rate[np.isnan(rate)] = 0.0

    def formulate(self, model):
        model.b -= self.rate * model.cell_area
 

class FiniteVolumeGroundwaterFlow:
    def __init__(
        self,
        cellsize,
        thickness,
        k,
        constant_head,
        packages,
    ):
        self.cellsize = cellsize
        self.cell_area = cellsize * cellsize
        self.thickness = thickness.ravel()
        self.k = k.ravel()
        self.active = self.k > 0
        self.inactive = ~self.active
        self.shape = k.shape
        self.nlayer, self.ny, self.nx = self.shape
        self.n = np.product(self.shape)
        self.nodes = np.arange(self.n).reshape(self.shape)
        self.constant_head = constant_head.ravel()
        self.packages = packages

    def _horizontal_connectivity(self):
        left = self.nodes[:, :, :-1].ravel()
        right = self.nodes[:, :, 1:].ravel()
        front = self.nodes[:, :-1, :].ravel()
        back = self.nodes[:, 1:, :].ravel()
        i = np.concatenate((left, right, front, back))
        j = np.concatenate((right, left, back, front))
        # Compute transmissivities
        T_i = self.k[i] * self.thickness[i]
        T_j = self.k[j] * self.thickness[j]
        T_product = T_i * T_j
        width = self.cellsize
        length = self.cellsize * 0.5
        # Compute harmonic mean conductance
        conductance = np.zeros_like(T_product)
        nonzero = T_product > 0
        conductance[nonzero] = (
            width * T_product[nonzero] / (T_i[nonzero] * length + T_j[nonzero] * length)
        )
        return i, j, conductance

    def _vertical_connectivity(self):
        upper = self.nodes[:-1, :, :].ravel()
        lower = self.nodes[1:, :, :].ravel()
        i = np.concatenate((upper, lower))
        j = np.concatenate((lower, upper))
        c_i = 0.5 * self.thickness[i] / self.k[i]
        c_j = 0.5 * self.thickness[j] / self.k[j]
        conductance = self.cell_area / (c_i + c_j)
        return i, j, conductance

    def _build_matrix(self):
        # Cell to cell connections
        h_i, h_j, h_cond = self._horizontal_connectivity()
        v_i, v_j, v_cond = self._vertical_connectivity()
        i = np.concatenate((h_i, v_i))
        j = np.concatenate((h_j, v_j))
        v = np.concatenate((h_cond, v_cond))

        # Remove constant head cells and inactive cells -- this results in an
        # unsymmetric matrix.
        is_variable = np.isnan(self.constant_head) & self.active
        i_variable = is_variable[i]
        # Also include nodes to make room for the diagonal entries.
        nodes = self.nodes.ravel()
        i = np.concatenate((i[i_variable], nodes))
        j = np.concatenate((j[i_variable], nodes))
        v = np.concatenate((v[i_variable], np.zeros(self.n)))
        self.A = scipy.sparse.csr_matrix((v, (i, j)), shape=(self.n, self.n))

    def solve(self):
        self._build_matrix()
        self.hcof = np.zeros(self.n)
        self.b = np.zeros(self.n)

        for package in self.packages:
            package.formulate(self)

        # Compute the diagonal values of the matrix.
        diagonal = self.hcof - np.squeeze(np.asarray(self.A.sum(axis=1)))
        # Set constant head values in the diagonal and rhs.
        is_constant = ~np.isnan(self.constant_head)
        diagonal[is_constant] = 1.0
        self.b[is_constant] = self.constant_head[is_constant]
        # Set inactive cells in the diagonal and rhs; 0.0 is a dummy value
        # which we'll replace by NaN later on, since introducing NaN here
        # pollutes them everywhere in the solve step.
        diagonal[self.inactive] = 1.0
        self.b[self.inactive] = 0.0
        self.A.setdiag(diagonal)
        # Do a direct linear solve
        self.h = scipy.sparse.linalg.spsolve(self.A, self.b)
        self.h[self.inactive] = np.nan

    @property
    def head(self):
        return self.h.reshape(self.shape)

We can add another method instead of to_mf6 to the ModelConverter:

    def to_finite_volume(self):
        return FiniteVolumeGroundwaterFlow(
            cellsize=self.cellsize,
            thickness=self.thickness,
            k=self.k,
            constant_head=self.constant_head,
            packages = (
                Recharge(self.top_recharge),
                Recharge(self.recharge),
                HeadBoundary(self.top_head, self.top_cond),
                HeadBoundary(self.bound_head, self.bound_cond),
                Well(self.bound_q),
            )
        )

And check it using this model instead:

(It's identical to the MODFLOW 6 result, as you'd expect:)

[Huite]

I think my current suggestion would be something along these lines.

We introduce a testing_utils namespace. Ideally this lives only in the testing directory, so that people cannot use it from TimML: it's fragile when used outside of a testing context with simple models and will not give the same results unless something simulating a finite boundary is included in the TimML model.

This testing_utils contains:

• a model_converter.py module
• and a finite_volume.py module

The converter contains two methods:

• to_finite_volume()
• to_mf6()

Then, we create a test-suite which tests simple models (as shown in the comments above). These simple models generally contain just a constant head around the domain, with some element in the middle. We test these models during continuous integration using the Python finite volume implementation, since this only requires numpy and scipy which are required anyway.

We keep the to_mf6 method so if a test fails or something unexpected happens, it can be reproduced locally using MODFLOW 6. This has the advantage that we can say that we do indeed cross-check correctness with MODFLOW 6 (and not just some guy's (me) finite volume implementation).

That's my thoughts so far, might be nice to discuss it during the next meeting.