Bug Fixes and Adjoint Test.

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):

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
 

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)
+
+