diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py index 7bcd1660..e472ca6a 100644 --- a/SimPEG/forward/Richards.py +++ b/SimPEG/forward/Richards.py @@ -1,6 +1,6 @@ from SimPEG.forward import Problem import numpy as np -from SimPEG.utils import sdiag, mkvc, setKwargs +from SimPEG.utils import sdiag, mkvc, setKwargs, Solver from SimPEG.inverse import NewtonRoot @@ -71,10 +71,10 @@ class RichardsProblem(Problem): @property def dataType(self): """Choose how your data is collected, must be 'saturation' or 'pressureHead'.""" - assert value in ['saturation','pressureHead'], "dataType must be 'saturation' or 'pressureHead'." return self._dataType @dataType.setter def dataType(self, value): + assert value in ['saturation','pressureHead'], "dataType must be 'saturation' or 'pressureHead'." self._dataType = value @@ -83,6 +83,7 @@ class RichardsProblem(Problem): Problem.__init__(self, mesh) self.empirical = empirical self.mesh.setCellGradBC('dirichlet') + self.dataType = 'pressureHead' self.doNewton = False # This also sets the rootFinder algorithm. setKwargs(self, **kwargs) @@ -142,6 +143,8 @@ class RichardsProblem(Problem): B = DdiagGh1*diagAVk2_AVdiagK2*dKa1 + Dz*diagAVk2_AVdiagK2*dKa1 + return Asub, Adiag, B + def getResidual(self, hn, h): """ Where h is the proposed value for the next time iterate (h_{n+1}) @@ -196,10 +199,12 @@ class RichardsProblem(Problem): JvC[ii] = Adiaginv.solve(B*v - Asub*JvC[ii-1]) if self.dataType is 'pressureHead': - Jv = P.Q*np.concatenate(JvC) + Jv = self.P*np.concatenate(JvC) elif self.dataType is 'saturation': dT = self.empirical.moistureContentDeriv(np.concatenate(Hs[1:])) - Jv = P.Q*dT*np.concatenate(JvC) + Jv = self.P*dT*np.concatenate(JvC) + + return Jv def Jt(self, m, v, u=None): if u is None: @@ -207,23 +212,24 @@ class RichardsProblem(Problem): Hs = u if self.dataType is 'pressureHead': - QTz = P.Q.T*z; + PTv = self.P.T*v; elif self.dataType is 'saturation': dT = self.empirical.moistureContentDeriv(np.concatenate(Hs[1:])) - QTz = dT.T*P.Q.T*z + PTv = dT.T*self.P.T*v # This is done via backward substitution. minus = 0 - BJtz = 0 - for ii in range(numel(Hs)-1,-1,-1): + BJtv = 0 + for ii in range(len(Hs)-1,0,-1): Asub, Adiag, B = self.diagsJacobian(Hs[ii-1], Hs[ii]) - #select the correct part of z - # TODO: This might be easier by reshaping the array - zpart = range((ii-2)*Adiag.shape[0], (ii-1)*Adiag.shape[0]) + #select the correct part of v + vpart = range((ii-1)*Adiag.shape[0], (ii)*Adiag.shape[0]) AdiaginvT = Solver(Adiag.T) - JTzC = AdiaginvT.solve(QTz[zpart] - minus) - minus = Asub.T*JTzC # this is now the super diagonal. - BJtz = BJtz + B.T*JTzC + JTvC = AdiaginvT.solve(PTv[vpart] - minus) + minus = Asub.T*JTvC # this is now the super diagonal. + BJtv = BJtv + B.T*JTvC + + return BJtv class Haverkamp(object): diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index 662a7884..3b744f76 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -189,7 +189,7 @@ def Rosenbrock(x, return_g=True, return_H=True): out += (H,) return out if len(out) > 1 else out[0] -def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.9, eps=1e-10): +def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.85, eps=1e-10): """ Basic derivative check diff --git a/SimPEG/tests/test_Richards.py b/SimPEG/tests/test_Richards.py index 8f0be566..e64d2759 100644 --- a/SimPEG/tests/test_Richards.py +++ b/SimPEG/tests/test_Richards.py @@ -1,4 +1,5 @@ import numpy as np +import scipy.sparse as sp import unittest from SimPEG.mesh import TensorMesh from TestUtils import OrderTest, checkDerivative @@ -6,6 +7,8 @@ from scipy.sparse.linalg import dsolve from SimPEG.forward import Richards +TOL = 1E-8 + class RichardsTests(unittest.TestCase): def setUp(self): @@ -13,16 +16,18 @@ class RichardsTests(unittest.TestCase): 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) - prob = Richards.RichardsProblem(M,E) - prob.timeStep = 1 - prob.boundaryConditions = np.array([-61.5,-20.7]) - prob.doNewton = True - prob.method = 'mixed' + bc = np.array([-61.5,-20.7]) + h = np.zeros(M.nC) + bc[0] + prob = Richards.RichardsProblem(M,E, timeStep=30, timeEnd=360, 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.P = P - h = np.zeros(M.nC) + prob.boundaryConditions[0] self.h0 = h self.M = M + self.Ks = Ks self.prob = prob def test_VanGenuchten_moistureContent(self): @@ -71,11 +76,44 @@ class RichardsTests(unittest.TestCase): def test_Richards_getResidual_Newton(self): self.prob.doNewton = True - checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False) + passed = checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False) + self.assertTrue(passed,True) def test_Richards_getResidual_Picard(self): self.prob.doNewton = False - checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False, expectedOrder=1) + passed = checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False, expectedOrder=1) + self.assertTrue(passed,True) + + 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_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) + +