mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-27 19:48:52 +08:00
SimPEG.FLOW
This commit is contained in:
Rowan Cockett
@@ -9,6 +9,7 @@ env:
|
||||
- TEST_DIR=tests/mesh
|
||||
- TEST_DIR=tests/base
|
||||
- TEST_DIR=tests/examples
|
||||
- TEST_DIR=tests/flow
|
||||
|
||||
# Setup anaconda
|
||||
before_install:
|
||||
|
||||
@@ -38,8 +38,6 @@ def run(plotIt=True):
|
||||
AinvrM = SolverLU(ArM)
|
||||
phirM = AinvrM*rhsrM
|
||||
|
||||
if not plotIt: return
|
||||
|
||||
#Step4: Making Figure
|
||||
fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
|
||||
label = ["(a)", "(b)"]
|
||||
@@ -71,7 +69,7 @@ def run(plotIt=True):
|
||||
axes[i].set_ylabel(" ")
|
||||
axes[i].set_xlabel("x")
|
||||
|
||||
plt.show()
|
||||
if plotIt: plt.show()
|
||||
|
||||
if __name__ == '__main__':
|
||||
run()
|
||||
|
||||
@@ -0,0 +1,52 @@
|
||||
from SimPEG import *
|
||||
from SimPEG.FLOW import Richards
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
def run(plotIt=True):
|
||||
M = Mesh.TensorMesh([np.ones(40)])
|
||||
M.setCellGradBC('dirichlet')
|
||||
params = Richards.Empirical.HaverkampParams().celia1990
|
||||
params['Ks'] = np.log(params['Ks'])
|
||||
E = Richards.Empirical.Haverkamp(M, **params)
|
||||
|
||||
bc = np.array([-61.5,-20.7])
|
||||
h = np.zeros(M.nC) + bc[0]
|
||||
|
||||
|
||||
def getFields(timeStep,method):
|
||||
timeSteps = np.ones(360/timeStep)*timeStep
|
||||
prob = Richards.RichardsProblem(M, mapping=E, timeSteps=timeSteps,
|
||||
boundaryConditions=bc, initialConditions=h,
|
||||
doNewton=False, method=method)
|
||||
return prob.fields(params['Ks'])
|
||||
|
||||
Hs_M10 = getFields(10., 'mixed')
|
||||
Hs_M30 = getFields(30., 'mixed')
|
||||
Hs_M120= getFields(120.,'mixed')
|
||||
Hs_H10 = getFields(10., 'head')
|
||||
Hs_H30 = getFields(30., 'head')
|
||||
Hs_H120= getFields(120.,'head')
|
||||
|
||||
plt.figure(figsize=(13,5))
|
||||
plt.subplot(121)
|
||||
plt.plot(40-M.gridCC, Hs_M10[-1],'b-')
|
||||
plt.plot(40-M.gridCC, Hs_M30[-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_H10[-1],'b-')
|
||||
plt.plot(40-M.gridCC, Hs_H30[-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 plotIt:plt.show()
|
||||
|
||||
if __name__ == '__main__':
|
||||
run()
|
||||
@@ -0,0 +1 @@
|
||||
import Celia1990
|
||||
@@ -0,0 +1,578 @@
|
||||
from SimPEG import Mesh, Maps, Utils, np
|
||||
|
||||
|
||||
class NonLinearMap(object):
|
||||
"""
|
||||
SimPEG NonLinearMap
|
||||
|
||||
"""
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
|
||||
counter = None #: A SimPEG.Utils.Counter object
|
||||
mesh = None #: A SimPEG Mesh
|
||||
|
||||
def __init__(self, mesh):
|
||||
self.mesh = mesh
|
||||
|
||||
def _transform(self, u, m):
|
||||
"""
|
||||
:param numpy.array u: fields
|
||||
:param numpy.array m: model
|
||||
:rtype: numpy.array
|
||||
:return: transformed model
|
||||
|
||||
The *transform* changes the model into the physical property.
|
||||
|
||||
"""
|
||||
return m
|
||||
|
||||
def derivU(self, u, m):
|
||||
"""
|
||||
:param numpy.array u: fields
|
||||
:param numpy.array m: model
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: derivative of transformed model
|
||||
|
||||
The *transform* changes the model into the physical property.
|
||||
The *transformDerivU* provides the derivative of the *transform* with respect to the fields.
|
||||
"""
|
||||
raise NotImplementedError('The transformDerivU is not implemented.')
|
||||
|
||||
|
||||
def derivM(self, u, m):
|
||||
"""
|
||||
:param numpy.array u: fields
|
||||
:param numpy.array m: model
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: derivative of transformed model
|
||||
|
||||
The *transform* changes the model into the physical property.
|
||||
The *transformDerivU* provides the derivative of the *transform* with respect to the model.
|
||||
"""
|
||||
raise NotImplementedError('The transformDerivM is not implemented.')
|
||||
|
||||
@property
|
||||
def nP(self):
|
||||
"""Number of parameters in the model."""
|
||||
return self.mesh.nC
|
||||
|
||||
def example(self):
|
||||
raise NotImplementedError('The example is not implemented.')
|
||||
|
||||
def test(self, m=None):
|
||||
raise NotImplementedError('The test is not implemented.')
|
||||
|
||||
|
||||
class RichardsMap(object):
|
||||
"""docstring for RichardsMap"""
|
||||
|
||||
mesh = None #: SimPEG mesh
|
||||
|
||||
@property
|
||||
def thetaModel(self):
|
||||
"""Model for moisture content"""
|
||||
return self._thetaModel
|
||||
|
||||
@property
|
||||
def kModel(self):
|
||||
"""Model for hydraulic conductivity"""
|
||||
return self._kModel
|
||||
|
||||
def __init__(self, mesh, thetaModel, kModel):
|
||||
self.mesh = mesh
|
||||
assert isinstance(thetaModel, NonLinearMap)
|
||||
assert isinstance(kModel, NonLinearMap)
|
||||
|
||||
self._thetaModel = thetaModel
|
||||
self._kModel = kModel
|
||||
|
||||
def theta(self, u, m):
|
||||
return self.thetaModel.transform(u, m)
|
||||
|
||||
def thetaDerivM(self, u, m):
|
||||
return self.thetaModel.transformDerivM(u, m)
|
||||
|
||||
def thetaDerivU(self, u, m):
|
||||
return self.thetaModel.transformDerivU(u, m)
|
||||
|
||||
def k(self, u, m):
|
||||
return self.kModel.transform(u, m)
|
||||
|
||||
def kDerivM(self, u, m):
|
||||
return self.kModel.transformDerivM(u, m)
|
||||
|
||||
def kDerivU(self, u, m):
|
||||
return self.kModel.transformDerivU(u, m)
|
||||
|
||||
def plot(self, m):
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
m = m[0]
|
||||
h = np.linspace(-100, 20, 1000)
|
||||
ax = plt.subplot(121)
|
||||
ax.plot(self.theta(h, m), h)
|
||||
ax = plt.subplot(122)
|
||||
ax.semilogx(self.k(h, m), h)
|
||||
|
||||
def _assertMatchesPair(self, pair):
|
||||
assert isinstance(self, pair), "Mapping object must be an instance of a %s class."%(pair.__name__)
|
||||
|
||||
|
||||
|
||||
def _ModelProperty(name, models, doc=None, default=None):
|
||||
|
||||
def fget(self):
|
||||
model = models[0]
|
||||
if getattr(self, model, None) is not None:
|
||||
MOD = getattr(self, model)
|
||||
return getattr(MOD, name, default)
|
||||
return default
|
||||
|
||||
def fset(self, value):
|
||||
for model in models:
|
||||
if getattr(self, model, None) is not None:
|
||||
MOD = getattr(self, model)
|
||||
setattr(MOD, name, value)
|
||||
|
||||
return property(fget, fset=fset, doc=doc)
|
||||
|
||||
|
||||
class HaverkampParams(object):
|
||||
"""Holds some default parameterizations for the Haverkamp model."""
|
||||
def __init__(self): pass
|
||||
@property
|
||||
def celia1990(self):
|
||||
"""
|
||||
Parameters used in:
|
||||
|
||||
Celia, Michael A., Efthimios T. Bouloutas, and Rebecca L. Zarba.
|
||||
"A general mass-conservative numerical solution for the unsaturated flow equation."
|
||||
Water Resources Research 26.7 (1990): 1483-1496.
|
||||
|
||||
"""
|
||||
return {'alpha':1.611e+06, 'beta':3.96,
|
||||
'theta_r':0.075, 'theta_s':0.287,
|
||||
'Ks':9.44e-03, 'A':1.175e+06,
|
||||
'gamma':4.74}
|
||||
|
||||
|
||||
class _haverkamp_theta(NonLinearMap):
|
||||
|
||||
theta_s = 0.430
|
||||
theta_r = 0.078
|
||||
alpha = 0.036
|
||||
beta = 3.960
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
NonLinearMap.__init__(self, mesh)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
def setModel(self, m):
|
||||
self._currentModel = m
|
||||
|
||||
def transform(self, u, m):
|
||||
self.setModel(m)
|
||||
f = (self.alpha*(self.theta_s - self.theta_r )/
|
||||
(self.alpha + abs(u)**self.beta) + self.theta_r)
|
||||
if Utils.isScalar(self.theta_s):
|
||||
f[u >= 0] = self.theta_s
|
||||
else:
|
||||
f[u >= 0] = self.theta_s[u >= 0]
|
||||
return f
|
||||
|
||||
def transformDerivM(self, u, m):
|
||||
self.setModel(m)
|
||||
|
||||
def transformDerivU(self, u, m):
|
||||
self.setModel(m)
|
||||
g = (self.alpha*((self.theta_s - self.theta_r)/
|
||||
(self.alpha + abs(u)**self.beta)**2)
|
||||
*(-self.beta*abs(u)**(self.beta-1)*np.sign(u)))
|
||||
g[u >= 0] = 0
|
||||
g = Utils.sdiag(g)
|
||||
return g
|
||||
|
||||
|
||||
class _haverkamp_k(NonLinearMap):
|
||||
|
||||
A = 1.175e+06
|
||||
gamma = 4.74
|
||||
Ks = np.log(24.96)
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
NonLinearMap.__init__(self, mesh)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
def setModel(self, m):
|
||||
self._currentModel = m
|
||||
#TODO: Fix me!
|
||||
self.Ks = m
|
||||
|
||||
def transform(self, u, m):
|
||||
self.setModel(m)
|
||||
f = np.exp(self.Ks)*self.A/(self.A+abs(u)**self.gamma)
|
||||
if Utils.isScalar(self.Ks):
|
||||
f[u >= 0] = np.exp(self.Ks)
|
||||
else:
|
||||
f[u >= 0] = np.exp(self.Ks[u >= 0])
|
||||
return f
|
||||
|
||||
def transformDerivM(self, u, m):
|
||||
self.setModel(m)
|
||||
#A
|
||||
# dA = np.exp(self.Ks)/(self.A+abs(u)**self.gamma) - np.exp(self.Ks)*self.A/(self.A+abs(u)**self.gamma)**2
|
||||
#gamma
|
||||
# dgamma = -(self.A*np.exp(self.Ks)*np.log(abs(u))*abs(u)**self.gamma)/(self.A + abs(u)**self.gamma)**2
|
||||
|
||||
# This assumes that the the model is Ks
|
||||
return Utils.sdiag(self.transform(u, m))
|
||||
|
||||
def transformDerivU(self, u, m):
|
||||
self.setModel(m)
|
||||
g = -(np.exp(self.Ks)*self.A*self.gamma*abs(u)**(self.gamma-1)*np.sign(u))/((self.A+abs(u)**self.gamma)**2)
|
||||
g[u >= 0] = 0
|
||||
g = Utils.sdiag(g)
|
||||
return g
|
||||
|
||||
class Haverkamp(RichardsMap):
|
||||
"""Haverkamp Model"""
|
||||
|
||||
alpha = _ModelProperty('alpha', ['thetaModel'], default=1.6110e+06)
|
||||
beta = _ModelProperty('beta', ['thetaModel'], default=3.96)
|
||||
theta_r = _ModelProperty('theta_r', ['thetaModel'], default=0.075)
|
||||
theta_s = _ModelProperty('theta_s', ['thetaModel'], default=0.287)
|
||||
|
||||
Ks = _ModelProperty('Ks', ['kModel'], default=np.log(24.96))
|
||||
A = _ModelProperty('A', ['kModel'], default=1.1750e+06)
|
||||
gamma = _ModelProperty('gamma', ['kModel'], default=4.74)
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
RichardsMap.__init__(self, mesh,
|
||||
_haverkamp_theta(mesh),
|
||||
_haverkamp_k(mesh))
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
|
||||
|
||||
|
||||
class _vangenuchten_theta(NonLinearMap):
|
||||
|
||||
theta_s = 0.430
|
||||
theta_r = 0.078
|
||||
alpha = 0.036
|
||||
n = 1.560
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
NonLinearMap.__init__(self, mesh)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
def setModel(self, m):
|
||||
self._currentModel = m
|
||||
|
||||
def transform(self, u, m):
|
||||
self.setModel(m)
|
||||
m = 1 - 1.0/self.n
|
||||
f = (( self.theta_s - self.theta_r )/
|
||||
((1+abs(self.alpha*u)**self.n)**m) + self.theta_r)
|
||||
if Utils.isScalar(self.theta_s):
|
||||
f[u >= 0] = self.theta_s
|
||||
else:
|
||||
f[u >= 0] = self.theta_s[u >= 0]
|
||||
|
||||
return f
|
||||
|
||||
def transformDerivM(self, u, m):
|
||||
self.setModel(m)
|
||||
|
||||
def transformDerivU(self, u, m):
|
||||
g = -self.alpha*self.n*abs(self.alpha*u)**(self.n - 1)*np.sign(self.alpha*u)*(1./self.n - 1)*(self.theta_r - self.theta_s)*(abs(self.alpha*u)**self.n + 1)**(1./self.n - 2)
|
||||
g[u >= 0] = 0
|
||||
g = Utils.sdiag(g)
|
||||
return g
|
||||
|
||||
|
||||
class _vangenuchten_k(NonLinearMap):
|
||||
|
||||
I = 0.500
|
||||
alpha = 0.036
|
||||
n = 1.560
|
||||
Ks = np.log(24.96)
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
NonLinearMap.__init__(self, mesh)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
def setModel(self, m):
|
||||
self._currentModel = m
|
||||
#TODO: Fix me!
|
||||
self.Ks = m
|
||||
|
||||
def transform(self, u, m):
|
||||
self.setModel(m)
|
||||
|
||||
alpha = self.alpha
|
||||
I = self.I
|
||||
n = self.n
|
||||
Ks = self.Ks
|
||||
m = 1.0 - 1.0/n
|
||||
|
||||
theta_e = 1.0/((1.0+abs(alpha*u)**n)**m)
|
||||
f = np.exp(Ks)*theta_e**I* ( ( 1.0 - ( 1.0 - theta_e**(1.0/m) )**m )**2 )
|
||||
if Utils.isScalar(self.Ks):
|
||||
f[u >= 0] = np.exp(self.Ks)
|
||||
else:
|
||||
f[u >= 0] = np.exp(self.Ks[u >= 0])
|
||||
return f
|
||||
|
||||
def transformDerivM(self, u, m):
|
||||
self.setModel(m)
|
||||
# #alpha
|
||||
# # dA = I*u*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*u)**n + 1)**(1.0/n - 2) - (2*u*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/(((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n));
|
||||
# #n
|
||||
# # dn = 2*np.exp(Ks)*((np.log(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n))/n**2 + ((1.0/n - 1)*(((np.log(abs(alpha*u)**n + 1)*(abs(alpha*u)**n + 1)**(1.0/n - 1))/n**2 - abs(alpha*u)**n*np.log(abs(alpha*u))*(1.0/n - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/((1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)) - np.log((abs(alpha*u)**n + 1)**(1.0/n - 1))/(n**2*(1.0/n - 1)**2*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))))/(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n))*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1) - I*np.exp(Ks)*((np.log(abs(alpha*u)**n + 1)*(abs(alpha*u)**n + 1)**(1.0/n - 1))/n**2 - abs(alpha*u)**n*np.log(abs(alpha*u))*(1.0/n - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2;
|
||||
# #I
|
||||
# # dI = np.exp(Ks)*np.log((abs(alpha*u)**n + 1)**(1.0/n - 1))*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2;
|
||||
return Utils.sdiag(self.transform(u, m)) # This assumes that the the model is Ks
|
||||
|
||||
def transformDerivU(self, u, m):
|
||||
self.setModel(m)
|
||||
alpha = self.alpha
|
||||
I = self.I
|
||||
n = self.n
|
||||
Ks = self.Ks
|
||||
m = 1.0 - 1.0/n
|
||||
|
||||
g = I*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1.0)*np.sign(alpha*u)*(1.0/n - 1.0)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*u)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/(((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n))
|
||||
g[u >= 0] = 0
|
||||
g = Utils.sdiag(g)
|
||||
return g
|
||||
|
||||
class VanGenuchten(RichardsMap):
|
||||
"""vanGenuchten Model"""
|
||||
|
||||
theta_r = _ModelProperty('theta_r', ['thetaModel'], default=0.075)
|
||||
theta_s = _ModelProperty('theta_s', ['thetaModel'], default=0.287)
|
||||
|
||||
alpha = _ModelProperty('alpha', ['thetaModel', 'kModel'], default=0.036)
|
||||
n = _ModelProperty('n', ['thetaModel', 'kModel'], default=1.560)
|
||||
|
||||
Ks = _ModelProperty('Ks', ['kModel'], default=np.log(24.96))
|
||||
I = _ModelProperty('I', ['kModel'], default=0.500)
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
RichardsMap.__init__(self, mesh,
|
||||
_vangenuchten_theta(mesh),
|
||||
_vangenuchten_k(mesh))
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
|
||||
class VanGenuchtenParams(object):
|
||||
"""
|
||||
The RETC code for quantifying the hydraulic functions of unsaturated soils,
|
||||
Van Genuchten, M Th, Leij, F J, Yates, S R
|
||||
|
||||
Table 3: Average values for selected soil water retention and hydraulic
|
||||
conductivity parameters for 11 major soil textural groups
|
||||
according to Rawls et al. [1982]
|
||||
|
||||
"""
|
||||
def __init__(self): pass
|
||||
@property
|
||||
def sand(self):
|
||||
return {"theta_r": 0.020, "theta_s": 0.417, "alpha": 0.138*100., "n": 1.592, "Ks": 504.0/100./24./60./60.}
|
||||
@property
|
||||
def loamySand(self):
|
||||
return {"theta_r": 0.035, "theta_s": 0.401, "alpha": 0.115*100., "n": 1.474, "Ks": 146.6/100./24./60./60.}
|
||||
@property
|
||||
def sandyLoam(self):
|
||||
return {"theta_r": 0.041, "theta_s": 0.412, "alpha": 0.068*100., "n": 1.322, "Ks": 62.16/100./24./60./60.}
|
||||
@property
|
||||
def loam(self):
|
||||
return {"theta_r": 0.027, "theta_s": 0.434, "alpha": 0.090*100., "n": 1.220, "Ks": 16.32/100./24./60./60.}
|
||||
@property
|
||||
def siltLoam(self):
|
||||
return {"theta_r": 0.015, "theta_s": 0.486, "alpha": 0.048*100., "n": 1.211, "Ks": 31.68/100./24./60./60.}
|
||||
@property
|
||||
def sandyClayLoam(self):
|
||||
return {"theta_r": 0.068, "theta_s": 0.330, "alpha": 0.036*100., "n": 1.250, "Ks": 10.32/100./24./60./60.}
|
||||
@property
|
||||
def clayLoam(self):
|
||||
return {"theta_r": 0.075, "theta_s": 0.390, "alpha": 0.039*100., "n": 1.194, "Ks": 5.52/100./24./60./60.}
|
||||
@property
|
||||
def siltyClayLoam(self):
|
||||
return {"theta_r": 0.040, "theta_s": 0.432, "alpha": 0.031*100., "n": 1.151, "Ks": 3.60/100./24./60./60.}
|
||||
@property
|
||||
def sandyClay(self):
|
||||
return {"theta_r": 0.109, "theta_s": 0.321, "alpha": 0.034*100., "n": 1.168, "Ks": 2.88/100./24./60./60.}
|
||||
@property
|
||||
def siltyClay(self):
|
||||
return {"theta_r": 0.056, "theta_s": 0.423, "alpha": 0.029*100., "n": 1.127, "Ks": 2.16/100./24./60./60.}
|
||||
@property
|
||||
def clay(self):
|
||||
return {"theta_r": 0.090, "theta_s": 0.385, "alpha": 0.027*100., "n": 1.131, "Ks": 1.44/100./24./60./60.}
|
||||
|
||||
|
||||
# From: INDIRECT METHODS FOR ESTIMATING THE HYDRAULIC PROPERTIES OF UNSATURATED SOILS
|
||||
# @property
|
||||
# def siltLoamGE3(self):
|
||||
# """Soil Index: 3310"""
|
||||
# return {"theta_r": 0.139, "theta_s": 0.394, "alpha": 0.00414, "n": 2.15}
|
||||
# @property
|
||||
# def yoloLightClayK_WC(self):
|
||||
# """Soil Index: None"""
|
||||
# return {"theta_r": 0.205, "theta_s": 0.499, "alpha": 0.02793, "n": 1.71}
|
||||
# @property
|
||||
# def yoloLightClayK_H(self):
|
||||
# """Soil Index: None"""
|
||||
# return {"theta_r": 0.205, "theta_s": 0.499, "alpha": 0.02793, "n": 1.71}
|
||||
# @property
|
||||
# def hygieneSandstone(self):
|
||||
# """Soil Index: 4130"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.256, "alpha": 0.00562, "n": 3.27}
|
||||
# @property
|
||||
# def lambcrgClay(self):
|
||||
# """Soil Index: 1003"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.502, "alpha": 0.140, "n": 1.93}
|
||||
# @property
|
||||
# def beitNetofaClaySoil(self):
|
||||
# """Soil Index: 1006"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.447, "alpha": 0.00156, "n": 1.17}
|
||||
# @property
|
||||
# def shiohotSiltyClay(self):
|
||||
# """Soil Index: 1101"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.456, "alpha": 183, "n":1.17}
|
||||
# @property
|
||||
# def siltColumbia(self):
|
||||
# """Soil Index: 2001"""
|
||||
# return {"theta_r": 0.146, "theta_s": 0.397, "alpha": 0.0145, "n": 1.85}
|
||||
# @property
|
||||
# def siltMontCenis(self):
|
||||
# """Soil Index: 2002"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.425, "alpha": 0.0103, "n": 1.34}
|
||||
# @property
|
||||
# def slateDust(self):
|
||||
# """Soil Index: 2004"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.498, "alpha": 0.00981, "n": 6.75}
|
||||
# @property
|
||||
# def weldSiltyClayLoam(self):
|
||||
# """Soil Index: 3001"""
|
||||
# return {"theta_r": 0.159, "theta_s": 0.496, "alpha": 0.0136, "n": 5.45}
|
||||
# @property
|
||||
# def rideauClayLoam_Wetting(self):
|
||||
# """Soil Index: 3101a"""
|
||||
# return {"theta_r": 0.279, "theta_s": 0.419, "alpha": 0.0661, "n": 1.89}
|
||||
# @property
|
||||
# def rideauClayLoam_Drying(self):
|
||||
# """Soil Index: 3101b"""
|
||||
# return {"theta_r": 0.290, "theta_s": 0.419, "alpha": 0.0177, "n": 3.18}
|
||||
# @property
|
||||
# def caribouSiltLoam_Drying(self):
|
||||
# """Soil Index: 3301a"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.451, "alpha": 0.00845, "n": 1.29}
|
||||
# @property
|
||||
# def caribouSiltLoam_Wetting(self):
|
||||
# """Soil Index: 3301b"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.450, "alpha": 0.140, "n": 1.09}
|
||||
# @property
|
||||
# def grenvilleSiltLoam_Wetting(self):
|
||||
# """Soil Index: 3302a"""
|
||||
# return {"theta_r": 0.013, "theta_s": 0523, "alpha": 0.0630, "n": 1.24}
|
||||
# @property
|
||||
# def grenvilleSiltLoam_Drying(self):
|
||||
# """Soil Index: 3302c"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.488, "alpha": 0.0112, "n": 1.23}
|
||||
# @property
|
||||
# def touchetSiltLoam(self):
|
||||
# """Soil Index: 3304"""
|
||||
# return {"theta_r": 0.183, "theta_s": 0.498, "alpha": 0.0104, "n": 5.78}
|
||||
# @property
|
||||
# def gilatLoam(self):
|
||||
# """Soil Index: 3402a"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.454, "alpha": 0.0291, "n": 1.47}
|
||||
# @property
|
||||
# def pachapaLoam(self):
|
||||
# """Soil Index: 3403"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.472, "alpha": 0.00829, "n": 1.62}
|
||||
# @property
|
||||
# def adelantoLoam(self):
|
||||
# """Soil Index: 3404"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.444, "alpha": 0.00710, "n": 1.26}
|
||||
# @property
|
||||
# def indioLoam(self):
|
||||
# """Soil Index: 3405a"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.507, "alpha": 0.00847, "n": 1.60}
|
||||
# @property
|
||||
# def guclphLoam(self):
|
||||
# """Soil Index: 3407a"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.563, "alpha": 0.0275, "n": 1.27}
|
||||
# @property
|
||||
# def guclphLoam(self):
|
||||
# """Soil Index: 3407b"""
|
||||
# return {"theta_r": 0.236, "theta_s": 0.435, "alpha": 0.0271, "n": 262}
|
||||
# @property
|
||||
# def rubiconSandyLoam(self):
|
||||
# """Soil Index: 3501a"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.393, "alpha": 0.00972, "n": 2.18}
|
||||
# @property
|
||||
# def rubiconSandyLoam(self):
|
||||
# """Soil Index: 350lb"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.433, "alpha": 0.147, "n": 1.28}
|
||||
# @property
|
||||
# def pachapaFmeSandyClay(self):
|
||||
# """Soil Index: 3503a"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.340, "alpha": 0.0194, "n": 1.45}
|
||||
# @property
|
||||
# def gilatSandyLoam(self):
|
||||
# """Soil Index: 3504"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.432, "alpha": 0.0103, "n": 1.48}
|
||||
# @property
|
||||
# def plainfieldSand_210to250(self):
|
||||
# """Soil Index: 4101a"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.351, "alpha": 0.0236, "n": 12.30}
|
||||
# @property
|
||||
# def plainfieldSand_210to250(self):
|
||||
# """Soil Index: 4101b"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.312, "alpha": 0.0387, "n": 4.48}
|
||||
# @property
|
||||
# def plainfieldSand_177to210(self):
|
||||
# """Soil Index: 4102a"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.361, "alpha": 0.0207, "n": 10.0}
|
||||
# @property
|
||||
# def plainfieldSand_177to210(self):
|
||||
# """Soil Index: 4102b"""
|
||||
# return {"theta_r": 0.022, "theta_s": 0.309, "alpha": 0.0328, "n": 6.23}
|
||||
# @property
|
||||
# def plainfieldSand_149to177(self):
|
||||
# """Soil Index: 4103a"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.387, "alpha": 0.0173, "n": 7.80}
|
||||
# @property
|
||||
# def plainfieldSand_149to177(self):
|
||||
# """Soil Index: 4103b"""
|
||||
# return {"theta_r": 0.025, "theta_s": 0.321, "alpha": 0.0272, "n": 6.69}
|
||||
# @property
|
||||
# def plainfieldSand_l25to149(self):
|
||||
# """Soil Index: 4104a"""
|
||||
# return {"theta_r": 0.000, "theta_s": 03770, "alpha": 0.0145, "n": 10.60}
|
||||
# @property
|
||||
# def plainfieldSand_125to149(self):
|
||||
# """Soil Index: 4104b"""
|
||||
# return {"theta_r": 0.000, "theta_s": 0.342, "alpha": 0.0230, "n": 5.18}
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
import matplotlib.pyplot as plt
|
||||
M = Mesh.TensorMesh([10])
|
||||
VGparams = VanGenuchtenParams()
|
||||
leg = []
|
||||
for p in dir(VGparams):
|
||||
if p[0] == '_': continue
|
||||
leg += [p]
|
||||
params = getattr(VGparams, p)
|
||||
model = VanGenuchten(M, **params)
|
||||
ks = np.log(np.r_[params['Ks']])
|
||||
model.plot(ks)
|
||||
|
||||
plt.legend(leg)
|
||||
|
||||
plt.show()
|
||||
@@ -0,0 +1,304 @@
|
||||
from SimPEG import *
|
||||
from Empirical import RichardsMap
|
||||
import time
|
||||
|
||||
|
||||
class RichardsRx(Survey.BaseTimeRx):
|
||||
"""Richards Receiver Object"""
|
||||
|
||||
knownRxTypes = ['saturation','pressureHead']
|
||||
|
||||
def projectFields(self, U, m, mapping, mesh, timeMesh):
|
||||
|
||||
if self.rxType == 'pressureHead':
|
||||
u = np.concatenate(U)
|
||||
elif self.rxType == 'saturation':
|
||||
u = np.concatenate([mapping.theta(ui, m) for ui in U])
|
||||
|
||||
return self.getP(mesh, timeMesh) * u
|
||||
|
||||
def projectFieldsDeriv(self, U, m, mapping, mesh, timeMesh):
|
||||
|
||||
P = self.getP(mesh, timeMesh)
|
||||
if self.rxType == 'pressureHead':
|
||||
return P
|
||||
elif self.rxType == 'saturation':
|
||||
#TODO: if m is a parameter in the theta
|
||||
# distribution, we may need to do
|
||||
# some more chain rule here.
|
||||
dT = sp.block_diag([mapping.thetaDerivU(ui, m) for ui in U])
|
||||
return P*dT
|
||||
|
||||
|
||||
class RichardsSurvey(Survey.BaseSurvey):
|
||||
"""docstring for RichardsSurvey"""
|
||||
|
||||
rxList = None
|
||||
|
||||
def __init__(self, rxList, **kwargs):
|
||||
self.rxList = rxList
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
|
||||
@property
|
||||
def nD(self):
|
||||
return np.array([rx.nD for rx in self.rxList]).sum()
|
||||
|
||||
@Utils.count
|
||||
@Utils.requires('prob')
|
||||
def dpred(self, m, u=None):
|
||||
"""
|
||||
Create the projected data from a model.
|
||||
The field, u, (if provided) will be used for the predicted data
|
||||
instead of recalculating the fields (which may be expensive!).
|
||||
|
||||
.. math::
|
||||
d_\\text{pred} = P(u(m), m)
|
||||
|
||||
Where P is a projection of the fields onto the data space.
|
||||
"""
|
||||
if u is None: u = self.prob.fields(m)
|
||||
return Utils.mkvc(self.projectFields(u, m))
|
||||
|
||||
@Utils.requires('prob')
|
||||
def projectFields(self, U, m):
|
||||
Ds = range(len(self.rxList))
|
||||
for ii, rx in enumerate(self.rxList):
|
||||
Ds[ii] = rx.projectFields(U, m,
|
||||
self.prob.mapping,
|
||||
self.prob.mesh,
|
||||
self.prob.timeMesh)
|
||||
|
||||
return np.concatenate(Ds)
|
||||
|
||||
@Utils.requires('prob')
|
||||
def projectFieldsDeriv(self, U, m):
|
||||
"""The Derivative with respect to the fields."""
|
||||
Ds = range(len(self.rxList))
|
||||
for ii, rx in enumerate(self.rxList):
|
||||
Ds[ii] = rx.projectFieldsDeriv(U, m,
|
||||
self.prob.mapping,
|
||||
self.prob.mesh,
|
||||
self.prob.timeMesh)
|
||||
|
||||
return sp.vstack(Ds)
|
||||
|
||||
class RichardsProblem(Problem.BaseTimeProblem):
|
||||
"""docstring for RichardsProblem"""
|
||||
|
||||
boundaryConditions = None
|
||||
initialConditions = None
|
||||
|
||||
surveyPair = RichardsSurvey
|
||||
mapPair = RichardsMap
|
||||
|
||||
debug=True
|
||||
|
||||
Solver = Solver
|
||||
solverOpts = {}
|
||||
|
||||
def __init__(self, mesh, mapping=None, **kwargs):
|
||||
Problem.BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
|
||||
|
||||
def getBoundaryConditions(self, ii, u_ii):
|
||||
if type(self.boundaryConditions) is np.ndarray:
|
||||
return self.boundaryConditions
|
||||
|
||||
time = self.timeMesh.vectorCCx[ii]
|
||||
|
||||
return self.boundaryConditions(time, u_ii)
|
||||
|
||||
@property
|
||||
def method(self):
|
||||
"""Method must be either 'mixed' or 'head'. See notes in Celia et al., 1990."""
|
||||
return getattr(self, '_method', 'mixed')
|
||||
@method.setter
|
||||
def method(self, value):
|
||||
assert value in ['mixed','head'], "method must be 'mixed' or 'head'."
|
||||
self._method = value
|
||||
|
||||
# Setting doNewton will clear the rootFinder, which will be reinitialized when called
|
||||
doNewton = Utils.dependentProperty('_doNewton', False, ['_rootFinder'],
|
||||
"Do a Newton iteration. If False, a Picard iteration will be completed.")
|
||||
|
||||
maxIterRootFinder = Utils.dependentProperty('_maxIterRootFinder', 30, ['_rootFinder'],
|
||||
"Maximum iterations for rootFinder iteration.")
|
||||
tolRootFinder = Utils.dependentProperty('_tolRootFinder', 1e-4, ['_rootFinder'],
|
||||
"Maximum iterations for rootFinder iteration.")
|
||||
|
||||
@property
|
||||
def rootFinder(self):
|
||||
"""Root-finding Algorithm"""
|
||||
if getattr(self, '_rootFinder', None) is None:
|
||||
self._rootFinder = Optimization.NewtonRoot(doLS=self.doNewton, maxIter=self.maxIterRootFinder, tol=self.tolRootFinder, Solver=self.Solver)
|
||||
return self._rootFinder
|
||||
|
||||
@Utils.timeIt
|
||||
def fields(self, m):
|
||||
tic = time.time()
|
||||
u = range(self.nT+1)
|
||||
u[0] = self.initialConditions
|
||||
for ii, dt in enumerate(self.timeSteps):
|
||||
bc = self.getBoundaryConditions(ii, u[ii])
|
||||
u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, dt, bc, return_g=return_g), u[ii])
|
||||
if self.debug: print "Solving Fields (%4d/%d - %3.1f%% Done) %d Iterations, %4.2f seconds"%(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic)
|
||||
return u
|
||||
|
||||
@Utils.timeIt
|
||||
def diagsJacobian(self, m, hn, hn1, dt, bc):
|
||||
|
||||
DIV = self.mesh.faceDiv
|
||||
GRAD = self.mesh.cellGrad
|
||||
BC = self.mesh.cellGradBC
|
||||
AV = self.mesh.aveF2CC.T
|
||||
if self.mesh.dim == 1:
|
||||
Dz = self.mesh.faceDivx
|
||||
elif self.mesh.dim == 2:
|
||||
Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.vnF[0]), self.mesh.faceDivy),format='csr')
|
||||
elif self.mesh.dim == 3:
|
||||
Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.vnF[0]+self.mesh.vnF[1]), self.mesh.faceDivz),format='csr')
|
||||
|
||||
dT = self.mapping.thetaDerivU(hn, m)
|
||||
dT1 = self.mapping.thetaDerivU(hn1, m)
|
||||
K1 = self.mapping.k(hn1, m)
|
||||
dK1 = self.mapping.kDerivU(hn1, m)
|
||||
dKm1 = self.mapping.kDerivM(hn1, m)
|
||||
|
||||
# Compute part of the derivative of:
|
||||
#
|
||||
# DIV*diag(GRAD*hn1+BC*bc)*(AV*(1.0/K))^-1
|
||||
|
||||
DdiagGh1 = DIV*Utils.sdiag(GRAD*hn1+BC*bc)
|
||||
diagAVk2_AVdiagK2 = Utils.sdiag((AV*(1./K1))**(-2)) * AV*Utils.sdiag(K1**(-2))
|
||||
|
||||
# The matrix that we are computing has the form:
|
||||
#
|
||||
# - - - - - -
|
||||
# | Adiag | | h1 | | b1 |
|
||||
# | Asub Adiag | | h2 | | b2 |
|
||||
# | Asub Adiag | | h3 | = | b3 |
|
||||
# | ... ... | | .. | | .. |
|
||||
# | Asub Adiag | | hn | | bn |
|
||||
# - - - - - -
|
||||
|
||||
Asub = (-1.0/dt)*dT
|
||||
|
||||
Adiag = (
|
||||
(1.0/dt)*dT1
|
||||
-DdiagGh1*diagAVk2_AVdiagK2*dK1
|
||||
-DIV*Utils.sdiag(1./(AV*(1./K1)))*GRAD
|
||||
-Dz*diagAVk2_AVdiagK2*dK1
|
||||
)
|
||||
|
||||
B = DdiagGh1*diagAVk2_AVdiagK2*dKm1 + Dz*diagAVk2_AVdiagK2*dKm1
|
||||
|
||||
return Asub, Adiag, B
|
||||
|
||||
@Utils.timeIt
|
||||
def getResidual(self, m, hn, h, dt, bc, return_g=True):
|
||||
"""
|
||||
Where h is the proposed value for the next time iterate (h_{n+1})
|
||||
"""
|
||||
DIV = self.mesh.faceDiv
|
||||
GRAD = self.mesh.cellGrad
|
||||
BC = self.mesh.cellGradBC
|
||||
AV = self.mesh.aveF2CC.T
|
||||
if self.mesh.dim == 1:
|
||||
Dz = self.mesh.faceDivx
|
||||
elif self.mesh.dim == 2:
|
||||
Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.vnF[0]), self.mesh.faceDivy),format='csr')
|
||||
elif self.mesh.dim == 3:
|
||||
Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.vnF[0]+self.mesh.vnF[1]), self.mesh.faceDivz),format='csr')
|
||||
|
||||
T = self.mapping.theta(h, m)
|
||||
dT = self.mapping.thetaDerivU(h, m)
|
||||
Tn = self.mapping.theta(hn, m)
|
||||
K = self.mapping.k(h, m)
|
||||
dK = self.mapping.kDerivU(h, m)
|
||||
|
||||
aveK = 1./(AV*(1./K))
|
||||
|
||||
RHS = DIV*Utils.sdiag(aveK)*(GRAD*h+BC*bc) + Dz*aveK
|
||||
if self.method == 'mixed':
|
||||
r = (T-Tn)/dt - RHS
|
||||
elif self.method == 'head':
|
||||
r = dT*(h - hn)/dt - RHS
|
||||
|
||||
if not return_g: return r
|
||||
|
||||
J = dT/dt - DIV*Utils.sdiag(aveK)*GRAD
|
||||
if self.doNewton:
|
||||
DDharmAve = Utils.sdiag(aveK**2)*AV*Utils.sdiag(K**(-2)) * dK
|
||||
J = J - DIV*Utils.sdiag(GRAD*h + BC*bc)*DDharmAve - Dz*DDharmAve
|
||||
|
||||
return r, J
|
||||
|
||||
@Utils.timeIt
|
||||
def Jfull(self, m, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
nn = len(u)-1
|
||||
Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
|
||||
for ii in range(nn):
|
||||
dt = self.timeSteps[ii]
|
||||
bc = self.getBoundaryConditions(ii, u[ii])
|
||||
Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, u[ii], u[ii+1], dt, bc)
|
||||
Ad = sp.block_diag(Adiags)
|
||||
zRight = Utils.spzeros((len(Asubs)-1)*Asubs[0].shape[0],Adiags[0].shape[1])
|
||||
zTop = Utils.spzeros(Adiags[0].shape[0], len(Adiags)*Adiags[0].shape[1])
|
||||
As = sp.vstack((zTop,sp.hstack((sp.block_diag(Asubs[1:]),zRight))))
|
||||
A = As + Ad
|
||||
B = np.array(sp.vstack(Bs).todense())
|
||||
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
P = self.survey.projectFieldsDeriv(u, m)
|
||||
AinvB = Ainv * B
|
||||
z = np.zeros((self.mesh.nC, B.shape[1]))
|
||||
zAinvB = np.vstack((z, AinvB))
|
||||
J = P * zAinvB
|
||||
return J
|
||||
|
||||
@Utils.timeIt
|
||||
def Jvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
JvC = range(len(u)-1) # Cell to hold each row of the long vector.
|
||||
|
||||
# This is done via forward substitution.
|
||||
bc = self.getBoundaryConditions(0, u[0])
|
||||
temp, Adiag, B = self.diagsJacobian(m, u[0], u[1], self.timeSteps[0], bc)
|
||||
Adiaginv = self.Solver(Adiag, **self.solverOpts)
|
||||
JvC[0] = Adiaginv * (B*v)
|
||||
|
||||
for ii in range(1,len(u)-1):
|
||||
bc = self.getBoundaryConditions(ii, u[ii])
|
||||
Asub, Adiag, B = self.diagsJacobian(m, u[ii], u[ii+1], self.timeSteps[ii], bc)
|
||||
Adiaginv = self.Solver(Adiag, **self.solverOpts)
|
||||
JvC[ii] = Adiaginv * (B*v - Asub*JvC[ii-1])
|
||||
|
||||
P = self.survey.projectFieldsDeriv(u, m)
|
||||
return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC)
|
||||
|
||||
@Utils.timeIt
|
||||
def Jtvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.field(m)
|
||||
|
||||
P = self.survey.projectFieldsDeriv(u, m)
|
||||
PTv = P.T*v
|
||||
|
||||
# This is done via backward substitution.
|
||||
minus = 0
|
||||
BJtv = 0
|
||||
for ii in range(len(u)-1,0,-1):
|
||||
bc = self.getBoundaryConditions(ii-1, u[ii-1])
|
||||
Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii], self.timeSteps[ii-1], bc)
|
||||
#select the correct part of v
|
||||
vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0])
|
||||
AdiaginvT = self.Solver(Adiag.T, **self.solverOpts)
|
||||
JTvC = AdiaginvT * (PTv[vpart] - minus)
|
||||
minus = Asub.T*JTvC # this is now the super diagonal.
|
||||
BJtv = BJtv + B.T*JTvC
|
||||
|
||||
return BJtv
|
||||
@@ -0,0 +1,2 @@
|
||||
import Empirical
|
||||
from RichardsProblem import *
|
||||
@@ -0,0 +1 @@
|
||||
import Richards
|
||||
@@ -0,0 +1,47 @@
|
||||
.. _api_Richards:
|
||||
|
||||
|
||||
Richards Equation
|
||||
*****************
|
||||
|
||||
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 [Celia et al., 1990]
|
||||
|
||||
.. 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 [Celia et al., 1990] 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 Celia et al. (1990):
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG.FLOW.Examples import Celia1990
|
||||
Celia1990.run()
|
||||
|
||||
Richards
|
||||
========
|
||||
|
||||
.. automodule:: simpegFLOW.Richards.Empirical
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -79,6 +79,14 @@ Utility Codes
|
||||
api_Tests
|
||||
|
||||
|
||||
Packages
|
||||
********
|
||||
|
||||
.. toctree::
|
||||
:maxdepth: 3
|
||||
|
||||
flow/index
|
||||
|
||||
Developer's Documentation
|
||||
*************************
|
||||
|
||||
|
||||
@@ -0,0 +1,11 @@
|
||||
if __name__ == '__main__':
|
||||
import os
|
||||
import glob
|
||||
import unittest
|
||||
test_file_strings = glob.glob('test_*.py')
|
||||
module_strings = [str[0:len(str)-3] for str in test_file_strings]
|
||||
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
|
||||
in module_strings]
|
||||
testSuite = unittest.TestSuite(suites)
|
||||
|
||||
unittest.TextTestRunner(verbosity=2).run(testSuite)
|
||||
@@ -0,0 +1,286 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG.Tests import OrderTest, checkDerivative
|
||||
from scipy.sparse.linalg import dsolve
|
||||
from SimPEG.FLOW import Richards
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
Solver = MumpsSolver
|
||||
except Exception, e:
|
||||
pass
|
||||
|
||||
|
||||
TOL = 1E-8
|
||||
|
||||
class TestModels(unittest.TestCase):
|
||||
|
||||
def test_BaseHaverkamp_Theta(self):
|
||||
mesh = Mesh.TensorMesh([50])
|
||||
hav = Richards.Empirical._haverkamp_theta(mesh)
|
||||
m = np.random.randn(50)
|
||||
def wrapper(u):
|
||||
return hav.transform(u, m), hav.transformDerivU(u, m)
|
||||
passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_vangenuchten_theta(self):
|
||||
mesh = Mesh.TensorMesh([50])
|
||||
hav = Richards.Empirical._vangenuchten_theta(mesh)
|
||||
m = np.random.randn(50)
|
||||
def wrapper(u):
|
||||
return hav.transform(u, m), hav.transformDerivU(u, m)
|
||||
passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_BaseHaverkamp_k(self):
|
||||
mesh = Mesh.TensorMesh([50])
|
||||
hav = Richards.Empirical._haverkamp_k(mesh)
|
||||
m = np.random.randn(50)
|
||||
def wrapper(u):
|
||||
return hav.transform(u, m), hav.transformDerivU(u, m)
|
||||
passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
hav = Richards.Empirical._haverkamp_k(mesh)
|
||||
u = np.random.randn(50)
|
||||
def wrapper(m):
|
||||
return hav.transform(u, m), hav.transformDerivM(u, m)
|
||||
passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_vangenuchten_k(self):
|
||||
mesh = Mesh.TensorMesh([50])
|
||||
hav = Richards.Empirical._vangenuchten_k(mesh)
|
||||
m = np.random.randn(50)
|
||||
def wrapper(u):
|
||||
return hav.transform(u, m), hav.transformDerivU(u, m)
|
||||
passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
hav = Richards.Empirical._vangenuchten_k(mesh)
|
||||
u = np.random.randn(50)
|
||||
def wrapper(m):
|
||||
return hav.transform(u, m), hav.transformDerivM(u, m)
|
||||
passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
|
||||
|
||||
class RichardsTests1D(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
M = Mesh.TensorMesh([np.ones(20)])
|
||||
M.setCellGradBC('dirichlet')
|
||||
|
||||
params = Richards.Empirical.HaverkampParams().celia1990
|
||||
params['Ks'] = np.log(params['Ks'])
|
||||
E = Richards.Empirical.Haverkamp(M, **params)
|
||||
|
||||
bc = np.array([-61.5,-20.7])
|
||||
h = np.zeros(M.nC) + bc[0]
|
||||
|
||||
prob = Richards.RichardsProblem(M, mapping=E, timeSteps=[(40,3),(60,3)], tolRootFinder=1e-6, debug=False,
|
||||
boundaryConditions=bc, initialConditions=h,
|
||||
doNewton=False, method='mixed')
|
||||
prob.Solver = Solver
|
||||
|
||||
locs = np.r_[5.,10,15]
|
||||
times = prob.times[3:5]
|
||||
rxSat = Richards.RichardsRx(locs, times, 'saturation')
|
||||
rxPre = Richards.RichardsRx(locs, times, 'pressureHead')
|
||||
survey = Richards.RichardsSurvey([rxSat, rxPre])
|
||||
|
||||
prob.pair(survey)
|
||||
|
||||
self.h0 = h
|
||||
self.M = M
|
||||
self.Ks = params['Ks']
|
||||
self.prob = prob
|
||||
self.survey = survey
|
||||
|
||||
def test_Richards_getResidual_Newton(self):
|
||||
self.prob.doNewton = True
|
||||
m = self.Ks
|
||||
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0], self.prob.boundaryConditions), self.h0, plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Richards_getResidual_Picard(self):
|
||||
self.prob.doNewton = False
|
||||
m = self.Ks
|
||||
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0], self.prob.boundaryConditions), self.h0, plotIt=False, expectedOrder=1)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Adjoint(self):
|
||||
v = np.random.rand(self.survey.nD)
|
||||
z = np.random.rand(self.M.nC)
|
||||
Hs = self.prob.fields(self.Ks)
|
||||
vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
|
||||
zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
|
||||
tol = TOL*(10**int(np.log10(np.abs(zJv))))
|
||||
passed = np.abs(vJz - zJv) < tol
|
||||
print 'Richards Adjoint Test - PressureHead'
|
||||
print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Sensitivity(self):
|
||||
mTrue = self.Ks*np.ones(self.M.nC)
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda v: self.prob.Jvec(m, v)]
|
||||
print 'Testing Richards Derivative'
|
||||
passed = checkDerivative(derChk, mTrue, num=4, plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
|
||||
def test_Sensitivity_full(self):
|
||||
mTrue = self.Ks*np.ones(self.M.nC)
|
||||
J = self.prob.Jfull(mTrue)
|
||||
derChk = lambda m: [self.survey.dpred(m), J]
|
||||
print 'Testing Richards Derivative FULL'
|
||||
passed = checkDerivative(derChk, mTrue, num=4, plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
|
||||
class RichardsTests2D(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
M = Mesh.TensorMesh([np.ones(8),np.ones(30)])
|
||||
|
||||
M.setCellGradBC(['neumann','dirichlet'])
|
||||
|
||||
params = Richards.Empirical.HaverkampParams().celia1990
|
||||
params['Ks'] = np.log(params['Ks'])
|
||||
E = Richards.Empirical.Haverkamp(M, **params)
|
||||
|
||||
bc = np.array([-61.5,-20.7])
|
||||
bc = np.r_[np.zeros(M.nCy*2),np.ones(M.nCx)*bc[0],np.ones(M.nCx)*bc[1]]
|
||||
h = np.zeros(M.nC) + bc[0]
|
||||
prob = Richards.RichardsProblem(M,E, timeSteps=[(40,3),(60,3)], boundaryConditions=bc, initialConditions=h, doNewton=False, method='mixed', tolRootFinder=1e-6, debug=False)
|
||||
prob.Solver = Solver
|
||||
|
||||
locs = Utils.ndgrid(np.array([5,7.]),np.array([5,15,25.]))
|
||||
times = prob.times[3:5]
|
||||
rxSat = Richards.RichardsRx(locs, times, 'saturation')
|
||||
rxPre = Richards.RichardsRx(locs, times, 'pressureHead')
|
||||
survey = Richards.RichardsSurvey([rxSat, rxPre])
|
||||
|
||||
prob.pair(survey)
|
||||
|
||||
self.h0 = h
|
||||
self.M = M
|
||||
self.Ks = params['Ks']
|
||||
self.prob = prob
|
||||
self.survey = survey
|
||||
|
||||
def test_Richards_getResidual_Newton(self):
|
||||
self.prob.doNewton = True
|
||||
m = self.Ks
|
||||
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0], self.prob.boundaryConditions), self.h0, plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Richards_getResidual_Picard(self):
|
||||
self.prob.doNewton = False
|
||||
m = self.Ks
|
||||
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0], self.prob.boundaryConditions), self.h0, plotIt=False, expectedOrder=1)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Adjoint(self):
|
||||
v = np.random.rand(self.survey.nD)
|
||||
z = np.random.rand(self.M.nC)
|
||||
Hs = self.prob.fields(self.Ks)
|
||||
vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
|
||||
zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
|
||||
tol = TOL*(10**int(np.log10(np.abs(zJv))))
|
||||
passed = np.abs(vJz - zJv) < tol
|
||||
print '2D: Richards Adjoint Test - PressureHead'
|
||||
print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Sensitivity(self):
|
||||
mTrue = self.Ks*np.ones(self.M.nC)
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda v: self.prob.Jvec(m, v)]
|
||||
print '2D: Testing Richards Derivative'
|
||||
passed = checkDerivative(derChk, mTrue, num=4, plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Sensitivity_full(self):
|
||||
mTrue = self.Ks*np.ones(self.M.nC)
|
||||
J = self.prob.Jfull(mTrue)
|
||||
derChk = lambda m: [self.survey.dpred(m), J]
|
||||
print '2D: Testing Richards Derivative FULL'
|
||||
passed = checkDerivative(derChk, mTrue, num=4, plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
|
||||
|
||||
class RichardsTests3D(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
M = Mesh.TensorMesh([np.ones(8),np.ones(20),np.ones(10)])
|
||||
|
||||
M.setCellGradBC(['neumann','neumann','dirichlet'])
|
||||
|
||||
params = Richards.Empirical.HaverkampParams().celia1990
|
||||
params['Ks'] = np.log(params['Ks'])
|
||||
E = Richards.Empirical.Haverkamp(M, **params)
|
||||
|
||||
bc = np.array([-61.5,-20.7])
|
||||
bc = np.r_[np.zeros(M.nCy*M.nCz*2),np.zeros(M.nCx*M.nCz*2),np.ones(M.nCx*M.nCy)*bc[0],np.ones(M.nCx*M.nCy)*bc[1]]
|
||||
h = np.zeros(M.nC) + bc[0]
|
||||
prob = Richards.RichardsProblem(M,E, timeSteps=[(40,3),(60,3)], boundaryConditions=bc, initialConditions=h, doNewton=False, method='mixed', tolRootFinder=1e-6, debug=False)
|
||||
prob.Solver = Solver
|
||||
|
||||
locs = Utils.ndgrid(np.r_[5,7.],np.r_[5,15.],np.r_[6,8.])
|
||||
times = prob.times[3:5]
|
||||
rxSat = Richards.RichardsRx(locs, times, 'saturation')
|
||||
rxPre = Richards.RichardsRx(locs, times, 'pressureHead')
|
||||
survey = Richards.RichardsSurvey([rxSat, rxPre])
|
||||
|
||||
prob.pair(survey)
|
||||
|
||||
self.h0 = h
|
||||
self.M = M
|
||||
self.Ks = params['Ks']
|
||||
self.prob = prob
|
||||
self.survey = survey
|
||||
|
||||
def test_Richards_getResidual_Newton(self):
|
||||
self.prob.doNewton = True
|
||||
m = self.Ks
|
||||
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0], self.prob.boundaryConditions), self.h0, plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Richards_getResidual_Picard(self):
|
||||
self.prob.doNewton = False
|
||||
m = self.Ks
|
||||
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0], self.prob.boundaryConditions), self.h0, plotIt=False, expectedOrder=1)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Adjoint(self):
|
||||
v = np.random.rand(self.survey.nD)
|
||||
z = np.random.rand(self.M.nC)
|
||||
Hs = self.prob.fields(self.Ks)
|
||||
vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
|
||||
zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
|
||||
tol = TOL*(10**int(np.log10(np.abs(zJv))))
|
||||
passed = np.abs(vJz - zJv) < tol
|
||||
print '3D: Richards Adjoint Test - PressureHead'
|
||||
print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Sensitivity(self):
|
||||
mTrue = self.Ks*np.ones(self.M.nC)
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda v: self.prob.Jvec(m, v)]
|
||||
print '3D: Testing Richards Derivative'
|
||||
passed = checkDerivative(derChk, mTrue, num=4, plotIt=False)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
# def test_Sensitivity_full(self):
|
||||
# mTrue = self.Ks*np.ones(self.M.nC)
|
||||
# J = self.prob.Jfull(mTrue)
|
||||
# derChk = lambda m: [self.survey.dpred(m), J]
|
||||
# print '3D: Testing Richards Derivative FULL'
|
||||
# passed = checkDerivative(derChk, mTrue, num=4, plotIt=False)
|
||||
# self.assertTrue(passed,True)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,12 @@
|
||||
import unittest
|
||||
import sys
|
||||
from SimPEG.FLOW.Examples import Celia1990
|
||||
import numpy as np
|
||||
|
||||
class TestCelia1990(unittest.TestCase):
|
||||
def test_running(self):
|
||||
Celia1990.run(plotIt=False)
|
||||
self.assertTrue(True)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
Reference in New Issue
Block a user