Add support for neumann, robin boundary conditions in laplace_interpolate #1470

Huite · 2025-03-03T21:01:31Z

Supersedes #1106

Here's some basic illustration:

# %%
import numpy as np
import xarray as xr

import imod

# %%

da = xr.DataArray(np.full((200, 100), np.nan), dims=("y", "x"))
da.data[:, 0] = -10.0
da.data[:, -1] = 10.0

# %%

rate = xr.ones_like(da) * 0.01
systemdim = xr.DataArray([0.0, 0.0], dims=("system",))
head = systemdim * xr.full_like(da, np.nan)
cond = systemdim * xr.full_like(da, np.nan)
head[0, :, 25] = 20.0
head[0, :, 75] = 20.0
cond[0, :, 25] = 1.0
cond[0, :, 75] = 1.0
filled = imod.prepare.laplace_interpolate(
    source=da,
    neumann_value=rate,
    robin_coefficient=cond,
    robin_value=head,
    direct=True,
)
filled.plot()
# %%
filled.isel(y=100).plot()
# %%

filled = imod.prepare.laplace_interpolate(source=da, direct=True)
filled.plot()

# %%

rate = xr.ones_like(da) * 0.01
filled = imod.prepare.laplace_interpolate(source=da, neumann_value=rate, direct=True)
filled.isel(y=50).plot()

# %%

head = xr.full_like(da, np.nan)
cond = xr.full_like(da, np.nan)
head.data[:, 50] = 0.0
cond.data[:, 50] = 1.0
filled = imod.prepare.laplace_interpolate(
    source=da,
    neumann_value=rate,
    robin_coefficient=cond,
    robin_value=head,
    direct=True,
)
filled.isel(y=50).plot()
# %%

I haven't figured out how to elegantly broadcast the Ugrid 2d face to face connectivity, but it shouldn't be too hard.
The broadcasting logic might be nice for xugrid. In principle having access to neumann and robin boundary conditions might be worthwhile for xugrid too, but it's arguably a bit niche; mostly groundwater related (which is why imod is a good fit).

What might be worth adding is spatial distances, i.e. a cell to cell conductance. Although this assumes all interpolation dimension are commensurate; for x and y this is true, for "layer" this is not.
Maybe interpolating over more than three dimensions should be disallowed anyway?

3D interpolation is still easy; just swap dims to get a z coordinate, and you'd have a commensurate coordinate for x and y.
Similarly for xugrid, I suppose...