simpeg/tests/flow/test_Richards.py

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

np.random.seed(0)

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,f=Hs))
        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=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,f=Hs))
        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=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=3, 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,f=Hs))
        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=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()