From 0514226ac2f3a24853c820a47bf89da0a41e9dbb Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 13 Nov 2013 15:07:40 -0800 Subject: [PATCH 01/18] Empirical Relations for Richards. --- SimPEG/forward/Richards.py | 219 ++++++++++++++++++++++++++++++++++ SimPEG/forward/__init__.py | 1 + SimPEG/tests/test_Richards.py | 66 ++++++++++ docs/api_Problem.rst | 7 ++ 4 files changed, 293 insertions(+) create mode 100644 SimPEG/forward/Richards.py create mode 100644 SimPEG/tests/test_Richards.py diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py new file mode 100644 index 00000000..a86f9697 --- /dev/null +++ b/SimPEG/forward/Richards.py @@ -0,0 +1,219 @@ +from SimPEG.forward import Problem +import numpy as np +from SimPEG.utils import sdiag, mkvc, setKwargs + + +class RichardsProblem(Problem): + """docstring for RichardsProblem""" + + dt = None + + + _method = 'mixed' + @property + def method(self): + """ + + Method must be either 'mixed' or 'head'. + + There are two different forms of Richards equation that differ + on how they deal with the non-linearity in the time-stepping term. + + The most fundamental form, referred to as the + 'mixed'-form of Richards Equation [Celia et al., 1990] + + .. math:: + + \\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. + This formulation of Richards equation is called the + 'mixed'-form because the equation is parameterized in psi + but the time-stepping is in terms of theta. + + As noted in [Celia et al., 1990] the 'head'-based form of Richards + 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 + \quad \psi \in \Omega + + However, it can be shown that this does not conserve mass in the discrete formulation. + + + """ + return self._method + @method.setter + def method(self, value): + assert value in ['mixed','head'] + self._method = value + + _doNewton = False + @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): + self._doNewton = value + + + def __init__(self, mesh, empirical): + Problem.__init__(self, mesh) + self.empirical = empirical + + + def getResidual(self, h, hn): + DIV = self.mesh.faceDiv + GRAD = self.mesh.cellGrad + BC = self.mesh.BC + AV = self.mesh.AV + Dz = self.mesh.Dz + + T = self.empirical.moistureContent(h) + dT = self.empirical.moistureContentDeriv(h) + Tn = self.empirical.moistureContent(hn) + K = self.empirical.hydraulicConductivity(h) + dK = self.empirical.hydraulicConductivityDeriv(h) + + aveK = 1./(AV*(1./K)); + + RHS = DIV*sdiag(aveK)*(GRAD*h+BC*P.bc) + Dz*aveK + if self.method is 'mixed': + r = (T-Tn)/dt - RHS + elif self.method is 'head': + r = dT*(h - hn)/dt - RHS + + J = + dT/dt - DIV*sdiag(aveK)*GRAD + if self.doNewton: + DDharmAve = sdiag(aveK**2)*AV*sdiag(K**(-2)) * dK + J = J - DIV*sdiag(GRAD*h + BC*P.bc)*DDharmAve - Dz*DDharmAve + + return r, J + +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): + setKwargs(self, **kwargs) + + 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 = 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 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 = sdiag(g) + #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 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): + setKwargs(self, **kwargs) + + def moistureContent(self, h): + m = 1 - 1/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 = sdiag(g) + return g + + def hydraulicConductivity(self, h): + alpha = self.alpha + I = self.I + n = self.n + Ks = self.Ks + m = 1 - 1/n + + theta_e = 1/((1+abs(alpha*h)**n)**m) + f = np.exp(Ks)*theta_e**I* ( ( 1 - ( 1 - theta_e**(1/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 hydraulicConductivityDeriv(self, h): + alpha = self.alpha + I = self.I + n = self.n + Ks = self.Ks + m = 1 - 1/n + + g = I*alpha*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2*(abs(alpha*h)**n + 1)**(1/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/(((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n)) + g[h >= 0] = 0 + g = sdiag(g) + + #alpha + # dA = I*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2*(abs(alpha*h)**n + 1)**(1/n - 2) - (2*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/(((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n)); + #n + # dn = 2*np.exp(Ks)*((np.log(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n))/n**2 + ((1/n - 1)*(((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1/n - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/((1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)) - np.log((abs(alpha*h)**n + 1)**(1/n - 1))/(n**2*(1/n - 1)**2*((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))))/(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n))*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1) - I*np.exp(Ks)*((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1/n - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2; + #I + # dI = np.exp(Ks)*np.log((abs(alpha*h)**n + 1)**(1/n - 1))*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2; + + return g diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py index 33c9a6b1..ec933f0a 100644 --- a/SimPEG/forward/__init__.py +++ b/SimPEG/forward/__init__.py @@ -2,3 +2,4 @@ from Problem import * import DCProblem from LinearProblem import LinearProblem import ModelTransforms +import Richards diff --git a/SimPEG/tests/test_Richards.py b/SimPEG/tests/test_Richards.py new file mode 100644 index 00000000..9c634e33 --- /dev/null +++ b/SimPEG/tests/test_Richards.py @@ -0,0 +1,66 @@ +import numpy as np +import unittest +from SimPEG.mesh import TensorMesh +from TestUtils import OrderTest, checkDerivative +from scipy.sparse.linalg import dsolve +from SimPEG.forward import Richards + + +class RichardsTests(unittest.TestCase): + + def setUp(self): + pass + # a = np.array([1, 1, 1]) + # b = np.array([1, 2]) + # c = np.array([1, 4]) + # self.mesh2 = TensorMesh([a, b], np.array([3, 5])) + # self.mesh3 = TensorMesh([a, b, c]) + + def test_VanGenuchten_moistureContent(self): + 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): + 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): + 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): + 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): + 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): + 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) + + + +if __name__ == '__main__': + unittest.main() diff --git a/docs/api_Problem.rst b/docs/api_Problem.rst index 068e30da..7fea7b5c 100644 --- a/docs/api_Problem.rst +++ b/docs/api_Problem.rst @@ -26,3 +26,10 @@ Linear Problem :members: :undoc-members: +Richards Problem +**************** + +.. automodule:: SimPEG.forward.Richards + :members: + :undoc-members: + From 81dea1a5dd5c30ac88b62b0cf9161fd6b09d55f1 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 14 Nov 2013 11:26:55 -0800 Subject: [PATCH 02/18] Added NewtonRoot finding algorithm. --- SimPEG/inverse/Optimize.py | 82 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 82 insertions(+) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 37c8b296..e47d61ae 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -643,6 +643,81 @@ class SteepestDescent(Minimize, Remember): def findSearchDirection(self): return -self.g + +class NewtonRoot(object): + """ + Newton Method - Root Finding + + root = newtonRoot(fun,x); + + Where fun is the function that returns the function value as well as the + gradient. + + For iterative solving of dh = -J\\r, use O.solveTol = TOL. For direct + solves, use SOLVETOL = 0 (default) + + Rowan Cockett + 16-May-2013 16:29:51 + University of British Columbia + rcockett@eos.ubc.ca + + """ + + tol = 1.000e-06; + solveTol = 0; # Default direct solve. + maxIter = 20; + stepDcr = 0.5; + maxLS = 30; + comments = False; + + def __init__(self, **kwargs): + setKwargs(self, **kwargs) + + def root(self, fun, x): + if self.comments: print 'Newton Method:\n' + + self._iter = 0 + while True: + + [r,J] = fun(x); + if self.solveTol == 0: + Jinv = Solver(J) + dh = - Jinv.solve(r) + else: + raise NotImplementedError('Iterative solve on NewtonRoot is not yet implemented.') + # M = @(x) tril(J)\(diag(J).*(triu(J)\x)); + # [dh, ~] = bicgstab(J,-r,O.solveTol,500,M); + + muLS = 1 + LScnt = 1 + + if self.comments: print '\tLinesearch:\n' + # Enter Linesearch + while True: + xt = x + muLS*dh + rt, Jt = fun(xt) # TODO: get rid of Jt + if self.comments: + print '\t\tResid: %e\n'%norm(rt) + if norm(rt) <= norm(r) or norm(rt) < self.tol: + break + + muLS = muLS*self.stepDcr + LScnt = LScnt + 1 + print '.' + if LScnt > self.maxLS: + print 'Newton Method: Line search break.' + root = NaN + return + + x = xt + self._iter += 1 + if norm(rt) < self.tol or self._iter > self.maxIter: + break + + return x + + + if __name__ == '__main__': from SimPEG.tests import Rosenbrock, checkDerivative import matplotlib.pyplot as plt @@ -657,3 +732,10 @@ if __name__ == '__main__': print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) xOpt = SteepestDescent(maxIter=30, maxIterLS=15,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock, x0) print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) + + + print 'test the newtonRoot finding.' + fun = lambda x: (np.sin(x), sdiag(np.cos(x))) + x = np.array([np.pi-0.3, np.pi+0.1, 0]) + pnt = NewtonRoot(comments=False).root(fun,x) + print pnt From 1e5969ec6a35811b142e97ee4d8b5c1564fdbc5e Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 14 Nov 2013 11:31:02 -0800 Subject: [PATCH 03/18] getResidual for Richards equation working. --- SimPEG/forward/Richards.py | 48 ++++++++++++++++++++++++++++------- SimPEG/tests/test_Richards.py | 24 +++++++++++++----- 2 files changed, 57 insertions(+), 15 deletions(-) diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py index a86f9697..08eb1242 100644 --- a/SimPEG/forward/Richards.py +++ b/SimPEG/forward/Richards.py @@ -6,7 +6,8 @@ from SimPEG.utils import sdiag, mkvc, setKwargs class RichardsProblem(Problem): """docstring for RichardsProblem""" - dt = None + timeStep = None + boundaryConditions = None _method = 'mixed' @@ -47,7 +48,7 @@ class RichardsProblem(Problem): return self._method @method.setter def method(self, value): - assert value in ['mixed','head'] + assert value in ['mixed','head'], "method must be 'mixed' or 'head'." self._method = value _doNewton = False @@ -57,20 +58,29 @@ class RichardsProblem(Problem): return self._doNewton @doNewton.setter def doNewton(self, value): + assert type(value) is bool, 'doNewton must be a boolean.' self._doNewton = value def __init__(self, mesh, empirical): Problem.__init__(self, mesh) self.empirical = empirical + self.mesh.setCellGradBC('dirichlet') - def getResidual(self, h, hn): + + def getResidual(self, hn, h): + """ + 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.BC - AV = self.mesh.AV - Dz = self.mesh.Dz + BC = self.mesh.cellGradBC + AV = self.mesh.aveCC2F + Dz = self.mesh.faceDiv #TODO: fix this for more than one dimension. + + bc = self.boundaryConditions + dt = self.timeStep T = self.empirical.moistureContent(h) dT = self.empirical.moistureContentDeriv(h) @@ -80,16 +90,16 @@ class RichardsProblem(Problem): aveK = 1./(AV*(1./K)); - RHS = DIV*sdiag(aveK)*(GRAD*h+BC*P.bc) + Dz*aveK + RHS = DIV*sdiag(aveK)*(GRAD*h+BC*bc) + Dz*aveK if self.method is 'mixed': r = (T-Tn)/dt - RHS elif self.method is 'head': r = dT*(h - hn)/dt - RHS - J = + dT/dt - DIV*sdiag(aveK)*GRAD + J = dT/dt - DIV*sdiag(aveK)*GRAD if self.doNewton: DDharmAve = sdiag(aveK**2)*AV*sdiag(K**(-2)) * dK - J = J - DIV*sdiag(GRAD*h + BC*P.bc)*DDharmAve - Dz*DDharmAve + J = J - DIV*sdiag(GRAD*h + BC*bc)*DDharmAve - Dz*DDharmAve return r, J @@ -217,3 +227,23 @@ class VanGenuchten(object): # dI = np.exp(Ks)*np.log((abs(alpha*h)**n + 1)**(1/n - 1))*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2; return g + + +if __name__ == '__main__': + from SimPEG.mesh import TensorMesh + from SimPEG.tests import checkDerivative + M = TensorMesh([np.ones(40)]) + Ks = 9.4400e-03 + E = 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) + + prob = RichardsProblem(M,E) + prob.timeStep = 1 + prob.boundaryConditions = np.array([-61.5,-20.7]) + prob.doNewton = True + prob.method = 'mixed' + + h = np.zeros(M.nC) + prob.boundaryConditions[0] + + checkDerivative(lambda hn1: prob.getResidual(h,hn1), h) + + diff --git a/SimPEG/tests/test_Richards.py b/SimPEG/tests/test_Richards.py index 9c634e33..36e044e4 100644 --- a/SimPEG/tests/test_Richards.py +++ b/SimPEG/tests/test_Richards.py @@ -9,12 +9,21 @@ from SimPEG.forward import Richards class RichardsTests(unittest.TestCase): def setUp(self): - pass - # a = np.array([1, 1, 1]) - # b = np.array([1, 2]) - # c = np.array([1, 4]) - # self.mesh2 = TensorMesh([a, b], np.array([3, 5])) - # self.mesh3 = TensorMesh([a, b, c]) + M = 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) + + prob = Richards.RichardsProblem(M,E) + prob.timeStep = 1 + prob.boundaryConditions = np.array([-61.5,-20.7]) + prob.doNewton = True + prob.method = 'mixed' + + h = np.zeros(M.nC) + prob.boundaryConditions[0] + + self.h0 = h + self.M = M + self.prob = prob def test_VanGenuchten_moistureContent(self): vanG = Richards.VanGenuchten() @@ -60,6 +69,9 @@ class RichardsTests(unittest.TestCase): passed = checkDerivative(wrapper, np.random.randn(n), plotIt=False) self.assertTrue(passed,True) + def test_Richards_getResidual(self): + checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False) + if __name__ == '__main__': From a5032e07b1bb8392b34b6d7c4a4f485aa7e9c0ec Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 14 Nov 2013 12:08:02 -0800 Subject: [PATCH 04/18] Minor updates to Utils documentation. --- docs/api_Solver.rst | 9 --------- docs/api_Utils.rst | 28 ++++++++++++++++++++++++++++ docs/index.rst | 1 - 3 files changed, 28 insertions(+), 10 deletions(-) delete mode 100644 docs/api_Solver.rst diff --git a/docs/api_Solver.rst b/docs/api_Solver.rst deleted file mode 100644 index db75860b..00000000 --- a/docs/api_Solver.rst +++ /dev/null @@ -1,9 +0,0 @@ -.. _api_Solver: - -Solver -****** - -.. automodule:: SimPEG.utils.Solver - :members: - :undoc-members: - diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst index 91d87d29..07903e3c 100644 --- a/docs/api_Utils.rst +++ b/docs/api_Utils.rst @@ -1,24 +1,52 @@ .. _api_Utils: + +Solver +****** + +.. automodule:: SimPEG.utils.Solver + :members: + :undoc-members: + + Utilities ********* +.. automodule:: SimPEG.utils + :members: + :undoc-members: + +Matrix Utilities +**************** + .. automodule:: SimPEG.utils.matutils :members: :undoc-members: +Sparse Utilities +**************** + .. automodule:: SimPEG.utils.sputils :members: :undoc-members: +LOM Utilities +************* + .. automodule:: SimPEG.utils.lomutils :members: :undoc-members: +Model Builder Utilities +*********************** + .. automodule:: SimPEG.utils.ModelBuilder :members: :undoc-members: +Interpolation Utilities +*********************** + .. automodule:: SimPEG.utils.interputils :members: :undoc-members: diff --git a/docs/index.rst b/docs/index.rst index 12db6f97..fb75f312 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -56,7 +56,6 @@ Utility Codes .. toctree:: :maxdepth: 2 - api_Solver api_Utils From 61e175cd12a7188b63d1f4ea287fdb07ac663ec1 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 14 Nov 2013 15:29:25 -0800 Subject: [PATCH 05/18] Fields calculation in Richards. Calculation of the diagonals of the Jacobian. --- SimPEG/forward/Richards.py | 121 +++++++++++++++++++++++++++++----- SimPEG/inverse/Optimize.py | 23 ++++--- SimPEG/tests/TestUtils.py | 8 +-- SimPEG/tests/test_Richards.py | 7 +- 4 files changed, 128 insertions(+), 31 deletions(-) diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py index 08eb1242..02318050 100644 --- a/SimPEG/forward/Richards.py +++ b/SimPEG/forward/Richards.py @@ -1,14 +1,21 @@ from SimPEG.forward import Problem import numpy as np from SimPEG.utils import sdiag, mkvc, setKwargs +from SimPEG.inverse import NewtonRoot class RichardsProblem(Problem): """docstring for RichardsProblem""" timeStep = None + timeEnd = None boundaryConditions = None + initialConditions = None + @property + def numIts(self): + """The number of iterations in the time domain problem.""" + return int(self.timeEnd/self.timeStep) _method = 'mixed' @property @@ -51,7 +58,6 @@ class RichardsProblem(Problem): assert value in ['mixed','head'], "method must be 'mixed' or 'head'." self._method = value - _doNewton = False @property def doNewton(self): """Do a Newton iteration. If False, a Picard iteration will be completed.""" @@ -59,15 +65,82 @@ class RichardsProblem(Problem): @doNewton.setter def doNewton(self, value): assert type(value) is bool, 'doNewton must be a boolean.' + self.rootFinder = NewtonRoot(doLS=value) self._doNewton = value + @property + def dataType(self): + """Choose how your data is collected, must be 'saturation' or 'pressureHead'.""" + assert value in ['saturation','pressureHead'], "dataType must be 'saturation' or 'pressureHead'." + return self._dataType + @dataType.setter + def dataType(self, value): + self._dataType = value - def __init__(self, mesh, empirical): + + + def __init__(self, mesh, empirical, **kwargs): Problem.__init__(self, mesh) self.empirical = empirical self.mesh.setCellGradBC('dirichlet') + self.doNewton = False # This also sets the rootFinder algorithm. + setKwargs(self, **kwargs) + + def field(self, m): + self.empirical.setModel(m) + Hs = range(self.numIts+1) + Hs[0] = self.initialConditions + for ii in range(self.numIts): + hn = Hs[ii] + hn1 = self.rootFinder.root(lambda hn1: self.getResidual(hn,hn1), hn) + Hs[ii+1] = hn1 + return Hs + def diagsJacobian(self, hn, hn1): + + DIV = self.mesh.faceDiv + GRAD = self.mesh.cellGrad + BC = self.mesh.cellGradBC + AV = self.mesh.aveCC2F + Dz = self.mesh.faceDiv #TODO: fix this for more than one dimension. + + bc = self.boundaryConditions + dt = self.timeStep + + dT = self.empirical.moistureContentDeriv(hn) + dT1 = self.empirical.moistureContentDeriv(hn1) + K1 = self.empirical.hydraulicConductivity(hn1) + dK1 = self.empirical.hydraulicConductivityDeriv(hn1) + dKa1 = self.empirical.hydraulicConductivityModelDeriv(hn1) + + # Compute part of the derivative of: + # + # DIV*diag(GRAD*hn1+BC*bc)*(AV*(1/K))^-1 + + DdiagGh1 = DIV*sdiag(GRAD*hn1+BC*bc) + diagAVk2_AVdiagK2 = sdiag((AV*(1./K1))**(-2)) * AV*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/dt)*dT + + Adiag = ( + (1/dt)*dT1 + -DdiagGh1*diagAVk2_AVdiagK2*dK1 + -DIV*sdiag(1./(AV*(1./K1)))*GRAD + -Dz*diagAVk2_AVdiagK2*dK1 + ) + + B = DdiagGh1*diagAVk2_AVdiagK2*dKa1 + Dz*diagAVk2_AVdiagK2*dKa1 def getResidual(self, hn, h): """ @@ -119,6 +192,9 @@ class Haverkamp(object): def __init__(self, **kwargs): 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) @@ -141,14 +217,17 @@ class Haverkamp(object): f[h >= 0] = np.exp(self.Ks) return f - 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 = sdiag(g) + 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 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 = sdiag(g) return g @@ -180,6 +259,9 @@ class VanGenuchten(object): def __init__(self, **kwargs): setKwargs(self, **kwargs) + def setModel(self, m): + self.Ks = m + def moistureContent(self, h): m = 1 - 1/self.n; f = (( self.theta_s - self.theta_r )/ @@ -208,6 +290,15 @@ class VanGenuchten(object): 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/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2*(abs(alpha*h)**n + 1)**(1/n - 2) - (2*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/(((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n)); + #n + # dn = 2*np.exp(Ks)*((np.log(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n))/n**2 + ((1/n - 1)*(((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1/n - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/((1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)) - np.log((abs(alpha*h)**n + 1)**(1/n - 1))/(n**2*(1/n - 1)**2*((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))))/(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n))*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1) - I*np.exp(Ks)*((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1/n - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2; + #I + # dI = np.exp(Ks)*np.log((abs(alpha*h)**n + 1)**(1/n - 1))*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2; + return sdiag(self.hydraulicConductivity(h)) # This assumes that the the model is Ks + def hydraulicConductivityDeriv(self, h): alpha = self.alpha I = self.I @@ -218,14 +309,6 @@ class VanGenuchten(object): g = I*alpha*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2*(abs(alpha*h)**n + 1)**(1/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/(((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n)) g[h >= 0] = 0 g = sdiag(g) - - #alpha - # dA = I*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2*(abs(alpha*h)**n + 1)**(1/n - 2) - (2*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/(((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n)); - #n - # dn = 2*np.exp(Ks)*((np.log(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n))/n**2 + ((1/n - 1)*(((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1/n - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/((1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)) - np.log((abs(alpha*h)**n + 1)**(1/n - 1))/(n**2*(1/n - 1)**2*((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))))/(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n))*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1) - I*np.exp(Ks)*((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1/n - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2; - #I - # dI = np.exp(Ks)*np.log((abs(alpha*h)**n + 1)**(1/n - 1))*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2; - return g @@ -244,6 +327,12 @@ if __name__ == '__main__': h = np.zeros(M.nC) + prob.boundaryConditions[0] - checkDerivative(lambda hn1: prob.getResidual(h,hn1), h) - + numIts = 10 + Hs = range(numIts+1) + Hs[0] = h + for ii in range(numIts): + hn = Hs[ii] + hn1 = NewtonRoot().root(lambda hn1: prob.getResidual(hn,hn1), hn) + Hs[ii+1] = hn1 + M.plotImage(hn1,showIt=True) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index e47d61ae..10493c68 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -663,12 +663,13 @@ class NewtonRoot(object): """ - tol = 1.000e-06; - solveTol = 0; # Default direct solve. - maxIter = 20; - stepDcr = 0.5; - maxLS = 30; - comments = False; + tol = 1.000e-06 + solveTol = 0 # Default direct solve. + maxIter = 20 + stepDcr = 0.5 + maxLS = 30 + comments = False + doLS = True def __init__(self, **kwargs): setKwargs(self, **kwargs) @@ -688,14 +689,14 @@ class NewtonRoot(object): # M = @(x) tril(J)\(diag(J).*(triu(J)\x)); # [dh, ~] = bicgstab(J,-r,O.solveTol,500,M); - muLS = 1 + muLS = 1. LScnt = 1 + xt = x + dh + rt, Jt = fun(xt) # TODO: get rid of Jt if self.comments: print '\tLinesearch:\n' # Enter Linesearch - while True: - xt = x + muLS*dh - rt, Jt = fun(xt) # TODO: get rid of Jt + while True and self.doLS: if self.comments: print '\t\tResid: %e\n'%norm(rt) if norm(rt) <= norm(r) or norm(rt) < self.tol: @@ -708,6 +709,8 @@ class NewtonRoot(object): print 'Newton Method: Line search break.' root = NaN return + xt = x + muLS*dh + rt, Jt = fun(xt) # TODO: get rid of Jt x = xt self._iter += 1 diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index 1cc2bbae..662a7884 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -189,7 +189,7 @@ def Rosenbrock(x, return_g=True, return_H=True): out += (H,) return out if len(out) > 1 else out[0] -def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): +def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.9, eps=1e-10): """ Basic derivative check @@ -201,6 +201,9 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): :param int num: number of times to reduce step length, h :param bool plotIt: if you would like to plot :param numpy.array dx: step direction + :param int expectedOrder: The order that you expect the derivative to yield. + :param float tolerance: The tolerance on the expected order. + :param float eps: What is zero? :rtype: bool :return: did you pass the test?! @@ -243,9 +246,6 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): order1 = np.log10(E1[:-1]/E1[1:]) print "%d\t%1.2e\t%1.3e\t\t%1.3e\t\t%1.3f" % (i, t[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1]) - tolerance = 0.9 - expectedOrder = 2 - eps = 1e-10 order0 = order0[E0[1:] > eps] order1 = order1[E1[1:] > eps] belowTol = order1.size == 0 and order0.size > 0 diff --git a/SimPEG/tests/test_Richards.py b/SimPEG/tests/test_Richards.py index 36e044e4..8f0be566 100644 --- a/SimPEG/tests/test_Richards.py +++ b/SimPEG/tests/test_Richards.py @@ -69,9 +69,14 @@ class RichardsTests(unittest.TestCase): passed = checkDerivative(wrapper, np.random.randn(n), plotIt=False) self.assertTrue(passed,True) - def test_Richards_getResidual(self): + def test_Richards_getResidual_Newton(self): + self.prob.doNewton = True checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False) + def test_Richards_getResidual_Picard(self): + self.prob.doNewton = False + checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False, expectedOrder=1) + if __name__ == '__main__': From 8f4e1174c816b052f4a999fc4ffb98ef9c9f45d1 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 14 Nov 2013 15:53:19 -0800 Subject: [PATCH 06/18] Added J and Jt, untested as of yet. --- SimPEG/forward/Richards.py | 50 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py index 02318050..7bcd1660 100644 --- a/SimPEG/forward/Richards.py +++ b/SimPEG/forward/Richards.py @@ -176,6 +176,56 @@ class RichardsProblem(Problem): return r, J + def J(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(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(Hs[ii],Hs[ii+1]) + Adiaginv = Solver(Adiag) + JvC[ii] = Adiaginv.solve(B*v - Asub*JvC[ii-1]) + + if self.dataType is 'pressureHead': + Jv = P.Q*np.concatenate(JvC) + elif self.dataType is 'saturation': + dT = self.empirical.moistureContentDeriv(np.concatenate(Hs[1:])) + Jv = P.Q*dT*np.concatenate(JvC) + + def Jt(self, m, v, u=None): + if u is None: + u = self.field(m) + Hs = u + + if self.dataType is 'pressureHead': + QTz = P.Q.T*z; + elif self.dataType is 'saturation': + dT = self.empirical.moistureContentDeriv(np.concatenate(Hs[1:])) + QTz = dT.T*P.Q.T*z + + # This is done via backward substitution. + minus = 0 + BJtz = 0 + for ii in range(numel(Hs)-1,-1,-1): + Asub, Adiag, B = self.diagsJacobian(Hs[ii-1], Hs[ii]) + #select the correct part of z + # TODO: This might be easier by reshaping the array + zpart = range((ii-2)*Adiag.shape[0], (ii-1)*Adiag.shape[0]) + AdiaginvT = Solver(Adiag.T) + JTzC = AdiaginvT.solve(QTz[zpart] - minus) + minus = Asub.T*JTzC # this is now the super diagonal. + BJtz = BJtz + B.T*JTzC + + class Haverkamp(object): """docstring for Haverkamp""" From bc49c6058fec3b0fdf2b9bce301f84bf313a9f57 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 14 Nov 2013 18:01:53 -0800 Subject: [PATCH 07/18] Bug Fixes and Adjoint Test. --- SimPEG/forward/Richards.py | 34 +++++++++++++--------- SimPEG/tests/TestUtils.py | 2 +- SimPEG/tests/test_Richards.py | 54 +++++++++++++++++++++++++++++------ 3 files changed, 67 insertions(+), 23 deletions(-) diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py index 7bcd1660..e472ca6a 100644 --- a/SimPEG/forward/Richards.py +++ b/SimPEG/forward/Richards.py @@ -1,6 +1,6 @@ from SimPEG.forward import Problem import numpy as np -from SimPEG.utils import sdiag, mkvc, setKwargs +from SimPEG.utils import sdiag, mkvc, setKwargs, Solver from SimPEG.inverse import NewtonRoot @@ -71,10 +71,10 @@ class RichardsProblem(Problem): @property def dataType(self): """Choose how your data is collected, must be 'saturation' or 'pressureHead'.""" - assert value in ['saturation','pressureHead'], "dataType must be 'saturation' or 'pressureHead'." return self._dataType @dataType.setter def dataType(self, value): + assert value in ['saturation','pressureHead'], "dataType must be 'saturation' or 'pressureHead'." self._dataType = value @@ -83,6 +83,7 @@ class RichardsProblem(Problem): Problem.__init__(self, mesh) self.empirical = empirical self.mesh.setCellGradBC('dirichlet') + self.dataType = 'pressureHead' self.doNewton = False # This also sets the rootFinder algorithm. setKwargs(self, **kwargs) @@ -142,6 +143,8 @@ class RichardsProblem(Problem): B = DdiagGh1*diagAVk2_AVdiagK2*dKa1 + Dz*diagAVk2_AVdiagK2*dKa1 + return Asub, Adiag, B + def getResidual(self, hn, h): """ Where h is the proposed value for the next time iterate (h_{n+1}) @@ -196,10 +199,12 @@ class RichardsProblem(Problem): JvC[ii] = Adiaginv.solve(B*v - Asub*JvC[ii-1]) if self.dataType is 'pressureHead': - Jv = P.Q*np.concatenate(JvC) + Jv = self.P*np.concatenate(JvC) elif self.dataType is 'saturation': dT = self.empirical.moistureContentDeriv(np.concatenate(Hs[1:])) - Jv = P.Q*dT*np.concatenate(JvC) + Jv = self.P*dT*np.concatenate(JvC) + + return Jv def Jt(self, m, v, u=None): if u is None: @@ -207,23 +212,24 @@ class RichardsProblem(Problem): Hs = u if self.dataType is 'pressureHead': - QTz = P.Q.T*z; + PTv = self.P.T*v; elif self.dataType is 'saturation': dT = self.empirical.moistureContentDeriv(np.concatenate(Hs[1:])) - QTz = dT.T*P.Q.T*z + PTv = dT.T*self.P.T*v # This is done via backward substitution. minus = 0 - BJtz = 0 - for ii in range(numel(Hs)-1,-1,-1): + BJtv = 0 + for ii in range(len(Hs)-1,0,-1): Asub, Adiag, B = self.diagsJacobian(Hs[ii-1], Hs[ii]) - #select the correct part of z - # TODO: This might be easier by reshaping the array - zpart = range((ii-2)*Adiag.shape[0], (ii-1)*Adiag.shape[0]) + #select the correct part of v + vpart = range((ii-1)*Adiag.shape[0], (ii)*Adiag.shape[0]) AdiaginvT = Solver(Adiag.T) - JTzC = AdiaginvT.solve(QTz[zpart] - minus) - minus = Asub.T*JTzC # this is now the super diagonal. - BJtz = BJtz + B.T*JTzC + JTvC = AdiaginvT.solve(PTv[vpart] - minus) + minus = Asub.T*JTvC # this is now the super diagonal. + BJtv = BJtv + B.T*JTvC + + return BJtv class Haverkamp(object): diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index 662a7884..3b744f76 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -189,7 +189,7 @@ def Rosenbrock(x, return_g=True, return_H=True): out += (H,) return out if len(out) > 1 else out[0] -def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.9, eps=1e-10): +def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.85, eps=1e-10): """ Basic derivative check diff --git a/SimPEG/tests/test_Richards.py b/SimPEG/tests/test_Richards.py index 8f0be566..e64d2759 100644 --- a/SimPEG/tests/test_Richards.py +++ b/SimPEG/tests/test_Richards.py @@ -1,4 +1,5 @@ import numpy as np +import scipy.sparse as sp import unittest from SimPEG.mesh import TensorMesh from TestUtils import OrderTest, checkDerivative @@ -6,6 +7,8 @@ from scipy.sparse.linalg import dsolve from SimPEG.forward import Richards +TOL = 1E-8 + class RichardsTests(unittest.TestCase): def setUp(self): @@ -13,16 +16,18 @@ class RichardsTests(unittest.TestCase): 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) - prob = Richards.RichardsProblem(M,E) - prob.timeStep = 1 - prob.boundaryConditions = np.array([-61.5,-20.7]) - prob.doNewton = True - prob.method = 'mixed' + bc = np.array([-61.5,-20.7]) + h = np.zeros(M.nC) + bc[0] + prob = Richards.RichardsProblem(M,E, timeStep=30, timeEnd=360, 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.P = P - h = np.zeros(M.nC) + prob.boundaryConditions[0] self.h0 = h self.M = M + self.Ks = Ks self.prob = prob def test_VanGenuchten_moistureContent(self): @@ -71,11 +76,44 @@ class RichardsTests(unittest.TestCase): def test_Richards_getResidual_Newton(self): self.prob.doNewton = True - checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False) + 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 - checkDerivative(lambda hn1: self.prob.getResidual(self.h0,hn1), self.h0, plotIt=False, expectedOrder=1) + 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) + + From 7663d127072d0f84e410430c1c386698ac5d3ed6 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 14 Nov 2013 18:28:28 -0800 Subject: [PATCH 08/18] Removed Synthetic Problem class, and added method to Problem. --- SimPEG/forward/DCProblem.py | 20 +++++++------------- SimPEG/forward/Problem.py | 18 +++--------------- SimPEG/forward/Richards.py | 6 +----- SimPEG/tests/test_Richards.py | 2 +- SimPEG/tests/test_forward_DCproblem.py | 13 ++++--------- 5 files changed, 16 insertions(+), 43 deletions(-) diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index 34046a52..074cdb8b 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -1,5 +1,5 @@ from SimPEG.mesh import TensorMesh -from SimPEG.forward import Problem, SyntheticProblem, ModelTransforms +from SimPEG.forward import Problem, ModelTransforms from SimPEG.tests import checkDerivative from SimPEG.utils import ModelBuilder, sdiag, mkvc from SimPEG import Solver @@ -201,23 +201,17 @@ if __name__ == '__main__': P = Q.T # Create some data - class syntheticDCProblem(DCProblem, SyntheticProblem): - pass + problem = DCProblem(mesh) + problem.P = P + problem.RHS = q + dobs, Wd = problem.createSyntheticData(mSynth, std=0.05) - synthetic = syntheticDCProblem(mesh); - synthetic.P = P - synthetic.RHS = q - dobs, Wd = synthetic.createData(mSynth, std=0.05) - - u = synthetic.field(mSynth) - u = synthetic.reshapeFields(u) + u = problem.field(mSynth) + u = problem.reshapeFields(u) mesh.plotImage(u[:,10]) # plt.show() # Now set up the problem to do some minimization - problem = DCProblem(mesh) - problem.P = P - problem.RHS = q problem.dobs = dobs problem.std = dobs*0 + 0.05 m0 = mesh.gridCC[:,0]*0+sig2 diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index cf22baae..35d41cb8 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -218,27 +218,15 @@ class Problem(object): """ return sp.eye(m.size) - - - -class SyntheticProblem(object): - """ - Has helpful functions when dealing with synthetic problems - - To use this class, inherit to your problem:: - - class mySyntheticExample(Problem, SyntheticProblem): - pass - """ - def createData(self, m, std=0.05): + def createSyntheticData(self, m, std=0.05): """ + Create synthetic data given a model, and a standard deviation. + :param numpy.array m: geophysical model :param numpy.array std: standard deviation :rtype: numpy.array, numpy.array :return: dobs, Wd - Create synthetic data given a model, and a standard deviation. - Returns the observed data with random Gaussian noise and Wd which is the same size as data, and can be used to weight the inversion. """ diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py index e472ca6a..d6d268be 100644 --- a/SimPEG/forward/Richards.py +++ b/SimPEG/forward/Richards.py @@ -78,7 +78,6 @@ class RichardsProblem(Problem): self._dataType = value - def __init__(self, mesh, empirical, **kwargs): Problem.__init__(self, mesh) self.empirical = empirical @@ -92,12 +91,9 @@ class RichardsProblem(Problem): Hs = range(self.numIts+1) Hs[0] = self.initialConditions for ii in range(self.numIts): - hn = Hs[ii] - hn1 = self.rootFinder.root(lambda hn1: self.getResidual(hn,hn1), hn) - Hs[ii+1] = hn1 + Hs[ii+1] = self.rootFinder.root(lambda hn1: self.getResidual(Hs[ii],hn1), Hs[ii]) return Hs - def diagsJacobian(self, hn, hn1): DIV = self.mesh.faceDiv diff --git a/SimPEG/tests/test_Richards.py b/SimPEG/tests/test_Richards.py index e64d2759..86f1a7e7 100644 --- a/SimPEG/tests/test_Richards.py +++ b/SimPEG/tests/test_Richards.py @@ -85,7 +85,7 @@ class RichardsTests(unittest.TestCase): self.assertTrue(passed,True) def test_Adjoint_PressureHead(self): - # self.prob.dataType = 'pressureHead' + self.prob.dataType = 'pressureHead' Ks = self.Ks v = np.random.rand(self.prob.P.shape[0]) z = np.random.rand(self.M.nC) diff --git a/SimPEG/tests/test_forward_DCproblem.py b/SimPEG/tests/test_forward_DCproblem.py index c6e6f9c2..847d1055 100644 --- a/SimPEG/tests/test_forward_DCproblem.py +++ b/SimPEG/tests/test_forward_DCproblem.py @@ -2,7 +2,7 @@ import numpy as np import unittest from SimPEG.mesh import TensorMesh from SimPEG.utils import ModelBuilder, sdiag -from SimPEG.forward import Problem, SyntheticProblem +from SimPEG.forward import Problem from SimPEG.forward.DCProblem import * from TestUtils import checkDerivative from scipy.sparse.linalg import dsolve @@ -40,18 +40,13 @@ class DCProblemTests(unittest.TestCase): P = Q.T # Create some data - class syntheticDCProblem(DCProblem, SyntheticProblem): - pass - synthetic = syntheticDCProblem(mesh); - synthetic.P = P - synthetic.RHS = q - dobs, Wd = synthetic.createData(mSynth, std=0.05) - - # Now set up the problem to do some minimization problem = DCProblem(mesh) problem.P = P problem.RHS = q + dobs, Wd = problem.createSyntheticData(mSynth, std=0.05) + + # Now set up the problem to do some minimization problem.W = Wd problem.dobs = dobs problem.std = dobs*0 + 0.05 From e98e41f34c510d0fc52942df8cb11881a08467e6 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 15 Nov 2013 10:10:10 -0800 Subject: [PATCH 09/18] dpred --- SimPEG/forward/Richards.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py index d6d268be..780902ab 100644 --- a/SimPEG/forward/Richards.py +++ b/SimPEG/forward/Richards.py @@ -86,6 +86,20 @@ class RichardsProblem(Problem): self.doNewton = False # This also sets the rootFinder algorithm. setKwargs(self, **kwargs) + def dpred(self, m, u=None): + """ + Predicted data. + + .. math:: + d_\\text{pred} = Pu(m) + """ + if u is None: + u = self.field(m) + u = np.concatenate(u[1:]) + if self.dataType is 'saturation': + u = self.empirical.moistureContent(u) + return self.P*u + def field(self, m): self.empirical.setModel(m) Hs = range(self.numIts+1) From 278eec94ecb9ec410a8612b624813ceb60e6959e Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 15 Nov 2013 13:59:41 -0800 Subject: [PATCH 10/18] Ensure that (1/dt) gives a float not an int. Use (1.0/dt). Added fullJ calculation. Tested Sensitivity. --- SimPEG/forward/Richards.py | 102 +++++++++++++++++++++++----------- SimPEG/inverse/Inversion.py | 3 + SimPEG/tests/test_Richards.py | 21 +++++-- 3 files changed, 89 insertions(+), 37 deletions(-) diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py index 780902ab..99090d84 100644 --- a/SimPEG/forward/Richards.py +++ b/SimPEG/forward/Richards.py @@ -1,17 +1,25 @@ from SimPEG.forward import Problem import numpy as np -from SimPEG.utils import sdiag, mkvc, setKwargs, Solver +from SimPEG.utils import sdiag, spzeros, mkvc, setKwargs, Solver from SimPEG.inverse import NewtonRoot +import scipy.sparse as sp class RichardsProblem(Problem): """docstring for RichardsProblem""" - timeStep = None timeEnd = None boundaryConditions = None initialConditions = None + @property + def timeStep(self): + """The time between steps.""" + return self._timeStep + @timeStep.setter + def timeStep(self, value): + self._timeStep = float(value) + @property def numIts(self): """The number of iterations in the time domain problem.""" @@ -127,7 +135,7 @@ class RichardsProblem(Problem): # Compute part of the derivative of: # - # DIV*diag(GRAD*hn1+BC*bc)*(AV*(1/K))^-1 + # DIV*diag(GRAD*hn1+BC*bc)*(AV*(1.0/K))^-1 DdiagGh1 = DIV*sdiag(GRAD*hn1+BC*bc) diagAVk2_AVdiagK2 = sdiag((AV*(1./K1))**(-2)) * AV*sdiag(K1**(-2)) @@ -142,10 +150,10 @@ class RichardsProblem(Problem): # | Asub Adiag | | hn | | bn | # - - - - - - - Asub = (-1/dt)*dT + Asub = (-1.0/dt)*dT Adiag = ( - (1/dt)*dT1 + (1.0/dt)*dT1 -DdiagGh1*diagAVk2_AVdiagK2*dK1 -DIV*sdiag(1./(AV*(1./K1)))*GRAD -Dz*diagAVk2_AVdiagK2*dK1 @@ -189,6 +197,26 @@ class RichardsProblem(Problem): 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(Hs[ii],Hs[ii+1]) + Ad = sp.block_diag(Adiags) + zRight = spzeros((len(Asubs)-1)*Asubs[0].shape[0],Adiags[0].shape[1]) + zTop = 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 J(self, m, v, u=None): if u is None: u = self.field(m) @@ -329,7 +357,7 @@ class VanGenuchten(object): self.Ks = m def moistureContent(self, h): - m = 1 - 1/self.n; + 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 @@ -346,10 +374,10 @@ class VanGenuchten(object): I = self.I n = self.n Ks = self.Ks - m = 1 - 1/n + m = 1 - 1.0/n - theta_e = 1/((1+abs(alpha*h)**n)**m) - f = np.exp(Ks)*theta_e**I* ( ( 1 - ( 1 - theta_e**(1/m) )**m )**2 ) + theta_e = 1.0/((1+abs(alpha*h)**n)**m) + f = np.exp(Ks)*theta_e**I* ( ( 1 - ( 1 - 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: @@ -358,11 +386,11 @@ class VanGenuchten(object): def hydraulicConductivityModelDeriv(self, h): #alpha - # dA = I*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2*(abs(alpha*h)**n + 1)**(1/n - 2) - (2*h*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/(((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n)); + # 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/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n))/n**2 + ((1/n - 1)*(((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1/n - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/((1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)) - np.log((abs(alpha*h)**n + 1)**(1/n - 1))/(n**2*(1/n - 1)**2*((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))))/(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n))*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1) - I*np.exp(Ks)*((np.log(abs(alpha*h)**n + 1)*(abs(alpha*h)**n + 1)**(1/n - 1))/n**2 - abs(alpha*h)**n*np.log(abs(alpha*h))*(1/n - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2; + # 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/n - 1))*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2; + # 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 sdiag(self.hydraulicConductivity(h)) # This assumes that the the model is Ks def hydraulicConductivityDeriv(self, h): @@ -370,35 +398,43 @@ class VanGenuchten(object): I = self.I n = self.n Ks = self.Ks - m = 1 - 1/n + m = 1 - 1.0/n - g = I*alpha*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**(I - 1)*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)**2*(abs(alpha*h)**n + 1)**(1/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*h)**(n - 1)*np.sign(alpha*h)*(1/n - 1)*((abs(alpha*h)**n + 1)**(1/n - 1))**I*((1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1 - 1/n) - 1)*(abs(alpha*h)**n + 1)**(1/n - 2))/(((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1) + 1)*(1 - 1/((abs(alpha*h)**n + 1)**(1/n - 1))**(1/(1/n - 1)))**(1/n)) + g = I*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 - 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 = sdiag(g) return g if __name__ == '__main__': - from SimPEG.mesh import TensorMesh - from SimPEG.tests import checkDerivative - M = TensorMesh([np.ones(40)]) + import SimPEG + from SimPEG import mesh, inverse, regularization, tests + import scipy.sparse as sp + import numpy as np + from SimPEG.forward import Problem, Richards + M = mesh.TensorMesh([np.ones(40)]) Ks = 9.4400e-03 - E = 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) + 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=10, timeEnd=60, boundaryConditions=bc, initialConditions=h, doNewton=False, method='mixed') - prob = RichardsProblem(M,E) - prob.timeStep = 1 - prob.boundaryConditions = np.array([-61.5,-20.7]) - prob.doNewton = True - prob.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.P = P - h = np.zeros(M.nC) + prob.boundaryConditions[0] - - numIts = 10 - Hs = range(numIts+1) - Hs[0] = h - for ii in range(numIts): - hn = Hs[ii] - hn1 = NewtonRoot().root(lambda hn1: prob.getResidual(hn,hn1), hn) - Hs[ii+1] = hn1 - M.plotImage(hn1,showIt=True) + prob.dataType = 'pressureHead' + mTrue = np.ones(M.nC)*np.log(Ks) + stdev = 0.01 # The standard deviation for the noise + dobs = prob.createSyntheticData(mTrue,std=stdev)[0] + # p = plot(dobs.reshape((-1,4))) + prob.dobs = dobs + prob.std = dobs*0 + stdev + opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + reg = regularization.Regularization(mesh) + inv = inverse.Inversion(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/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index e3d500b1..871a7f06 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -38,6 +38,9 @@ class BaseInversion(object): eps = np.linalg.norm(mkvc(self.prob.dobs),2)*1e-5 self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps) return self._Wd + @Wd.setter + def Wd(self, value): + self._Wd = value @property def phi_d_target(self): diff --git a/SimPEG/tests/test_Richards.py b/SimPEG/tests/test_Richards.py index 86f1a7e7..4510aab0 100644 --- a/SimPEG/tests/test_Richards.py +++ b/SimPEG/tests/test_Richards.py @@ -1,7 +1,7 @@ import numpy as np import scipy.sparse as sp import unittest -from SimPEG.mesh import TensorMesh +from SimPEG import mesh, regularization, inverse from TestUtils import OrderTest, checkDerivative from scipy.sparse.linalg import dsolve from SimPEG.forward import Richards @@ -12,7 +12,7 @@ TOL = 1E-8 class RichardsTests(unittest.TestCase): def setUp(self): - M = TensorMesh([np.ones(40)]) + 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) @@ -113,8 +113,21 @@ class RichardsTests(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_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(mesh) + 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__': From fc435feb9e2ce03ff3c62611ce5f4fd05ca14c01 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 19 Nov 2013 15:59:37 -0800 Subject: [PATCH 11/18] BFGS implementation in the notebook. --- SimPEG/inverse/Optimize.py | 10 +- notebooks/BFGS.ipynb | 197 +++++++++++++++++++++++++++++++++++++ 2 files changed, 203 insertions(+), 4 deletions(-) create mode 100644 notebooks/BFGS.ipynb diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 37c8b296..8a9abaed 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -76,7 +76,7 @@ class IterationPrinters(object): itType = {"title": "itType", "value": lambda M: M._itType, "width": 8, "format": "%s"} aSet = {"title": "aSet", "value": lambda M: np.sum(M.activeSet(M.xc)), "width": 8, "format": "%d"} bSet = {"title": "bSet", "value": lambda M: np.sum(M.bindingSet(M.xc)), "width": 8, "format": "%d"} - comment = {"title": "Comment", "value": lambda M: M.projComment, "width": 7, "format": "%s"} + comment = {"title": "Comment", "value": lambda M: M.comment, "width": 12, "format": "%s"} beta = {"title": "beta", "value": lambda M: M.parent._beta, "width": 10, "format": "%1.2e"} phi_d = {"title": "phi_d", "value": lambda M: M.parent.phi_d, "width": 10, "format": "%1.2e"} @@ -106,6 +106,8 @@ class Minimize(object): debug = False debugLS = False + comment = '' + def __init__(self, **kwargs): self._id = int(np.random.rand()*1e6) # create a unique identifier to this program to be used in pubsub self.stoppers = [StoppingCriteria.tolerance_f, StoppingCriteria.moving_x, StoppingCriteria.tolerance_g, StoppingCriteria.norm_g, StoppingCriteria.iteration] @@ -525,7 +527,7 @@ class ProjectedGradient(Minimize, Remember): self.stopDoingPG = False self._itType = 'SD' - self.projComment = '' + self.comment = '' self.aSet_prev = self.activeSet(x0) @@ -600,7 +602,7 @@ class ProjectedGradient(Minimize, Remember): self.exploreCG = np.all(aSet == bSet) # explore conjugate gradient f_current_decrease = self.f_last - self.f - self.projComment = '' + self.comment = '' if self._iter < 1: # Note that this is reset on every CG iteration. self.f_decrease_max = -np.inf @@ -608,7 +610,7 @@ class ProjectedGradient(Minimize, Remember): self.f_decrease_max = max(self.f_decrease_max, f_current_decrease) self.stopDoingPG = f_current_decrease < 0.25 * self.f_decrease_max if self.stopDoingPG: - self.projComment = 'Stop SD' + self.comment = 'Stop SD' self.explorePG = False self.exploreCG = True # implement 3.8, MoreToraldo91 diff --git a/notebooks/BFGS.ipynb b/notebooks/BFGS.ipynb new file mode 100644 index 00000000..23de6e06 --- /dev/null +++ b/notebooks/BFGS.ipynb @@ -0,0 +1,197 @@ +{ + "metadata": { + "name": "" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import SimPEG\n", + "from SimPEG import Solver\n", + "from SimPEG.mesh import TensorMesh\n", + "from SimPEG.regularization import Regularization\n", + "import SimPEG.inverse as inverse\n", + "from SimPEG.inverse import Minimize, Remember, IterationPrinters\n", + "import numpy as np\n", + "import scipy.sparse as sp" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 1 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "FUN = SimPEG.tests.Rosenbrock\n", + "FUN = SimPEG.tests.getQuadratic(sp.csr_matrix(([100,1],([0,1],[0,1])),shape=(2,2)),np.array([-5,-5]),100)\n", + "\n", + "\n", + "class BFGS(Minimize, Remember):\n", + " name = 'BFGS'\n", + " nbfgs = 10\n", + " \n", + " @property\n", + " def H0(self):\n", + " \"\"\"\n", + " Approximate Hessian used in preconditioning the problem.\n", + " \n", + " Must be a SimPEG.Solver\n", + " \"\"\"\n", + " _H0 = getattr(self,'_H0',None)\n", + " if _H0 is None:\n", + " return Solver(sp.identity(self.xc.size).tocsc())\n", + " return _H0\n", + " @H0.setter\n", + " def H0(self, value):\n", + " self._H0 = value\n", + " \n", + " def _startup_BFGS(self,x0):\n", + " self._bfgscnt = -1\n", + " self._bfgsY = np.zeros((x0.size, self.nbfgs))\n", + " self._bfgsS = np.zeros((x0.size, self.nbfgs))\n", + " if not np.any([p is IterationPrinters.comment for p in self.printers]):\n", + " self.printers.append(IterationPrinters.comment)\n", + " \n", + " def bfgs(self,n,d):\n", + " nn = min(self._bfgsS.shape[1],n)\n", + " ktop = nn\n", + " d = self.bfgsrec(ktop,n,nn,self._bfgsS,self._bfgsY,d)\n", + " return d\n", + "\n", + " def bfgsrec(self,k,n,nn,S,Y,d):\n", + " \"\"\"BFGS recursion\"\"\"\n", + " if k < 0:\n", + " d = self.H0.solve(d)\n", + " else:\n", + " khat = mod(n-nn+k,nn)\n", + " gamma = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat])\n", + " d = d - gamma*Y[:,khat]\n", + " d = self.bfgsrec(k-1,n,nn,S,Y,d)\n", + " d = d + (gamma - np.vdot(Y[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]))*S[:,khat]\n", + " \n", + " return d\n", + " \n", + " def findSearchDirection(self):\n", + " return self.bfgs(self._bfgscnt,-self.g)\n", + " \n", + " def _doEndIteration_BFGS(self, xt):\n", + " if self._iter is 0: \n", + " self.g_last = self.g\n", + " return\n", + " \n", + " yy = self.g - self.g_last;\n", + " ss = self.xc - xt;\n", + " self.g_last = self.g\n", + " \n", + " if yy.dot(ss) > 0:\n", + " self._bfgscnt += 1\n", + " ktop = np.mod(self._bfgscnt,self.nbfgs)\n", + " self._bfgsY[:,ktop] = yy\n", + " self._bfgsS[:,ktop] = ss\n", + " self.comment = ''\n", + " else:\n", + " self.comment = 'Skip BFGS'\n", + " \n", + "\n", + "\n", + "x0 = np.array([1,0])\n", + "opt = BFGS()\n", + "xopt = opt.minimize(FUN,x0)\n", + "print xopt\n", + "opt = inverse.GaussNewton()\n", + "xopt = opt.minimize(FUN,x0)\n", + "print xopt\n", + "opt = inverse.SteepestDescent()\n", + "xopt = opt.minimize(FUN,x0)\n", + "print xopt" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "===================== BFGS =====================\n", + " # f |proj(x-g)-x| LS Comment \n", + "-----------------------------------------------\n", + " 0 1.45e+02 9.51e+01 0 \n", + " 1 1.14e+02 5.37e+01 6 \n", + " 2 1.04e+02 3.04e+01 6 \n", + " 3 8.83e+01 1.37e+01 0 \n", + " 4 8.76e+01 5.97e+00 0 Skip BFGS \n", + " 5 8.74e+01 2.61e+00 0 Skip BFGS \n", + " 6 8.74e+01 1.14e+00 0 Skip BFGS \n", + " 7 8.74e+01 5.01e-01 0 Skip BFGS \n", + " 8 8.74e+01 2.19e-01 0 Skip BFGS \n", + " 9 8.74e+01 9.60e-02 0 Skip BFGS \n", + "------------------------- STOP! -------------------------\n", + "1 : |fc-fOld| = 1.9437e-04 <= tolF*(1+|f0|) = 1.4600e+01\n", + "1 : |xc-x_last| = 1.2663e-03 <= tolX*(1+|x0|) = 2.0000e-01\n", + "1 : |proj(x-g)-x| = 9.5952e-02 <= tolG = 1.0000e-01\n", + "0 : |proj(x-g)-x| = 9.5952e-02 <= 1e3*eps = 1.0000e-02\n", + "0 : maxIter = 20 <= iter = 9\n", + "------------------------- DONE! -------------------------\n", + "[ 0.05095952 4.99977449]\n", + "=========== Gauss Newton ===========\n", + " # f |proj(x-g)-x| LS \n", + "-----------------------------------\n", + " 0 1.45e+02 9.51e+01 0 \n", + " 1 8.74e+01 4.44e-15 0 \n", + "------------------------- STOP! -------------------------\n", + "0 : |fc-fOld| = 5.7625e+01 <= tolF*(1+|f0|) = 1.4600e+01\n", + "0 : |xc-x_last| = 5.0894e+00 <= tolX*(1+|x0|) = 2.0000e-01\n", + "1 : |proj(x-g)-x| = 4.4409e-15 <= tolG = 1.0000e-01\n", + "1 : |proj(x-g)-x| = 4.4409e-15 <= 1e3*eps = 1.0000e-02\n", + "0 : maxIter = 20 <= iter = 1\n", + "------------------------- DONE! -------------------------\n", + "[ 0.05 5. ]\n", + "========= Steepest Descent =========\n", + " # f |proj(x-g)-x| LS \n", + "-----------------------------------\n", + " 0 1.45e+02 9.51e+01 0 \n", + " 1 1.14e+02 5.37e+01 6 \n", + " 2 1.04e+02 3.04e+01 6 \n", + " 3 1.00e+02 1.76e+01 6 \n", + " 4 9.88e+01 1.06e+01 6 \n", + " 5 9.82e+01 7.07e+00 6 \n", + " 6 9.80e+01 1.22e+01 5 \n", + " 7 9.73e+01 7.77e+00 6 \n", + " 8 9.68e+01 5.64e+00 6 \n", + " 9 9.65e+01 8.72e+00 5 \n", + " 10 9.60e+01 5.97e+00 6 \n", + " 11 9.58e+01 9.98e+00 5 \n", + " 12 9.53e+01 6.48e+00 6 \n", + " 13 9.53e+01 1.16e+01 5 \n", + " 14 9.46e+01 7.20e+00 6 \n", + " 15 9.43e+01 5.07e+00 6 \n", + " 16 9.41e+01 8.17e+00 5 \n", + " 17 9.37e+01 5.43e+00 6 \n", + " 18 9.36e+01 9.42e+00 5 \n", + " 19 9.32e+01 5.98e+00 6 \n", + " 20 9.29e+01 4.32e+00 6 \n", + "------------------------- STOP! -------------------------\n", + "1 : |fc-fOld| = 2.5913e-01 <= tolF*(1+|f0|) = 1.4600e+01\n", + "1 : |xc-x_last| = 9.3379e-02 <= tolX*(1+|x0|) = 2.0000e-01\n", + "0 : |proj(x-g)-x| = 4.3246e+00 <= tolG = 1.0000e-01\n", + "0 : |proj(x-g)-x| = 4.3246e+00 <= 1e3*eps = 1.0000e-02\n", + "1 : maxIter = 20 <= iter = 20\n", + "------------------------- DONE! -------------------------\n", + "[ 0.07777107 1.6849632 ]\n" + ] + } + ], + "prompt_number": 14 + } + ], + "metadata": {} + } + ] +} \ No newline at end of file From 4472b39d6e66ced02c7985b164d07ce8874f19e0 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 19 Nov 2013 16:00:36 -0800 Subject: [PATCH 12/18] Iterative Solvers in SimPEG.Solver --- SimPEG/tests/HTMLTestRunner.py | 3 +- SimPEG/tests/TestUtils.py | 8 ++-- SimPEG/utils/Solver.py | 78 ++++++++++++++++++++++++++++------ 3 files changed, 72 insertions(+), 17 deletions(-) diff --git a/SimPEG/tests/HTMLTestRunner.py b/SimPEG/tests/HTMLTestRunner.py index af384971..05ae09df 100644 --- a/SimPEG/tests/HTMLTestRunner.py +++ b/SimPEG/tests/HTMLTestRunner.py @@ -241,7 +241,8 @@ function showClassDetail(cid, count) { for (var i = 0; i < count; i++) { tid = id_list[i]; if (toHide) { - document.getElementById('div_'+tid).style.display = 'none' + var divTid = document.getElementById('div_'+tid); + if(divTid !== null){divTid.style.display = 'none';} document.getElementById(tid).className = 'hiddenRow'; } else { diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index 1cc2bbae..bfb150c3 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -276,16 +276,16 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): -def getQuadratic(A, b): +def getQuadratic(A, b, c=0): """ - Given A and b, this returns a quadratic, Q + Given A, b and c, this returns a quadratic, Q .. math:: - \mathbf{Q( x ) = 0.5 x A x + b x} + \mathbf{Q( x ) = 0.5 x A x + b x} + c """ def Quadratic(x, return_g=True, return_H=True): - f = 0.5 * x.dot( A.dot(x)) + b.dot( x ) + f = 0.5 * x.dot( A.dot(x)) + b.dot( x ) + c out = (f,) if return_g: g = A.dot(x) + b diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index a1194487..82c4b709 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -1,9 +1,10 @@ import numpy as np -import scipy.sparse as sparse +import scipy.sparse as sp import scipy.sparse.linalg as linalg -from SimPEG.utils import mkvc +from SimPEG.utils import mkvc, sdiag +import warnings -DEFAULTS = {'direct':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'} +DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'} try: import TriSolve @@ -45,13 +46,41 @@ class Solver(object): def __init__(self, A, doDirect=True, flag=None, options={}): assert type(doDirect) is bool, 'doDirect must be a boolean' assert flag in [None, 'L', 'U', 'D'], "flag must be set to None, 'L', 'U', or 'D'" - + assert type(options) is dict, 'options must be a dictionary object' self.A = A self.dsolve = None self.doDirect = doDirect self.flag = flag self.options = options + if doDirect: return + + # Now deal with iterative stuff only + if 'M' not in options: + warnings.warn("You should provide a preconditioner, M.", UserWarning) + return + M = options['M'] + if type(M) is sp.linalg.LinearOperator: + return + elif type(M) is tuple: + PreconditionerList = ['J','GS'] + assert type(M[0]) is str and M[0] in PreconditionerList, "M as a tuple must be (str, Matrix) where str is in ['J','GS']: e.g. ('J', WtW) where J stands for Jacobi, and WtW is a sparse matrix." + if M[0] is 'J': + Jacobi = sdiag(1.0/M[1].diagonal()) + options['M'] = Jacobi + elif M[0] is 'GS': + LL = sp.tril(M[1]) + UU = sp.triu(M[1]) + DD = sdiag(M[1].diagonal()) + Uinv = Solver(UU, flag='U') + Linv = Solver(LL, flag='L') + def GS(f): + return Uinv.solve(DD*Linv.solve(f)) + options['M'] = sp.linalg.LinearOperator( A.shape, GS, dtype=A.dtype ) + + else: + raise Exception('M must be a LinearOperator or a tuple') + def solve(self, b): """ @@ -118,8 +147,20 @@ class Solver(object): return X - def solveIter(self, b, M=None, iterSolver='CG'): - pass + def solveIter(self, b, backend=None, M=None, iterSolver='CG', tol=1e-6, maxIter=50): + if backend is None: backend = DEFAULTS['iter'] + + algorithms = {'CG':sp.linalg.cg} + assert iterSolver in algorithms, "iterSolver must be 'CG', or implement it yourself and add it here!" + alg = algorithms[iterSolver] + + if len(b.shape) == 1 or b.shape[1] == 1: + x, self.info = alg(self.A, b, M=M, tol=tol, maxiter=maxIter) + else: + x = np.empty_like(b) + for i in range(b.shape[1]): + x[:,i], self.info = alg(self.A, b[:,i], M=M, tol=tol, maxiter=maxIter) + return x def solveBackward(self, b, backend=None): """ @@ -132,9 +173,8 @@ class Solver(object): :return: x """ if backend is None: backend = DEFAULTS['backward'] - if type(self.A) is not sparse.csr.csr_matrix: - from scipy.sparse import csr_matrix - self.A = csr_matrix(self.A) + if type(self.A) is not sp.csr.csr_matrix: + self.A = sp.csr_matrix(self.A) vals = self.A.data rowptr = self.A.indptr colind = self.A.indices @@ -164,7 +204,7 @@ class Solver(object): :return: x """ if backend is None: backend = DEFAULTS['forward'] - if type(self.A) is not sparse.csr.csr_matrix: + if type(self.A) is not sp.csr.csr_matrix: from scipy.sparse import csr_matrix self.A = csr_matrix(self.A) vals = self.A.data @@ -240,13 +280,13 @@ if __name__ == '__main__': print np.linalg.norm(e-x,np.inf) - n = 6000 + n = 600 A_dense = np.random.random((n,n)) L = np.tril(np.dot(A_dense, A_dense)) # Positive definite is better conditioned. e = np.ones(n) b = np.dot(L, e) - A = sparse.csr_matrix(L) + A = sp.csr_matrix(L) pSolve = Solver(A,flag='L',options={'backend':'python'}); fSolve = Solver(A,flag='L',options={'backend':'fortran'}) tic = time() @@ -257,3 +297,17 @@ if __name__ == '__main__': x = fSolve.solve(b) toc = time() - tic print 'Error Forward Fortran = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc + + + + A = -D*D.T + A[0,0] *= 10 # remove the constant null space from the matrix + e = np.ones(M.nC) + b = A.dot(e) + + iSolve = Solver(A, doDirect=False,options={'M':('GS',A)}) + tic = time() + x = iSolve.solve(b) + toc = time() - tic + print x + print 'Error CG = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc, 'Info: ', iSolve.info From 500844ccb5f182f987355ef212b759be49f66747 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 19 Nov 2013 17:41:51 -0800 Subject: [PATCH 13/18] Brought BFGS into Optimize. updated InexactGaussNewton to use BFGS as a preconditioner. --- SimPEG/inverse/Optimize.py | 114 +++++++++++++++++++++++++++++++++++-- SimPEG/utils/Solver.py | 7 ++- notebooks/BFGS.ipynb | 105 +++++++++++----------------------- 3 files changed, 147 insertions(+), 79 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 8a9abaed..ff38d67b 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -588,7 +588,7 @@ class ProjectedGradient(Minimize, Remember): def reduceHess(v): # Z is tall and skinny return Z.T*(self.H*(Z*v)) - operator = sp.linalg.LinearOperator( (shape[1], shape[1]), reduceHess, dtype=float ) + operator = sp.linalg.LinearOperator( (shape[1], shape[1]), reduceHess, dtype=self.xc.dtype ) p, info = sp.linalg.cg(operator, -Z.T*self.g, tol=self.tolCG, maxiter=self.maxIterCG) p = Z*p # bring up to full size # aSet_after = self.activeSet(self.xc+p) @@ -622,21 +622,125 @@ class ProjectedGradient(Minimize, Remember): if self.debug: print 'doEndIteration.ProjGrad, f_decrease_max: ', self.f_decrease_max if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD + + +class BFGS(Minimize, Remember): + name = 'BFGS' + nbfgs = 10 + + @property + def bfgsH0(self): + """ + Approximate Hessian used in preconditioning the problem. + + Must be a SimPEG.Solver + """ + _bfgsH0 = getattr(self,'_bfgsH0',None) + if _bfgsH0 is None: + return Solver(sp.identity(self.xc.size).tocsc(), flag='D') + return _bfgsH0 + @bfgsH0.setter + def bfgsH0(self, value): + assert type(value) is Solver, 'bfgsH0 must be a SimPEG.Solver' + self._bfgsH0 = value + + def _startup_BFGS(self,x0): + self._bfgscnt = -1 + self._bfgsY = np.zeros((x0.size, self.nbfgs)) + self._bfgsS = np.zeros((x0.size, self.nbfgs)) + if not np.any([p is IterationPrinters.comment for p in self.printers]): + self.printers.append(IterationPrinters.comment) + + def bfgs(self, d): + n = self._bfgscnt + nn = ktop = min(self._bfgsS.shape[1],n) + return self.bfgsrec(ktop,n,nn,self._bfgsS,self._bfgsY,d) + + def bfgsrec(self,k,n,nn,S,Y,d): + """BFGS recursion""" + if k < 0: + d = self.bfgsH0.solve(d) + else: + khat = np.mod(n-nn+k,nn) + gamma = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]) + d = d - gamma*Y[:,khat] + d = self.bfgsrec(k-1,n,nn,S,Y,d) + d = d + (gamma - np.vdot(Y[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]))*S[:,khat] + return d + + def findSearchDirection(self): + return self.bfgs(-self.g) + + def _doEndIteration_BFGS(self, xt): + if self._iter is 0: + self.g_last = self.g + return + + yy = self.g - self.g_last; + ss = self.xc - xt; + self.g_last = self.g + + if yy.dot(ss) > 0: + self._bfgscnt += 1 + ktop = np.mod(self._bfgscnt,self.nbfgs) + self._bfgsY[:,ktop] = yy + self._bfgsS[:,ktop] = ss + self.comment = '' + else: + self.comment = 'Skip BFGS' + + class GaussNewton(Minimize, Remember): name = 'Gauss Newton' def findSearchDirection(self): return Solver(self.H).solve(-self.g) -class InexactGaussNewton(Minimize, Remember): +class InexactGaussNewton(BFGS, Minimize, Remember): + """ + Minimizes using CG as the inexact solver of + + .. math:: + + \mathbf{H p = -g} + + By default BFGS is used as the preconditioner. + + Use *nbfgs* to set the memory limitation of BFGS. + + To set the initial H0 to be used in BFGS, set *bfgsH0* to be a SimPEG.Solver + + """ + + def __init__(self, **kwargs): + Minimize.__init__(self, **kwargs) + name = 'Inexact Gauss Newton' maxIterCG = 10 - tolCG = 1e-5 + tolCG = 1e-3 + + @property + def approxHinv(self): + """ + The approximate Hessian inverse is used to precondition CG. + + Default uses BFGS, with an initial H0 of *bfgsH0*. + + Must be a scipy.sparse.linalg.LinearOperator + """ + _approxHinv = getattr(self,'_approxHinv',None) + if _approxHinv is None: + M = sp.linalg.LinearOperator( (self.xc.size, self.xc.size), self.bfgs, dtype=self.xc.dtype ) + return M + return _approxHinv + @approxHinv.setter + def approxHinv(self, value): + self._approxHinv = value def findSearchDirection(self): - # TODO: use BFGS as a preconditioner or gauss sidel of the WtW or solve WtW directly - p, info = sp.linalg.cg(self.H, -self.g, tol=self.tolCG, maxiter=self.maxIterCG) + Hinv = Solver(self.H, doDirect=False, options={'iterSolver': 'CG', 'M': self.approxHinv, 'tol': self.tolCG, 'maxIter': self.maxIterCG}) + p = Hinv.solve(-self.g) return p diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py index 82c4b709..16db872f 100644 --- a/SimPEG/utils/Solver.py +++ b/SimPEG/utils/Solver.py @@ -62,8 +62,11 @@ class Solver(object): M = options['M'] if type(M) is sp.linalg.LinearOperator: return - elif type(M) is tuple: - PreconditionerList = ['J','GS'] + PreconditionerList = ['J','GS'] + if type(M) is str: + assert M in PreconditionerList, "M must be in the known preconditioner list. ['J','GS']" + M = (M,A) # use A as the base for the preconditioner. + if type(M) is tuple: assert type(M[0]) is str and M[0] in PreconditionerList, "M as a tuple must be (str, Matrix) where str is in ['J','GS']: e.g. ('J', WtW) where J stands for Jacobi, and WtW is a sparse matrix." if M[0] is 'J': Jacobi = sdiag(1.0/M[1].diagonal()) diff --git a/notebooks/BFGS.ipynb b/notebooks/BFGS.ipynb index 23de6e06..c2f68c99 100644 --- a/notebooks/BFGS.ipynb +++ b/notebooks/BFGS.ipynb @@ -23,7 +23,7 @@ "language": "python", "metadata": {}, "outputs": [], - "prompt_number": 1 + "prompt_number": 2 }, { "cell_type": "code", @@ -32,77 +32,8 @@ "FUN = SimPEG.tests.Rosenbrock\n", "FUN = SimPEG.tests.getQuadratic(sp.csr_matrix(([100,1],([0,1],[0,1])),shape=(2,2)),np.array([-5,-5]),100)\n", "\n", - "\n", - "class BFGS(Minimize, Remember):\n", - " name = 'BFGS'\n", - " nbfgs = 10\n", - " \n", - " @property\n", - " def H0(self):\n", - " \"\"\"\n", - " Approximate Hessian used in preconditioning the problem.\n", - " \n", - " Must be a SimPEG.Solver\n", - " \"\"\"\n", - " _H0 = getattr(self,'_H0',None)\n", - " if _H0 is None:\n", - " return Solver(sp.identity(self.xc.size).tocsc())\n", - " return _H0\n", - " @H0.setter\n", - " def H0(self, value):\n", - " self._H0 = value\n", - " \n", - " def _startup_BFGS(self,x0):\n", - " self._bfgscnt = -1\n", - " self._bfgsY = np.zeros((x0.size, self.nbfgs))\n", - " self._bfgsS = np.zeros((x0.size, self.nbfgs))\n", - " if not np.any([p is IterationPrinters.comment for p in self.printers]):\n", - " self.printers.append(IterationPrinters.comment)\n", - " \n", - " def bfgs(self,n,d):\n", - " nn = min(self._bfgsS.shape[1],n)\n", - " ktop = nn\n", - " d = self.bfgsrec(ktop,n,nn,self._bfgsS,self._bfgsY,d)\n", - " return d\n", - "\n", - " def bfgsrec(self,k,n,nn,S,Y,d):\n", - " \"\"\"BFGS recursion\"\"\"\n", - " if k < 0:\n", - " d = self.H0.solve(d)\n", - " else:\n", - " khat = mod(n-nn+k,nn)\n", - " gamma = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat])\n", - " d = d - gamma*Y[:,khat]\n", - " d = self.bfgsrec(k-1,n,nn,S,Y,d)\n", - " d = d + (gamma - np.vdot(Y[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]))*S[:,khat]\n", - " \n", - " return d\n", - " \n", - " def findSearchDirection(self):\n", - " return self.bfgs(self._bfgscnt,-self.g)\n", - " \n", - " def _doEndIteration_BFGS(self, xt):\n", - " if self._iter is 0: \n", - " self.g_last = self.g\n", - " return\n", - " \n", - " yy = self.g - self.g_last;\n", - " ss = self.xc - xt;\n", - " self.g_last = self.g\n", - " \n", - " if yy.dot(ss) > 0:\n", - " self._bfgscnt += 1\n", - " ktop = np.mod(self._bfgscnt,self.nbfgs)\n", - " self._bfgsY[:,ktop] = yy\n", - " self._bfgsS[:,ktop] = ss\n", - " self.comment = ''\n", - " else:\n", - " self.comment = 'Skip BFGS'\n", - " \n", - "\n", - "\n", "x0 = np.array([1,0])\n", - "opt = BFGS()\n", + "opt = inverse.BFGS()\n", "xopt = opt.minimize(FUN,x0)\n", "print xopt\n", "opt = inverse.GaussNewton()\n", @@ -186,9 +117,39 @@ "------------------------- DONE! -------------------------\n", "[ 0.07777107 1.6849632 ]\n" ] + }, + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "/Users/rowan/git/simpeg/SimPEG/inverse/Optimize.py:664: RuntimeWarning: divide by zero encountered in remainder\n", + " khat = np.mod(n-nn+k,nn)\n" + ] } ], - "prompt_number": 14 + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "A = sp.identity(2)\n", + "S = Solver(A)\n", + "\n", + "assert type(S) is Solver" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 6 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] } ], "metadata": {} From 0afe42aa34fb82f14376354f4171be518620ce62 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 19 Nov 2013 17:49:18 -0800 Subject: [PATCH 14/18] Bug fix in timeStep (must return something, not an error.) Minor updates to Optimize and Inversion. --- SimPEG/forward/Richards.py | 4 ++-- SimPEG/inverse/Inversion.py | 2 ++ SimPEG/inverse/Optimize.py | 9 ++++++++- 3 files changed, 12 insertions(+), 3 deletions(-) diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py index 99090d84..b23680d7 100644 --- a/SimPEG/forward/Richards.py +++ b/SimPEG/forward/Richards.py @@ -15,10 +15,10 @@ class RichardsProblem(Problem): @property def timeStep(self): """The time between steps.""" - return self._timeStep + return getattr(self, '_timeStep', None) @timeStep.setter def timeStep(self, value): - self._timeStep = float(value) + self._timeStep = float(value) # Because integers suck. @property def numIts(self): diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 871a7f06..9ab84fbd 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -93,6 +93,8 @@ class BaseInversion(object): self.m = m0 self._iter = 0 self._beta = None + self.phi_d_last = np.nan + self.phi_m_last = np.nan def doEndIteration(self): """ diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index ca2061f5..bf8b63be 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -189,7 +189,9 @@ class Minimize(object): self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) if doPub: pub.sendMessage('Minimize.evalFunction', minimize=self, f=self.f, g=self.g, H=self.H) self.printIter() - if self.stoppingCriteria(): break + if self.stoppingCriteria(): + self.doEndIteration(self.xc) + break p = self.findSearchDirection() if doPub: pub.sendMessage('Minimize.searchDirection', minimize=self, p=p) p = self.scaleSearchDirection(p) @@ -273,6 +275,11 @@ class Minimize(object): parent.printIter function and call that. """ + + for method in [posible for posible in dir(self) if '_printIter' in posible]: + if self.debug: print 'printIter is calling self.'+method + getattr(self,method)(inLS) + if doPub and not inLS: pub.sendMessage('Minimize.printIter', minimize=self) pad = ' '*10 if inLS else '' printLine(self, self.printers if not inLS else self.printersLS, pad=pad) From 1c895f397d014677db5c8f45d13c14b31f08e2bc Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 20 Nov 2013 19:37:09 -0800 Subject: [PATCH 15/18] Set Mref in inversion. --- SimPEG/inverse/Inversion.py | 11 ++++++++++- SimPEG/regularization/Regularization.py | 5 ++--- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 9ab84fbd..e0b62866 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -29,6 +29,11 @@ class BaseInversion(object): self.opt.printers.insert(3,SimPEG.inverse.IterationPrinters.phi_m) self.opt.stoppers.append(SimPEG.inverse.StoppingCriteria.phi_d_target_Minimize) + if not hasattr(opt, '_bfgsH0'): # Check if it has been set by the user and the default is not being used. + print 'Setting bfgsH0 to the inverse of the modelObj2Deriv.' + opt.bfgsH0 = SimPEG.Solver(reg.modelObj2Deriv(),doDirect=True,options={'factorize':True}) # False, options={'M':'GS','maxIter':15} + + @property def Wd(self): """ @@ -90,6 +95,10 @@ class BaseInversion(object): if self.debug: print 'startup is calling self.'+method getattr(self,method)(m0) + if not hasattr(self.reg, '_mref'): + print 'Regularization has not set mref. SimPEG will set it to m0.' + self.reg.mref = m0 + self.m = m0 self._iter = 0 self._beta = None @@ -159,7 +168,7 @@ class BaseInversion(object): if return_H: def H_fun(v): phi_d2Deriv = self.dataObj2Deriv(m, v, u=u) - phi_m2Deriv = self.reg.modelObj2Deriv(m)*v + phi_m2Deriv = self.reg.modelObj2Deriv()*v return phi_d2Deriv + self._beta * phi_m2Deriv diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index 6f5970e6..72335661 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/regularization/Regularization.py @@ -7,7 +7,7 @@ class Regularization(object): @property def mref(self): if getattr(self, '_mref', None) is None: - self._mref = np.zeros(self.mesh.nC); + return np.zeros(self.mesh.nC); return self._mref @mref.setter def mref(self, value): @@ -104,8 +104,7 @@ class Regularization(object): return mobjDeriv - def modelObj2Deriv(self, m): - mresid = m - self.mref + def modelObj2Deriv(self): mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws From 808036ddc8c5997b4d9ebdfe6377a8b8b6f4837b Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 20 Nov 2013 19:37:48 -0800 Subject: [PATCH 16/18] Bug fixes in optimization. --- SimPEG/inverse/Optimize.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index bf8b63be..c2f95dd9 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -189,9 +189,7 @@ class Minimize(object): self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) if doPub: pub.sendMessage('Minimize.evalFunction', minimize=self, f=self.f, g=self.g, H=self.H) self.printIter() - if self.stoppingCriteria(): - self.doEndIteration(self.xc) - break + if self.stoppingCriteria(): break p = self.findSearchDirection() if doPub: pub.sendMessage('Minimize.searchDirection', minimize=self, p=p) p = self.scaleSearchDirection(p) @@ -668,7 +666,7 @@ class BFGS(Minimize, Remember): if k < 0: d = self.bfgsH0.solve(d) else: - khat = np.mod(n-nn+k,nn) + khat = 0 if nn is 0 else np.mod(n-nn+k,nn) gamma = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]) d = d - gamma*Y[:,khat] d = self.bfgsrec(k-1,n,nn,S,Y,d) From 95fd98c23eb5c49539af2e08b81ca92361e22602 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 21 Nov 2013 11:52:41 -0800 Subject: [PATCH 17/18] Remove Richards problem from SimPEG. This will be held in a different repo. --- SimPEG/forward/Richards.py | 440 ---------------------------------- SimPEG/tests/test_Richards.py | 134 ----------- docs/api_Problem.rst | 6 - 3 files changed, 580 deletions(-) delete mode 100644 SimPEG/forward/Richards.py delete mode 100644 SimPEG/tests/test_Richards.py diff --git a/SimPEG/forward/Richards.py b/SimPEG/forward/Richards.py deleted file mode 100644 index b23680d7..00000000 --- a/SimPEG/forward/Richards.py +++ /dev/null @@ -1,440 +0,0 @@ -from SimPEG.forward import Problem -import numpy as np -from SimPEG.utils import sdiag, spzeros, mkvc, setKwargs, Solver -from SimPEG.inverse import NewtonRoot -import scipy.sparse as sp - - -class RichardsProblem(Problem): - """docstring for RichardsProblem""" - - timeEnd = None - boundaryConditions = None - initialConditions = None - - @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) - - _method = 'mixed' - @property - def method(self): - """ - - Method must be either 'mixed' or 'head'. - - There are two different forms of Richards equation that differ - on how they deal with the non-linearity in the time-stepping term. - - The most fundamental form, referred to as the - 'mixed'-form of Richards Equation [Celia et al., 1990] - - .. math:: - - \\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. - This formulation of Richards equation is called the - 'mixed'-form because the equation is parameterized in psi - but the time-stepping is in terms of theta. - - As noted in [Celia et al., 1990] the 'head'-based form of Richards - 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 - \quad \psi \in \Omega - - However, it can be shown that this does not conserve mass in the discrete formulation. - - - """ - return self._method - @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): - assert type(value) is bool, 'doNewton must be a boolean.' - self.rootFinder = NewtonRoot(doLS=value) - self._doNewton = value - - @property - def dataType(self): - """Choose how your data is collected, must be 'saturation' or 'pressureHead'.""" - return self._dataType - @dataType.setter - def dataType(self, value): - assert value in ['saturation','pressureHead'], "dataType must be 'saturation' or 'pressureHead'." - self._dataType = value - - - def __init__(self, mesh, empirical, **kwargs): - Problem.__init__(self, mesh) - self.empirical = empirical - self.mesh.setCellGradBC('dirichlet') - self.dataType = 'pressureHead' - self.doNewton = False # This also sets the rootFinder algorithm. - setKwargs(self, **kwargs) - - def dpred(self, m, u=None): - """ - Predicted data. - - .. math:: - d_\\text{pred} = Pu(m) - """ - if u is None: - u = self.field(m) - u = np.concatenate(u[1:]) - if self.dataType is 'saturation': - u = self.empirical.moistureContent(u) - return self.P*u - - def field(self, m): - self.empirical.setModel(m) - Hs = range(self.numIts+1) - Hs[0] = self.initialConditions - for ii in range(self.numIts): - Hs[ii+1] = self.rootFinder.root(lambda hn1: self.getResidual(Hs[ii],hn1), Hs[ii]) - return Hs - - def diagsJacobian(self, hn, hn1): - - DIV = self.mesh.faceDiv - GRAD = self.mesh.cellGrad - BC = self.mesh.cellGradBC - AV = self.mesh.aveCC2F - Dz = self.mesh.faceDiv #TODO: fix this for more than one dimension. - - bc = self.boundaryConditions - dt = self.timeStep - - dT = self.empirical.moistureContentDeriv(hn) - dT1 = self.empirical.moistureContentDeriv(hn1) - K1 = self.empirical.hydraulicConductivity(hn1) - dK1 = self.empirical.hydraulicConductivityDeriv(hn1) - dKa1 = self.empirical.hydraulicConductivityModelDeriv(hn1) - - # Compute part of the derivative of: - # - # DIV*diag(GRAD*hn1+BC*bc)*(AV*(1.0/K))^-1 - - DdiagGh1 = DIV*sdiag(GRAD*hn1+BC*bc) - diagAVk2_AVdiagK2 = sdiag((AV*(1./K1))**(-2)) * AV*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*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, hn, h): - """ - 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 - Dz = self.mesh.faceDiv #TODO: fix this for more than one dimension. - - bc = self.boundaryConditions - dt = self.timeStep - - T = self.empirical.moistureContent(h) - dT = self.empirical.moistureContentDeriv(h) - Tn = self.empirical.moistureContent(hn) - K = self.empirical.hydraulicConductivity(h) - dK = self.empirical.hydraulicConductivityDeriv(h) - - aveK = 1./(AV*(1./K)); - - RHS = DIV*sdiag(aveK)*(GRAD*h+BC*bc) + Dz*aveK - if self.method is 'mixed': - r = (T-Tn)/dt - RHS - elif self.method is 'head': - r = dT*(h - hn)/dt - RHS - - J = dT/dt - DIV*sdiag(aveK)*GRAD - if self.doNewton: - DDharmAve = sdiag(aveK**2)*AV*sdiag(K**(-2)) * dK - J = J - DIV*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(Hs[ii],Hs[ii+1]) - Ad = sp.block_diag(Adiags) - zRight = spzeros((len(Asubs)-1)*Asubs[0].shape[0],Adiags[0].shape[1]) - zTop = 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 J(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(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(Hs[ii],Hs[ii+1]) - Adiaginv = Solver(Adiag) - JvC[ii] = Adiaginv.solve(B*v - Asub*JvC[ii-1]) - - if self.dataType is 'pressureHead': - Jv = self.P*np.concatenate(JvC) - elif self.dataType is 'saturation': - dT = self.empirical.moistureContentDeriv(np.concatenate(Hs[1:])) - Jv = self.P*dT*np.concatenate(JvC) - - return Jv - - def Jt(self, m, v, u=None): - if u is None: - u = self.field(m) - Hs = u - - if self.dataType is 'pressureHead': - PTv = self.P.T*v; - elif self.dataType is 'saturation': - dT = self.empirical.moistureContentDeriv(np.concatenate(Hs[1:])) - 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(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 - - -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): - 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 = 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 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 = 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): - 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 = sdiag(g) - return g - - def hydraulicConductivity(self, h): - alpha = self.alpha - I = self.I - n = self.n - Ks = self.Ks - m = 1 - 1.0/n - - theta_e = 1.0/((1+abs(alpha*h)**n)**m) - f = np.exp(Ks)*theta_e**I* ( ( 1 - ( 1 - 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 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 - 1.0/n - - g = I*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 - 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 = sdiag(g) - return g - - -if __name__ == '__main__': - import SimPEG - from SimPEG import mesh, inverse, regularization, tests - import scipy.sparse as sp - import numpy as np - from SimPEG.forward import Problem, Richards - 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] - prob = Richards.RichardsProblem(M,E, timeStep=10, timeEnd=60, 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.P = P - - prob.dataType = 'pressureHead' - mTrue = np.ones(M.nC)*np.log(Ks) - stdev = 0.01 # The standard deviation for the noise - dobs = prob.createSyntheticData(mTrue,std=stdev)[0] - # p = plot(dobs.reshape((-1,4))) - prob.dobs = dobs - prob.std = dobs*0 + stdev - opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) - reg = regularization.Regularization(mesh) - inv = inverse.Inversion(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/SimPEG/tests/test_Richards.py b/SimPEG/tests/test_Richards.py deleted file mode 100644 index 4510aab0..00000000 --- a/SimPEG/tests/test_Richards.py +++ /dev/null @@ -1,134 +0,0 @@ -import numpy as np -import scipy.sparse as sp -import unittest -from SimPEG import mesh, regularization, inverse -from TestUtils import OrderTest, checkDerivative -from scipy.sparse.linalg import dsolve -from SimPEG.forward import Richards - - -TOL = 1E-8 - -class RichardsTests(unittest.TestCase): - - def setUp(self): - 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] - prob = Richards.RichardsProblem(M,E, timeStep=30, timeEnd=360, 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.P = P - - - self.h0 = h - self.M = M - self.Ks = Ks - self.prob = prob - - def test_VanGenuchten_moistureContent(self): - 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): - 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): - 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): - 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): - 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): - 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_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(mesh) - 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() diff --git a/docs/api_Problem.rst b/docs/api_Problem.rst index 7fea7b5c..3d5a6d2c 100644 --- a/docs/api_Problem.rst +++ b/docs/api_Problem.rst @@ -26,10 +26,4 @@ Linear Problem :members: :undoc-members: -Richards Problem -**************** - -.. automodule:: SimPEG.forward.Richards - :members: - :undoc-members: From ef45abddbb0f0bec93f31028af16c4edea47b21e Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 21 Nov 2013 12:03:25 -0800 Subject: [PATCH 18/18] Test results and bug fix. --- SimPEG/forward/__init__.py | 1 - docs/api_TestResults.rst | 1681 ++++++++++++++++++++++-------------- 2 files changed, 1024 insertions(+), 658 deletions(-) diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py index ec933f0a..33c9a6b1 100644 --- a/SimPEG/forward/__init__.py +++ b/SimPEG/forward/__init__.py @@ -2,4 +2,3 @@ from Problem import * import DCProblem from LinearProblem import LinearProblem import ModelTransforms -import Richards diff --git a/docs/api_TestResults.rst b/docs/api_TestResults.rst index 9e8b2df3..fed77e4c 100644 --- a/docs/api_TestResults.rst +++ b/docs/api_TestResults.rst @@ -137,7 +137,8 @@ Test Results for (var i = 0; i < count; i++) { tid = id_list[i]; if (toHide) { - document.getElementById('div_'+tid).style.display = 'none' + var divTid = document.getElementById('div_'+tid); + if(divTid !== null){divTid.style.display = 'none';} document.getElementById(tid).className = 'hiddenRow'; } else { @@ -185,9 +186,9 @@ Test Results -->
-

