H process: Theis solution (Pumping well)

import os
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import ogstools as ot
from scipy.special import exp1

def calc_u(r, S, T, t):
    """Calculate and return the dimensionless time parameter, u."""
    return r**2 * S / 4 / T / t


def theis_drawdown(t, S, T, Q, r):
    """Calculate and return the drawdown s(r,t) for parameters S, T.

    This version uses the Theis equation, s(r,t) = Q * W(u) / (4.pi.T),
    where W(u) is the Well function for u = Sr^2 / (4Tt).
    S is the aquifer storage coefficient,
    T is the transmissivity (m2/day),
    r is the distance from the well (m), and
    Q is the pumping rate (m3/day).
    """
    u = calc_u(r, S, T, t)
    return Q / 4 / np.pi / T * exp1(u)

Q = 2000  # Pumping rate from well (m3/day)
r = 10  # Distance from well (m)

# Time grid, days.
t = np.array([1, 2, 4, 8, 12, 16, 20, 30, 40, 50, 60, 70, 80, 90, 100])

# Calculate some synthetic data to fit.
S, T = 0.0003, 1000
s = theis_drawdown(t, S, T, Q, r)

# Plot the data
fig, ax = plt.subplots(figsize=(16, 10))

ax.plot(t, s, label=f"r = {r} m")
ax.set(
    xlim=(0, 100),
    ylim=(1.6, 2.6),
    title="Theis: Analytical solution",
    xlabel=r"$t\;/\;\mathrm{days}$",
    ylabel=r"$s\;/\;\mathrm{m}$",
)
ax.legend()
ax.grid()
ot.plot.utils.update_font_sizes(ax)

# Recalculation from days in sec
Q = 0.016  # Pumping rate from well (m3/s)
t = 864000  # Time in s.

# Distance from well (m)
r = np.arange(1, 41, 1)

# Calculate some synthetic data to fit.
S = 0.001
T = 9.2903e-4
u = calc_u(r, S, T, t)
s = theis_drawdown(t, S, T, Q, r)
s = s - 5  # reference head

# Plot the data
fig, ax = plt.subplots(figsize=(16, 10))
ax.plot(r, s, label=f"t = {t} days")
ax.set(
    xlim=(1, 40),
    title="Theis: Analytical solution",
    xlabel=r"$r\;/\mathrm{m}$",
    ylabel=r"$hydraulic head\;/\;\mathrm{m}$",
)
ax.legend()
ax.grid()
ot.plot.utils.update_font_sizes(ax)

fig = ot.plot.contourf(mesh, "OGS5_pressure")
fig.get_axes()[0].set_title("Initial pressure")
ot.plot.utils.update_font_sizes(fig.get_axes())

out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
out_dir.mkdir(parents=True, exist_ok=True)

model = ot.Project(input_file="axisym_theis.prj", output_file="axisym_theis.prj")
model.run_model(logfile=out_dir / "out.log", args=f"-o {out_dir}")

%%capture --no-display

xaxis = np.column_stack((np.linspace(1.0, 40, 40), np.zeros((40, 2))))
pressure = ot.variables.pressure.replace(
    data_unit="m", output_unit="m", output_name="hydraulic head", symbol=""
)

ms_probe = ot.MeshSeries.extract_probe(ms, xaxis)

labels = [f"$t={np.round(x, 2)}s$" for x in ms_probe[1:].timevalues]
fig, ax = plt.subplots(figsize=(18, 12))

ot.plot.line(ms_probe[1:], "x", pressure, labels=labels, ax=ax)
ot.plot.line(
    ms_probe[ms_probe.closest_timestep(1728)],
    "x",
    "OGS5_pressure",
    marker="x",
    markersize=8,
    linestyle="",
    label="OGS, $t=1728 s$",
    ax=ax,
)
fig.tight_layout()

time = 864000
timestep = ms_probe.closest_timestep(time)

ana_sol = np.zeros_like(ms_probe.point_data[pressure])
ana_sol[timestep] = s

mask = np.nonzero(ana_sol)
abs_error = np.zeros_like(ms_probe.point_data[pressure])
abs_error[mask] = np.abs(ana_sol[mask] - ms_probe.point_data[pressure][mask])

rel_error = np.zeros_like(ms_probe.point_data[pressure])
rel_error[mask] = np.abs(abs_error[mask] / (ana_sol[mask]))

np.testing.assert_array_less(abs_error, 0.9)
np.testing.assert_array_less(rel_error, 0.07)

ms_probe.point_data[pressure.anasol.data_name] = ana_sol
ms_probe.point_data[pressure.abs_error.data_name] = abs_error
ms_probe.point_data[pressure.rel_error.data_name] = rel_error

fig, ax = plt.subplots(ncols=3, figsize=(40, 12))

ot.plot.line(
    ms_probe[timestep],
    "x",
    pressure,
    marker="x",
    markersize=8,
    linestyle="",
    label="numerical solution (ogs6)",
    ax=ax[0],
)
ot.plot.line(
    ms_probe[timestep],
    "x",
    pressure.anasol,
    marker="o",
    label="analytical solution",
    ax=ax[0],
    color="C1",
)
ot.plot.line(ms_probe[timestep], "x", pressure.abs_error, marker="x", ax=ax[1])
ot.plot.line(ms_probe[timestep], "x", pressure.rel_error, marker="x", ax=ax[2])

fig.suptitle("Theis: Comparison analytical and numerical solution", fontsize=38)
fig.tight_layout()