From 96e129aae240cf8939ce129fe77519c351636e21 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 20 Apr 2014 13:12:14 -0700 Subject: [PATCH] Time projections --- simpegFLOW/Richards/RichardsProblem.py | 83 ++++++++++------ simpegFLOW/Tests/test_Richards.py | 128 +++---------------------- 2 files changed, 67 insertions(+), 144 deletions(-) diff --git a/simpegFLOW/Richards/RichardsProblem.py b/simpegFLOW/Richards/RichardsProblem.py index 43485aeb..30715c55 100644 --- a/simpegFLOW/Richards/RichardsProblem.py +++ b/simpegFLOW/Richards/RichardsProblem.py @@ -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. diff --git a/simpegFLOW/Tests/test_Richards.py b/simpegFLOW/Tests/test_Richards.py index 5387a18e..bcda57bf 100644 --- a/simpegFLOW/Tests/test_Richards.py +++ b/simpegFLOW/Tests/test_Richards.py @@ -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):