mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-02 15:23:26 +08:00
VanG empirical relations
This commit is contained in:
rowanc1
+321
-127
@@ -1,4 +1,4 @@
|
||||
from SimPEG import Model, Utils, np
|
||||
from SimPEG import Mesh, Model, Utils, np
|
||||
|
||||
|
||||
class RichardsModel(object):
|
||||
@@ -42,19 +42,32 @@ class RichardsModel(object):
|
||||
def kDerivU(self, u, m):
|
||||
return self.kModel.transformDerivU(u, m)
|
||||
|
||||
def plot(self, m):
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
def _ModelProperty(name, model, doc=None, default=None):
|
||||
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 _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):
|
||||
if getattr(self, model, None) is not None:
|
||||
MOD = getattr(self, model)
|
||||
setattr(MOD, name, 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)
|
||||
|
||||
@@ -156,14 +169,14 @@ class _haverkamp_k(Model.BaseNonLinearModel):
|
||||
class Haverkamp(RichardsModel):
|
||||
"""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)
|
||||
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)
|
||||
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):
|
||||
RichardsModel.__init__(self, mesh,
|
||||
@@ -174,138 +187,319 @@ class Haverkamp(RichardsModel):
|
||||
|
||||
|
||||
|
||||
class _vangenuchten_theta(Model.BaseNonLinearModel):
|
||||
|
||||
# class Haverkamp(object):
|
||||
# """docstring for Haverkamp"""
|
||||
theta_s = 0.430
|
||||
theta_r = 0.078
|
||||
alpha = 0.036
|
||||
n = 1.560
|
||||
|
||||
# empiricalModelName = "VanGenuchten"
|
||||
def __init__(self, mesh, **kwargs):
|
||||
Model.BaseNonLinearModel.__init__(self, mesh)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
# theta_s = 0.430
|
||||
# theta_r = 0.078
|
||||
# alpha = 0.036
|
||||
# beta = 3.960
|
||||
# A = 1.175e+06
|
||||
# gamma = 4.74
|
||||
# Ks = np.log(24.96)
|
||||
def setModel(self, m):
|
||||
self._currentModel = m
|
||||
|
||||
# def __init__(self, **kwargs):
|
||||
# Utils.setKwargs(self, **kwargs)
|
||||
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)
|
||||
f[u > 0] = self.theta_s
|
||||
return f
|
||||
|
||||
# def setModel(self, m):
|
||||
# self.Ks = m
|
||||
def transformDerivM(self, u, m):
|
||||
self.setModel(m)
|
||||
|
||||
# def moistureContent(self, h):
|
||||
# f = (self.alpha*(self.theta_s - self.theta_r )/
|
||||
# (self.alpha + abs(h)**self.beta) + self.theta_r)
|
||||
# f[h > 0] = self.theta_s
|
||||
# return f
|
||||
|
||||
# def moistureContentDeriv(self, h):
|
||||
# g = (self.alpha*((self.theta_s - self.theta_r)/
|
||||
# (self.alpha + abs(h)**self.beta)**2)
|
||||
# *(-self.beta*abs(h)**(self.beta-1)*np.sign(h)));
|
||||
# g[h >= 0] = 0
|
||||
# g = Utils.sdiag(g)
|
||||
# return g
|
||||
|
||||
# def hydraulicConductivity(self, h):
|
||||
# f = np.exp(self.Ks)*self.A/(self.A+abs(h)**self.gamma)
|
||||
# if type(self.Ks) is np.ndarray and self.Ks.size > 1:
|
||||
# f[h >= 0] = np.exp(self.Ks[h >= 0])
|
||||
# else:
|
||||
# f[h >= 0] = np.exp(self.Ks)
|
||||
# return f
|
||||
|
||||
# def hydraulicConductivityModelDeriv(self, h):
|
||||
# #A
|
||||
# # dA = np.exp(self.Ks)/(self.A+abs(h)**self.gamma) - np.exp(self.Ks)*self.A/(self.A+abs(h)**self.gamma)**2;
|
||||
# #gamma
|
||||
# # dgamma = -(self.A*np.exp(self.Ks)*np.log(abs(h))*abs(h)**self.gamma)/(self.A + abs(h)**self.gamma)**2;
|
||||
# return Utils.sdiag(self.hydraulicConductivity(h)) # This assumes that the the model is Ks
|
||||
|
||||
# def hydraulicConductivityDeriv(self, h):
|
||||
# g = -(np.exp(self.Ks)*self.A*self.gamma*abs(h)**(self.gamma-1)*np.sign(h))/((self.A+abs(h)**self.gamma)**2)
|
||||
# g[h >= 0] = 0
|
||||
# g = Utils.sdiag(g)
|
||||
# return g
|
||||
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(object):
|
||||
# """
|
||||
class _vangenuchten_k(Model.BaseNonLinearModel):
|
||||
|
||||
# .. math::
|
||||
I = 0.500
|
||||
alpha = 0.036
|
||||
n = 1.560
|
||||
Ks = np.log(24.96)
|
||||
|
||||
# \\theta(h) = \\frac{\\alpha (\\theta_s - \\theta_r)}{\\alpha + |h|^\\beta} + \\theta_r
|
||||
def __init__(self, mesh, **kwargs):
|
||||
Model.BaseNonLinearModel.__init__(self, mesh)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
# Where parameters alpha, beta, gamma, A are constants in the media;
|
||||
# theta_r and theta_s are the residual and saturated moisture
|
||||
# contents; and K_s is the saturated hydraulic conductivity.
|
||||
def setModel(self, m):
|
||||
self._currentModel = m
|
||||
#TODO: Fix me!
|
||||
self.Ks = m
|
||||
|
||||
# Celia1990
|
||||
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
|
||||
|
||||
# empiricalModelName = "VanGenuchten"
|
||||
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 type(self.Ks) is np.ndarray and self.Ks.size > 1:
|
||||
f[u >= 0] = np.exp(self.Ks[u >= 0])
|
||||
else:
|
||||
f[u >= 0] = np.exp(self.Ks)
|
||||
return f
|
||||
|
||||
# theta_s = 0.430
|
||||
# theta_r = 0.078
|
||||
# alpha = 0.036
|
||||
# n = 1.560
|
||||
# beta = 3.960
|
||||
# I = 0.500
|
||||
# Ks = np.log(24.96)
|
||||
|
||||
# def __init__(self, **kwargs):
|
||||
# Utils.setKwargs(self, **kwargs)
|
||||
|
||||
# def setModel(self, m):
|
||||
# self.Ks = m
|
||||
|
||||
# def moistureContent(self, h):
|
||||
# m = 1 - 1.0/self.n;
|
||||
# f = (( self.theta_s - self.theta_r )/
|
||||
# ((1+abs(self.alpha*h)**self.n)**m) + self.theta_r)
|
||||
# f[h > 0] = self.theta_s
|
||||
# return f
|
||||
|
||||
# def moistureContentDeriv(self, h):
|
||||
# g = -self.alpha*self.n*abs(self.alpha*h)**(self.n - 1)*np.sign(self.alpha*h)*(1./self.n - 1)*(self.theta_r - self.theta_s)*(abs(self.alpha*h)**self.n + 1)**(1./self.n - 2)
|
||||
# g[h > 0] = 0
|
||||
# g = Utils.sdiag(g)
|
||||
# return g
|
||||
|
||||
# def hydraulicConductivity(self, h):
|
||||
# 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*h)**n)**m)
|
||||
# f = np.exp(Ks)*theta_e**I* ( ( 1.0 - ( 1.0 - theta_e**(1.0/m) )**m )**2 )
|
||||
# if type(self.Ks) is np.ndarray and self.Ks.size > 1:
|
||||
# f[h >= 0] = np.exp(self.Ks[h >= 0])
|
||||
# else:
|
||||
# f[h >= 0] = np.exp(self.Ks)
|
||||
# return f
|
||||
|
||||
# def hydraulicConductivityModelDeriv(self, h):
|
||||
def transformDerivM(self, u, m):
|
||||
self.setModel(m)
|
||||
# #alpha
|
||||
# # dA = I*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1.0/n - 1)*((abs(alpha*h)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*h)**n + 1)**(1.0/n - 2) - (2*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1.0/n - 1)*((abs(alpha*h)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*h)**n + 1)**(1.0/n - 2))/(((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)*(1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n));
|
||||
# # 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*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))*(1 - 1.0/((abs(alpha*h)**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*h)**n + 1)*(abs(alpha*h)**n + 1)**(1.0/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1.0/n - 1)*(abs(alpha*h)**n + 1)**(1.0/n - 2))/((1.0/n - 1)*((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)) - np.log((abs(alpha*h)**n + 1)**(1.0/n - 1))/(n**2*(1.0/n - 1)**2*((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))))/(1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n))*((abs(alpha*h)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*h)**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*h)**n + 1)*(abs(alpha*h)**n + 1)**(1.0/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1.0/n - 1)*(abs(alpha*h)**n + 1)**(1.0/n - 2))*((abs(alpha*h)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2;
|
||||
# # 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*h)**n + 1)**(1.0/n - 1))*((abs(alpha*h)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2;
|
||||
# return Utils.sdiag(self.hydraulicConductivity(h)) # This assumes that the the model is Ks
|
||||
# # 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 hydraulicConductivityDeriv(self, h):
|
||||
# alpha = self.alpha
|
||||
# I = self.I
|
||||
# n = self.n
|
||||
# Ks = self.Ks
|
||||
# m = 1.0 - 1.0/n
|
||||
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*h)**(n - 1.0)*np.sign(alpha*h)*(1.0/n - 1.0)*((abs(alpha*h)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*h)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1.0/n - 1)*((abs(alpha*h)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*h)**n + 1)**(1.0/n - 2))/(((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)*(1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n))
|
||||
# g[h >= 0] = 0
|
||||
# g = Utils.sdiag(g)
|
||||
# return g
|
||||
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(RichardsModel):
|
||||
"""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):
|
||||
RichardsModel.__init__(self, mesh,
|
||||
_haverkamp_theta(mesh),
|
||||
_haverkamp_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()
|
||||
|
||||
@@ -69,6 +69,9 @@ class RichardsProblem(Problem.BaseProblem):
|
||||
surveyPair = RichardsSurvey
|
||||
modelPair = RichardsModel
|
||||
|
||||
Solver = Solver
|
||||
solverOpts = {}
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, model, **kwargs)
|
||||
|
||||
@@ -102,7 +105,7 @@ class RichardsProblem(Problem.BaseProblem):
|
||||
def rootFinder(self):
|
||||
"""Root-finding Algorithm"""
|
||||
if getattr(self, '_rootFinder', None) is None:
|
||||
self._rootFinder = Optimization.NewtonRoot(doLS=self.doNewton)
|
||||
self._rootFinder = Optimization.NewtonRoot(doLS=self.doNewton, Solver=self.Solver)
|
||||
return self._rootFinder
|
||||
|
||||
def fields(self, m):
|
||||
@@ -121,9 +124,9 @@ class RichardsProblem(Problem.BaseProblem):
|
||||
if self.mesh.dim == 1:
|
||||
Dz = self.mesh.faceDivx
|
||||
elif self.mesh.dim == 2:
|
||||
Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.nFv[0]), self.mesh.faceDivy),format='csr')
|
||||
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.nFv[0]+self.mesh.nFv[1]), self.mesh.faceDivz),format='csr')
|
||||
Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.vnF[0]+self.mesh.vnF[1]), self.mesh.faceDivz),format='csr')
|
||||
|
||||
bc = self.boundaryConditions
|
||||
dt = self.timeStep
|
||||
@@ -175,9 +178,9 @@ class RichardsProblem(Problem.BaseProblem):
|
||||
if self.mesh.dim == 1:
|
||||
Dz = self.mesh.faceDivx
|
||||
elif self.mesh.dim == 2:
|
||||
Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.nFv[0]), self.mesh.faceDivy),format='csr')
|
||||
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.nFv[0]+self.mesh.nFv[1]), self.mesh.faceDivz),format='csr')
|
||||
Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.vnF[0]+self.mesh.vnF[1]), self.mesh.faceDivz),format='csr')
|
||||
|
||||
bc = self.boundaryConditions
|
||||
dt = self.timeStep
|
||||
@@ -220,7 +223,7 @@ class RichardsProblem(Problem.BaseProblem):
|
||||
A = As + Ad
|
||||
B = np.array(sp.vstack(Bs).todense())
|
||||
|
||||
Ainv = Solver(A)
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
P = self.survey.projectFieldsDeriv(u, m)
|
||||
J = P * Ainv.solve(B)
|
||||
return J
|
||||
@@ -233,7 +236,7 @@ class RichardsProblem(Problem.BaseProblem):
|
||||
|
||||
# This is done via forward substitution.
|
||||
temp, Adiag, B = self.diagsJacobian(m, u[0],u[1])
|
||||
Adiaginv = Solver(Adiag)
|
||||
Adiaginv = self.Solver(Adiag, **self.solverOpts)
|
||||
JvC[0] = Adiaginv.solve(B*v)
|
||||
|
||||
# M = @(x) tril(Adiag)\(diag(Adiag).*(triu(Adiag)\x));
|
||||
@@ -241,7 +244,7 @@ class RichardsProblem(Problem.BaseProblem):
|
||||
|
||||
for ii in range(1,len(u)-1):
|
||||
Asub, Adiag, B = self.diagsJacobian(m, u[ii],u[ii+1])
|
||||
Adiaginv = Solver(Adiag)
|
||||
Adiaginv = self.Solver(Adiag, **self.solverOpts)
|
||||
JvC[ii] = Adiaginv.solve(B*v - Asub*JvC[ii-1])
|
||||
|
||||
P = self.survey.projectFieldsDeriv(u, m)
|
||||
@@ -261,7 +264,7 @@ class RichardsProblem(Problem.BaseProblem):
|
||||
Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii])
|
||||
#select the correct part of v
|
||||
vpart = range((ii-1)*Adiag.shape[0], (ii)*Adiag.shape[0])
|
||||
AdiaginvT = Solver(Adiag.T)
|
||||
AdiaginvT = self.Solver(Adiag.T, **self.solverOpts)
|
||||
JTvC = AdiaginvT.solve(PTv[vpart] - minus)
|
||||
minus = Asub.T*JTvC # this is now the super diagonal.
|
||||
BJtv = BJtv + B.T*JTvC
|
||||
|
||||
@@ -17,6 +17,14 @@ class TestModels(unittest.TestCase):
|
||||
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])
|
||||
@@ -34,6 +42,22 @@ class TestModels(unittest.TestCase):
|
||||
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)
|
||||
|
||||
# def test_Haverkamp_hydraulicConductivity(self):
|
||||
# print 'Haverkamp_hydraulicConductivity'
|
||||
# hav = Richards.Haverkamp()
|
||||
|
||||
Reference in New Issue
Block a user