Add rhs argument to enable "poisson interpolation" in laplace_interpo…

[Huite]

This would allow us to do fancy downscaling, basically getting rid of respighi: instead of MODFLOW, we can directly call this function.

[JoerivanEngelen]

Some screenshots for context: The poisson interpolation allows interpolating ditches with recharge.

This would looks something like this for ditches:

Where purple line is the ditch level (boundary condition), blue line is the laplacian interpolation, and red line the poisson interpolation. This is an example for negative recharge AFAIK (was still work in progress and the sign got swapped accidentily, but it was the best screenshot I have to explain the context somewhat).

After fixing the sign errors this turned into:

[sonarqubecloud]

Quality Gate passed

Issues
1 New issue
0 Accepted issues

Measures
0 Security Hotspots
0.0% Coverage on New Code
0.0% Duplication on New Code

See analysis details on SonarCloud

[Huite]

After some reflection: the respighi scaling method requires us to use "weak" boundaries as well (e.g. RIV or GHB).

With only hard boundaries, if you downscale at 25.0 m, the value at 25.0 m gets burned into the grid directly. This means that if the ditch distance is <50 m or so, you will only get the boundary conditions (like ditch stages). This would mean that for e.g. the peaty areas of South-Holland, the scaling essentially fails.

It's relatively easy to add a robin boundary (a GHB) as optional arguments for the system of equations.

This allows the interpolated head to deviate somewhat from the provided robin values.

[Huite]

Few quick notes:

#1239 removed the pcg module.

However, the same method can be easily applied.
We should add robin_coeficient and robin_value for the robin boundary conditions, as well as neumann_value for the neumann boundary condition.

In terms of the _interpolate function:

def _interpolate(
    arr: np.ndarray,
    connectivity: sparse.csr_matrix,
    direct: bool,
    delta: float,
    relax: float,
    rtol: float,
    atol: float,
    maxiter: int,
    robin_coefficient: FloatArray | None = None,
    robin_value: FloatArray | None = None,
    neumann_value: FloatArray | None = None,
):
    ar1d = arr.ravel()
    unknown = np.isnan(ar1d)
    known = ~unknown

    if robin_coefficient is None:
        robin_coefficient = 0.0
    if robin_value is None:
        robin_value = 0.0
    if neumann_value is None:
        neumann_value = 0.0

    # Set up system of equations.
    matrix = connectivity.copy()
    matrix.setdiag(-matrix.sum(axis=1).A[:, 0] - robin_coefficient
    rhs = -matrix[:, known].dot(ar1d[known]) - (robin_coefficient * robin_value)

We can make it a bit nicer and support multiple boundary conditions in a single cell by allowing an additional dimension for these, either by using np.add.at or by simply looping over the extra dimension.

It would also be good to support non-equidistant and unstructured grids (by taking the centroid-to-centroid distance into account for the coefficient matrix), as xugrid can do. For unstructured grids, it's a simple case of taking the connectivity from the Ugrid2d topology instead of generating it for the cartesian neighbors.

Might be easiest to scrap this PR and start from current master.

[JoerivanEngelen]

Will be followed up upon in #1470