Note
Click here to download the full example code
FLOW: Richards: 1D: InversionΒΆ
The example shows an inversion of Richards equation in 1D with a heterogeneous hydraulic conductivity function.
The haverkamp model is used with the same parameters as Celia1990 the boundary and initial conditions are also the same. The simulation domain is 40cm deep and is run for an hour with an exponentially increasing time step that has a maximum of one minute. The general setup of the experiment is an infiltration front that advances downward through the model over time.
The model chosen is the saturated hydraulic conductivity inside the hydraulic conductivity function (using haverkamp). The initial model is chosen to be the background (1e-3 cm/s). The saturation data has 2% random Gaussian noise added.
The figure shows the recovered saturated hydraulic conductivity next to the true model. The other two figures show the saturation field for the entire simulation for the true and recovered models.
Rowan Cockett - 21/12/2016
Out:
SimPEG.InvProblem will set Regularization.mref to m0.
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver and solverOpts as the problem***
model has any nan: 0
============================ Inexact Gauss Newton ============================
# beta phi_d phi_m f |proj(x-g)-x| LS Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
0 1.73e+05 9.78e+04 0.00e+00 9.78e+04 1.46e+04 0
1 1.73e+05 9.12e+04 1.32e-02 9.35e+04 4.83e+03 0
2 1.73e+05 8.89e+04 2.33e-02 9.30e+04 1.46e+03 0 Skip BFGS
3 4.33e+04 8.83e+04 2.69e-02 8.95e+04 1.10e+04 0 Skip BFGS
4 4.33e+04 7.84e+04 1.10e-01 8.32e+04 8.46e+03 0
5 4.33e+04 7.07e+04 2.04e-01 7.95e+04 5.85e+03 0 Skip BFGS
6 1.08e+04 6.56e+04 2.80e-01 6.87e+04 1.25e+04 0 Skip BFGS
7 1.08e+04 5.12e+04 5.57e-01 5.72e+04 1.18e+04 0
8 1.08e+04 3.89e+04 8.44e-01 4.81e+04 1.19e+04 0 Skip BFGS
9 2.71e+03 2.94e+04 1.11e+00 3.24e+04 1.57e+04 0
10 2.71e+03 1.68e+04 1.58e+00 2.11e+04 1.46e+04 0
NewtonRoot stopped by maxIters (30). norm: 8.1784e-04
11 2.71e+03 7.13e+03 1.93e+00 1.23e+04 1.46e+04 0
12 6.76e+02 2.00e+03 2.25e+00 3.52e+03 7.74e+03 0
13 6.76e+02 8.63e+02 2.54e+00 2.58e+03 1.32e+03 0 Skip BFGS
14 6.76e+02 5.92e+02 2.77e+00 2.46e+03 1.47e+03 0 Skip BFGS
15 1.69e+02 5.50e+02 2.80e+00 1.02e+03 5.17e+02 0
16 1.69e+02 3.09e+02 3.47e+00 8.96e+02 5.11e+02 0
17 1.69e+02 2.93e+02 3.55e+00 8.92e+02 2.15e+02 0 Skip BFGS
18 4.23e+01 2.94e+02 3.54e+00 4.43e+02 2.37e+02 0
19 4.23e+01 2.41e+02 4.09e+00 4.14e+02 4.35e+02 0
20 4.23e+01 2.41e+02 4.07e+00 4.13e+02 1.35e+02 0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.3195e+00 <= tolF*(1+|f0|) = 9.7784e+03
1 : |xc-x_last| = 1.5320e-01 <= tolX*(1+|x0|) = 4.4688e+00
0 : |proj(x-g)-x| = 1.3482e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 1.3482e+02 <= 1e3*eps = 1.0000e-02
1 : maxIter = 20 <= iter = 20
------------------------- DONE! -------------------------
/Users/josephcapriotti/codes/simpeg/examples/09-flow/plot_inv_flow_richards_1D.py:154: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
plt.show()
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import discretize
from SimPEG import maps
from SimPEG import regularization
from SimPEG import data_misfit
from SimPEG import optimization
from SimPEG import inverse_problem
from SimPEG import directives
from SimPEG import inversion
from SimPEG.flow import richards
def run(plotIt=True):
M = discretize.TensorMesh([np.ones(40)], x0="N")
M.setCellGradBC("dirichlet")
# We will use the haverkamp empirical model with parameters from Celia1990
k_fun, theta_fun = richards.empirical.haverkamp(
M,
A=1.1750e06,
gamma=4.74,
alpha=1.6110e06,
theta_s=0.287,
theta_r=0.075,
beta=3.96,
)
# Here we are making saturated hydraulic conductivity
# an exponential mapping to the model (defined below)
k_fun.KsMap = maps.ExpMap(nP=M.nC)
# Setup the boundary and initial conditions
bc = np.array([-61.5, -20.7])
h = np.zeros(M.nC) + bc[0]
prob = richards.SimulationNDCellCentered(
M,
hydraulic_conductivity=k_fun,
water_retention=theta_fun,
boundary_conditions=bc,
initial_conditions=h,
do_newton=False,
method="mixed",
debug=False,
)
prob.time_steps = [(5, 25, 1.1), (60, 40)]
# Create the survey
locs = -np.arange(2, 38, 4.0).reshape(-1, 1)
times = np.arange(30, prob.time_mesh.vectorCCx[-1], 60)
rxSat = richards.receivers.Saturation(locs, times)
survey = richards.Survey([rxSat])
prob.survey = survey
# Create a simple model for Ks
Ks = 1e-3
mtrue = np.ones(M.nC) * np.log(Ks)
mtrue[15:20] = np.log(5e-2)
mtrue[20:35] = np.log(3e-3)
mtrue[35:40] = np.log(1e-2)
m0 = np.ones(M.nC) * np.log(Ks)
# Create some synthetic data and fields
relative = 0.02 # The standard deviation for the noise
Hs = prob.fields(mtrue)
data = prob.make_synthetic_data(
mtrue, relative_error=relative, f=Hs, add_noise=True
)
# Setup a pretty standard inversion
reg = regularization.Tikhonov(M, alpha_s=1e-1)
dmis = data_misfit.L2DataMisfit(simulation=prob, data=data)
opt = optimization.InexactGaussNewton(maxIter=20, maxIterCG=10)
invProb = inverse_problem.BaseInvProblem(dmis, reg, opt)
beta = directives.BetaSchedule(coolingFactor=4)
betaest = directives.BetaEstimate_ByEig(beta0_ratio=1e2)
target = directives.TargetMisfit()
dir_list = [beta, betaest, target]
inv = inversion.BaseInversion(invProb, directiveList=dir_list)
mopt = inv.run(m0)
Hs_opt = prob.fields(mopt)
if plotIt:
plt.figure(figsize=(14, 9))
ax = plt.subplot(121)
plt.semilogx(np.exp(np.c_[mopt, mtrue]), M.gridCC)
plt.xlabel("Saturated Hydraulic Conductivity, $K_s$")
plt.ylabel("Depth, cm")
plt.semilogx([10 ** -3.9] * len(locs), locs, "ro")
plt.legend(("$m_{rec}$", "$m_{true}$", "Data locations"), loc=4)
ax = plt.subplot(222)
mesh2d = discretize.TensorMesh([prob.time_mesh.hx / 60, prob.mesh.hx], "0N")
sats = [theta_fun(_) for _ in Hs]
clr = mesh2d.plotImage(np.c_[sats][1:, :], ax=ax)
cmap0 = matplotlib.cm.RdYlBu_r
clr[0].set_cmap(cmap0)
c = plt.colorbar(clr[0])
c.set_label("Saturation $\\theta$")
plt.xlabel("Time, minutes")
plt.ylabel("Depth, cm")
plt.title("True saturation over time")
ax = plt.subplot(224)
mesh2d = discretize.TensorMesh([prob.time_mesh.hx / 60, prob.mesh.hx], "0N")
sats = [theta_fun(_) for _ in Hs_opt]
clr = mesh2d.plotImage(np.c_[sats][1:, :], ax=ax)
cmap0 = matplotlib.cm.RdYlBu_r
clr[0].set_cmap(cmap0)
c = plt.colorbar(clr[0])
c.set_label("Saturation $\\theta$")
plt.xlabel("Time, minutes")
plt.ylabel("Depth, cm")
plt.title("Recovered saturation over time")
plt.tight_layout()
if __name__ == "__main__":
run()
plt.show()
Total running time of the script: ( 6 minutes 23.112 seconds)