diff --git a/.travis.yml b/.travis.yml index 08ab3851..6a8a995f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -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: diff --git a/SimPEG/Examples/DCfwd.py b/SimPEG/Examples/DCfwd.py index 3cd5f66a..509de6bc 100644 --- a/SimPEG/Examples/DCfwd.py +++ b/SimPEG/Examples/DCfwd.py @@ -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() diff --git a/SimPEG/FLOW/Examples/Celia1990.py b/SimPEG/FLOW/Examples/Celia1990.py new file mode 100644 index 00000000..c866c329 --- /dev/null +++ b/SimPEG/FLOW/Examples/Celia1990.py @@ -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() diff --git a/SimPEG/FLOW/Examples/__init__.py b/SimPEG/FLOW/Examples/__init__.py new file mode 100644 index 00000000..7a894e11 --- /dev/null +++ b/SimPEG/FLOW/Examples/__init__.py @@ -0,0 +1 @@ +import Celia1990 diff --git a/SimPEG/FLOW/Richards/Empirical.py b/SimPEG/FLOW/Richards/Empirical.py new file mode 100644 index 00000000..1bb163b5 --- /dev/null +++ b/SimPEG/FLOW/Richards/Empirical.py @@ -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() diff --git a/SimPEG/FLOW/Richards/RichardsProblem.py b/SimPEG/FLOW/Richards/RichardsProblem.py new file mode 100644 index 00000000..7c61c221 --- /dev/null +++ b/SimPEG/FLOW/Richards/RichardsProblem.py @@ -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 diff --git a/SimPEG/FLOW/Richards/__init__.py b/SimPEG/FLOW/Richards/__init__.py new file mode 100644 index 00000000..6e25d292 --- /dev/null +++ b/SimPEG/FLOW/Richards/__init__.py @@ -0,0 +1,2 @@ +import Empirical +from RichardsProblem import * diff --git a/SimPEG/FLOW/__init__.py b/SimPEG/FLOW/__init__.py new file mode 100644 index 00000000..622d2f2f --- /dev/null +++ b/SimPEG/FLOW/__init__.py @@ -0,0 +1 @@ +import Richards diff --git a/docs/flow/index.rst b/docs/flow/index.rst new file mode 100644 index 00000000..3c8083fc --- /dev/null +++ b/docs/flow/index.rst @@ -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: diff --git a/docs/index.rst b/docs/index.rst index 628a079b..c66b2741 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -79,6 +79,14 @@ Utility Codes api_Tests +Packages +******** + +.. toctree:: + :maxdepth: 3 + + flow/index + Developer's Documentation ************************* diff --git a/tests/flow/__init__.py b/tests/flow/__init__.py new file mode 100644 index 00000000..38d84328 --- /dev/null +++ b/tests/flow/__init__.py @@ -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) diff --git a/tests/flow/test_Richards.py b/tests/flow/test_Richards.py new file mode 100644 index 00000000..b079e371 --- /dev/null +++ b/tests/flow/test_Richards.py @@ -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() diff --git a/tests/flow/test_examples.py b/tests/flow/test_examples.py new file mode 100644 index 00000000..bc519e6d --- /dev/null +++ b/tests/flow/test_examples.py @@ -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()