From 715bed21ce24bdd4453dc5f43dc591cc47c3f622 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 25 Feb 2014 11:35:13 -0800 Subject: [PATCH] Initial commit of richards equation code. Forward working. Inverse untested. --- docs/api_Richards.rst | 8 +- simpegFLOW/Richards/BaseRichards.py | 256 +++++++++++++++++++++++ simpegFLOW/Richards/RichardsProblem.py | 278 +++++++++++++++++++++++++ simpegFLOW/Richards/__init__.py | 3 +- simpegFLOW/Tests/test_Richards.py | 213 ++++++++++++++++++- 5 files changed, 752 insertions(+), 6 deletions(-) create mode 100644 simpegFLOW/Richards/BaseRichards.py create mode 100644 simpegFLOW/Richards/RichardsProblem.py diff --git a/docs/api_Richards.rst b/docs/api_Richards.rst index 084ae86e..9bb787e3 100644 --- a/docs/api_Richards.rst +++ b/docs/api_Richards.rst @@ -1,5 +1,9 @@ .. _api_Richards: + +Richards Equation +***************** + There are two different forms of Richards equation that differ on how they deal with the non-linearity in the time-stepping term. @@ -8,7 +12,7 @@ The most fundamental form, referred to as the .. math:: - \\frac{\partial \\theta(\psi)}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0 + \frac{\partial \theta(\psi)}{\partial t} - \nabla \cdot k(\psi) \nabla \psi - \frac{\partial k(\psi)}{\partial z} = 0 \quad \psi \in \Omega where theta is water content, and psi is pressure head. @@ -21,7 +25,7 @@ equation can be written in the continuous form as: .. math:: - \\frac{\partial \\theta}{\partial \psi}\\frac{\partial \psi}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0 + \frac{\partial \theta}{\partial \psi}\frac{\partial \psi}{\partial t} - \nabla \cdot k(\psi) \nabla \psi - \frac{\partial k(\psi)}{\partial z} = 0 \quad \psi \in \Omega However, it can be shown that this does not conserve mass in the discrete formulation. diff --git a/simpegFLOW/Richards/BaseRichards.py b/simpegFLOW/Richards/BaseRichards.py new file mode 100644 index 00000000..04359512 --- /dev/null +++ b/simpegFLOW/Richards/BaseRichards.py @@ -0,0 +1,256 @@ +from SimPEG import Model, Utils, np + + +class RichardsModel(object): + """docstring for RichardsModel""" + + mesh = None #: SimPEG mesh + + @property + def thetaModel(self): + """Model for moisture content""" + return self._thetaModel + + @property + def kModel(self): + """Model for hydraulic conductivity""" + return self._kModel + + def __init__(self, mesh, thetaModel, kModel): + self.mesh = mesh + assert isinstance(thetaModel, Model.BaseNonLinearModel) + assert isinstance(kModel, Model.BaseNonLinearModel) + + self._thetaModel = thetaModel + self._kModel = kModel + + def theta(self, u, m): + return self.thetaModel.transform(u, m) + + def thetaDerivM(self, u, m): + return self.thetaModel.transformDerivM(u, m) + + def thetaDerivU(self, u, m): + return self.thetaModel.transformDerivU(u, m) + + def k(self, u, m): + return self.kModel.transform(u, m) + + def kDerivM(self, u, m): + return self.kModel.transformDerivM(u, m) + + def kDerivU(self, u, m): + return self.kModel.transformDerivU(u, m) + + +class BaseHaverkamp_theta(Model.BaseNonLinearModel): + + theta_s = 0.430 + theta_r = 0.078 + alpha = 0.036 + beta = 3.960 + + def __init__(self, mesh, **kwargs): + Model.BaseNonLinearModel.__init__(self, mesh) + Utils.setKwargs(self, **kwargs) + + def setModel(self, m): + self._currentModel = m + + def transform(self, u, m): + self.setModel(m) + f = (self.alpha*(self.theta_s - self.theta_r )/ + (self.alpha + abs(u)**self.beta) + self.theta_r) + f[u >= 0] = self.theta_s + return f + + def transformDerivM(self, u, m): + self.setModel(m) + + def transformDerivU(self, u, m): + self.setModel(m) + g = (self.alpha*((self.theta_s - self.theta_r)/ + (self.alpha + abs(u)**self.beta)**2) + *(-self.beta*abs(u)**(self.beta-1)*np.sign(u))) + g[u >= 0] = 0 + g = Utils.sdiag(g) + return g + + +class BaseHaverkamp_k(Model.BaseNonLinearModel): + + A = 1.175e+06 + gamma = 4.74 + Ks = np.log(24.96) + + def __init__(self, mesh, **kwargs): + Model.BaseNonLinearModel.__init__(self, mesh) + Utils.setKwargs(self, **kwargs) + + def setModel(self, m): + self._currentModel = m + #TODO: Fix me! + self.Ks = m + + def transform(self, u, m): + self.setModel(m) + f = np.exp(self.Ks)*self.A/(self.A+abs(u)**self.gamma) + if type(self.Ks) is np.ndarray and self.Ks.size > 1: + f[u >= 0] = np.exp(self.Ks[u >= 0]) + else: + f[u >= 0] = np.exp(self.Ks) + return f + + def transformDerivM(self, u, m): + self.setModel(m) + #A + # dA = np.exp(self.Ks)/(self.A+abs(u)**self.gamma) - np.exp(self.Ks)*self.A/(self.A+abs(u)**self.gamma)**2 + #gamma + # dgamma = -(self.A*np.exp(self.Ks)*np.log(abs(u))*abs(u)**self.gamma)/(self.A + abs(u)**self.gamma)**2 + + # This assumes that the the model is Ks + return Utils.sdiag(self.transform(u, m)) + + def transformDerivU(self, u, m): + self.setModel(m) + g = -(np.exp(self.Ks)*self.A*self.gamma*abs(u)**(self.gamma-1)*np.sign(u))/((self.A+abs(u)**self.gamma)**2) + g[u >= 0] = 0 + g = Utils.sdiag(g) + return g + + + +# class Haverkamp(object): +# """docstring for Haverkamp""" + +# empiricalModelName = "VanGenuchten" + +# theta_s = 0.430 +# theta_r = 0.078 +# alpha = 0.036 +# beta = 3.960 +# A = 1.175e+06 +# gamma = 4.74 +# Ks = np.log(24.96) + +# def __init__(self, **kwargs): +# Utils.setKwargs(self, **kwargs) + +# def setModel(self, m): +# self.Ks = m + +# def moistureContent(self, h): +# f = (self.alpha*(self.theta_s - self.theta_r )/ +# (self.alpha + abs(h)**self.beta) + self.theta_r) +# f[h > 0] = self.theta_s +# return f + +# def moistureContentDeriv(self, h): +# g = (self.alpha*((self.theta_s - self.theta_r)/ +# (self.alpha + abs(h)**self.beta)**2) +# *(-self.beta*abs(h)**(self.beta-1)*np.sign(h))); +# g[h >= 0] = 0 +# g = Utils.sdiag(g) +# return g + +# def hydraulicConductivity(self, h): +# f = np.exp(self.Ks)*self.A/(self.A+abs(h)**self.gamma) +# if type(self.Ks) is np.ndarray and self.Ks.size > 1: +# f[h >= 0] = np.exp(self.Ks[h >= 0]) +# else: +# f[h >= 0] = np.exp(self.Ks) +# return f + +# def hydraulicConductivityModelDeriv(self, h): +# #A +# # dA = np.exp(self.Ks)/(self.A+abs(h)**self.gamma) - np.exp(self.Ks)*self.A/(self.A+abs(h)**self.gamma)**2; +# #gamma +# # dgamma = -(self.A*np.exp(self.Ks)*np.log(abs(h))*abs(h)**self.gamma)/(self.A + abs(h)**self.gamma)**2; +# return Utils.sdiag(self.hydraulicConductivity(h)) # This assumes that the the model is Ks + +# def hydraulicConductivityDeriv(self, h): +# g = -(np.exp(self.Ks)*self.A*self.gamma*abs(h)**(self.gamma-1)*np.sign(h))/((self.A+abs(h)**self.gamma)**2) +# g[h >= 0] = 0 +# g = Utils.sdiag(g) +# return g + + +# class VanGenuchten(object): +# """ + +# .. math:: + +# \\theta(h) = \\frac{\\alpha (\\theta_s - \\theta_r)}{\\alpha + |h|^\\beta} + \\theta_r + +# Where parameters alpha, beta, gamma, A are constants in the media; +# theta_r and theta_s are the residual and saturated moisture +# contents; and K_s is the saturated hydraulic conductivity. + +# Celia1990 + +# """ + +# empiricalModelName = "VanGenuchten" + +# theta_s = 0.430 +# theta_r = 0.078 +# alpha = 0.036 +# n = 1.560 +# beta = 3.960 +# I = 0.500 +# Ks = np.log(24.96) + +# def __init__(self, **kwargs): +# Utils.setKwargs(self, **kwargs) + +# def setModel(self, m): +# self.Ks = m + +# def moistureContent(self, h): +# m = 1 - 1.0/self.n; +# f = (( self.theta_s - self.theta_r )/ +# ((1+abs(self.alpha*h)**self.n)**m) + self.theta_r) +# f[h > 0] = self.theta_s +# return f + +# def moistureContentDeriv(self, h): +# g = -self.alpha*self.n*abs(self.alpha*h)**(self.n - 1)*np.sign(self.alpha*h)*(1./self.n - 1)*(self.theta_r - self.theta_s)*(abs(self.alpha*h)**self.n + 1)**(1./self.n - 2) +# g[h > 0] = 0 +# g = Utils.sdiag(g) +# return g + +# def hydraulicConductivity(self, h): +# alpha = self.alpha +# I = self.I +# n = self.n +# Ks = self.Ks +# m = 1.0 - 1.0/n + +# theta_e = 1.0/((1.0+abs(alpha*h)**n)**m) +# f = np.exp(Ks)*theta_e**I* ( ( 1.0 - ( 1.0 - theta_e**(1.0/m) )**m )**2 ) +# if type(self.Ks) is np.ndarray and self.Ks.size > 1: +# f[h >= 0] = np.exp(self.Ks[h >= 0]) +# else: +# f[h >= 0] = np.exp(self.Ks) +# return f + +# def hydraulicConductivityModelDeriv(self, h): +# #alpha +# # dA = I*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1.0/n - 1)*((abs(alpha*h)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*h)**n + 1)**(1.0/n - 2) - (2*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1.0/n - 1)*((abs(alpha*h)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*h)**n + 1)**(1.0/n - 2))/(((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)*(1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n)); +# #n +# # dn = 2*np.exp(Ks)*((np.log(1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))*(1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n))/n**2 + ((1.0/n - 1)*(((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1.0/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1.0/n - 1)*(abs(alpha*h)**n + 1)**(1.0/n - 2))/((1.0/n - 1)*((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)) - np.log((abs(alpha*h)**n + 1)**(1.0/n - 1))/(n**2*(1.0/n - 1)**2*((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))))/(1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n))*((abs(alpha*h)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1) - I*np.exp(Ks)*((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1.0/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1.0/n - 1)*(abs(alpha*h)**n + 1)**(1.0/n - 2))*((abs(alpha*h)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2; +# #I +# # dI = np.exp(Ks)*np.log((abs(alpha*h)**n + 1)**(1.0/n - 1))*((abs(alpha*h)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2; +# return Utils.sdiag(self.hydraulicConductivity(h)) # This assumes that the the model is Ks + +# def hydraulicConductivityDeriv(self, h): +# alpha = self.alpha +# I = self.I +# n = self.n +# Ks = self.Ks +# m = 1.0 - 1.0/n + +# g = I*alpha*n*np.exp(Ks)*abs(alpha*h)**(n - 1.0)*np.sign(alpha*h)*(1.0/n - 1.0)*((abs(alpha*h)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*h)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1.0/n - 1)*((abs(alpha*h)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*h)**n + 1)**(1.0/n - 2))/(((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)*(1 - 1.0/((abs(alpha*h)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n)) +# g[h >= 0] = 0 +# g = Utils.sdiag(g) +# return g diff --git a/simpegFLOW/Richards/RichardsProblem.py b/simpegFLOW/Richards/RichardsProblem.py new file mode 100644 index 00000000..870d857b --- /dev/null +++ b/simpegFLOW/Richards/RichardsProblem.py @@ -0,0 +1,278 @@ +from SimPEG import * +from BaseRichards import RichardsModel + +class RichardsData(Data.BaseData): + """docstring for RichardsData""" + + P = None + + def __init__(self, **kwargs): + Data.BaseData.__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 projectFields(self, u): + u = np.concatenate(u[1:]) + if self.dataType == 'saturation': + #TODO: Fix this: + u = self.prob.model.theta(MODEL, u) + return self.P*u + + +class RichardsProblem(Problem.BaseProblem): + """docstring for RichardsProblem""" + + timeEnd = None + boundaryConditions = None + initialConditions = None + + dataPair = RichardsData + modelPair = RichardsModel + + def __init__(self, mesh, model, **kwargs): + self.doNewton = False # This also sets the rootFinder algorithm. + Problem.BaseProblem.__init__(self, mesh, model, **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) + + @property + def method(self): + """Method must be either 'mixed' or 'head'. See notes in Celia et al., 1990.""" + return getattr(self, '_method', 'mixed') + @method.setter + def method(self, value): + assert value in ['mixed','head'], "method must be 'mixed' or 'head'." + self._method = value + + @property + def doNewton(self): + """Do a Newton iteration. If False, a Picard iteration will be completed.""" + return self._doNewton + @doNewton.setter + def doNewton(self, value): + value = bool(value) + self.rootFinder = Optimization.NewtonRoot(doLS=value) + self._doNewton = value + + def fields(self, m): + Hs = range(self.numIts+1) + Hs[0] = self.initialConditions + for ii in range(self.numIts): + Hs[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, Hs[ii], hn1m, return_g=return_g), Hs[ii]) + return Hs + + def diagsJacobian(self, m, hn, hn1): + + DIV = self.mesh.faceDiv + GRAD = self.mesh.cellGrad + BC = self.mesh.cellGradBC + AV = self.mesh.aveCC2F + if self.mesh.dim == 1: + Dz = self.mesh.faceDivx + elif self.mesh.dim == 2: + Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.nFv[0]), self.mesh.faceDivy),format='csr') + elif self.mesh.dim == 3: + Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.nFv[0]+self.mesh.nFv[1]), self.mesh.faceDivz),format='csr') + + bc = self.boundaryConditions + dt = self.timeStep + + dT = self.model.thetaDerivU(hn, m) + dT1 = self.model.thetaDerivU(hn1, m) + K1 = self.model.k(hn1, m) + dK1 = self.model.kDerivU(hn1, m) + dKa1 = self.model.kDerivM(hn1, m) + + # Compute part of the derivative of: + # + # DIV*diag(GRAD*hn1+BC*bc)*(AV*(1.0/K))^-1 + + DdiagGh1 = DIV*Utils.sdiag(GRAD*hn1+BC*bc) + diagAVk2_AVdiagK2 = Utils.sdiag((AV*(1./K1))**(-2)) * AV*Utils.sdiag(K1**(-2)) + + # The matrix that we are computing has the form: + # + # - - - - - - + # | Adiag | | h1 | | b1 | + # | Asub Adiag | | h2 | | b2 | + # | Asub Adiag | | h3 | = | b3 | + # | ... ... | | .. | | .. | + # | Asub Adiag | | hn | | bn | + # - - - - - - + + Asub = (-1.0/dt)*dT + + Adiag = ( + (1.0/dt)*dT1 + -DdiagGh1*diagAVk2_AVdiagK2*dK1 + -DIV*Utils.sdiag(1./(AV*(1./K1)))*GRAD + -Dz*diagAVk2_AVdiagK2*dK1 + ) + + B = DdiagGh1*diagAVk2_AVdiagK2*dKa1 + Dz*diagAVk2_AVdiagK2*dKa1 + + return Asub, Adiag, B + + def getResidual(self, m, hn, h, return_g=True): + """ + Where h is the proposed value for the next time iterate (h_{n+1}) + """ + DIV = self.mesh.faceDiv + GRAD = self.mesh.cellGrad + BC = self.mesh.cellGradBC + AV = self.mesh.aveCC2F + if self.mesh.dim == 1: + Dz = self.mesh.faceDivx + elif self.mesh.dim == 2: + Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.nFv[0]), self.mesh.faceDivy),format='csr') + elif self.mesh.dim == 3: + Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.nFv[0]+self.mesh.nFv[1]), self.mesh.faceDivz),format='csr') + + bc = self.boundaryConditions + dt = self.timeStep + + T = self.model.theta(h, m) + dT = self.model.thetaDerivU(h, m) + Tn = self.model.theta(hn, m) + K = self.model.k(h, m) + dK = self.model.kDerivU(h, m) + + aveK = 1./(AV*(1./K)); + + RHS = DIV*Utils.sdiag(aveK)*(GRAD*h+BC*bc) + Dz*aveK + if self.method == 'mixed': + r = (T-Tn)/dt - RHS + elif self.method == 'head': + r = dT*(h - hn)/dt - RHS + + if not return_g: return r + + J = dT/dt - DIV*Utils.sdiag(aveK)*GRAD + if self.doNewton: + DDharmAve = Utils.sdiag(aveK**2)*AV*Utils.sdiag(K**(-2)) * dK + J = J - DIV*Utils.sdiag(GRAD*h + BC*bc)*DDharmAve - Dz*DDharmAve + + return r, J + + def fullJ(self, m, u=None): + if u is None: + u = self.field(m) + Hs = u + nn = len(Hs)-1 + Asubs, Adiags, Bs = range(nn), range(nn), range(nn) + for ii in range(nn): + Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, Hs[ii],Hs[ii+1]) + Ad = sp.block_diag(Adiags) + zRight = Utils.spzeros((len(Asubs)-1)*Asubs[0].shape[0],Adiags[0].shape[1]) + zTop = Utils.spzeros(Adiags[0].shape[0], len(Adiags)*Adiags[0].shape[1]) + As = sp.vstack((zTop,sp.hstack((sp.block_diag(Asubs[1:]),zRight)))) + A = As + Ad + B = np.array(sp.vstack(Bs).todense()) + + Ainv = Solver(A) + J = Ainv.solve(B) + return J + + + def Jvec(self, m, v, u=None): + if u is None: + u = self.field(m) + Hs = u + JvC = range(len(Hs)-1) # Cell to hold each row of the long vector. + + # This is done via forward substitution. + temp, Adiag, B = self.diagsJacobian(m, Hs[0],Hs[1]) + Adiaginv = Solver(Adiag) + 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(Hs)-1): + Asub, Adiag, B = self.diagsJacobian(m, Hs[ii],Hs[ii+1]) + Adiaginv = Solver(Adiag) + JvC[ii] = Adiaginv.solve(B*v - Asub*JvC[ii-1]) + + if self.dataType == 'pressureHead': + Jv = self.P*np.concatenate(JvC) + elif self.dataType == 'saturation': + dT = self.model.thetaDerivU(np.concatenate(Hs[1:]), m) + Jv = self.P*dT*np.concatenate(JvC) + + return Jv + + def Jtvec(self, m, v, u=None): + if u is None: + u = self.field(m) + Hs = u + + if self.dataType == 'pressureHead': + PTv = self.P.T*v; + elif self.dataType == 'saturation': + dT = self.model.thetaDerivU(np.concatenate(Hs[1:]), m) + PTv = dT.T*self.P.T*v + + # This is done via backward substitution. + minus = 0 + BJtv = 0 + for ii in range(len(Hs)-1,0,-1): + Asub, Adiag, B = self.diagsJacobian(m, Hs[ii-1], Hs[ii]) + #select the correct part of v + vpart = range((ii-1)*Adiag.shape[0], (ii)*Adiag.shape[0]) + AdiaginvT = Solver(Adiag.T) + JTvC = AdiaginvT.solve(PTv[vpart] - minus) + minus = Asub.T*JTvC # this is now the super diagonal. + BJtv = BJtv + B.T*JTvC + + return BJtv + + + +if __name__ == '__main__': + from SimPEG import * + import Richards + import matplotlib.pyplot as plt + M = Mesh.TensorMesh([np.ones(40)]) + Ks = 9.4400e-03 + E = Richards.Haverkamp(Ks=np.log(Ks), A=1.1750e+06, gamma=4.74, alpha=1.6110e+06, theta_s=0.287, theta_r=0.075, beta=3.96) + bc = np.array([-61.5,-20.7]) + h = np.zeros(M.nC) + bc[0] + + # data = R + prob = Richards.RichardsProblem(M,E, timeStep=10, timeEnd=100, boundaryConditions=bc, initialConditions=h, doNewton=False, method='mixed') + + q = sp.csr_matrix((np.ones(4),(np.arange(4),np.array([20, 30, 35, 38]))), shape=(4,M.nCx)) + P = sp.kron(sp.identity(prob.numIts),q) + + prob.dataType = 'pressureHead' + mTrue = np.ones(M.nC)*np.log(Ks) + stdev = 0.01 # The standard deviation for the noise + data = prob.createSyntheticData(mTrue,std=stdev, P=P) + p = plt.plot(data.dobs.reshape((-1,4))) + plt.show() + # opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + # reg = Regularization.Tikhonov(model) + # inv = Inversion.BaseInversion(prob, reg, opt, beta0=1e4) + # derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)] + # print inv.dataObj(mTrue*0+np.log(1e-5)) + # print inv.dataObj(mTrue) + # tests.checkDerivative(derChk, mTrue, plotIt=False) + diff --git a/simpegFLOW/Richards/__init__.py b/simpegFLOW/Richards/__init__.py index 85657732..bf9832f7 100644 --- a/simpegFLOW/Richards/__init__.py +++ b/simpegFLOW/Richards/__init__.py @@ -1 +1,2 @@ -#blank! +from BaseRichards import * +from RichardsProblem import * diff --git a/simpegFLOW/Tests/test_Richards.py b/simpegFLOW/Tests/test_Richards.py index 72c1cbf5..08befb01 100644 --- a/simpegFLOW/Tests/test_Richards.py +++ b/simpegFLOW/Tests/test_Richards.py @@ -2,12 +2,11 @@ import unittest from SimPEG import * from SimPEG.Tests.TestUtils import OrderTest, checkDerivative from scipy.sparse.linalg import dsolve -import simpegFLOW.Richards +import simpegFLOW.Richards as Richards TOL = 1E-8 -class EmpiricalRelations(unittest.TestCase): - +class TestModels(unittest.TestCase): def test_BaseHaverkamp_Theta(self): mesh = Mesh.TensorMesh([50]) @@ -52,5 +51,213 @@ class EmpiricalRelations(unittest.TestCase): # 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): + +# def setUp(self): +# M = Mesh.TensorMesh([np.ones(20)]) +# M.setCellGradBC('dirichlet') + +# Ks = 9.4400e-03 +# E = Richards.Haverkamp(Ks=np.log(Ks), A=1.1750e+06, gamma=4.74, alpha=1.6110e+06, theta_s=0.287, theta_r=0.075, beta=3.96) + +# bc = np.array([-61.5,-20.7]) +# h = np.zeros(M.nC) + bc[0] +# prob = Richards.RichardsProblem(M,E, timeStep=60, timeEnd=180, 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) +# prob.P = P + +# self.h0 = h +# self.M = M +# self.Ks = Ks +# self.prob = prob + +# def test_Richards_getResidual_Newton(self): +# self.prob.doNewton = True +# passed = checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False) +# self.assertTrue(passed,True) + +# def test_Richards_getResidual_Picard(self): +# self.prob.doNewton = False +# passed = checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False, expectedOrder=1) +# self.assertTrue(passed,True) + +# def test_Adjoint_PressureHead(self): +# self.prob.dataType = 'pressureHead' +# Ks = self.Ks +# v = np.random.rand(self.prob.P.shape[0]) +# z = np.random.rand(self.M.nC) +# Hs = self.prob.field(np.log(Ks)) +# vJz = v.dot(self.prob.J(np.log(Ks),z,u=Hs)) +# zJv = z.dot(self.prob.Jt(np.log(Ks),v,u=Hs)) +# tol = TOL*(10**int(np.log10(zJv))) +# passed = np.abs(vJz - zJv) < tol +# print '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_Adjoint_Saturation(self): +# self.prob.dataType = 'saturation' +# Ks = self.Ks +# v = np.random.rand(self.prob.P.shape[0]) +# z = np.random.rand(self.M.nC) +# Hs = self.prob.field(np.log(Ks)) +# vJz = v.dot(self.prob.J(np.log(Ks),z,u=Hs)) +# zJv = z.dot(self.prob.Jt(np.log(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) +# self.assertTrue(passed,True) + +# def test_Sensitivity(self): +# self.prob.dataType = 'pressureHead' +# mTrue = np.ones(self.M.nC)*np.log(self.Ks) +# stdev = 0.01 # The standard deviation for the noise +# dobs = self.prob.createSyntheticData(mTrue,std=stdev)[0] +# self.prob.dobs = dobs +# self.prob.std = dobs*0 + stdev +# Hs = self.prob.field(mTrue) +# opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) +# reg = regularization.Regularization(self.M) +# inv = inverse.Inversion(self.prob, reg, opt, beta0=1e4) +# derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)] +# print 'Testing Richards Derivative' +# passed = checkDerivative(derChk, mTrue, num=5, plotIt=False) +# self.assertTrue(passed,True) + + + +# class RichardsTests2D(object): + +# def setUp(self): +# M = mesh.TensorMesh([np.ones(8),np.ones(30)]) +# Ks = 9.4400e-03 +# E = Richards.Haverkamp(Ks=np.log(Ks), A=1.1750e+06, gamma=4.74, alpha=1.6110e+06, theta_s=0.287, theta_r=0.075, beta=3.96) + +# bc = np.array([-61.5,-20.7]) +# 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, timeStep=60, timeEnd=180, boundaryConditions=bc, initialConditions=h, doNewton=False, method='mixed') + + +# XY = utils.ndgrid(np.array([5,7.]),np.array([5,15,25.])) +# q = M.getInterpolationMat(XY,'CC') +# P = sp.kron(sp.identity(prob.numIts),q) +# prob.P = P + +# self.h0 = h +# self.M = M +# self.Ks = Ks +# self.prob = prob + +# def test_Richards_getResidual_Newton(self): +# self.prob.doNewton = True +# passed = checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False) +# self.assertTrue(passed,True) + +# def test_Richards_getResidual_Picard(self): +# self.prob.doNewton = False +# passed = checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False, expectedOrder=1) +# self.assertTrue(passed,True) + +# def test_Adjoint_PressureHead(self): +# self.prob.dataType = 'pressureHead' +# Ks = self.Ks +# v = np.random.rand(self.prob.P.shape[0]) +# z = np.random.rand(self.M.nC) +# Hs = self.prob.field(np.log(Ks)) +# vJz = v.dot(self.prob.J(np.log(Ks),z,u=Hs)) +# zJv = z.dot(self.prob.Jt(np.log(Ks),v,u=Hs)) +# tol = TOL*(10**int(np.log10(zJv))) +# passed = np.abs(vJz - zJv) < tol +# print '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_Adjoint_Saturation(self): +# self.prob.dataType = 'saturation' +# Ks = self.Ks +# v = np.random.rand(self.prob.P.shape[0]) +# z = np.random.rand(self.M.nC) +# Hs = self.prob.field(np.log(Ks)) +# vJz = v.dot(self.prob.J(np.log(Ks),z,u=Hs)) +# zJv = z.dot(self.prob.Jt(np.log(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) +# self.assertTrue(passed,True) + +# def test_Sensitivity(self): +# self.prob.dataType = 'pressureHead' +# mTrue = np.ones(self.M.nC)*np.log(self.Ks) +# stdev = 0.01 # The standard deviation for the noise +# dobs = self.prob.createSyntheticData(mTrue,std=stdev)[0] +# self.prob.dobs = dobs +# self.prob.std = dobs*0 + stdev +# Hs = self.prob.field(mTrue) +# opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) +# reg = regularization.Regularization(self.M) +# inv = inverse.Inversion(self.prob, reg, opt, beta0=1e4) +# derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)] +# print 'Testing Richards Derivative' +# passed = checkDerivative(derChk, mTrue, num=5, plotIt=False) +# self.assertTrue(passed,True) + if __name__ == '__main__': unittest.main()