Start Time: 2013-11-12 13:34:53

-

Duration: 0:00:31.732498

-

Status: Pass 117

+

Start Time: 2013-11-21 11:59:36

+

Duration: 0:00:32.074343

+

Status: Pass 124

SimPEG Test Report was automatically generated.

@@ -356,7 +357,27 @@ Test Results
test_adjoint
- pass + + + + + pass + + + + + @@ -374,18 +395,19 @@ Test Results 
             
-    pt3.2: ==================== checkDerivative ====================
+    pt3.2: Setting bfgsH0 to the inverse of the modelObj2Deriv.
+    ==================== checkDerivative ====================
     iter	h		|J0-Jt|		|J0+h*dJ'*dx-Jt|	Order
     ---------------------------------------------------------
-    0	1.00e-01	5.517e+01		4.463e+01		nan
-    1	1.00e-02	1.516e+00		4.609e-01		1.986
-    2	1.00e-03	1.101e-01		4.624e-03		1.999
-    3	1.00e-04	1.059e-02		4.626e-05		2.000
-    4	1.00e-05	1.055e-03		4.626e-07		2.000
-    5	1.00e-06	1.055e-04		4.623e-09		2.000
-    6	1.00e-07	1.055e-05		4.504e-11		2.011
+    0	1.00e-01	6.680e+01		6.164e+01		nan
+    1	1.00e-02	1.154e+00		6.381e-01		1.985
+    2	1.00e-03	5.798e-02		6.401e-03		1.999
+    3	1.00e-04	5.222e-03		6.403e-05		2.000
+    4	1.00e-05	5.164e-04		6.403e-07		2.000
+    5	1.00e-06	5.159e-05		6.403e-09		2.000
+    6	1.00e-07	5.158e-06		6.377e-11		2.002
     ========================= PASS! =========================
