From 52c9cf83c94064f079e2608899dcc3c3c7863cd3 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 13 Nov 2013 19:34:11 -0800 Subject: [PATCH] Tests on Cell Grad (bug fixes for non-uniform mesh). and aveCC2F with extrapolation. --- SimPEG/__init__.py | 5 +++ SimPEG/mesh/DiffOperators.py | 80 +++++++++++++++++++++++++++++----- SimPEG/tests/test_operators.py | 53 +++++++++++++++++++++- 3 files changed, 126 insertions(+), 12 deletions(-) diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 7f059a74..5e399ebd 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -4,3 +4,8 @@ import mesh import inverse import forward import regularization + + +import scipy.version as _v +if _v.version < '0.13.0': + print 'Warning: upgrade your scipy to 0.13.0' diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 27fe3cc2..2c94a0cf 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -99,7 +99,7 @@ def ddxCellGradBC(n, bc): """ bc = checkBC(bc) - ij = (np.array([0, n+1]),np.array([0, 1])) + ij = (np.array([0, n]),np.array([0, 1])) vals = np.zeros(2) # Set the first side @@ -233,17 +233,19 @@ class DiffOperators(object): for i, bc_i in enumerate(BC): BC[i] = checkBC(bc_i) - self._cellGrad = None # ensure we create a new gradient next time we call it - self._cellGradBC = BC + # ensure we create a new gradient next time we call it + self._cellGrad = None + self._cellGradBC = None + self._cellGradBC_list = BC return BC - _cellGradBC = 'neumann' + _cellGradBC_list = 'neumann' def cellGrad(): doc = "The cell centered Gradient, takes you to cell faces." def fget(self): if(self._cellGrad is None): - BC = self.setCellGradBC(self._cellGradBC) + BC = self.setCellGradBC(self._cellGradBC_list) n = self.n if(self.dim == 1): G = ddxCellGrad(n[0], BC[0]) @@ -258,13 +260,40 @@ class DiffOperators(object): G = sp.vstack((G1, G2, G3), format="csr") # Compute areas of cell faces & volumes S = self.area - V = self.vol - self._cellGrad = sdiag(S)*G*sdiag(1/V) + V = self.aveCC2F*self.vol # Average volume between adjacent cells + self._cellGrad = sdiag(S/V)*G return self._cellGrad return locals() _cellGrad = None cellGrad = property(**cellGrad()) + def cellGradBC(): + doc = "The cell centered Gradient boundary condition matrix" + + def fget(self): + if(self._cellGradBC is None): + BC = self.setCellGradBC(self._cellGradBC_list) + n = self.n + if(self.dim == 1): + G = ddxCellGradBC(n[0], BC[0]) + elif(self.dim == 2): + G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0])) + G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0])) + G = sp.vstack((G1, G2), format="csr") + elif(self.dim == 3): + G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0])) + G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0])) + G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0])) + G = sp.vstack((G1, G2, G3), format="csr") + # Compute areas of cell faces & volumes + S = self.area + V = self.vol + self._cellGradBC = sdiag(S)*G*sdiag(1/V[[0,-1]]) + return self._cellGradBC + return locals() + _cellGradBC = None + cellGradBC = property(**cellGradBC()) + def cellGradx(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." @@ -377,16 +406,47 @@ class DiffOperators(object): self._aveF2CC = av(n[0]) elif(self.dim == 2): self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])), - sp.kron(av(n[1]), speye(n[0]))), format="csr") + sp.kron(av(n[1]), speye(n[0]))), format="csr") elif(self.dim == 3): self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") + kron3(speye(n[2]), av(n[1]), speye(n[0])), + kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") return self._aveF2CC return locals() _aveF2CC = None aveF2CC = property(**aveF2CC()) + def aveCC2F(): + doc = "Construct the averaging operator on cell cell centers to faces." + + def fget(self): + if(self._aveCC2F is None): + n = self.n + if(self.dim == 1): + Av = av(n[0]).T + Av = sdiag(1/Av.sum(axis=1))*Av + self._aveCC2F = Av + elif(self.dim == 2): + Av1 = av(n[0]).T + Av1 = sdiag(1/Av1.sum(axis=1))*Av1 + Av2 = av(n[1]).T + Av2 = sdiag(1/Av2.sum(axis=1))*Av2 + Av = sp.vstack((sp.kron(speye(n[1]), Av1), sp.kron(Av2, speye(n[0]))), format="csr") + self._aveCC2F = Av + elif(self.dim == 3): + Av1 = av(n[0]).T + Av1 = sdiag(1/Av1.sum(axis=1))*Av1 + Av2 = av(n[1]).T + Av2 = sdiag(1/Av2.sum(axis=1))*Av2 + Av3 = av(n[2]).T + Av3 = sdiag(1/Av3.sum(axis=1))*Av3 + Av = sp.vstack((kron3(speye(n[2]), speye(n[1]), Av1), kron3(speye(n[2]), Av2, speye(n[0])), kron3(Av3, speye(n[1]), speye(n[0]))), format="csr") + self._aveCC2F = Av + return self._aveCC2F + return locals() + _aveCC2F = None + aveCC2F = property(**aveCC2F()) + def aveE2CC(): doc = "Construct the averaging operator on cell edges to cell centers." diff --git a/SimPEG/tests/test_operators.py b/SimPEG/tests/test_operators.py index d67a382b..d1c38443 100644 --- a/SimPEG/tests/test_operators.py +++ b/SimPEG/tests/test_operators.py @@ -49,6 +49,34 @@ class TestCurl(OrderTest): self.orderTest() +class TestCellGrad1D_InhomogeneousDirichlet(OrderTest): + name = "Cell Grad 1D - Dirichlet" + meshTypes = ['uniformTensorMesh'] + meshDimension = 1 + expectedOrders = 1 # because of the averaging involved in the ghost point. u_b = (u_n + u_g)/2 + meshSizes = [8, 16, 32, 64] + + def getError(self): + #Test function + fx = lambda x: -2*np.pi*np.sin(2*np.pi*x) + sol = lambda x: np.cos(2*np.pi*x) + + + xc = sol(self.M.gridCC) + + gradX_anal = fx(self.M.gridFx) + + bc = np.array([1,1]) + self.M.setCellGradBC('dirichlet') + gradX = self.M.cellGrad.dot(xc) + self.M.cellGradBC*bc + + err = np.linalg.norm((gradX-gradX_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + class TestCellGrad2D_Dirichlet(OrderTest): name = "Cell Grad 2D - Dirichlet" meshTypes = ['uniformTensorMesh'] @@ -81,7 +109,7 @@ class TestCellGrad3D_Dirichlet(OrderTest): name = "Cell Grad 3D - Dirichlet" meshTypes = ['uniformTensorMesh'] meshDimension = 3 - meshSizes = [8, 16, 32, 64] + meshSizes = [8, 16, 32] def getError(self): #Test function @@ -137,7 +165,7 @@ class TestCellGrad3D_Neumann(OrderTest): name = "Cell Grad 3D - Neumann" meshTypes = ['uniformTensorMesh'] meshDimension = 3 - meshSizes = [8, 16, 32, 64] + meshSizes = [8, 16, 32] def getError(self): #Test function @@ -310,6 +338,16 @@ class TestAveraging2D(OrderTest): self.getAve = lambda M: M.aveF2CC self.orderTest() + def test_orderCC2F(self): + self.name = "Averaging 2D: CC2F" + fun = lambda x, y: (np.cos(x)+np.sin(y)) + self.getHere = lambda M: call2(fun, M.gridCC) + self.getThere = lambda M: np.r_[call2(fun, M.gridFx), call2(fun, M.gridFy)] + self.getAve = lambda M: M.aveCC2F + self.expectedOrders = 1 + self.orderTest() + self.expectedOrders = 2 + def test_orderE2CC(self): self.name = "Averaging 2D: E2CC" @@ -371,6 +409,17 @@ class TestAveraging3D(OrderTest): self.getAve = lambda M: M.aveE2CC self.orderTest() + def test_orderCC2F(self): + self.name = "Averaging 3D: CC2F" + fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z)) + self.getHere = lambda M: call3(fun, M.gridCC) + self.getThere = lambda M: np.r_[call3(fun, M.gridFx), call3(fun, M.gridFy), call3(fun, M.gridFz)] + self.getAve = lambda M: M.aveCC2F + self.expectedOrders = 1 + self.orderTest() + self.expectedOrders = 2 + + if __name__ == '__main__': unittest.main()