From 2831b66ea9a8b4f5b7f7b5f1e451e0efd04187c9 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 20 Apr 2014 11:33:59 -0700 Subject: [PATCH] Variable Time Stepping --- simpegFLOW/Richards/RichardsProblem.py | 36 +++++++------------------- simpegFLOW/Tests/test_Richards.py | 8 +++--- 2 files changed, 14 insertions(+), 30 deletions(-) diff --git a/simpegFLOW/Richards/RichardsProblem.py b/simpegFLOW/Richards/RichardsProblem.py index 58a6f0fb..43485aeb 100644 --- a/simpegFLOW/Richards/RichardsProblem.py +++ b/simpegFLOW/Richards/RichardsProblem.py @@ -59,10 +59,9 @@ class RichardsSurvey(Survey.BaseSurvey): return self.P*dT -class RichardsProblem(Problem.BaseProblem): +class RichardsProblem(Problem.BaseTimeProblem): """docstring for RichardsProblem""" - timeEnd = None boundaryConditions = None initialConditions = None @@ -73,20 +72,7 @@ class RichardsProblem(Problem.BaseProblem): solverOpts = {} def __init__(self, mesh, mapping=None, **kwargs): - Problem.BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs) - - @property - def timeStep(self): - """The time between steps.""" - return getattr(self, '_timeStep', None) - @timeStep.setter - def timeStep(self, value): - self._timeStep = float(value) # Because integers suck. - - @property - def numIts(self): - """The number of iterations in the time domain problem.""" - return int(self.timeEnd/self.timeStep) + Problem.BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs) @property def method(self): @@ -112,13 +98,13 @@ class RichardsProblem(Problem.BaseProblem): return self._rootFinder def fields(self, m): - u = range(self.numIts+1) + u = range(self.nT+1) u[0] = self.initialConditions - for ii in range(self.numIts): - u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, return_g=return_g), u[ii]) + for ii, dt in enumerate(self.timeSteps): + u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, dt, return_g=return_g), u[ii]) return u - def diagsJacobian(self, m, hn, hn1): + def diagsJacobian(self, m, hn, hn1, dt): DIV = self.mesh.faceDiv GRAD = self.mesh.cellGrad @@ -132,7 +118,6 @@ class RichardsProblem(Problem.BaseProblem): Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.vnF[0]+self.mesh.vnF[1]), self.mesh.faceDivz),format='csr') bc = self.boundaryConditions - dt = self.timeStep dT = self.mapping.thetaDerivU(hn, m) dT1 = self.mapping.thetaDerivU(hn1, m) @@ -170,7 +155,7 @@ class RichardsProblem(Problem.BaseProblem): return Asub, Adiag, B - def getResidual(self, m, hn, h, return_g=True): + def getResidual(self, m, hn, h, dt, return_g=True): """ Where h is the proposed value for the next time iterate (h_{n+1}) """ @@ -186,7 +171,6 @@ class RichardsProblem(Problem.BaseProblem): Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.vnF[0]+self.mesh.vnF[1]), self.mesh.faceDivz),format='csr') bc = self.boundaryConditions - dt = self.timeStep T = self.mapping.theta(h, m) dT = self.mapping.thetaDerivU(h, m) @@ -238,7 +222,7 @@ class RichardsProblem(Problem.BaseProblem): JvC = range(len(u)-1) # Cell to hold each row of the long vector. # This is done via forward substitution. - temp, Adiag, B = self.diagsJacobian(m, u[0],u[1]) + temp, Adiag, B = self.diagsJacobian(m, u[0], u[1], self.timeSteps[0]) Adiaginv = self.Solver(Adiag, **self.solverOpts) JvC[0] = Adiaginv.solve(B*v) @@ -246,7 +230,7 @@ class RichardsProblem(Problem.BaseProblem): # 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]) + 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]) @@ -264,7 +248,7 @@ class RichardsProblem(Problem.BaseProblem): minus = 0 BJtv = 0 for ii in range(len(u)-1,0,-1): - Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii]) + 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]) AdiaginvT = self.Solver(Adiag.T, **self.solverOpts) diff --git a/simpegFLOW/Tests/test_Richards.py b/simpegFLOW/Tests/test_Richards.py index 10dcf5e7..5387a18e 100644 --- a/simpegFLOW/Tests/test_Richards.py +++ b/simpegFLOW/Tests/test_Richards.py @@ -138,12 +138,12 @@ class RichardsTests1D(unittest.TestCase): bc = np.array([-61.5,-20.7]) h = np.zeros(M.nC) + bc[0] - prob = Richards.RichardsProblem(M, mapping=E, timeStep=60, timeEnd=180, + prob = Richards.RichardsProblem(M, mapping=E, timeSteps=[(40,3),(60,3)], 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.numIts),q) + P = sp.kron(sp.identity(prob.nT),q) survey = Richards.RichardsSurvey(P=P) prob.pair(survey) @@ -157,13 +157,13 @@ 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.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.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):