-    Awesome, Rowan, just awesome.
+    You get a gold star!
     
     
     
@@ -411,18 +433,19 @@ Test Results
             
             
             
-    pt3.3: ==================== checkDerivative ====================
+    pt3.3: Setting bfgsH0 to the inverse of the modelObj2Deriv.
+    ==================== checkDerivative ====================
     iter	h		|J0-Jt|		|J0+h*dJ'*dx-Jt|	Order
     ---------------------------------------------------------
-    0	1.00e-01	1.180e-01		6.488e-03		nan
-    1	1.00e-02	1.175e-02		6.401e-05		2.006
-    2	1.00e-03	1.175e-03		6.390e-07		2.001
-    3	1.00e-04	1.175e-04		6.389e-09		2.000
-    4	1.00e-05	1.175e-05		6.389e-11		2.000
-    5	1.00e-06	1.175e-06		6.392e-13		2.000
-    6	1.00e-07	1.175e-07		7.710e-15		1.919
+    0	1.00e-01	1.203e-01		1.785e-02		nan
+    1	1.00e-02	1.175e-02		1.805e-04		1.995
+    2	1.00e-03	1.173e-03		1.805e-06		2.000
+    3	1.00e-04	1.173e-04		1.805e-08		2.000
+    4	1.00e-05	1.173e-05		1.805e-10		2.000
+    5	1.00e-06	1.173e-06		1.805e-12		2.000
+    6	1.00e-07	1.173e-07		1.915e-14		1.975
     ========================= PASS! =========================
