From 56ee22d1cea65b1cb1bbe56e800a6431a2e53f1a Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 13 Nov 2013 15:09:13 -0800 Subject: [PATCH 1/7] =?UTF-8?q?Bug=20Fix=20in=20TensorView.plotImage=20(di?= =?UTF-8?q?dn=E2=80=99t=20work=20in=202D)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- SimPEG/mesh/TensorView.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py index 0b9ff7b3..3ae91259 100644 --- a/SimPEG/mesh/TensorView.py +++ b/SimPEG/mesh/TensorView.py @@ -89,18 +89,18 @@ class TensorView(object): # Determine the subplot number: 131, 121 numPlots = 130 if plotAll else len(imageType)/2*10+100 pltNum = 1 - fx, fy, fz = self.r(I,'F','F','M') + fxyz = self.r(I,'F','F','M') if plotAll or 'Fx' in imageType: ax_x = plt.subplot(numPlots+pltNum) - self.plotImage(fx, imageType='Fx', ax=ax_x, **options) + self.plotImage(fxyz[0], imageType='Fx', ax=ax_x, **options) pltNum +=1 if plotAll or 'Fy' in imageType: ax_y = plt.subplot(numPlots+pltNum) - self.plotImage(fy, imageType='Fy', ax=ax_y, **options) + self.plotImage(fxyz[1], imageType='Fy', ax=ax_y, **options) pltNum +=1 if plotAll or 'Fz' in imageType: ax_z = plt.subplot(numPlots+pltNum) - self.plotImage(fz, imageType='Fz', ax=ax_z, **options) + self.plotImage(fxyz[2], imageType='Fz', ax=ax_z, **options) pltNum +=1 return else: From 4c82ce7dc27ba595321232e5c303dcc41d439f99 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 13 Nov 2013 15:14:33 -0800 Subject: [PATCH 2/7] Documentation for cellGradient. Boundary Conditions. Testing in 2D & 3D. --- SimPEG/mesh/DiffOperators.py | 86 +++++++++++++++++++++++- SimPEG/tests/test_operators.py | 117 ++++++++++++++++++++++++++++++++- 2 files changed, 199 insertions(+), 4 deletions(-) diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index fa3f42e0..27fe3cc2 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -4,7 +4,13 @@ from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av def checkBC(bc): - """ Checks if boundary condition 'bc' is valid. """ + """ + + Checks if boundary condition 'bc' is valid. + + Each bc must be either 'dirichlet' or 'neumann' + + """ if(type(bc) is str): bc = [bc, bc] assert type(bc) is list, 'bc must be a list' @@ -17,7 +23,33 @@ def checkBC(bc): def ddxCellGrad(n, bc): - """Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1""" + """ + Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1 + + For Cell-Centered **Dirichlet**, use a ghost point:: + + (u_1 - u_g)/hf = grad + + u_g u_1 u_2 + * | * | * ... + ^ + 0 + + u_g = - u_1 + grad = 2*u1/dx + negitive on the other side. + + For Cell-Centered **Neumann**, use a ghost point:: + + (u_1 - u_g)/hf = 0 + + u_g u_1 u_2 + * | * | * ... + + u_g = u_1 + grad = 0; put a zero in. + + """ bc = checkBC(bc) D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr") @@ -33,6 +65,56 @@ def ddxCellGrad(n, bc): D[-1, -1] = 0 return D +def ddxCellGradBC(n, bc): + """ + + Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1 + + For Cell-Centered **Dirichlet**, use a ghost point:: + + (u_1 - u_g)/hf = grad + + u_g u_1 u_2 + * | * | * ... + ^ + u_b + + We know the value at the boundary (u_b):: + + (u_g+u_1)/2 = u_b (the average) + u_g = 2*u_b - u_1 + + So plug in to gradient: + + (u_1 - (2*u_b - u_1))/hf = grad + 2*(u_1-u_b)/hf = grad + + Separate, because BC are known (and can move to RHS later):: + + ( 2/hf )*u_1 + ( -2/hf )*u_b = grad + + ( ^ ) JUST RETURN THIS + + + """ + bc = checkBC(bc) + + ij = (np.array([0, n+1]),np.array([0, 1])) + vals = np.zeros(2) + + # Set the first side + if(bc[0] == 'dirichlet'): + vals[0] = -2 + elif(bc[0] == 'neumann'): + vals[0] = 0 + # Set the second side + if(bc[1] == 'dirichlet'): + vals[1] = 2 + elif(bc[1] == 'neumann'): + vals[1] = 0 + D = sp.csr_matrix((vals, ij), shape=(n+1,2)) + return D + class DiffOperators(object): """ diff --git a/SimPEG/tests/test_operators.py b/SimPEG/tests/test_operators.py index 26233626..d67a382b 100644 --- a/SimPEG/tests/test_operators.py +++ b/SimPEG/tests/test_operators.py @@ -1,6 +1,7 @@ import numpy as np import unittest from TestUtils import OrderTest +import matplotlib.pyplot as plt MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM'] call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) @@ -48,8 +49,120 @@ class TestCurl(OrderTest): self.orderTest() -class TestFaceDiv(OrderTest): - name = "Face Divergence" +class TestCellGrad2D_Dirichlet(OrderTest): + name = "Cell Grad 2D - Dirichlet" + meshTypes = ['uniformTensorMesh'] + meshDimension = 2 + meshSizes = [8, 16, 32, 64] + + def getError(self): + #Test function + fx = lambda x, y: 2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y) + fy = lambda x, y: 2*np.pi*np.cos(2*np.pi*y)*np.sin(2*np.pi*x) + sol = lambda x, y: np.sin(2*np.pi*x)*np.sin(2*np.pi*y) + + xc = call2(sol, self.M.gridCC) + + Fc = cartF2(self.M, fx, fy) + gradX_anal = self.M.projectFaceVector(Fc) + + self.M.setCellGradBC('dirichlet') + gradX = self.M.cellGrad.dot(xc) + + err = np.linalg.norm((gradX-gradX_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestCellGrad3D_Dirichlet(OrderTest): + name = "Cell Grad 3D - Dirichlet" + meshTypes = ['uniformTensorMesh'] + meshDimension = 3 + meshSizes = [8, 16, 32, 64] + + def getError(self): + #Test function + fx = lambda x, y, z: 2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)*np.sin(2*np.pi*z) + fy = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)*np.sin(2*np.pi*z) + fz = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)*np.sin(2*np.pi*y)*np.cos(2*np.pi*z) + sol = lambda x, y, z: np.sin(2*np.pi*x)*np.sin(2*np.pi*y)*np.sin(2*np.pi*z) + + xc = call3(sol, self.M.gridCC) + + Fc = cartF3(self.M, fx, fy, fz) + gradX_anal = self.M.projectFaceVector(Fc) + + self.M.setCellGradBC('dirichlet') + gradX = self.M.cellGrad.dot(xc) + + err = np.linalg.norm((gradX-gradX_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + +class TestCellGrad2D_Neumann(OrderTest): + name = "Cell Grad 2D - Neumann" + meshTypes = ['uniformTensorMesh'] + meshDimension = 2 + meshSizes = [8, 16, 32, 64] + + def getError(self): + #Test function + fx = lambda x, y: -2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y) + fy = lambda x, y: -2*np.pi*np.sin(2*np.pi*y)*np.cos(2*np.pi*x) + sol = lambda x, y: np.cos(2*np.pi*x)*np.cos(2*np.pi*y) + + xc = call2(sol, self.M.gridCC) + + Fc = cartF2(self.M, fx, fy) + gradX_anal = self.M.projectFaceVector(Fc) + + self.M.setCellGradBC('neumann') + gradX = self.M.cellGrad.dot(xc) + + err = np.linalg.norm((gradX-gradX_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestCellGrad3D_Neumann(OrderTest): + name = "Cell Grad 3D - Neumann" + meshTypes = ['uniformTensorMesh'] + meshDimension = 3 + meshSizes = [8, 16, 32, 64] + + def getError(self): + #Test function + fx = lambda x, y, z: -2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)*np.cos(2*np.pi*z) + fy = lambda x, y, z: -2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)*np.cos(2*np.pi*z) + fz = lambda x, y, z: -2*np.pi*np.cos(2*np.pi*x)*np.cos(2*np.pi*y)*np.sin(2*np.pi*z) + sol = lambda x, y, z: np.cos(2*np.pi*x)*np.cos(2*np.pi*y)*np.cos(2*np.pi*z) + + xc = call3(sol, self.M.gridCC) + + Fc = cartF3(self.M, fx, fy, fz) + gradX_anal = self.M.projectFaceVector(Fc) + + self.M.setCellGradBC('neumann') + gradX = self.M.cellGrad.dot(xc) + + err = np.linalg.norm((gradX-gradX_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + +class TestFaceDiv3D(OrderTest): + name = "Face Divergence 3D" meshTypes = MESHTYPES meshSizes = [8, 16, 32] From 6fbe3a616ba61c3a44eb1e9ac4cf2c56e111ad07 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 13 Nov 2013 15:52:56 -0800 Subject: [PATCH 3/7] Added shortcut to making a TensorMesh on a unit cube: mesh = TensorMesh([10, 12, 15]) --- SimPEG/mesh/TensorMesh.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 90f8fd10..30eebc13 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -26,18 +26,25 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): .. plot:: examples/mesh/plot_TensorMesh.py + For a quick tensor mesh on a (10x12x15) unit cube:: + + mesh = TensorMesh([10, 12, 15]) + """ _meshType = 'TENSOR' def __init__(self, h, x0=None): - super(TensorMesh, self).__init__(np.array([x.size for x in h]), x0) - - assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)" - for i, h_i in enumerate(h): + if type(h_i) in [int, long, float]: + # This gives you something over the unit cube. + h_i = np.ones(int(h_i))/int(h_i) + h[i] = h_i assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i) assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) + BaseMesh.__init__(self, np.array([x.size for x in h]), x0) + assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)" + # Ensure h contains 1D vectors self._h = [mkvc(x.astype(float)) for x in h] From 52c9cf83c94064f079e2608899dcc3c3c7863cd3 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 13 Nov 2013 19:34:11 -0800 Subject: [PATCH 4/7] 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() From 358411792dac59a44556b240e2daf0942fbdbc68 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 13 Nov 2013 19:37:14 -0800 Subject: [PATCH 5/7] Cleaned up averaging. Put it in an AveExtrap function. --- SimPEG/mesh/DiffOperators.py | 25 +++++++------------------ SimPEG/utils/__init__.py | 2 +- SimPEG/utils/sputils.py | 7 ++++++- 3 files changed, 14 insertions(+), 20 deletions(-) diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 2c94a0cf..7176b8fd 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -1,6 +1,6 @@ import numpy as np from scipy import sparse as sp -from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av +from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av, avExtrap def checkBC(bc): @@ -423,25 +423,14 @@ class DiffOperators(object): 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 + self._aveCC2F = avExtrap(n[0]) 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 + self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])), + sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr") 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 + self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])), + kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])), + kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr") return self._aveCC2F return locals() _aveCC2F = None diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index d2170721..75e7c360 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -4,7 +4,7 @@ import lomutils import interputils import ModelBuilder from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind -from sputils import spzeros, kron3, speye, sdiag, ddx, av +from sputils import spzeros, kron3, speye, sdiag, ddx, av, avExtrap from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird from interputils import interpmat from ipythonUtils import easyAnimate as animate diff --git a/SimPEG/utils/sputils.py b/SimPEG/utils/sputils.py index e066322b..06eef80d 100644 --- a/SimPEG/utils/sputils.py +++ b/SimPEG/utils/sputils.py @@ -29,5 +29,10 @@ def ddx(n): def av(n): - """Define 1D averaging operator from cell-centres to nodes.""" + """Define 1D averaging operator from nodes to cell-centers.""" return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") + +def avExtrap(n): + """Define 1D averaging operator from cell-centers to nodes.""" + Av = sp.spdiags((0.5*np.ones((n, 1))*[1, 1]).T, [-1, 0], n+1, n, format="csr") + sp.csr_matrix(([0.5,0.5],([0,n],[0,n-1])),shape=(n+1,n)) + return Av From dc07ea02fddebf876c529ff2ab480fdf07e9d4a3 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 13 Nov 2013 19:44:27 -0800 Subject: [PATCH 6/7] Make BC work in 2D and 3D --- SimPEG/mesh/DiffOperators.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 7176b8fd..1f453d58 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -287,8 +287,8 @@ class DiffOperators(object): 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]]) + V = self.aveCC2F*self.vol # Average volume between adjacent cells + self._cellGradBC = sdiag(S/V)*G return self._cellGradBC return locals() _cellGradBC = None From c6c0d3db4fb80fab1a5adfa4dbb53743d53d7e95 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 13 Nov 2013 23:00:51 -0800 Subject: [PATCH 7/7] Bug fix in cellGradx/y/z --- SimPEG/mesh/DiffOperators.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 1f453d58..19a3ab34 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -308,9 +308,9 @@ class DiffOperators(object): elif(self.dim == 3): G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC)) # Compute areas of cell faces & volumes - S = self.r(self.area, 'F','Fx', 'V') - V = self.vol - self._cellGradx = sdiag(S)*G1*sdiag(1/V) + V = self.aveCC2F*self.vol + L = self.r(self.area/V, 'F','Fx', 'V') + self._cellGradx = sdiag(L)*G1 return self._cellGradx return locals() cellGradx = property(**cellGradx()) @@ -328,9 +328,9 @@ class DiffOperators(object): elif(self.dim == 3): G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0])) # Compute areas of cell faces & volumes - S = self.r(self.area, 'F','Fy', 'V') - V = self.vol - self._cellGrady = sdiag(S)*G2*sdiag(1/V) + V = self.aveCC2F*self.vol + L = self.r(self.area/V, 'F','Fy', 'V') + self._cellGrady = sdiag(L)*G2 return self._cellGrady return locals() cellGrady = property(**cellGrady()) @@ -345,9 +345,9 @@ class DiffOperators(object): n = self.n G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0])) # Compute areas of cell faces & volumes - S = self.r(self.area, 'F','Fz', 'V') - V = self.vol - self._cellGradz = sdiag(S)*G3*sdiag(1/V) + V = self.aveCC2F*self.vol + L = self.r(self.area/V, 'F','Fz', 'V') + self._cellGradz = sdiag(L)*G3 return self._cellGradz return locals() cellGradz = property(**cellGradz())