<s> import unittest <EOL> from SimPEG import * <EOL> from SimPEG . Tests import OrderTest , checkDerivative <EOL> from scipy . sparse . linalg import dsolve <EOL> from SimPEG . FLOW import Richards <EOL> try : <EOL> from pymatsolver import MumpsSolver <EOL> Solver = MumpsSolver <EOL> except Exception , e : <EOL> pass <EOL> TOL = <NUM_LIT> <EOL> np . random . seed ( <NUM_LIT:0> ) <EOL> class TestModels ( unittest . TestCase ) : <EOL> def test_BaseHaverkamp_Theta ( self ) : <EOL> mesh = Mesh . TensorMesh ( [ <NUM_LIT:50> ] ) <EOL> hav = Richards . Empirical . _haverkamp_theta ( mesh ) <EOL> m = np . random . randn ( <NUM_LIT:50> ) <EOL> def wrapper ( u ) : <EOL> return hav . transform ( u , m ) , hav . transformDerivU ( u , m ) <EOL> passed = checkDerivative ( wrapper , np . random . randn ( <NUM_LIT:50> ) , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_vangenuchten_theta ( self ) : <EOL> mesh = Mesh . TensorMesh ( [ <NUM_LIT:50> ] ) <EOL> hav = Richards . Empirical . _vangenuchten_theta ( mesh ) <EOL> m = np . random . randn ( <NUM_LIT:50> ) <EOL> def wrapper ( u ) : <EOL> return hav . transform ( u , m ) , hav . transformDerivU ( u , m ) <EOL> passed = checkDerivative ( wrapper , np . random . randn ( <NUM_LIT:50> ) , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_BaseHaverkamp_k ( self ) : <EOL> mesh = Mesh . TensorMesh ( [ <NUM_LIT:50> ] ) <EOL> hav = Richards . Empirical . _haverkamp_k ( mesh ) <EOL> m = np . random . randn ( <NUM_LIT:50> ) <EOL> def wrapper ( u ) : <EOL> return hav . transform ( u , m ) , hav . transformDerivU ( u , m ) <EOL> passed = checkDerivative ( wrapper , np . random . randn ( <NUM_LIT:50> ) , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> hav = Richards . Empirical . _haverkamp_k ( mesh ) <EOL> u = np . random . randn ( <NUM_LIT:50> ) <EOL> def wrapper ( m ) : <EOL> return hav . transform ( u , m ) , hav . transformDerivM ( u , m ) <EOL> passed = checkDerivative ( wrapper , np . random . randn ( <NUM_LIT:50> ) , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_vangenuchten_k ( self ) : <EOL> mesh = Mesh . TensorMesh ( [ <NUM_LIT:50> ] ) <EOL> hav = Richards . Empirical . _vangenuchten_k ( mesh ) <EOL> m = np . random . randn ( <NUM_LIT:50> ) <EOL> def wrapper ( u ) : <EOL> return hav . transform ( u , m ) , hav . transformDerivU ( u , m ) <EOL> passed = checkDerivative ( wrapper , np . random . randn ( <NUM_LIT:50> ) , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> hav = Richards . Empirical . _vangenuchten_k ( mesh ) <EOL> u = np . random . randn ( <NUM_LIT:50> ) <EOL> def wrapper ( m ) : <EOL> return hav . transform ( u , m ) , hav . transformDerivM ( u , m ) <EOL> passed = checkDerivative ( wrapper , np . random . randn ( <NUM_LIT:50> ) , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> class RichardsTests1D ( unittest . TestCase ) : <EOL> def setUp ( self ) : <EOL> M = Mesh . TensorMesh ( [ np . ones ( <NUM_LIT:20> ) ] ) <EOL> M . setCellGradBC ( '<STR_LIT>' ) <EOL> params = Richards . Empirical . HaverkampParams ( ) . celia1990 <EOL> params [ '<STR_LIT>' ] = np . log ( params [ '<STR_LIT>' ] ) <EOL> E = Richards . Empirical . Haverkamp ( M , ** params ) <EOL> bc = np . array ( [ - <NUM_LIT> , - <NUM_LIT> ] ) <EOL> h = np . zeros ( M . nC ) + bc [ <NUM_LIT:0> ] <EOL> prob = Richards . RichardsProblem ( M , mapping = E , timeSteps = [ ( <NUM_LIT> , <NUM_LIT:3> ) , ( <NUM_LIT> , <NUM_LIT:3> ) ] , tolRootFinder = <NUM_LIT> , debug = False , <EOL> boundaryConditions = bc , initialConditions = h , <EOL> doNewton = False , method = '<STR_LIT>' ) <EOL> prob . Solver = Solver <EOL> locs = np . r_ [ <NUM_LIT> , <NUM_LIT:10> , <NUM_LIT:15> ] <EOL> times = prob . times [ <NUM_LIT:3> : <NUM_LIT:5> ] <EOL> rxSat = Richards . RichardsRx ( locs , times , '<STR_LIT>' ) <EOL> rxPre = Richards . RichardsRx ( locs , times , '<STR_LIT>' ) <EOL> survey = Richards . RichardsSurvey ( [ rxSat , rxPre ] ) <EOL> prob . pair ( survey ) <EOL> self . h0 = h <EOL> self . M = M <EOL> self . Ks = params [ '<STR_LIT>' ] <EOL> self . prob = prob <EOL> self . survey = survey <EOL> def test_Richards_getResidual_Newton ( self ) : <EOL> self . prob . doNewton = True <EOL> m = self . Ks <EOL> passed = checkDerivative ( lambda hn1 : self . prob . getResidual ( m , self . h0 , hn1 , self . prob . timeSteps [ <NUM_LIT:0> ] , self . prob . boundaryConditions ) , self . h0 , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_Richards_getResidual_Picard ( self ) : <EOL> self . prob . doNewton = False <EOL> m = self . Ks <EOL> passed = checkDerivative ( lambda hn1 : self . prob . getResidual ( m , self . h0 , hn1 , self . prob . timeSteps [ <NUM_LIT:0> ] , self . prob . boundaryConditions ) , self . h0 , plotIt = False , expectedOrder = <NUM_LIT:1> ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_Adjoint ( self ) : <EOL> v = np . random . rand ( self . survey . nD ) <EOL> z = np . random . rand ( self . M . nC ) <EOL> Hs = self . prob . fields ( self . Ks ) <EOL> vJz = v . dot ( self . prob . Jvec ( self . Ks , z , u = Hs ) ) <EOL> zJv = z . dot ( self . prob . Jtvec ( self . Ks , v , u = Hs ) ) <EOL> tol = TOL * ( <NUM_LIT:10> ** int ( np . log10 ( np . abs ( zJv ) ) ) ) <EOL> passed = np . abs ( vJz - zJv ) < tol <EOL> print '<STR_LIT>' <EOL> print '<STR_LIT>' % ( vJz , zJv , np . abs ( vJz - zJv ) , tol ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_Sensitivity ( self ) : <EOL> mTrue = self . Ks * np . ones ( self . M . nC ) <EOL> derChk = lambda m : [ self . survey . dpred ( m ) , lambda v : self . prob . Jvec ( m , v ) ] <EOL> print '<STR_LIT>' <EOL> passed = checkDerivative ( derChk , mTrue , num = <NUM_LIT:4> , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_Sensitivity_full ( self ) : <EOL> mTrue = self . Ks * np . ones ( self . M . nC ) <EOL> J = self . prob . Jfull ( mTrue ) <EOL> derChk = lambda m : [ self . survey . dpred ( m ) , J ] <EOL> print '<STR_LIT>' <EOL> passed = checkDerivative ( derChk , mTrue , num = <NUM_LIT:4> , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> class RichardsTests2D ( unittest . TestCase ) : <EOL> def setUp ( self ) : <EOL> M = Mesh . TensorMesh ( [ np . ones ( <NUM_LIT:8> ) , np . ones ( <NUM_LIT:30> ) ] ) <EOL> M . setCellGradBC ( [ '<STR_LIT>' , '<STR_LIT>' ] ) <EOL> params = Richards . Empirical . HaverkampParams ( ) . celia1990 <EOL> params [ '<STR_LIT>' ] = np . log ( params [ '<STR_LIT>' ] ) <EOL> E = Richards . Empirical . Haverkamp ( M , ** params ) <EOL> bc = np . array ( [ - <NUM_LIT> , - <NUM_LIT> ] ) <EOL> bc = np . r_ [ np . zeros ( M . nCy * <NUM_LIT:2> ) , np . ones ( M . nCx ) * bc [ <NUM_LIT:0> ] , np . ones ( M . nCx ) * bc [ <NUM_LIT:1> ] ] <EOL> h = np . zeros ( M . nC ) + bc [ <NUM_LIT:0> ] <EOL> prob = Richards . RichardsProblem ( M , E , timeSteps = [ ( <NUM_LIT> , <NUM_LIT:3> ) , ( <NUM_LIT> , <NUM_LIT:3> ) ] , boundaryConditions = bc , initialConditions = h , doNewton = False , method = '<STR_LIT>' , tolRootFinder = <NUM_LIT> , debug = False ) <EOL> prob . Solver = Solver <EOL> locs = Utils . ndgrid ( np . array ( [ <NUM_LIT:5> , <NUM_LIT> ] ) , np . array ( [ <NUM_LIT:5> , <NUM_LIT:15> , <NUM_LIT> ] ) ) <EOL> times = prob . times [ <NUM_LIT:3> : <NUM_LIT:5> ] <EOL> rxSat = Richards . RichardsRx ( locs , times , '<STR_LIT>' ) <EOL> rxPre = Richards . RichardsRx ( locs , times , '<STR_LIT>' ) <EOL> survey = Richards . RichardsSurvey ( [ rxSat , rxPre ] ) <EOL> prob . pair ( survey ) <EOL> self . h0 = h <EOL> self . M = M <EOL> self . Ks = params [ '<STR_LIT>' ] <EOL> self . prob = prob <EOL> self . survey = survey <EOL> def test_Richards_getResidual_Newton ( self ) : <EOL> self . prob . doNewton = True <EOL> m = self . Ks <EOL> passed = checkDerivative ( lambda hn1 : self . prob . getResidual ( m , self . h0 , hn1 , self . prob . timeSteps [ <NUM_LIT:0> ] , self . prob . boundaryConditions ) , self . h0 , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_Richards_getResidual_Picard ( self ) : <EOL> self . prob . doNewton = False <EOL> m = self . Ks <EOL> passed = checkDerivative ( lambda hn1 : self . prob . getResidual ( m , self . h0 , hn1 , self . prob . timeSteps [ <NUM_LIT:0> ] , self . prob . boundaryConditions ) , self . h0 , plotIt = False , expectedOrder = <NUM_LIT:1> ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_Adjoint ( self ) : <EOL> v = np . random . rand ( self . survey . nD ) <EOL> z = np . random . rand ( self . M . nC ) <EOL> Hs = self . prob . fields ( self . Ks ) <EOL> vJz = v . dot ( self . prob . Jvec ( self . Ks , z , u = Hs ) ) <EOL> zJv = z . dot ( self . prob . Jtvec ( self . Ks , v , u = Hs ) ) <EOL> tol = TOL * ( <NUM_LIT:10> ** int ( np . log10 ( np . abs ( zJv ) ) ) ) <EOL> passed = np . abs ( vJz - zJv ) < tol <EOL> print '<STR_LIT>' <EOL> print '<STR_LIT>' % ( vJz , zJv , np . abs ( vJz - zJv ) , tol ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_Sensitivity ( self ) : <EOL> mTrue = self . Ks * np . ones ( self . M . nC ) <EOL> derChk = lambda m : [ self . survey . dpred ( m ) , lambda v : self . prob . Jvec ( m , v ) ] <EOL> print '<STR_LIT>' <EOL> passed = checkDerivative ( derChk , mTrue , num = <NUM_LIT:3> , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_Sensitivity_full ( self ) : <EOL> mTrue = self . Ks * np . ones ( self . M . nC ) <EOL> J = self . prob . Jfull ( mTrue ) <EOL> derChk = lambda m : [ self . survey . dpred ( m ) , J ] <EOL> print '<STR_LIT>' <EOL> passed = checkDerivative ( derChk , mTrue , num = <NUM_LIT:4> , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> class RichardsTests3D ( unittest . TestCase ) : <EOL> def setUp ( self ) : <EOL> M = Mesh . TensorMesh ( [ np . ones ( <NUM_LIT:8> ) , np . ones ( <NUM_LIT:20> ) , np . ones ( <NUM_LIT:10> ) ] ) <EOL> M . setCellGradBC ( [ '<STR_LIT>' , '<STR_LIT>' , '<STR_LIT>' ] ) <EOL> params = Richards . Empirical . HaverkampParams ( ) . celia1990 <EOL> params [ '<STR_LIT>' ] = np . log ( params [ '<STR_LIT>' ] ) <EOL> E = Richards . Empirical . Haverkamp ( M , ** params ) <EOL> bc = np . array ( [ - <NUM_LIT> , - <NUM_LIT> ] ) <EOL> bc = np . r_ [ np . zeros ( M . nCy * M . nCz * <NUM_LIT:2> ) , np . zeros ( M . nCx * M . nCz * <NUM_LIT:2> ) , np . ones ( M . nCx * M . nCy ) * bc [ <NUM_LIT:0> ] , np . ones ( M . nCx * M . nCy ) * bc [ <NUM_LIT:1> ] ] <EOL> h = np . zeros ( M . nC ) + bc [ <NUM_LIT:0> ] <EOL> prob = Richards . RichardsProblem ( M , E , timeSteps = [ ( <NUM_LIT> , <NUM_LIT:3> ) , ( <NUM_LIT> , <NUM_LIT:3> ) ] , boundaryConditions = bc , initialConditions = h , doNewton = False , method = '<STR_LIT>' , tolRootFinder = <NUM_LIT> , debug = False ) <EOL> prob . Solver = Solver <EOL> locs = Utils . ndgrid ( np . r_ [ <NUM_LIT:5> , <NUM_LIT> ] , np . r_ [ <NUM_LIT:5> , <NUM_LIT> ] , np . r_ [ <NUM_LIT:6> , <NUM_LIT> ] ) <EOL> times = prob . times [ <NUM_LIT:3> : <NUM_LIT:5> ] <EOL> rxSat = Richards . RichardsRx ( locs , times , '<STR_LIT>' ) <EOL> rxPre = Richards . RichardsRx ( locs , times , '<STR_LIT>' ) <EOL> survey = Richards . RichardsSurvey ( [ rxSat , rxPre ] ) <EOL> prob . pair ( survey ) <EOL> self . h0 = h <EOL> self . M = M <EOL> self . Ks = params [ '<STR_LIT>' ] <EOL> self . prob = prob <EOL> self . survey = survey <EOL> def test_Richards_getResidual_Newton ( self ) : <EOL> self . prob . doNewton = True <EOL> m = self . Ks <EOL> passed = checkDerivative ( lambda hn1 : self . prob . getResidual ( m , self . h0 , hn1 , self . prob . timeSteps [ <NUM_LIT:0> ] , self . prob . boundaryConditions ) , self . h0 , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_Richards_getResidual_Picard ( self ) : <EOL> self . prob . doNewton = False <EOL> m = self . Ks <EOL> passed = checkDerivative ( lambda hn1 : self . prob . getResidual ( m , self . h0 , hn1 , self . prob . timeSteps [ <NUM_LIT:0> ] , self . prob . boundaryConditions ) , self . h0 , plotIt = False , expectedOrder = <NUM_LIT:1> ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_Adjoint ( self ) : <EOL> v = np . random . rand ( self . survey . nD ) <EOL> z = np . random . rand ( self . M . nC ) <EOL> Hs = self . prob . fields ( self . Ks ) <EOL> vJz = v . dot ( self . prob . Jvec ( self . Ks , z , u = Hs ) ) <EOL> zJv = z . dot ( self . prob . Jtvec ( self . Ks , v , u = Hs ) ) <EOL> tol = TOL * ( <NUM_LIT:10> ** int ( np . log10 ( np . abs ( zJv ) ) ) ) <EOL> passed = np . abs ( vJz - zJv ) < tol <EOL> print '<STR_LIT>' <EOL> print '<STR_LIT>' % ( vJz , zJv , np . abs ( vJz - zJv ) , tol ) <EOL> self . assertTrue ( passed , True ) <EOL> def test_Sensitivity ( self ) : <EOL> mTrue = self . Ks * np . ones ( self . M . nC ) <EOL> derChk = lambda m : [ self . survey . dpred ( m ) , lambda v : self . prob . Jvec ( m , v ) ] <EOL> print '<STR_LIT>' <EOL> passed = checkDerivative ( derChk , mTrue , num = <NUM_LIT:4> , plotIt = False ) <EOL> self . assertTrue ( passed , True ) <EOL> if __name__ == '<STR_LIT:__main__>' : <EOL> unittest . main ( ) </s>