-    Go Test Go!
+    The test be workin!
     
     
     
@@ -448,18 +471,19 @@ Test Results
             
             
             
-    pt3.4: ==================== checkDerivative ====================
+    pt3.4: Setting bfgsH0 to the inverse of the modelObj2Deriv.
+    ==================== checkDerivative ====================
     iter	h		|J0-Jt|		|J0+h*dJ'*dx-Jt|	Order
     ---------------------------------------------------------
-    0	1.00e-01	6.898e+00		7.074e+00		nan
-    1	1.00e-02	5.316e-02		7.074e-02		2.000
-    2	1.00e-03	1.051e-03		7.074e-04		2.000
-    3	1.00e-04	1.688e-04		7.074e-06		2.000
-    4	1.00e-05	1.752e-05		7.074e-08		2.000
-    5	1.00e-06	1.758e-06		7.074e-10		2.000
-    6	1.00e-07	1.759e-07		7.076e-12		2.000
+    0	1.00e-01	6.895e+00		6.894e+00		nan
+    1	1.00e-02	6.903e-02		6.894e-02		2.000
+    2	1.00e-03	6.982e-04		6.894e-04		2.000
+    3	1.00e-04	7.769e-06		6.894e-06		2.000
+    4	1.00e-05	1.564e-07		6.894e-08		2.000
+    5	1.00e-06	9.431e-09		6.894e-10		2.000
+    6	1.00e-07	8.810e-10		6.895e-12		2.000
     ========================= PASS! =========================
