wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-08 11:42:49 +08:00
Code Issues Packages Projects Releases Wiki Activity

update the tests to test the 3D version of the equations.

Browse Source
This commit is contained in:
Rowan Cockett
2015-01-29 14:58:01 -08:00
parent b2d57e8892
commit 5c675a7992
1 changed files with 80 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegFLOW/Tests/test_Richards.py
+80
View File
@@ -3,6 +3,12 @@ from SimPEG import *
from SimPEG.Tests.TestUtils import OrderTest, checkDerivative
from scipy.sparse.linalg import dsolve
from simpegFLOW import Richards
try:
from pymatsolver import MumpsSolver
Solver = MumpsSolver
except Exception, e:
pass
TOL = 1E-8
@@ -76,6 +82,7 @@ class RichardsTests1D(unittest.TestCase):
prob = Richards.RichardsProblem(M, mapping=E, timeSteps=[(40,3),(60,3)],
boundaryConditions=bc, initialConditions=h,
doNewton=False, method='mixed')
prob.Solver = Solver
locs = np.r_[5.,10,15]
times = prob.times[3:5]
@@ -147,6 +154,7 @@ class RichardsTests2D(unittest.TestCase):
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')
prob.Solver = Solver
locs = Utils.ndgrid(np.array([5,7.]),np.array([5,15,25.]))
times = prob.times[3:5]
@@ -202,5 +210,77 @@ class RichardsTests2D(unittest.TestCase):
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')
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,u=Hs))
zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=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()