derivative tests

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:

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

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

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

simpegFLOW/Richards/__init__.py

@@ -1,2 +1,2 @@
-from BaseRichards import *
+import Empirical
 from RichardsProblem import *

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