-    Testing is important.
+    Go Test Go!
     
     
     
@@ -498,15 +522,15 @@ Test Results
     ==================== checkDerivative ====================
     iter	h		|J0-Jt|		|J0+h*dJ'*dx-Jt|	Order
     ---------------------------------------------------------
-    0	1.00e-01	2.762e-01		4.105e-17		nan
-    1	1.00e-02	2.762e-02		2.532e-17		0.210
-    2	1.00e-03	2.762e-03		4.066e-17		-0.206
-    3	1.00e-04	2.762e-04		4.297e-17		-0.024
-    4	1.00e-05	2.762e-05		5.494e-17		-0.107
-    5	1.00e-06	2.762e-06		4.696e-17		0.068
-    6	1.00e-07	2.762e-07		2.834e-17		0.219
+    0	1.00e-01	1.936e-01		7.633e-17		nan
+    1	1.00e-02	1.936e-02		7.268e-17		0.021
+    2	1.00e-03	1.936e-03		6.545e-17		0.046
+    3	1.00e-04	1.936e-04		8.778e-17		-0.127
+    4	1.00e-05	1.936e-05		6.714e-17		0.116
+    5	1.00e-06	1.936e-06		5.514e-17		0.086
+    6	1.00e-07	1.936e-07		8.359e-17		-0.181
     ========================= PASS! =========================
