From 9252ed7980d5bde45b5afe3f99fbea1f451f1d7b Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 9 Apr 2014 10:34:47 -0600 Subject: [PATCH] VanG empirical relations --- simpegFLOW/Richards/Empirical.py | 448 ++++++++++++++++++------- simpegFLOW/Richards/RichardsProblem.py | 21 +- simpegFLOW/Tests/test_Richards.py | 24 ++ 3 files changed, 357 insertions(+), 136 deletions(-) diff --git a/simpegFLOW/Richards/Empirical.py b/simpegFLOW/Richards/Empirical.py index 0688eff5..8565c05c 100644 --- a/simpegFLOW/Richards/Empirical.py +++ b/simpegFLOW/Richards/Empirical.py @@ -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() diff --git a/simpegFLOW/Richards/RichardsProblem.py b/simpegFLOW/Richards/RichardsProblem.py index b77b52d5..7ea788f6 100644 --- a/simpegFLOW/Richards/RichardsProblem.py +++ b/simpegFLOW/Richards/RichardsProblem.py @@ -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 diff --git a/simpegFLOW/Tests/test_Richards.py b/simpegFLOW/Tests/test_Richards.py index 84c31c2e..399e3e94 100644 --- a/simpegFLOW/Tests/test_Richards.py +++ b/simpegFLOW/Tests/test_Richards.py @@ -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()