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

Time projections

Browse Source
This commit is contained in:
rowanc1
2014-04-20 13:12:14 -07:00
parent 2831b66ea9
commit 96e129aae2
2 changed files with 67 additions and 144 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegFLOW/Richards/RichardsProblem.py
+53 -30
View File
@@ -1,22 +1,44 @@
from SimPEG import *
from Empirical import RichardsMap
class RichardsRx(Survey.BaseTimeRx):
"""Richards Receiver Object"""
knownRxTypes = ['saturation','pressureHead']
def projectFields(self, u, m, mapping, mesh, timeMesh):
if self.rxType == 'saturation':
u = mapping.theta(u, m)
return self.getP(mesh, timeMesh) * u
def projectFieldsDeriv(self, u, m, mapping, mesh, timeMesh):
P = self.getP(mesh, timeMesh)
if self.rxType == 'pressureHead':
return P
elif self.rxType == 'saturation':
#TODO: if m is a parameter in the theta
# distribution, we may need to do
# some more chain rule here.
dT = mapping.thetaDerivU(u, m)
return P*dT
class RichardsSurvey(Survey.BaseSurvey):
"""docstring for RichardsSurvey"""
P = None
rxList = None
def __init__(self, **kwargs):
def __init__(self, rxList, **kwargs):
self.rxList = rxList
Survey.BaseSurvey.__init__(self, **kwargs)
@property
def dataType(self):
"""Choose how your data is collected, must be 'saturation' or 'pressureHead'."""
return getattr(self, '_dataType', 'pressureHead')
@dataType.setter
def dataType(self, value):
assert value in ['saturation','pressureHead'], "dataType must be 'saturation' or 'pressureHead'."
self._dataType = value
def nD(self):
return np.array([rx.nD for rx in self.rxList]).sum()
@Utils.count
@Utils.requires('prob')
@@ -27,7 +49,7 @@ class RichardsSurvey(Survey.BaseSurvey):
instead of recalculating the fields (which may be expensive!).
.. math::
d_\\text{pred} = P(u(m))
d_\\text{pred} = P(u(m), m)
Where P is a projection of the fields onto the data space.
"""
@@ -37,27 +59,31 @@ class RichardsSurvey(Survey.BaseSurvey):
@Utils.requires('prob')
def projectFields(self, U, m):
u = np.concatenate(U[1:])
u = np.concatenate(U)
if self.dataType == 'saturation':
u = self.prob.model.theta(u, m)
return self.P*u
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.projectFields(u, m,
self.prob.mapping,
self.prob.mesh,
self.prob.timeMesh)
return np.concatenate(Ds)
@Utils.requires('prob')
def projectFieldsDeriv(self, U, m):
"""The Derivative with respect to the fields."""
u = np.concatenate(U[1:])
u = np.concatenate(U)
if self.dataType == 'pressureHead':
return self.P
elif self.dataType == 'saturation':
#TODO: if m is a parameter in the theta
# distribution, we may need to do
# some more chain rule here.
dT = self.mapping.thetaDerivU(u, m)
return self.P*dT
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.projectFieldsDeriv(u, m,
self.prob.mapping,
self.prob.mesh,
self.prob.timeMesh)
return sp.vstack(Ds)
class RichardsProblem(Problem.BaseTimeProblem):
"""docstring for RichardsProblem"""
@@ -197,7 +223,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
def Jfull(self, m, u=None):
if u is None:
u = self.field(m)
u = self.fields(m)
nn = len(u)-1
Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
@@ -217,7 +243,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
def Jvec(self, m, v, u=None):
if u is None:
u = self.field(m)
u = self.fields(m)
JvC = range(len(u)-1) # Cell to hold each row of the long vector.
@@ -226,16 +252,13 @@ class RichardsProblem(Problem.BaseTimeProblem):
Adiaginv = self.Solver(Adiag, **self.solverOpts)
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(u)-1):
Asub, Adiag, B = self.diagsJacobian(m, u[ii], u[ii+1], self.timeSteps[ii])
Adiaginv = self.Solver(Adiag, **self.solverOpts)
JvC[ii] = Adiaginv.solve(B*v - Asub*JvC[ii-1])
P = self.survey.projectFieldsDeriv(u, m)
return P * np.concatenate(JvC)
return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC)
def Jtvec(self, m, v, u=None):
if u is None:
@@ -250,7 +273,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
for ii in range(len(u)-1,0,-1):
Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii], self.timeSteps[ii-1])
#select the correct part of v
vpart = range((ii-1)*Adiag.shape[0], (ii)*Adiag.shape[0])
vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0])
AdiaginvT = self.Solver(Adiag.T, **self.solverOpts)
JTvC = AdiaginvT.solve(PTv[vpart] - minus)
minus = Asub.T*JTvC # this is now the super diagonal.
simpegFLOW/Tests/test_Richards.py
+14 -114
View File
@@ -58,72 +58,6 @@ class TestModels(unittest.TestCase):
passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
self.assertTrue(passed,True)
# def test_Haverkamp_hydraulicConductivity(self):
# print 'Haverkamp_hydraulicConductivity'
# hav = Richards.Haverkamp()
# def wrapper(x):
# return hav.hydraulicConductivity(x), hav.hydraulicConductivityDeriv(x)
# passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
# self.assertTrue(passed,True)
# def test_Haverkamp_hydraulicConductivity_FullKs(self):
# print 'Haverkamp_hydraulicConductivity_FullKs'
# n = 50
# hav = Richards.Haverkamp(Ks=np.random.rand(n))
# def wrapper(x):
# return hav.hydraulicConductivity(x), hav.hydraulicConductivityDeriv(x)
# passed = checkDerivative(wrapper, np.random.randn(n), plotIt=False)
# self.assertTrue(passed,True)
# def test_VanGenuchten_moistureContent(self):
# print 'VanGenuchten_moistureContent'
# vanG = Richards.VanGenuchten()
# def wrapper(x):
# return vanG.moistureContent(x), vanG.moistureContentDeriv(x)
# passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
# self.assertTrue(passed,True)
# def test_VanGenuchten_hydraulicConductivity(self):
# print 'VanGenuchten_hydraulicConductivity'
# hav = Richards.VanGenuchten()
# def wrapper(x):
# return hav.hydraulicConductivity(x), hav.hydraulicConductivityDeriv(x)
# passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
# self.assertTrue(passed,True)
# def test_VanGenuchten_hydraulicConductivity_FullKs(self):
# print 'VanGenuchten_hydraulicConductivity_FullKs'
# n = 50
# hav = Richards.VanGenuchten(Ks=np.random.rand(n))
# def wrapper(x):
# return hav.hydraulicConductivity(x), hav.hydraulicConductivityDeriv(x)
# passed = checkDerivative(wrapper, np.random.randn(n), plotIt=False)
# self.assertTrue(passed,True)
# def test_Haverkamp_moistureContent(self):
# print 'Haverkamp_moistureContent'
# hav = Richards.Haverkamp()
# def wrapper(x):
# return hav.moistureContent(x), hav.moistureContentDeriv(x)
# passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
# self.assertTrue(passed,True)
# def test_Haverkamp_hydraulicConductivity(self):
# print 'Haverkamp_hydraulicConductivity'
# hav = Richards.Haverkamp()
# def wrapper(x):
# return hav.hydraulicConductivity(x), hav.hydraulicConductivityDeriv(x)
# passed = checkDerivative(wrapper, np.random.randn(50), plotIt=False)
# self.assertTrue(passed,True)
# def test_Haverkamp_hydraulicConductivity_FullKs(self):
# print 'Haverkamp_hydraulicConductivity_FullKs'
# n = 50
# hav = Richards.Haverkamp(Ks=np.random.rand(n))
# def wrapper(x):
# return hav.hydraulicConductivity(x), hav.hydraulicConductivityDeriv(x)
# passed = checkDerivative(wrapper, np.random.randn(n), plotIt=False)
# self.assertTrue(passed,True)
class RichardsTests1D(unittest.TestCase):
@@ -142,9 +76,11 @@ class RichardsTests1D(unittest.TestCase):
boundaryConditions=bc, initialConditions=h,
doNewton=False, method='mixed')
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.nT),q)
survey = Richards.RichardsSurvey(P=P)
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)
@@ -157,18 +93,17 @@ class RichardsTests1D(unittest.TestCase):
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.h0, plotIt=False)
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0]), 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.h0, plotIt=False, expectedOrder=1)
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0]), self.h0, plotIt=False, expectedOrder=1)
self.assertTrue(passed,True)
def test_Adjoint_PressureHead(self):
self.prob.dataType = 'pressureHead'
v = np.random.rand(self.survey.P.shape[0])
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))
@@ -179,48 +114,13 @@ class RichardsTests1D(unittest.TestCase):
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'
v = np.random.rand(self.survey.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)
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_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
survey = self.prob.createSyntheticSurvey(mTrue, std=stdev, P=self.survey.P)
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
reg = Regularization.Tikhonov(self.M)
obj = ObjFunction.BaseObjFunction(survey, 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
survey = self.prob.createSyntheticSurvey(mTrue, std=stdev, P=self.survey.P)
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
reg = Regularization.Tikhonov(self.M)
obj = ObjFunction.BaseObjFunction(survey, 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)
# class RichardsTests2D(object):