-    You deserve a pat on the back!
+    The test be workin!
     
     
     
@@ -535,15 +559,15 @@ Test Results
     pt4.2: ==================== checkDerivative ====================
     iter	h		|J0-Jt|		|J0+h*dJ'*dx-Jt|	Order
     ---------------------------------------------------------
-    0	1.00e-01	5.835e-01		2.875e-01		nan
-    1	1.00e-02	3.247e-02		2.875e-03		2.000
-    2	1.00e-03	2.989e-03		2.875e-05		2.000
-    3	1.00e-04	2.963e-04		2.875e-07		2.000
-    4	1.00e-05	2.960e-05		2.875e-09		2.000
-    5	1.00e-06	2.960e-06		2.875e-11		2.000
-    6	1.00e-07	2.960e-07		2.879e-13		1.999
+    0	1.00e-01	3.925e-01		6.586e-02		nan
+    1	1.00e-02	4.517e-02		6.586e-04		2.000
+    2	1.00e-03	4.577e-03		6.586e-06		2.000
+    3	1.00e-04	4.582e-04		6.586e-08		2.000
+    4	1.00e-05	4.583e-05		6.586e-10		2.000
+    5	1.00e-06	4.583e-06		6.587e-12		2.000
+    6	1.00e-07	4.583e-07		6.827e-14		1.984
     ========================= PASS! =========================
-    Awesome, Rowan, just awesome.
+    Not just a pretty face Rowan
     
     
     
@@ -583,22 +607,22 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  2.14e-01   |
-      16  |  6.31e-02   |   3.3843    |  1.7589
-      32  |  1.72e-02   |   3.6712    |  1.8763
+       8  |  2.91e-01   |
+      16  |  6.70e-02   |   4.3512    |  2.1214
+      32  |  1.48e-02   |   4.5346    |  2.1810
     ---------------------------------------------
-    Well done Rowan!
+    Not just a pretty face Rowan
     
     
     randomTensorMesh:  Interpolation 1D: CC
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  3.35e-01   |
-      16  |  1.18e-01   |   2.8469    |  1.8016
-      32  |  2.27e-02   |   5.1806    |  3.9741
+       8  |  4.81e-01   |
+      16  |  9.23e-02   |   5.2190    |  1.9811
+      32  |  1.61e-02   |   5.7473    |  3.1176
     ---------------------------------------------
-    Happy little convergence test!
+    You get a gold star!
     
     
     
@@ -629,22 +653,22 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  2.98e-01   |
-      16  |  6.69e-02   |   4.4487    |  2.1534
-      32  |  1.55e-02   |   4.3044    |  2.1058
+       8  |  2.46e-01   |
+      16  |  6.96e-02   |   3.5355    |  1.8219
+      32  |  1.88e-02   |   3.6930    |  1.8848
     ---------------------------------------------
-    Well done Rowan!
+    Awesome, Rowan, just awesome.
     
     
     randomTensorMesh:  Interpolation 1D: N
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.16e-01   |
-      16  |  1.60e-01   |   4.4731    |  1.8851
-      32  |  3.72e-02   |   4.3056    |  2.5981
+       8  |  6.95e-01   |
+      16  |  1.12e-01   |   6.2065    |  2.7830
+      32  |  4.60e-02   |   2.4348    |  2.0359
     ---------------------------------------------
-    Not just a pretty face Rowan
+    The test be workin!
     
     
     
@@ -684,22 +708,22 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.48e-02   |
-      16  |  1.79e-02   |   4.1840    |  2.0649
-      32  |  4.38e-03   |   4.0830    |  2.0296
-      64  |  1.16e-03   |   3.7604    |  1.9109
+       8  |  7.59e-02   |
+      16  |  1.90e-02   |   3.9894    |  1.9962
+      32  |  4.72e-03   |   4.0321    |  2.0115
+      64  |  1.20e-03   |   3.9485    |  1.9813
     ---------------------------------------------
