forked from simpeg/simpeg
/
plot_richards_celia1990.py
107 lines (87 loc) · 3.22 KB
/
plot_richards_celia1990.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
"""
FLOW: Richards: 1D: Celia1990
=============================
There are two different forms of Richards equation that differ
on how they deal with the non-linearity in the time-stepping term.
The most fundamental form, referred to as the
'mixed'-form of Richards Equation Celia1990_
.. math::
\\frac{\partial \\theta(\psi)}{\partial t} -
\\nabla \cdot k(\psi) \\nabla \psi -
\\frac{\partial k(\psi)}{\partial z} = 0
\quad \psi \in \Omega
where \\\\(\\\\theta\\\\) is water content, and \\\\(\\\\psi\\\\)
is pressure head. This formulation of Richards equation is called the
'mixed'-form because the equation is parameterized in \\\\(\\\\psi\\\\)
but the time-stepping is in terms of \\\\(\\\\theta\\\\).
As noted in Celia1990_ the 'head'-based form of Richards
equation can be written in the continuous form as:
.. math::
\\frac{\partial \\theta}{\partial \psi}
\\frac{\partial \psi}{\partial t} -
\\nabla \cdot k(\psi) \\nabla \psi -
\\frac{\partial k(\psi)}{\partial z} = 0
\quad \psi \in \Omega
However, it can be shown that this does not conserve mass in the
discrete formulation.
Here we reproduce the results from Celia1990_ demonstrating the
head-based formulation and the mixed-formulation.
.. _Celia1990: http://www.webpages.uidaho.edu/ch/papers/Celia.pdf
"""
import matplotlib.pyplot as plt
import numpy as np
from SimPEG import Mesh, Maps
from SimPEG.FLOW import Richards
def run(plotIt=True):
M = Mesh.TensorMesh([np.ones(40)])
M.setCellGradBC('dirichlet')
params = Richards.Empirical.HaverkampParams().celia1990
k_fun, theta_fun = Richards.Empirical.haverkamp(M, **params)
k_fun.KsMap = Maps.IdentityMap(nP=M.nC)
bc = np.array([-61.5, -20.7])
h = np.zeros(M.nC) + bc[0]
def getFields(timeStep, method):
timeSteps = np.ones(int(360/timeStep))*timeStep
prob = Richards.RichardsProblem(
M,
hydraulic_conductivity=k_fun,
water_retention=theta_fun,
boundary_conditions=bc, initial_conditions=h,
do_newton=False, method=method
)
prob.timeSteps = timeSteps
return prob.fields(params['Ks'] * np.ones(M.nC))
Hs_M010 = getFields(10, 'mixed')
Hs_M030 = getFields(30, 'mixed')
Hs_M120 = getFields(120, 'mixed')
Hs_H010 = getFields(10, 'head')
Hs_H030 = getFields(30, 'head')
Hs_H120 = getFields(120, 'head')
if not plotIt:
return
plt.figure(figsize=(13, 5))
plt.subplot(121)
plt.plot(40-M.gridCC, Hs_M010[-1], 'b-')
plt.plot(40-M.gridCC, Hs_M030[-1], 'r-')
plt.plot(40-M.gridCC, Hs_M120[-1], 'k-')
plt.ylim([-70, -10])
plt.title('Mixed Method')
plt.xlabel('Depth, cm')
plt.ylabel('Pressure Head, cm')
plt.legend(
('$\Delta t$ = 10 sec', '$\Delta t$ = 30 sec', '$\Delta t$ = 120 sec')
)
plt.subplot(122)
plt.plot(40-M.gridCC, Hs_H010[-1], 'b-')
plt.plot(40-M.gridCC, Hs_H030[-1], 'r-')
plt.plot(40-M.gridCC, Hs_H120[-1], 'k-')
plt.ylim([-70, -10])
plt.title('Head-Based Method')
plt.xlabel('Depth, cm')
plt.ylabel('Pressure Head, cm')
plt.legend((
'$\Delta t$ = 10 sec', '$\Delta t$ = 30 sec', '$\Delta t$ = 120 sec'
))
if __name__ == '__main__':
run()
plt.show()