diff --git a/docs/api_Richards.rst b/docs/api_Richards.rst index 9bb787e3..d1d51809 100644 --- a/docs/api_Richards.rst +++ b/docs/api_Richards.rst @@ -31,12 +31,14 @@ equation can be written in the continuous form as: However, it can be shown that this does not conserve mass in the discrete formulation. +Here we reproduce the results from Ceilia et al. (1990): + +.. plot:: examples/richards/comparisonToCeiliaEtAl1990.py Richards ======== -.. automodule:: simpegFLOW.Richards +.. automodule:: simpegFLOW.Richards.Empirical :show-inheritance: :members: :undoc-members: - :inherited-members: diff --git a/docs/examples/richards/comparisonToCeiliaEtAl1990.py b/docs/examples/richards/comparisonToCeiliaEtAl1990.py new file mode 100644 index 00000000..22abb2f2 --- /dev/null +++ b/docs/examples/richards/comparisonToCeiliaEtAl1990.py @@ -0,0 +1,45 @@ +from SimPEG import * +from simpegFLOW import Richards +import matplotlib.pyplot as plt + +M = Mesh.TensorMesh([np.ones(40)]) +M.setCellGradBC('dirichlet') +params = Richards.Empirical.HaverkampParams().celia1990 +model = Richards.Empirical.Haverkamp(M, **params) + +bc = np.array([-61.5,-20.7]) +h = np.zeros(M.nC) + bc[0] + +def getFields(timeStep,method): + prob = Richards.RichardsProblem(M,model, timeStep=timeStep, timeEnd=360, + 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')) +plt.show() diff --git a/simpegFLOW/Richards/BaseRichards.py b/simpegFLOW/Richards/Empirical.py similarity index 83% rename from simpegFLOW/Richards/BaseRichards.py rename to simpegFLOW/Richards/Empirical.py index 04359512..0688eff5 100644 --- a/simpegFLOW/Richards/BaseRichards.py +++ b/simpegFLOW/Richards/Empirical.py @@ -43,7 +43,42 @@ class RichardsModel(object): return self.kModel.transformDerivU(u, m) -class BaseHaverkamp_theta(Model.BaseNonLinearModel): +def _ModelProperty(name, model, doc=None, default=None): + + def fget(self): + 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) + + 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':np.log(9.44e-03), 'A':1.175e+06, + 'gamma':4.74} + + +class _haverkamp_theta(Model.BaseNonLinearModel): theta_s = 0.430 theta_r = 0.078 @@ -77,7 +112,7 @@ class BaseHaverkamp_theta(Model.BaseNonLinearModel): return g -class BaseHaverkamp_k(Model.BaseNonLinearModel): +class _haverkamp_k(Model.BaseNonLinearModel): A = 1.175e+06 gamma = 4.74 @@ -118,6 +153,26 @@ class BaseHaverkamp_k(Model.BaseNonLinearModel): g = Utils.sdiag(g) return g +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) + + 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, + _haverkamp_theta(mesh), + _haverkamp_k(mesh)) + Utils.setKwargs(self, **kwargs) + + + # class Haverkamp(object): diff --git a/simpegFLOW/Richards/RichardsProblem.py b/simpegFLOW/Richards/RichardsProblem.py index 870d857b..5016f2f6 100644 --- a/simpegFLOW/Richards/RichardsProblem.py +++ b/simpegFLOW/Richards/RichardsProblem.py @@ -18,13 +18,40 @@ class RichardsData(Data.BaseData): assert value in ['saturation','pressureHead'], "dataType must be 'saturation' or 'pressureHead'." self._dataType = value - def projectFields(self, u): - u = np.concatenate(u[1:]) + @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)) + + 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)) + + def projectFields(self, U, m): + + u = np.concatenate(U[1:]) + if self.dataType == 'saturation': - #TODO: Fix this: - u = self.prob.model.theta(MODEL, u) + u = self.prob.model.theta(u, m) return self.P*u + def projectFieldsDerivU(self, U, m): + + u = np.concatenate(U[1:]) + + if self.dataType == 'pressureHead': + return self.P + elif self.dataType == 'saturation': + dT = self.model.thetaDerivU(u, m) + return self.P*dT + class RichardsProblem(Problem.BaseProblem): """docstring for RichardsProblem""" @@ -73,11 +100,11 @@ class RichardsProblem(Problem.BaseProblem): self._doNewton = value def fields(self, m): - Hs = range(self.numIts+1) - Hs[0] = self.initialConditions + u = range(self.numIts+1) + u[0] = self.initialConditions for ii in range(self.numIts): - Hs[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, Hs[ii], hn1m, return_g=return_g), Hs[ii]) - return Hs + u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, return_g=return_g), u[ii]) + return u def diagsJacobian(self, m, hn, hn1): @@ -172,14 +199,14 @@ class RichardsProblem(Problem.BaseProblem): return r, J - def fullJ(self, m, u=None): + def Jfull(self, m, u=None): if u is None: u = self.field(m) - Hs = u - nn = len(Hs)-1 + + nn = len(u)-1 Asubs, Adiags, Bs = range(nn), range(nn), range(nn) for ii in range(nn): - Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, Hs[ii],Hs[ii+1]) + Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, u[ii],u[ii+1]) 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]) @@ -195,46 +222,38 @@ class RichardsProblem(Problem.BaseProblem): def Jvec(self, m, v, u=None): if u is None: u = self.field(m) - Hs = u - JvC = range(len(Hs)-1) # Cell to hold each row of the long vector. + + JvC = range(len(u)-1) # Cell to hold each row of the long vector. # This is done via forward substitution. - temp, Adiag, B = self.diagsJacobian(m, Hs[0],Hs[1]) + temp, Adiag, B = self.diagsJacobian(m, u[0],u[1]) Adiaginv = Solver(Adiag) JvC[0] = Adiaginv.solve(B*v) # M = @(x) tril(Adiag)\(diag(Adiag).*(triu(Adiag)\x)); # JvC{1} = bicgstab(Adiag,(B*v),tolbcg,500,M); - for ii in range(1,len(Hs)-1): - Asub, Adiag, B = self.diagsJacobian(m, Hs[ii],Hs[ii+1]) + for ii in range(1,len(u)-1): + Asub, Adiag, B = self.diagsJacobian(m, u[ii],u[ii+1]) Adiaginv = Solver(Adiag) JvC[ii] = Adiaginv.solve(B*v - Asub*JvC[ii-1]) - if self.dataType == 'pressureHead': - Jv = self.P*np.concatenate(JvC) - elif self.dataType == 'saturation': - dT = self.model.thetaDerivU(np.concatenate(Hs[1:]), m) - Jv = self.P*dT*np.concatenate(JvC) - - return Jv + P = self.data.projectFieldsDerivU(u, m) + return P * np.concatenate(JvC) def Jtvec(self, m, v, u=None): if u is None: u = self.field(m) - Hs = u - if self.dataType == 'pressureHead': - PTv = self.P.T*v; - elif self.dataType == 'saturation': - dT = self.model.thetaDerivU(np.concatenate(Hs[1:]), m) - PTv = dT.T*self.P.T*v + + P = self.data.projectFieldsDerivU(u, m) + PTv = P.T*v # This is done via backward substitution. minus = 0 BJtv = 0 - for ii in range(len(Hs)-1,0,-1): - Asub, Adiag, B = self.diagsJacobian(m, Hs[ii-1], Hs[ii]) + for ii in range(len(u)-1,0,-1): + 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) @@ -243,36 +262,3 @@ class RichardsProblem(Problem.BaseProblem): BJtv = BJtv + B.T*JTvC return BJtv - - - -if __name__ == '__main__': - from SimPEG import * - import Richards - import matplotlib.pyplot as plt - M = Mesh.TensorMesh([np.ones(40)]) - Ks = 9.4400e-03 - E = Richards.Haverkamp(Ks=np.log(Ks), A=1.1750e+06, gamma=4.74, alpha=1.6110e+06, theta_s=0.287, theta_r=0.075, beta=3.96) - bc = np.array([-61.5,-20.7]) - h = np.zeros(M.nC) + bc[0] - - # data = R - prob = Richards.RichardsProblem(M,E, timeStep=10, timeEnd=100, boundaryConditions=bc, initialConditions=h, doNewton=False, method='mixed') - - q = sp.csr_matrix((np.ones(4),(np.arange(4),np.array([20, 30, 35, 38]))), shape=(4,M.nCx)) - P = sp.kron(sp.identity(prob.numIts),q) - - prob.dataType = 'pressureHead' - mTrue = np.ones(M.nC)*np.log(Ks) - stdev = 0.01 # The standard deviation for the noise - data = prob.createSyntheticData(mTrue,std=stdev, P=P) - p = plt.plot(data.dobs.reshape((-1,4))) - plt.show() - # opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) - # reg = Regularization.Tikhonov(model) - # inv = Inversion.BaseInversion(prob, reg, opt, beta0=1e4) - # derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)] - # print inv.dataObj(mTrue*0+np.log(1e-5)) - # print inv.dataObj(mTrue) - # tests.checkDerivative(derChk, mTrue, plotIt=False) - diff --git a/simpegFLOW/Richards/__init__.py b/simpegFLOW/Richards/__init__.py index bf9832f7..6e25d292 100644 --- a/simpegFLOW/Richards/__init__.py +++ b/simpegFLOW/Richards/__init__.py @@ -1,2 +1,2 @@ -from BaseRichards import * +import Empirical from RichardsProblem import * diff --git a/simpegFLOW/Tests/test_Richards.py b/simpegFLOW/Tests/test_Richards.py index 08befb01..6dcde6c5 100644 --- a/simpegFLOW/Tests/test_Richards.py +++ b/simpegFLOW/Tests/test_Richards.py @@ -2,7 +2,7 @@ import unittest from SimPEG import * from SimPEG.Tests.TestUtils import OrderTest, checkDerivative from scipy.sparse.linalg import dsolve -import simpegFLOW.Richards as Richards +from simpegFLOW import Richards TOL = 1E-8 @@ -10,7 +10,7 @@ class TestModels(unittest.TestCase): def test_BaseHaverkamp_Theta(self): mesh = Mesh.TensorMesh([50]) - hav = Richards.BaseHaverkamp_theta(mesh) + hav = Richards.Empirical._haverkamp_theta(mesh) m = np.random.randn(50) def wrapper(u): return hav.transform(u, m), hav.transformDerivU(u, m) @@ -20,14 +20,14 @@ class TestModels(unittest.TestCase): def test_BaseHaverkamp_k(self): mesh = Mesh.TensorMesh([50]) - hav = Richards.BaseHaverkamp_k(mesh) + 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.BaseHaverkamp_k(mesh) + hav = Richards.Empirical._haverkamp_k(mesh) u = np.random.randn(50) def wrapper(m): return hav.transform(u, m), hav.transformDerivM(u, m) @@ -102,82 +102,99 @@ class TestModels(unittest.TestCase): # self.assertTrue(passed,True) -# class RichardsTests1D(unittest.TestCase): +class RichardsTests1D(unittest.TestCase): -# def setUp(self): -# M = Mesh.TensorMesh([np.ones(20)]) -# M.setCellGradBC('dirichlet') + def setUp(self): + M = Mesh.TensorMesh([np.ones(20)]) + M.setCellGradBC('dirichlet') -# Ks = 9.4400e-03 -# E = Richards.Haverkamp(Ks=np.log(Ks), A=1.1750e+06, gamma=4.74, alpha=1.6110e+06, theta_s=0.287, theta_r=0.075, beta=3.96) + params = Richards.Empirical.HaverkampParams().celia1990 + model = Richards.Empirical.Haverkamp(M, **params) -# bc = np.array([-61.5,-20.7]) -# h = np.zeros(M.nC) + bc[0] -# prob = Richards.RichardsProblem(M,E, timeStep=60, timeEnd=180, boundaryConditions=bc, initialConditions=h, doNewton=False, method='mixed') + bc = np.array([-61.5,-20.7]) + h = np.zeros(M.nC) + bc[0] -# q = sp.csr_matrix((np.ones(3),(np.arange(3),np.array([5,10,15]))),shape=(3,M.nC)) -# P = sp.kron(sp.identity(prob.numIts),q) -# prob.P = P + prob = Richards.RichardsProblem(M,model, timeStep=60, timeEnd=180, + boundaryConditions=bc, initialConditions=h, + doNewton=False, method='mixed') -# self.h0 = h -# self.M = M -# self.Ks = Ks -# self.prob = prob + q = sp.csr_matrix((np.ones(3),(np.arange(3),np.array([5,10,15]))),shape=(3,M.nC)) + P = sp.kron(sp.identity(prob.numIts),q) + data = Richards.RichardsData(P=P) -# def test_Richards_getResidual_Newton(self): -# self.prob.doNewton = True -# passed = checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False) -# self.assertTrue(passed,True) + prob.pair(data) -# def test_Richards_getResidual_Picard(self): -# self.prob.doNewton = False -# passed = checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False, expectedOrder=1) -# self.assertTrue(passed,True) + self.h0 = h + self.M = M + self.Ks = params['Ks'] + self.prob = prob + self.data = data -# def test_Adjoint_PressureHead(self): -# self.prob.dataType = 'pressureHead' -# Ks = self.Ks -# v = np.random.rand(self.prob.P.shape[0]) -# z = np.random.rand(self.M.nC) -# Hs = self.prob.field(np.log(Ks)) -# vJz = v.dot(self.prob.J(np.log(Ks),z,u=Hs)) -# zJv = z.dot(self.prob.Jt(np.log(Ks),v,u=Hs)) -# tol = TOL*(10**int(np.log10(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_Richards_getResidual_Newton(self): + self.prob.doNewton = True + m = self.Ks + passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0,hn1), 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.h0, plotIt=False, expectedOrder=1) + self.assertTrue(passed,True) -# def test_Adjoint_Saturation(self): -# self.prob.dataType = 'saturation' -# Ks = self.Ks -# v = np.random.rand(self.prob.P.shape[0]) -# z = np.random.rand(self.M.nC) -# Hs = self.prob.field(np.log(Ks)) -# vJz = v.dot(self.prob.J(np.log(Ks),z,u=Hs)) -# zJv = z.dot(self.prob.Jt(np.log(Ks),v,u=Hs)) -# tol = TOL*(10**int(np.log10(zJv))) -# passed = np.abs(vJz - zJv) < tol -# print 'Richards Adjoint Test - Saturation' -# print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol) -# self.assertTrue(passed,True) + def test_Adjoint_PressureHead(self): + self.prob.dataType = 'pressureHead' + v = np.random.rand(self.data.P.shape[0]) + 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(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): -# self.prob.dataType = 'pressureHead' -# mTrue = np.ones(self.M.nC)*np.log(self.Ks) -# stdev = 0.01 # The standard deviation for the noise -# dobs = self.prob.createSyntheticData(mTrue,std=stdev)[0] -# self.prob.dobs = dobs -# self.prob.std = dobs*0 + stdev -# Hs = self.prob.field(mTrue) -# opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) -# reg = regularization.Regularization(self.M) -# inv = inverse.Inversion(self.prob, reg, opt, beta0=1e4) -# derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)] -# print 'Testing Richards Derivative' -# passed = checkDerivative(derChk, mTrue, num=5, plotIt=False) -# self.assertTrue(passed,True) + def test_Adjoint_Saturation(self): + self.prob.dataType = 'saturation' + v = np.random.rand(self.data.P.shape[0]) + 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(zJv))) + passed = np.abs(vJz - zJv) < tol + print 'Richards Adjoint Test - Saturation' + print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol) + self.assertTrue(passed,True) + + def test_SensitivityPressureHead(self): + self.prob.dataType = 'pressureHead' + self.prob.unpair() + mTrue = np.ones(self.M.nC)*self.Ks + stdev = 0.01 # The standard deviation for the noise + data = self.prob.createSyntheticData(mTrue, std=stdev, P=self.data.P) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + reg = Regularization.Tikhonov(Model.BaseModel(self.M)) + obj = ObjFunction.BaseObjFunction(data, reg) + derChk = lambda m: [obj.dataObj(m), obj.dataObjDeriv(m)] + print 'Testing Richards Derivative - Pressure Head' + passed = checkDerivative(derChk, mTrue, num=5, plotIt=False) + self.assertTrue(passed,True) + + def test_SensitivitySaturation(self): + self.prob.unpair() + self.prob.dataType = 'saturation' + mTrue = np.ones(self.M.nC)*self.Ks + stdev = 0.01 # The standard deviation for the noise + data = self.prob.createSyntheticData(mTrue, std=stdev, P=self.data.P) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + reg = Regularization.Tikhonov(Model.BaseModel(self.M)) + obj = ObjFunction.BaseObjFunction(data, reg) + derChk = lambda m: [obj.dataObj(m), obj.dataObjDeriv(m)] + print 'Testing Richards Derivative - Saturation' + passed = checkDerivative(derChk, mTrue, num=5, plotIt=False) + self.assertTrue(passed,True) @@ -245,7 +262,7 @@ class TestModels(unittest.TestCase): # def test_Sensitivity(self): # self.prob.dataType = 'pressureHead' -# mTrue = np.ones(self.M.nC)*np.log(self.Ks) +# mTrue = np.ones(self.M.nC)*self.Ks # stdev = 0.01 # The standard deviation for the noise # dobs = self.prob.createSyntheticData(mTrue,std=stdev)[0] # self.prob.dobs = dobs