-    Testing is important.
+    Happy little convergence test!
     
     
     randomTensorMesh:  Interpolation 2D: CC
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  4.76e-02   |
-      16  |  6.08e-02   |   0.7830    |  -1.0896
-      32  |  1.14e-02   |   5.3162    |  1.7097
-      64  |  1.19e-03   |   9.6340    |  3.5674
+       8  |  1.52e-01   |
+      16  |  2.90e-02   |   5.2472    |  1.9020
+      32  |  9.09e-03   |   3.1868    |  1.4098
+      64  |  2.68e-03   |   3.3935    |  2.1529
     ---------------------------------------------
     Well done Rowan!
     
@@ -732,24 +756,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  6.98e-02   |
-      16  |  1.85e-02   |   3.7771    |  1.9173
-      32  |  4.76e-03   |   3.8827    |  1.9571
-      64  |  1.20e-03   |   3.9676    |  1.9883
+       8  |  6.57e-02   |
+      16  |  1.88e-02   |   3.4893    |  1.8030
+      32  |  4.79e-03   |   3.9343    |  1.9761
+      64  |  1.17e-03   |   4.1029    |  2.0366
     ---------------------------------------------
-    That was easy!
+    Go Test Go!
     
     
     randomTensorMesh:  Interpolation 2D: Ex
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  2.61e-01   |
-      16  |  3.89e-02   |   6.7177    |  2.4144
-      32  |  1.41e-02   |   2.7587    |  1.3921
-      64  |  4.59e-03   |   3.0696    |  1.9255
+       8  |  1.16e-01   |
+      16  |  3.68e-02   |   3.1394    |  1.5218
+      32  |  8.77e-03   |   4.1974    |  3.1728
+      64  |  3.14e-03   |   2.7963    |  1.7619
     ---------------------------------------------
-    That was easy!
+    Testing is important.
     
     
     
@@ -780,24 +804,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  6.97e-02   |
-      16  |  1.87e-02   |   3.7164    |  1.8939
-      32  |  4.74e-03   |   3.9516    |  1.9825
-      64  |  1.14e-03   |   4.1658    |  2.0586
+       8  |  7.04e-02   |
+      16  |  1.88e-02   |   3.7494    |  1.9067
+      32  |  4.55e-03   |   4.1251    |  2.0444
+      64  |  1.16e-03   |   3.9315    |  1.9751
     ---------------------------------------------
-    That was easy!
+    You get a gold star!
     
     
     randomTensorMesh:  Interpolation 2D: Ey
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  1.82e-01   |
-      16  |  3.85e-02   |   4.7305    |  2.2574
-      32  |  9.33e-03   |   4.1268    |  2.3056
-      64  |  3.29e-03   |   2.8348    |  1.1080
+       8  |  1.99e-01   |
+      16  |  8.02e-02   |   2.4773    |  4.5721
+      32  |  1.30e-02   |   6.1949    |  1.7850
+      64  |  3.95e-03   |   3.2812    |  1.9890
     ---------------------------------------------
-    Happy little convergence test!
+    Once upon a time, a happy little test passed.
     
     
     
@@ -828,24 +852,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.48e-02   |
-      16  |  1.79e-02   |   4.1840    |  2.0649
-      32  |  4.38e-03   |   4.0830    |  2.0296
-      64  |  1.16e-03   |   3.7604    |  1.9109
+       8  |  7.59e-02   |
+      16  |  1.90e-02   |   3.9894    |  1.9962
+      32  |  4.72e-03   |   4.0321    |  2.0115
+      64  |  1.20e-03   |   3.9485    |  1.9813
     ---------------------------------------------
-    You are awesome.
+    Yay passed!
     
     
     randomTensorMesh:  Interpolation 2D: Fx
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  9.78e-02   |
-      16  |  4.87e-02   |   2.0098    |  1.1486
-      32  |  6.39e-03   |   7.6229    |  2.5269
-      64  |  3.06e-03   |   2.0873    |  1.0490
+       8  |  5.31e-02   |
+      16  |  4.42e-02   |   1.2007    |  0.2070
+      32  |  1.05e-02   |   4.2086    |  2.5806
+      64  |  1.80e-03   |   5.8391    |  2.4462
     ---------------------------------------------
-    You are awesome.
+    The test be workin!
     
     
     
@@ -876,24 +900,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.59e-02   |
-      16  |  1.90e-02   |   3.9914    |  1.9969
-      32  |  4.62e-03   |   4.1142    |  2.0406
-      64  |  1.15e-03   |   4.0248    |  2.0089
+       8  |  7.61e-02   |
+      16  |  1.92e-02   |   3.9697    |  1.9890
+      32  |  4.76e-03   |   4.0235    |  2.0084
+      64  |  1.15e-03   |   4.1348    |  2.0478
     ---------------------------------------------
-    Happy little convergence test!
+    Not just a pretty face Rowan
     
     
     randomTensorMesh:  Interpolation 2D: Fy
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  1.34e-01   |
-      16  |  2.73e-02   |   4.9171    |  2.4681
-      32  |  1.10e-02   |   2.4900    |  1.2825
-      64  |  3.06e-03   |   3.5774    |  1.9074
+       8  |  1.37e-01   |
+      16  |  3.12e-02   |   4.3873    |  2.7317
+      32  |  7.83e-03   |   3.9895    |  2.3562
+      64  |  2.26e-03   |   3.4691    |  1.3948
     ---------------------------------------------
-    Testing is important.
+    Happy little convergence test!
     
     
     
@@ -924,24 +948,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  6.98e-02   |
-      16  |  1.85e-02   |   3.7771    |  1.9173
-      32  |  4.76e-03   |   3.8827    |  1.9571
-      64  |  1.20e-03   |   3.9676    |  1.9883
+       8  |  6.57e-02   |
+      16  |  1.88e-02   |   3.4893    |  1.8030
+      32  |  4.79e-03   |   3.9343    |  1.9761
+      64  |  1.17e-03   |   4.1029    |  2.0366
     ---------------------------------------------
-    You deserve a pat on the back!
+    Well done Rowan!
     
     
     randomTensorMesh:  Interpolation 2D: N
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  2.05e-01   |
-      16  |  4.37e-02   |   4.6889    |  2.7894
-      32  |  1.30e-02   |   3.3576    |  1.7086
-      64  |  3.69e-03   |   3.5290    |  1.8138
+       8  |  2.16e-01   |
+      16  |  3.36e-02   |   6.4220    |  1.9825
+      32  |  2.07e-02   |   1.6177    |  1.1560
+      64  |  3.05e-03   |   6.8009    |  2.5482
     ---------------------------------------------
-    Go Test Go!
+    The test be workin!
     
     
     
@@ -981,24 +1005,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.56e-02   |
-      16  |  1.87e-02   |   4.0396    |  2.0142
-      32  |  4.52e-03   |   4.1409    |  2.0499
-      64  |  1.19e-03   |   3.8144    |  1.9315
+       8  |  7.38e-02   |
+      16  |  1.70e-02   |   4.3363    |  2.1165
+      32  |  4.67e-03   |   3.6472    |  1.8668
+      64  |  1.12e-03   |   4.1563    |  2.0553
     ---------------------------------------------
-    Not just a pretty face Rowan
+    Well done Rowan!
     
     
     randomTensorMesh:  Interpolation CC
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  1.03e-01   |
-      16  |  5.43e-02   |   1.8919    |  1.5786
-      32  |  1.33e-02   |   4.0699    |  1.4626
-      64  |  3.40e-03   |   3.9265    |  2.7539
+       8  |  2.39e-01   |
+      16  |  2.83e-02   |   8.4410    |  3.3650
+      32  |  8.35e-03   |   3.3902    |  1.4106
+      64  |  2.59e-03   |   3.2260    |  1.7120
     ---------------------------------------------
-    Awesome, Rowan, just awesome.
+    You are awesome.
     
     
     
@@ -1029,24 +1053,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.04e-02   |
-      16  |  1.86e-02   |   3.7784    |  1.9178
-      32  |  4.78e-03   |   3.8981    |  1.9628
-      64  |  1.20e-03   |   3.9741    |  1.9906
+       8  |  6.97e-02   |
+      16  |  1.88e-02   |   3.7094    |  1.8912
+      32  |  4.68e-03   |   4.0208    |  2.0075
+      64  |  1.17e-03   |   3.9957    |  1.9985
     ---------------------------------------------
-    Well done Rowan!
+    You deserve a pat on the back!
     
     
     randomTensorMesh:  Interpolation Ex
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  2.52e-01   |
-      16  |  2.58e-02   |   9.7516    |  5.3251
-      32  |  1.69e-02   |   1.5285    |  0.4736
-      64  |  3.83e-03   |   4.4134    |  2.6664
+       8  |  2.13e-01   |
+      16  |  5.12e-02   |   4.1677    |  1.3954
+      32  |  1.99e-02   |   2.5661    |  1.5524
+      64  |  3.48e-03   |   5.7313    |  2.1382
     ---------------------------------------------
-    You deserve a pat on the back!
+    You get a gold star!
     
     
     
@@ -1077,24 +1101,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.03e-02   |
-      16  |  1.83e-02   |   3.8430    |  1.9422
-      32  |  4.72e-03   |   3.8783    |  1.9554
-      64  |  1.17e-03   |   4.0244    |  2.0088
+       8  |  7.04e-02   |
+      16  |  1.86e-02   |   3.7826    |  1.9194
+      32  |  4.79e-03   |   3.8833    |  1.9573
+      64  |  1.17e-03   |   4.1003    |  2.0357
     ---------------------------------------------
-    Yay passed!
+    That was easy!
     
     
     randomTensorMesh:  Interpolation Ey
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.70e-02   |
-      16  |  2.92e-02   |   2.6359    |  2.0836
-      32  |  1.66e-02   |   1.7612    |  0.7062
-      64  |  3.79e-03   |   4.3764    |  2.5121
+       8  |  2.43e-01   |
+      16  |  5.66e-02   |   4.2904    |  1.4531
+      32  |  1.55e-02   |   3.6614    |  2.1536
+      64  |  3.27e-03   |   4.7300    |  2.2808
     ---------------------------------------------
-    Well done Rowan!
+    That was easy!
     
     
     
@@ -1125,24 +1149,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.04e-02   |
-      16  |  1.87e-02   |   3.7587    |  1.9102
-      32  |  4.79e-03   |   3.9106    |  1.9674
-      64  |  1.17e-03   |   4.1061    |  2.0378
+       8  |  6.98e-02   |
+      16  |  1.81e-02   |   3.8612    |  1.9491
+      32  |  4.70e-03   |   3.8479    |  1.9441
+      64  |  1.20e-03   |   3.9268    |  1.9733
     ---------------------------------------------
-    You are awesome.
+    And then everyone was happy.
     
     
     randomTensorMesh:  Interpolation Ez
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  1.56e-01   |
-      16  |  3.70e-02   |   4.2142    |  1.7123
-      32  |  1.16e-02   |   3.1865    |  1.2409
-      64  |  4.18e-03   |   2.7767    |  1.8503
+       8  |  9.77e-02   |
+      16  |  4.93e-02   |   1.9820    |  2.0322
+      32  |  1.34e-02   |   3.6656    |  1.4108
+      64  |  4.90e-03   |   2.7436    |  1.4520
     ---------------------------------------------
-    That was easy!
+    Yay passed!
     
     
     
@@ -1173,22 +1197,22 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.56e-02   |
-      16  |  1.87e-02   |   4.0396    |  2.0142
-      32  |  4.52e-03   |   4.1409    |  2.0499
-      64  |  1.19e-03   |   3.8144    |  1.9315
+       8  |  7.38e-02   |
+      16  |  1.70e-02   |   4.3363    |  2.1165
+      32  |  4.67e-03   |   3.6472    |  1.8668
+      64  |  1.12e-03   |   4.1563    |  2.0553
     ---------------------------------------------
-    The test be workin!
+    Yay passed!
     
     
     randomTensorMesh:  Interpolation Fx
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  4.85e-02   |
-      16  |  3.26e-02   |   1.4891    |  0.4457
-      32  |  7.59e-03   |   4.2884    |  2.0210
-      64  |  2.21e-03   |   3.4371    |  2.3960
+       8  |  1.51e-01   |
+      16  |  2.31e-02   |   6.5343    |  3.0196
+      32  |  1.32e-02   |   1.7456    |  0.8096
+      64  |  2.24e-03   |   5.8991    |  2.1039
     ---------------------------------------------
     Go Test Go!
     
@@ -1221,24 +1245,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.55e-02   |
-      16  |  1.86e-02   |   4.0668    |  2.0239
-      32  |  4.25e-03   |   4.3625    |  2.1252
-      64  |  1.12e-03   |   3.7855    |  1.9205
+       8  |  7.61e-02   |
+      16  |  1.92e-02   |   3.9630    |  1.9866
+      32  |  4.81e-03   |   3.9962    |  1.9986
+      64  |  1.20e-03   |   4.0135    |  2.0049
     ---------------------------------------------
-    Once upon a time, a happy little test passed.
+    Well done Rowan!
     
     
     randomTensorMesh:  Interpolation Fy
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  1.26e-01   |
-      16  |  5.43e-02   |   2.3127    |  1.1973
-      32  |  6.95e-03   |   7.8090    |  3.7839
-      64  |  2.07e-03   |   3.3670    |  1.4675
+       8  |  1.80e-01   |
+      16  |  4.05e-02   |   4.4477    |  1.4993
+      32  |  1.59e-02   |   2.5460    |  1.1426
+      64  |  2.54e-03   |   6.2522    |  2.3141
     ---------------------------------------------
-    You get a gold star!
+    Not just a pretty face Rowan
     
     
     
@@ -1269,24 +1293,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.56e-02   |
-      16  |  1.87e-02   |   4.0457    |  2.0164
-      32  |  4.58e-03   |   4.0787    |  2.0281
-      64  |  1.19e-03   |   3.8390    |  1.9407
+       8  |  7.51e-02   |
+      16  |  1.82e-02   |   4.1301    |  2.0462
+      32  |  4.16e-03   |   4.3667    |  2.1266
+      64  |  1.13e-03   |   3.6758    |  1.8781
     ---------------------------------------------
-    You are awesome.
+    And then everyone was happy.
     
     
     randomTensorMesh:  Interpolation Fz
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  1.17e-01   |
-      16  |  2.69e-02   |   4.3374    |  2.1585
-      32  |  1.64e-02   |   1.6402    |  0.8257
-      64  |  3.86e-03   |   4.2449    |  1.7674
+       8  |  1.32e-01   |
+      16  |  5.78e-02   |   2.2902    |  1.6547
+      32  |  9.69e-03   |   5.9708    |  2.9133
+      64  |  4.41e-03   |   2.1983    |  1.1867
     ---------------------------------------------
-    Happy little convergence test!
+    Well done Rowan!
     
     
     
@@ -1317,24 +1341,24 @@ Test Results
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  7.04e-02   |
-      16  |  1.86e-02   |   3.7784    |  1.9178
-      32  |  4.78e-03   |   3.8981    |  1.9628
-      64  |  1.20e-03   |   3.9741    |  1.9906
+       8  |  6.97e-02   |
+      16  |  1.88e-02   |   3.7094    |  1.8912
+      32  |  4.68e-03   |   4.0208    |  2.0075
+      64  |  1.17e-03   |   3.9957    |  1.9985
     ---------------------------------------------
-    The test be workin!
+    Yay passed!
     
     
     randomTensorMesh:  Interpolation N
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       8  |  1.78e-01   |
-      16  |  3.66e-02   |   4.8776    |  3.1782
-      32  |  1.61e-02   |   2.2707    |  1.1549
-      64  |  3.17e-03   |   5.0740    |  2.0894
+       8  |  7.50e-02   |
+      16  |  7.65e-02   |   0.9805    |  -0.0299
+      32  |  1.75e-02   |   4.3794    |  2.0088
+      64  |  4.46e-03   |   3.9163    |  1.9951
     ---------------------------------------------
-    Happy little convergence test!
+    And then everyone was happy.
     
     
     
@@ -1426,7 +1450,7 @@ Test Results
       16  |  2.42e-03   |
       32  |  6.06e-04   |   4.0001    |  2.0000
     ---------------------------------------------
-    Happy little convergence test!
+    Testing is important.
     
     
     uniformLOM:  Edge Inner Product - Isotropic
@@ -1446,7 +1470,7 @@ Test Results
       16  |  2.81e-03   |
       32  |  7.11e-04   |   3.9432    |  1.9794
     ---------------------------------------------
-    You deserve a pat on the back!
+    Not just a pretty face Rowan
     
     
     
@@ -1480,7 +1504,7 @@ Test Results
       16  |  6.29e-04   |
       32  |  1.57e-04   |   3.9978    |  1.9992
     ---------------------------------------------
-    Testing is important.
+    And then everyone was happy.
     
     
     uniformLOM:  Face Inner Product - Isotropic
@@ -1490,7 +1514,7 @@ Test Results
       16  |  6.29e-04   |
       32  |  1.57e-04   |   3.9978    |  1.9992
     ---------------------------------------------
-    And then everyone was happy.
+    Not just a pretty face Rowan
     
     
     rotateLOM:  Face Inner Product - Isotropic
@@ -1500,7 +1524,7 @@ Test Results
       16  |  3.08e-04   |
       32  |  7.07e-05   |   4.3564    |  2.1231
     ---------------------------------------------
-    Not just a pretty face Rowan
+    Well done Rowan!
     
     
     
@@ -1534,7 +1558,7 @@ Test Results
       16  |  6.99e-03   |
       32  |  1.75e-03   |   3.9996    |  1.9998
     ---------------------------------------------
-    That was easy!
+    And then everyone was happy.
     
     
     uniformLOM:  Edge Inner Product - Anisotropic
@@ -1544,7 +1568,7 @@ Test Results
       16  |  6.99e-03   |
       32  |  1.75e-03   |   3.9996    |  1.9998
     ---------------------------------------------
-    You get a gold star!
+    Go Test Go!
     
     
     rotateLOM:  Edge Inner Product - Anisotropic
@@ -1554,7 +1578,7 @@ Test Results
       16  |  7.70e-03   |
       32  |  1.94e-03   |   3.9622    |  1.9863
     ---------------------------------------------
-    Happy little convergence test!
+    Go Test Go!
     
     
     
@@ -1588,7 +1612,7 @@ Test Results
       16  |  2.68e-03   |
       32  |  6.69e-04   |   3.9982    |  1.9993
     ---------------------------------------------
-    And then everyone was happy.
+    You deserve a pat on the back!
     
     
     uniformLOM:  Face Inner Product - Anisotropic
@@ -1598,7 +1622,7 @@ Test Results
       16  |  2.68e-03   |
       32  |  6.69e-04   |   3.9982    |  1.9993
     ---------------------------------------------
-    Yay passed!
+    You are awesome.
     
     
     rotateLOM:  Face Inner Product - Anisotropic
@@ -1608,7 +1632,7 @@ Test Results
       16  |  2.15e-03   |
       32  |  5.25e-04   |   4.0845    |  2.0302
     ---------------------------------------------
-    Yay passed!
+    The test be workin!
     
     
     
@@ -1642,7 +1666,7 @@ Test Results
       16  |  6.79e-03   |
       32  |  1.70e-03   |   3.9996    |  1.9998
     ---------------------------------------------
-    The test be workin!
+    Not just a pretty face Rowan
     
     
     uniformLOM:  Edge Inner Product - Full Tensor
@@ -1652,7 +1676,7 @@ Test Results
       16  |  6.79e-03   |
       32  |  1.70e-03   |   3.9996    |  1.9998
     ---------------------------------------------
-    And then everyone was happy.
+    Go Test Go!
     
     
     rotateLOM:  Edge Inner Product - Full Tensor
@@ -1662,7 +1686,7 @@ Test Results
       16  |  7.49e-03   |
       32  |  1.89e-03   |   3.9617    |  1.9861
     ---------------------------------------------
-    The test be workin!
+    Not just a pretty face Rowan
     
     
     
@@ -1696,7 +1720,7 @@ Test Results
       16  |  3.10e-03   |
       32  |  7.74e-04   |   3.9981    |  1.9993
     ---------------------------------------------
-    You get a gold star!
+    You are awesome.
     
     
     uniformLOM:  Face Inner Product - Full Tensor
@@ -1706,7 +1730,7 @@ Test Results
       16  |  3.10e-03   |
       32  |  7.74e-04   |   3.9981    |  1.9993
     ---------------------------------------------
-    The test be workin!
+    You are awesome.
     
     
     rotateLOM:  Face Inner Product - Full Tensor
@@ -1716,7 +1740,7 @@ Test Results
       16  |  2.54e-03   |
       32  |  6.23e-04   |   4.0741    |  2.0265
     ---------------------------------------------
-    You deserve a pat on the back!
+    Testing is important.
     
     
     
@@ -1763,7 +1787,7 @@ Test Results
       64  |  2.18e-02   |   4.0001    |  2.0000
      128  |  5.46e-03   |   4.0000    |  2.0000
     ---------------------------------------------
-    The test be workin!
+    Happy little convergence test!
     
     
     uniformLOM:  2D Edge Inner Product - Isotropic
@@ -1777,7 +1801,7 @@ Test Results
       64  |  2.18e-02   |   4.0001    |  2.0000
      128  |  5.46e-03   |   4.0000    |  2.0000
     ---------------------------------------------
-    You get a gold star!
+    You deserve a pat on the back!
     
     
     rotateLOM:  2D Edge Inner Product - Isotropic
@@ -1791,7 +1815,7 @@ Test Results
       64  |  2.00e-02   |   4.0155    |  2.0056
      128  |  5.00e-03   |   4.0038    |  2.0014
     ---------------------------------------------
-    You are awesome.
+    Go Test Go!
     
     
     
@@ -1829,7 +1853,7 @@ Test Results
       64  |  2.53e-02   |   4.0000    |  2.0000
      128  |  6.32e-03   |   4.0000    |  2.0000
     ---------------------------------------------
-    Happy little convergence test!
+    You get a gold star!
     
     
     uniformLOM:  2D Face Inner Product - Isotropic
@@ -1843,7 +1867,7 @@ Test Results
       64  |  2.53e-02   |   4.0000    |  2.0000
      128  |  6.32e-03   |   4.0000    |  2.0000
     ---------------------------------------------
-    Not just a pretty face Rowan
+    Yay passed!
     
     
     rotateLOM:  2D Face Inner Product - Isotropic
@@ -1857,7 +1881,7 @@ Test Results
       64  |  2.30e-02   |   4.0132    |  2.0048
      128  |  5.74e-03   |   4.0009    |  2.0003
     ---------------------------------------------
-    That was easy!
+    Well done Rowan!
     
     
     
@@ -1895,7 +1919,7 @@ Test Results
       64  |  1.30e-01   |   4.0000    |  2.0000
      128  |  3.24e-02   |   4.0000    |  2.0000
     ---------------------------------------------
-    That was easy!
+    And then everyone was happy.
     
     
     uniformLOM:  2D Face Inner Product - Anisotropic
@@ -1909,7 +1933,7 @@ Test Results
       64  |  1.30e-01   |   4.0000    |  2.0000
      128  |  3.24e-02   |   4.0000    |  2.0000
     ---------------------------------------------
-    Go Test Go!
+    Awesome, Rowan, just awesome.
     
     
     rotateLOM:  2D Face Inner Product - Anisotropic
@@ -1961,7 +1985,7 @@ Test Results
       64  |  3.66e-02   |   4.0000    |  2.0000
      128  |  9.14e-03   |   4.0000    |  2.0000
     ---------------------------------------------
-    Not just a pretty face Rowan
+    You get a gold star!
     
     
     uniformLOM:  2D Edge Inner Product - Anisotropic
@@ -1975,7 +1999,7 @@ Test Results
       64  |  3.66e-02   |   4.0000    |  2.0000
      128  |  9.14e-03   |   4.0000    |  2.0000
     ---------------------------------------------
-    Well done Rowan!
+    Awesome, Rowan, just awesome.
     
     
     rotateLOM:  2D Edge Inner Product - Anisotropic
@@ -1989,7 +2013,7 @@ Test Results
       64  |  2.23e-02   |   4.0739    |  2.0264
      128  |  5.54e-03   |   4.0255    |  2.0092
     ---------------------------------------------
-    Testing is important.
+    Well done Rowan!
     
     
     
@@ -2027,7 +2051,7 @@ Test Results
       64  |  1.52e-01   |   4.0000    |  2.0000
      128  |  3.80e-02   |   4.0000    |  2.0000
     ---------------------------------------------
-    Testing is important.
+    And then everyone was happy.
     
     
     uniformLOM:  2D Face Inner Product - Full Tensor
@@ -2041,7 +2065,7 @@ Test Results
       64  |  1.52e-01   |   4.0000    |  2.0000
      128  |  3.80e-02   |   4.0000    |  2.0000
     ---------------------------------------------
-    Awesome, Rowan, just awesome.
+    Go Test Go!
     
     
     rotateLOM:  2D Face Inner Product - Full Tensor
@@ -2055,7 +2079,7 @@ Test Results
       64  |  1.45e-01   |   4.0074    |  2.0027
      128  |  3.62e-02   |   3.9984    |  1.9994
     ---------------------------------------------
-    You deserve a pat on the back!
+    Once upon a time, a happy little test passed.
     
     
     
@@ -2093,7 +2117,7 @@ Test Results
       64  |  3.78e-03   |   4.0001    |  2.0000
      128  |  9.46e-04   |   4.0000    |  2.0000
     ---------------------------------------------
-    Not just a pretty face Rowan
+    Go Test Go!
     
     
     uniformLOM:  2D Edge Inner Product - Full Tensor
@@ -2107,7 +2131,7 @@ Test Results
       64  |  3.78e-03   |   4.0001    |  2.0000
      128  |  9.46e-04   |   4.0000    |  2.0000
     ---------------------------------------------
-    Testing is important.
+    You get a gold star!
     
     
     rotateLOM:  2D Edge Inner Product - Full Tensor
@@ -2121,7 +2145,7 @@ Test Results
       64  |  1.98e-02   |   3.8708    |  1.9526
      128  |  5.03e-03   |   3.9418    |  1.9789
     ---------------------------------------------
-    Happy little convergence test!
+    Testing is important.
     
     
     
@@ -2134,15 +2158,15 @@ Test Results
     
     
         test_operators.TestAveraging2D
-        5
-        5
+        6
+        6
         0
         0
-        Detail
+        Detail
     
     
     
-        
test_orderE2CC

+        
test_orderCC2F

         
     
         
@@ -2157,38 +2181,38 @@ Test Results
             
             
     pt11.1: 
-    uniformTensorMesh:  Averaging 2D: E2CC
+    uniformTensorMesh:  Averaging 2D: CC2F
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       4  |  6.87e-03   |
-       8  |  1.76e-03   |   3.8978    |  1.9627
-      16  |  4.45e-04   |   3.9561    |  1.9841
-      32  |  1.12e-04   |   3.9799    |  1.9927
+       4  |  1.25e-01   |
+       8  |  6.25e-02   |   1.9961    |  0.9972
+      16  |  3.12e-02   |   1.9990    |  0.9993
+      32  |  1.56e-02   |   1.9998    |  0.9998
     ---------------------------------------------
-    Yay passed!
+    You are awesome.
     
     
-    uniformLOM:  Averaging 2D: E2CC
+    uniformLOM:  Averaging 2D: CC2F
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       4  |  6.87e-03   |
-       8  |  1.76e-03   |   3.8978    |  1.9627
-      16  |  4.45e-04   |   3.9561    |  1.9841
-      32  |  1.12e-04   |   3.9799    |  1.9927
+       4  |  1.25e-01   |
+       8  |  6.25e-02   |   1.9961    |  0.9972
+      16  |  3.12e-02   |   1.9990    |  0.9993
+      32  |  1.56e-02   |   1.9998    |  0.9998
     ---------------------------------------------
-    That was easy!
+    And then everyone was happy.
     
     
-    rotateLOM:  Averaging 2D: E2CC
+    rotateLOM:  Averaging 2D: CC2F
     _____________________________________________
        h  |    error    | e(i-1)/e(i) |  order
     ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
-       4  |  6.88e-03   |
-       8  |  1.82e-03   |   3.7846    |  1.9202
-      16  |  4.66e-04   |   3.8970    |  1.9624
-      32  |  1.18e-04   |   3.9552    |  1.9837
+       4  |  1.28e-01   |
+       8  |  7.22e-02   |   1.7667    |  0.8210
+      16  |  3.81e-02   |   1.8958    |  0.9228
+      32  |  1.95e-02   |   1.9507    |  0.9640
     ---------------------------------------------
     Go Test Go!
     
@@ -2202,7 +2226,7 @@ Test Results
     
     
     
-        
test_orderF2CC

+        
test_orderE2CC

         
     
         
@@ -2217,6 +2241,66 @@ Test Results
             
             
     pt11.2: 
+    uniformTensorMesh:  Averaging 2D: E2CC
+    _____________________________________________
+       h  |    error    | e(i-1)/e(i) |  order
+    ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
+       4  |  6.87e-03   |
+       8  |  1.76e-03   |   3.8978    |  1.9627
+      16  |  4.45e-04   |   3.9561    |  1.9841
+      32  |  1.12e-04   |   3.9799    |  1.9927
+    ---------------------------------------------
+    That was easy!
+    
+    
+    uniformLOM:  Averaging 2D: E2CC
+    _____________________________________________
+       h  |    error    | e(i-1)/e(i) |  order
+    ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
+       4  |  6.87e-03   |
+       8  |  1.76e-03   |   3.8978    |  1.9627
+      16  |  4.45e-04   |   3.9561    |  1.9841
+      32  |  1.12e-04   |   3.9799    |  1.9927
+    ---------------------------------------------
+    And then everyone was happy.
+    
+    
+    rotateLOM:  Averaging 2D: E2CC
+    _____________________________________________
+       h  |    error    | e(i-1)/e(i) |  order
+    ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
+       4  |  6.88e-03   |
+       8  |  1.82e-03   |   3.7846    |  1.9202
+      16  |  4.66e-04   |   3.8970    |  1.9624
+      32  |  1.18e-04   |   3.9552    |  1.9837
+    ---------------------------------------------
+    That was easy!
+    
+    
+    
+            

+        
+        
+    
+        
+    
+    
+    
+        
test_orderF2CC

+        
+    
+        
+        
+            pass
+    
+        
-        
-    
-        
-    
-    
-    
-        
test_orderN2CC

-        
-    
-        
-        
-            pass
-    
-        
+        
+    
+        
+    
+    
+    
+        
test_orderN2F

+        
+    
+        
+        
+            pass
+    
+        
@@ -2691,7 +2775,7 @@ Test Results
     
     
     
-        
test_orderN2F

+        
test_orderN2E

         
     
         
@@ -2706,6 +2790,66 @@ Test Results
             
             
     pt12.5: 
+    uniformTensorMesh:  Averaging 3D: N2E
+    _____________________________________________
+       h  |    error    | e(i-1)/e(i) |  order
+    ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
+       4  |  1.88e-02   |
+       8  |  4.99e-03   |   3.7613    |  1.9112
+      16  |  1.29e-03   |   3.8779    |  1.9553
+      32  |  3.27e-04   |   3.9382    |  1.9775
+    ---------------------------------------------
+    Yay passed!
+    
+    
+    uniformLOM:  Averaging 3D: N2E
+    _____________________________________________
+       h  |    error    | e(i-1)/e(i) |  order
+    ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
+       4  |  1.88e-02   |
+       8  |  4.99e-03   |   3.7613    |  1.9112
+      16  |  1.29e-03   |   3.8779    |  1.9553
+      32  |  3.27e-04   |   3.9382    |  1.9775
+    ---------------------------------------------
+    Yay passed!
+    
+    
+    rotateLOM:  Averaging 3D: N2E
+    _____________________________________________
+       h  |    error    | e(i-1)/e(i) |  order
+    ~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~
+       4  |  2.41e-02   |
+       8  |  6.11e-03   |   3.9523    |  1.9827
+      16  |  1.70e-03   |   3.5908    |  1.8443
+      32  |  4.41e-04   |   3.8566    |  1.9473
+    ---------------------------------------------
+    Go Test Go!
+    
+    
+    
+            

+        
+        
+    
+        
+    
+    
+    
+        
test_orderN2F

+        
+    
+        
+        
+            pass
+    
+        
+        
+    
+        
+    
+    
+    
+        test_operators.TestCellGrad3D_Neumann
+        1
+        1
+        0
+        0
+        Detail
+    
+    
+    
+        
test_order

+        
+    
+        
+        
+            pass
+    
+        
+        
+    
+        
+    
+    
+    
+        test_operators.TestCurl
+        1
+        1
+        0
+        0
+        Detail
+    
+    
+    
+        
test_order

+        
+    
+        
+        
+            pass
+    
+        
+        
+    
+        
+    
+    
+    
+        test_operators.TestFaceDiv2D
+        1
+        1
+        0
+        0
+        Detail
+    
+    
+    
+        
test_order

+        
+    
+        
+        
+            pass
+    
+        
+        
+    
+        
+    
+    
+    
+        test_operators.TestFaceDiv3D
+        1
+        1
+        0
+        0
+        Detail
+    
+    
+    
+        
test_order

+        
+    
+        
+        
+            pass
+    
+        
+        
+    
+        
+    
+    
+    
+        test_operators.TestNodalGrad
+        1
+        1
+        0
+        0
+        Detail
+    
+    
+    
+        
test_order

+        
+    
+        
+        
+            pass
+    
+