From 082df8f1104b1240dd5a9d0ecf04f6eae9dd5198 Mon Sep 17 00:00:00 2001 From: SEOGI KANG Date: Wed, 17 Jul 2013 15:03:56 -0700 Subject: [PATCH 01/11] Testing differential operators (Div, Grad, Curl) --- SimPEG/getDiffop.py | 198 ++++++++++++++++++++++++++++++++++++++ SimPEG/sputils.py | 16 ++- SimPEG/tests/test_curl.py | 48 +++++++++ SimPEG/tests/test_div.py | 27 +++--- SimPEG/tests/test_grad.py | 46 +++++++++ 5 files changed, 319 insertions(+), 16 deletions(-) create mode 100644 SimPEG/getDiffop.py create mode 100644 SimPEG/tests/test_curl.py create mode 100644 SimPEG/tests/test_grad.py diff --git a/SimPEG/getDiffop.py b/SimPEG/getDiffop.py new file mode 100644 index 00000000..6debc348 --- /dev/null +++ b/SimPEG/getDiffop.py @@ -0,0 +1,198 @@ +import numpy as np +from scipy import sparse +from utils import mkvc +from sputils import ddx, sdiag, speye, kron3, spzeros, av + +def getvol(h): + """Construct cell volumes of the 3D model as 1d array.""" + # Cell sizes in each direction + h1 = h[0] + h2 = h[1] + h3 = h[2] + + # Compute cell volumes + v12 = h1.T*h2 + V = mkvc(v12.reshape(-1,1)*h3) + + return V + +def getarea(h): + """Construct face areas of the 3D model as 1d array.""" + # Cell sizes in each direction + h1 = h[0] + h2 = h[1] + h3 = h[2] + + # The number of cell centers in each direction + n1 = np.size(h1) + n2 = np.size(h2) + n3 = np.size(h3) + # Compute areas of cell faces + area1 = np.ones((n1+1,1))*mkvc(h2.T*h3) + area2 = h1.T*mkvc(np.ones((n2+1,1))*h3) + area3 = h1.T*mkvc(h2.T*np.ones(n3+1)) + area = np.concatenate((mkvc(area1), mkvc(area2), mkvc(area3)), axis=0) + + return area + +def getlength_e(h): + """Construct edge legnths of the 3D model as 1d array.""" + # Cell sizes in each direction + h1 = h[0] + h2 = h[1] + h3 = h[2] + + # The number of cell centers in each direction + n1 = np.size(h1) + n2 = np.size(h2) + n3 = np.size(h3) + # Compute areas of cell faces + l1 = h1.T*mkvc(np.ones((n2+1,1))*np.ones(n3+1)) + l2 = np.ones((n1+1,1))*mkvc(h2.T*np.ones(n3+1)) + l3 = np.ones((n1+1,1))*mkvc(np.ones((n2+1,1))*h3) + #l = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3))) + l = np.concatenate((mkvc(l1), mkvc(l2), mkvc(l3)), axis=0) + + return l + +def getDivMatrix(h): + """Construct the 3D divergence operator on Faces.""" + + # Cell sizes in each direction + h1 = h[0] + h2 = h[1] + h3 = h[2] + + # The number of cell centers in each direction + n1 = np.size(h1) + n2 = np.size(h2) + n3 = np.size(h3) + + # Compute areas of cell faces + S = getarea(h) + + # Compute cell volumes + V = getvol(h) + + # Compute divergence operator on faces + d1 = ddx(n1) + d2 = ddx(n2) + d3 = ddx(n3) + D1 = kron3(speye(n3), speye(n2), d1) + D2 = kron3(speye(n3), d2, speye(n1)) + D3 = kron3(d3, speye(n2), speye(n1)) + + D = sparse.hstack((D1, D2, D3), format="csr") + return sdiag(1/V)*D*sdiag(S) + +def getGradMatrix(h): + """Construct the 3D nodal gradient operator.""" + + # Cell sizes in each direction + h1 = h[0] + h2 = h[1] + h3 = h[2] + + # The number of cell centers in each direction + n1 = np.size(h1) + n2 = np.size(h2) + n3 = np.size(h3) + + # Compute lengths of cell edges + L = getlength_e(h) + + # Compute divergence operator on faces + d1 = ddx(n1) + d2 = ddx(n2) + d3 = ddx(n3) + D1 = kron3(speye(n3+1), speye(n2+1), d1) + D2 = kron3(speye(n3+1), d2, speye(n1+1)) + D3 = kron3(d3, speye(n2+1), speye(n1+1)) + + G = sparse.vstack((D1, D2, D3), format="csr") + return sdiag(1/L)*G + +def getCurlMatrix(h): + """Construct the 3D curl operator.""" + + # Cell sizes in each direction + h1 = h[0] + h2 = h[1] + h3 = h[2] + + # The number of cell centers in each direction + n1 = np.size(h1) + n2 = np.size(h2) + n3 = np.size(h3) + + # Compute lengths of cell edges + L = getlength_e(h) + + # Compute areas of cell faces + S = getarea(h) + + # Compute divergence operator on faces + d1 = ddx(n1) + d2 = ddx(n2) + d3 = ddx(n3) + + D32 = kron3(d3, speye(n2), speye(n1+1)) + D23 = kron3(speye(n3), d2, speye(n1+1)) + D31 = kron3(d3, speye(n2+1), speye(n1)) + D13 = kron3(speye(n3), speye(n2+1), d1) + D21 = kron3(speye(n3+1), d2, speye(n1)) + D12 = kron3(speye(n3+1), speye(n2), d1) + + O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) + O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) + O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) + + C = sparse.vstack((sparse.hstack((O1,-D32, D23)), + sparse.hstack((D31,O2, -D13)), + sparse.hstack((-D21,D12, O3))), format="csr") + + return sdiag(1/S)*(C*sdiag(L)) + +def getAverageMatrixF(h): + """Construct the 3D averaging operator on cell faces.""" + + # Cell sizes in each direction + h1 = h[0] + h2 = h[1] + h3 = h[2] + + # The number of cell centers in each direction + n1 = np.size(h1) + n2 = np.size(h2) + n3 = np.size(h3) + + av1 = av(n1) + av2 = av(n2) + av3 = av(n3) + + AvF = sparse.hstack(kron3(speye(n3), speye(n2), av1), + kron3(speye(n3), av2, speye(n3)), + kron3(av3, speye(n2), speye(n3)), format="csr") + return AvF + +def getAverageMatrixE(h): + """Construct the 3D averaging operator on cell edges.""" + + # Cell sizes in each direction + h1 = h[0] + h2 = h[1] + h3 = h[2] + + # The number of cell centers in each direction + n1 = np.size(h1) + n2 = np.size(h2) + n3 = np.size(h3) + + av1 = av(n1) + av2 = av(n2) + av3 = av(n3) + + AvE = sparse.hstack(kron3(av3, av2, speye(n1)), + kron3(av3, speye(n2), av1), + kron3(speye(n3), av2, av1), format="csr") + return AvE \ No newline at end of file diff --git a/SimPEG/sputils.py b/SimPEG/sputils.py index 204e71cc..241d9067 100644 --- a/SimPEG/sputils.py +++ b/SimPEG/sputils.py @@ -3,16 +3,24 @@ from scipy import sparse def ddx(n): """Define 1D derivatives""" - return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1) + return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1, format="csr") def sdiag(h): """Sparse diagonal matrix""" - return sparse.spdiags(h, 0, np.size(h), np.size(h)) + return sparse.spdiags(h, 0, np.size(h), np.size(h), format="csr") def speye(n): """Sparse identity""" - return sparse.identity(n) + return sparse.identity(n, format="csr") def kron3(A, B, C): """Two kron prods""" - return sparse.kron(sparse.kron(A, B), C) \ No newline at end of file + return sparse.kron(sparse.kron(A, B), C, format="csr") + +def spzeros(n1, n2): + """spzeros""" + return sparse.coo_matrix((n1, n2)).tocsr() + +def av(n): + """Define 1D averaging operator""" + return sparse.spdiags((0.5*np.ones((n+1,1))*[1,1]).T, [0,1], n, n+1, format="csr") \ No newline at end of file diff --git a/SimPEG/tests/test_curl.py b/SimPEG/tests/test_curl.py new file mode 100644 index 00000000..8175ab2f --- /dev/null +++ b/SimPEG/tests/test_curl.py @@ -0,0 +1,48 @@ +import numpy as np + +import sys +sys.path.append('../') +from TensorMesh import TensorMesh +from getDiffop import getCurlMatrix + + +err=0. +print '>> Test Curl operator' +for i in range(4): + icount=i+1 + nc = 2**icount + # Define the mesh + h1 = np.ones((1,nc))/nc + h2 = np.ones((1,nc))/nc + h3 = np.ones((1,nc))/nc + h = [h1, h2, h3] + x0 = np.zeros((3, 1)) + M = TensorMesh(h, x0) + #n = M.plotGrid() + + # Generate DIV matrix + CURL = getCurlMatrix(h) + #Test function + fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) + sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) + + Ex = fun(M.gridEx[:,1]) + Ey = fun(M.gridEy[:,2]) + Ez = fun(M.gridEz[:,0]) + E = np.concatenate((Ex,Ey,Ez)) + + Fx = sol(M.gridFx[:,2]) + Fy = sol(M.gridFy[:,0]) + Fz = sol(M.gridFz[:,1]) + curlE_anal = np.concatenate((Fx,Fy,Fz)) + + curlE = CURL*E + err = np.linalg.norm((curlE-curlE_anal), np.inf) + + if icount == 1: + print 'h | inf norm | error ratio' + print '---------------------------------------' + print '%6.4f | %8.2e |'% (h1[0,0], err) + else: + print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) + err_old = err \ No newline at end of file diff --git a/SimPEG/tests/test_div.py b/SimPEG/tests/test_div.py index bc0c1a9d..3ab8f6e2 100644 --- a/SimPEG/tests/test_div.py +++ b/SimPEG/tests/test_div.py @@ -5,16 +5,18 @@ sys.path.append('../') from TensorMesh import TensorMesh from getDIV import getDivMatrix, getarea, getvol -# Define the mesh + err=0. +print '>> Test face Divergence operator' for i in range(4): - icount=i+1; - nc = 2*icount; - h1 = np.pi/nc*np.ones((1,nc)) - h2 = np.pi/nc*np.ones((1,nc)) - h3 = np.pi/nc*np.ones((1,nc)) + icount=i+1 + nc = 2**icount + # Define the mesh + h1 = np.ones((1,nc))/nc + h2 = np.ones((1,nc))/nc + h3 = np.ones((1,nc))/nc h = [h1, h2, h3] - x0 = -np.pi/2*np.ones((3, 1)) + x0 = np.zeros((3, 1)) M = TensorMesh(h, x0) #n = M.plotGrid() @@ -34,12 +36,13 @@ for i in range(4): area = getarea(h) vol = getvol(h) - err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2) + #err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2) + err = np.linalg.norm((divF-divF_anal), np.inf) + if icount == 1: - err1 = err - print 'h | 2 norm | error ratio' + print 'h | inf norm | error ratio' print '---------------------------------------' print '%6.4f | %8.2e |'% (h1[0,0], err) else: - print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err1/err) - + print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) + err_old = err diff --git a/SimPEG/tests/test_grad.py b/SimPEG/tests/test_grad.py new file mode 100644 index 00000000..a433c28d --- /dev/null +++ b/SimPEG/tests/test_grad.py @@ -0,0 +1,46 @@ +import numpy as np + +import sys +sys.path.append('../') +from TensorMesh import TensorMesh +from getDiffop import getGradMatrix + + +err=0. +print '>> Test nodal Gradient operator' +for i in range(4): + icount=i+1 + nc = 2**icount + # Define the mesh + h1 = np.ones((1,nc))/nc + h2 = np.ones((1,nc))/nc + h3 = np.ones((1,nc))/nc + h = [h1, h2, h3] + x0 = np.zeros((3, 1)) + M = TensorMesh(h, x0) + #n = M.plotGrid() + + # Generate DIV matrix + GRAD = getGradMatrix(h) + #Test function + fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) + + phi = fun(M.gridN[:,0], M.gridN[:,1], M.gridN[:,2]) + gradE = GRAD*phi + + Ex = sol(M.gridEx[:,0]) + Ey = sol(M.gridEy[:,1]) + Ez = sol(M.gridEz[:,2]) + + gradE_anal = np.concatenate((Ex,Ey,Ez)) + err = np.linalg.norm((gradE-gradE_anal), np.inf) + + if icount == 1: + print 'h | inf norm | error ratio' + print '---------------------------------------' + print '%6.4f | %8.2e |'% (h1[0,0], err) + else: + print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) + err_old = err + From 3680295e5422576281ef361c50e347448ae70acf Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 17 Jul 2013 15:11:38 -0700 Subject: [PATCH 02/11] Deleted unused code. --- SimPEG/getDIV.py | 75 ------------------------------------------------ 1 file changed, 75 deletions(-) delete mode 100644 SimPEG/getDIV.py diff --git a/SimPEG/getDIV.py b/SimPEG/getDIV.py deleted file mode 100644 index 6c45c57c..00000000 --- a/SimPEG/getDIV.py +++ /dev/null @@ -1,75 +0,0 @@ -import numpy as np -from scipy import sparse -from utils import mkvc -from sputils import ddx, sdiag, speye, kron3 - -def getvol(h): - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # Compute cell volumes - v12 = h1.T*h2 - V = mkvc(v12.reshape(-1,1)*h3) - - return V - -def getarea(h): - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - # Compute areas of cell faces - area1 = np.ones((n1+1,1))*mkvc(h2.T*h3) - area2 = h1.T*mkvc(np.ones((n2+1,1))*h3) - area3 = h1.T*mkvc(h2.T*np.ones(n3+1)) - area = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3))) - - return area - -def getDivMatrix(h): - """Consturct the 3D divergence operator on Faces.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - # Compute areas of cell faces - #area1 = np.ones((n1+1,1))*mkvc(h2.T*h3) - #area2 = h1.T*mkvc(np.ones((n2+1,1))*h3) - #area3 = h1.T*mkvc(h2.T*np.ones(n3+1)) - #area = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3))) - area = getarea(h) - - S = sdiag(area) - - # Compute cell volumes - #v12 = h1.T*h2 - #V = mkvc(v12.reshape(-1,1)*h3) - V = getvol(h) - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - D1 = kron3(speye(n3), speye(n2), d1) - D2 = kron3(speye(n3), d2, speye(n1)) - D3 = kron3(d3, speye(n2), speye(n1)) - - D = sparse.hstack((sparse.hstack((D1, D2)), D3)) - return sdiag(1/V)*D*S - From ef557beb84ce3d6f39f216e011d2a0fdf0d24cc8 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 18 Jul 2013 10:41:58 -0700 Subject: [PATCH 03/11] Merged Seogi's code into the Tensor Mesh code base. This is under the DiffOperators class, which can be inherited with BaseMesh to any Mesh Object. Cell area and other dimension calculations are included in the TensorMesh class. Wrote unit tests for cell vol, area, and edges. --- SimPEG/DiffOperators.py | 156 +++++++++++++++++++++++++ SimPEG/TensorMesh.py | 79 ++++++++++++- SimPEG/__init__.py | 1 + SimPEG/getDiffop.py | 198 -------------------------------- SimPEG/tests/test_curl.py | 19 ++- SimPEG/tests/test_div.py | 19 ++- SimPEG/tests/test_grad.py | 21 ++-- SimPEG/tests/test_tensorMesh.py | 30 ++++- 8 files changed, 289 insertions(+), 234 deletions(-) create mode 100644 SimPEG/DiffOperators.py delete mode 100644 SimPEG/getDiffop.py diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py new file mode 100644 index 00000000..f240a44f --- /dev/null +++ b/SimPEG/DiffOperators.py @@ -0,0 +1,156 @@ +import numpy as np +from scipy import sparse +from sputils import ddx, sdiag, speye, kron3, spzeros, av + + +class DiffOperators(object): + """ + Class creates the differential operators that you need! + """ + def __init__(self): + raise Exception('You should use a Mesh class.') + + def DIV(): + doc = "Construct the 3D divergence operator on Faces." + + def fget(self): + if(self._DIV is None): + # The number of cell centers in each direction + n = [x.size for x in self.h] + # Compute divergence operator on faces + dd = [ddx(x) for x in n] + if(self.dim == 1): + D = dd[0] + elif(self.dim == 2): + D1 = sparse.kron(speye(n[1]), dd[0]) + D2 = sparse.kron(dd[1], speye(n[0])) + D = sparse.hstack((D1, D2), format="csr") + elif(self.dim == 3): + D1 = kron3(speye(n[2]), speye(n[1]), dd[0]) + D2 = kron3(speye(n[2]), dd[1], speye(n[0])) + D3 = kron3(dd[2], speye(n[1]), speye(n[0])) + D = sparse.hstack((D1, D2, D3), format="csr") + # Compute areas of cell faces + S = self.area + # Compute cell volumes + V = self.vol + self._DIV = sdiag(1/V)*D*sdiag(S) + + return self._DIV + return locals() + _DIV = None + DIV = property(**DIV()) + + def GRAD(): + doc = "Construct the 3D nodal gradient operator." + + def fget(self): + if(self._GRAD is None): + # The number of cell centers in each direction + n1 = np.size(self.hx) + n2 = np.size(self.hy) + n3 = np.size(self.hz) + + # Compute lengths of cell edges + L = self.edge + + # Compute divergence operator on faces + d1 = ddx(n1) + d2 = ddx(n2) + d3 = ddx(n3) + D1 = kron3(speye(n3+1), speye(n2+1), d1) + D2 = kron3(speye(n3+1), d2, speye(n1+1)) + D3 = kron3(d3, speye(n2+1), speye(n1+1)) + + G = sparse.vstack((D1, D2, D3), format="csr") + self._GRAD = sdiag(1/L)*G + return self._GRAD + return locals() + _GRAD = None + GRAD = property(**GRAD()) + + def CURL(): + doc = "Construct the 3D curl operator." + + def fget(self): + if(self._CURL is None): + # The number of cell centers in each direction + n1 = np.size(self.hx) + n2 = np.size(self.hy) + n3 = np.size(self.hz) + + # Compute lengths of cell edges + L = self.edge + + # Compute areas of cell faces + S = self.area + + # Compute divergence operator on faces + d1 = ddx(n1) + d2 = ddx(n2) + d3 = ddx(n3) + + D32 = kron3(d3, speye(n2), speye(n1+1)) + D23 = kron3(speye(n3), d2, speye(n1+1)) + D31 = kron3(d3, speye(n2+1), speye(n1)) + D13 = kron3(speye(n3), speye(n2+1), d1) + D21 = kron3(speye(n3+1), d2, speye(n1)) + D12 = kron3(speye(n3+1), speye(n2), d1) + + O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) + O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) + O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) + + C = sparse.vstack((sparse.hstack((O1, -D32, D23)), + sparse.hstack((D31, O2, -D13)), + sparse.hstack((-D21, D12, O3))), format="csr") + + self._CURL = sdiag(1/S)*(C*sdiag(L)) + return self._CURL + return locals() + _CURL = None + CURL = property(**CURL()) + + def AVE_F(): + doc = "Construct the 3D averaging operator on cell faces." + + def fget(self): + if(self._AVE_F is None): + # The number of cell centers in each direction + n1 = np.size(self.hx) + n2 = np.size(self.hy) + n3 = np.size(self.hz) + + av1 = av(n1) + av2 = av(n2) + av3 = av(n3) + + self._AVE_F = sparse.hstack(kron3(speye(n3), speye(n2), av1), + kron3(speye(n3), av2, speye(n3)), + kron3(av3, speye(n2), speye(n3)), format="csr") + return self._AVE_F + return locals() + _AVE_F = None + AVE_F = property(**AVE_F()) + + def AVE_E(): + doc = "Construct the 3D averaging operator on cell edges." + + def fget(self): + if(self._AVE_E is None): + # The number of cell centers in each direction + n1 = np.size(self.hx) + n2 = np.size(self.hy) + n3 = np.size(self.hz) + + av1 = av(n1) + av2 = av(n2) + av3 = av(n3) + + self._AVE_E = sparse.hstack(kron3(av3, av2, speye(n1)), + kron3(av3, speye(n2), av1), + kron3(speye(n3), av2, av1), format="csr") + return self._AVE_E + return locals() + _AVE_E = None + AVE_E = property(**AVE_E()) diff --git a/SimPEG/TensorMesh.py b/SimPEG/TensorMesh.py index 5f5e44b3..023362ae 100644 --- a/SimPEG/TensorMesh.py +++ b/SimPEG/TensorMesh.py @@ -1,10 +1,11 @@ import numpy as np from BaseMesh import BaseMesh from TensorView import TensorView -from utils import ndgrid +from DiffOperators import DiffOperators +from utils import ndgrid, mkvc -class TensorMesh(BaseMesh, TensorView): +class TensorMesh(BaseMesh, TensorView, DiffOperators): """ TensorMesh is a mesh class that deals with tensor product meshes. @@ -186,6 +187,80 @@ class TensorMesh(BaseMesh, TensorView): def getCellNumbering(self): pass + # --------------- Geometries --------------------- + def vol(): + doc = "Construct cell volumes of the 3D model as 1d array." + + def fget(self): + if(self._vol is None): + vh = [mkvc(x, 2) for x in self.h] + # Compute cell volumes + if(self.dim == 1): + self._vol = mkvc(vh[0]) + elif(self.dim == 2): + # Cell sizes in each direction + self._vol = mkvc(vh[0]*vh[1].T) + elif(self.dim == 3): + # Cell sizes in each direction + v12 = vh[0]*vh[1].T + self._vol = mkvc(mkvc(v12, 2)*vh[2].T) + return self._vol + return locals() + _vol = None + vol = property(**vol()) + + def area(): + doc = "Construct face areas of the 3D model as 1d array." + + def fget(self): + if(self._area is None): + # Ensure that we are working with column vectors + vh = [mkvc(x, 2) for x in self.h] + # The number of cell centers in each direction + n = [x.size for x in self.h] + # Compute areas of cell faces + if(self.dim == 1): + self._area = np.ones((n[0]+1, 1)) + elif(self.dim == 2): + area1 = np.ones((n[0]+1, 1))*vh[1].T + area2 = vh[0]*np.ones((n[1]+1, 1)).T + self._area = np.r_[mkvc(area1), mkvc(area2)] + elif(self.dim == 3): + area1 = np.ones((n[0]+1, 1))*mkvc(vh[1]*vh[2].T) + area2 = vh[0]*mkvc(np.ones((n[1]+1, 1))*vh[2].T) + area3 = vh[0]*mkvc(vh[1]*np.ones((n[2]+1, 1)).T) + self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)] + return self._area + return locals() + _area = None + area = property(**area()) + + def edge(): + doc = "Construct edge legnths of the 3D model as 1d array." + + def fget(self): + if(self._area is None): + # Ensure that we are working with column vectors + vh = [mkvc(x, 2) for x in self.h] + # The number of cell centers in each direction + n = [x.size for x in self.h] + # Compute edge lengths + if(self.dim == 1): + self._edge = mkvc(vh[0]) + elif(self.dim == 2): + l1 = vh[0]*np.ones((n[1]+1, 1)).T + l2 = np.ones((n[0]+1, 1))*vh[1].T + self._edge = np.r_[mkvc(l1), mkvc(l2)] + elif(self.dim == 3): + l1 = vh[0]*mkvc(np.ones((n[1]+1, 1))*np.ones((n[2]+1, 1)).T) + l2 = np.ones((n[0]+1, 1))*mkvc(vh[1]*np.ones((n[2]+1, 1)).T) + l3 = np.ones((n[0]+1, 1))*mkvc(np.ones((n[1]+1, 1))*vh[2].T) + self._edge = np.r_[mkvc(l1), mkvc(l2), mkvc(l3)] + return self._edge + return locals() + _edge = None + edge = property(**edge()) + if __name__ == '__main__': print('Welcome to tensor mesh!') diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 2d1b1e99..bcbba681 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1 +1,2 @@ from TensorMesh import TensorMesh +import utils diff --git a/SimPEG/getDiffop.py b/SimPEG/getDiffop.py deleted file mode 100644 index 6debc348..00000000 --- a/SimPEG/getDiffop.py +++ /dev/null @@ -1,198 +0,0 @@ -import numpy as np -from scipy import sparse -from utils import mkvc -from sputils import ddx, sdiag, speye, kron3, spzeros, av - -def getvol(h): - """Construct cell volumes of the 3D model as 1d array.""" - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # Compute cell volumes - v12 = h1.T*h2 - V = mkvc(v12.reshape(-1,1)*h3) - - return V - -def getarea(h): - """Construct face areas of the 3D model as 1d array.""" - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - # Compute areas of cell faces - area1 = np.ones((n1+1,1))*mkvc(h2.T*h3) - area2 = h1.T*mkvc(np.ones((n2+1,1))*h3) - area3 = h1.T*mkvc(h2.T*np.ones(n3+1)) - area = np.concatenate((mkvc(area1), mkvc(area2), mkvc(area3)), axis=0) - - return area - -def getlength_e(h): - """Construct edge legnths of the 3D model as 1d array.""" - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - # Compute areas of cell faces - l1 = h1.T*mkvc(np.ones((n2+1,1))*np.ones(n3+1)) - l2 = np.ones((n1+1,1))*mkvc(h2.T*np.ones(n3+1)) - l3 = np.ones((n1+1,1))*mkvc(np.ones((n2+1,1))*h3) - #l = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3))) - l = np.concatenate((mkvc(l1), mkvc(l2), mkvc(l3)), axis=0) - - return l - -def getDivMatrix(h): - """Construct the 3D divergence operator on Faces.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - # Compute areas of cell faces - S = getarea(h) - - # Compute cell volumes - V = getvol(h) - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - D1 = kron3(speye(n3), speye(n2), d1) - D2 = kron3(speye(n3), d2, speye(n1)) - D3 = kron3(d3, speye(n2), speye(n1)) - - D = sparse.hstack((D1, D2, D3), format="csr") - return sdiag(1/V)*D*sdiag(S) - -def getGradMatrix(h): - """Construct the 3D nodal gradient operator.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - # Compute lengths of cell edges - L = getlength_e(h) - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - D1 = kron3(speye(n3+1), speye(n2+1), d1) - D2 = kron3(speye(n3+1), d2, speye(n1+1)) - D3 = kron3(d3, speye(n2+1), speye(n1+1)) - - G = sparse.vstack((D1, D2, D3), format="csr") - return sdiag(1/L)*G - -def getCurlMatrix(h): - """Construct the 3D curl operator.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - # Compute lengths of cell edges - L = getlength_e(h) - - # Compute areas of cell faces - S = getarea(h) - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - - D32 = kron3(d3, speye(n2), speye(n1+1)) - D23 = kron3(speye(n3), d2, speye(n1+1)) - D31 = kron3(d3, speye(n2+1), speye(n1)) - D13 = kron3(speye(n3), speye(n2+1), d1) - D21 = kron3(speye(n3+1), d2, speye(n1)) - D12 = kron3(speye(n3+1), speye(n2), d1) - - O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) - O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) - O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) - - C = sparse.vstack((sparse.hstack((O1,-D32, D23)), - sparse.hstack((D31,O2, -D13)), - sparse.hstack((-D21,D12, O3))), format="csr") - - return sdiag(1/S)*(C*sdiag(L)) - -def getAverageMatrixF(h): - """Construct the 3D averaging operator on cell faces.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - av1 = av(n1) - av2 = av(n2) - av3 = av(n3) - - AvF = sparse.hstack(kron3(speye(n3), speye(n2), av1), - kron3(speye(n3), av2, speye(n3)), - kron3(av3, speye(n2), speye(n3)), format="csr") - return AvF - -def getAverageMatrixE(h): - """Construct the 3D averaging operator on cell edges.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - av1 = av(n1) - av2 = av(n2) - av3 = av(n3) - - AvE = sparse.hstack(kron3(av3, av2, speye(n1)), - kron3(av3, speye(n2), av1), - kron3(speye(n3), av2, av1), format="csr") - return AvE \ No newline at end of file diff --git a/SimPEG/tests/test_curl.py b/SimPEG/tests/test_curl.py index 8175ab2f..e4d9812c 100644 --- a/SimPEG/tests/test_curl.py +++ b/SimPEG/tests/test_curl.py @@ -3,7 +3,6 @@ import numpy as np import sys sys.path.append('../') from TensorMesh import TensorMesh -from getDiffop import getCurlMatrix err=0. @@ -11,7 +10,7 @@ print '>> Test Curl operator' for i in range(4): icount=i+1 nc = 2**icount - # Define the mesh + # Define the mesh h1 = np.ones((1,nc))/nc h2 = np.ones((1,nc))/nc h3 = np.ones((1,nc))/nc @@ -21,28 +20,28 @@ for i in range(4): #n = M.plotGrid() # Generate DIV matrix - CURL = getCurlMatrix(h) + CURL = M.CURL #Test function fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) - + Ex = fun(M.gridEx[:,1]) Ey = fun(M.gridEy[:,2]) Ez = fun(M.gridEz[:,0]) - E = np.concatenate((Ex,Ey,Ez)) + E = np.concatenate((Ex,Ey,Ez)) Fx = sol(M.gridFx[:,2]) Fy = sol(M.gridFy[:,0]) Fz = sol(M.gridFz[:,1]) - curlE_anal = np.concatenate((Fx,Fy,Fz)) + curlE_anal = np.concatenate((Fx,Fy,Fz)) - curlE = CURL*E + curlE = CURL*E err = np.linalg.norm((curlE-curlE_anal), np.inf) if icount == 1: - print 'h | inf norm | error ratio' - print '---------------------------------------' + print 'h | inf norm | error ratio' + print '---------------------------------------' print '%6.4f | %8.2e |'% (h1[0,0], err) else: print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) - err_old = err \ No newline at end of file + err_old = err \ No newline at end of file diff --git a/SimPEG/tests/test_div.py b/SimPEG/tests/test_div.py index 3ab8f6e2..35f4e408 100644 --- a/SimPEG/tests/test_div.py +++ b/SimPEG/tests/test_div.py @@ -3,7 +3,6 @@ import numpy as np import sys sys.path.append('../') from TensorMesh import TensorMesh -from getDIV import getDivMatrix, getarea, getvol err=0. @@ -11,7 +10,7 @@ print '>> Test face Divergence operator' for i in range(4): icount=i+1 nc = 2**icount - # Define the mesh + # Define the mesh h1 = np.ones((1,nc))/nc h2 = np.ones((1,nc))/nc h3 = np.ones((1,nc))/nc @@ -21,28 +20,26 @@ for i in range(4): #n = M.plotGrid() # Generate DIV matrix - DIV = getDivMatrix(h) - + DIV = M.DIV + #Test function fun = lambda x: np.sin(x) Fx = fun(M.gridFx[:,0]) Fy = fun(M.gridFy[:,1]) Fz = fun(M.gridFz[:,2]) - + F = np.concatenate((Fx,Fy,Fz)) divF = DIV*F sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) divF_anal = sol(M.gridCC[:,0], M.gridCC[:,1], M.gridCC[:,2]) - - area = getarea(h) - vol = getvol(h) + #err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2) err = np.linalg.norm((divF-divF_anal), np.inf) if icount == 1: - print 'h | inf norm | error ratio' - print '---------------------------------------' + print 'h | inf norm | error ratio' + print '---------------------------------------' print '%6.4f | %8.2e |'% (h1[0,0], err) else: print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) - err_old = err + err_old = err diff --git a/SimPEG/tests/test_grad.py b/SimPEG/tests/test_grad.py index a433c28d..245d9111 100644 --- a/SimPEG/tests/test_grad.py +++ b/SimPEG/tests/test_grad.py @@ -3,15 +3,14 @@ import numpy as np import sys sys.path.append('../') from TensorMesh import TensorMesh -from getDiffop import getGradMatrix err=0. -print '>> Test nodal Gradient operator' +print '>> Test nodal Gradient operator' for i in range(4): icount=i+1 nc = 2**icount - # Define the mesh + # Define the mesh h1 = np.ones((1,nc))/nc h2 = np.ones((1,nc))/nc h3 = np.ones((1,nc))/nc @@ -21,11 +20,11 @@ for i in range(4): #n = M.plotGrid() # Generate DIV matrix - GRAD = getGradMatrix(h) + GRAD = M.GRAD #Test function - fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) - + phi = fun(M.gridN[:,0], M.gridN[:,1], M.gridN[:,2]) gradE = GRAD*phi @@ -33,14 +32,14 @@ for i in range(4): Ey = sol(M.gridEy[:,1]) Ez = sol(M.gridEz[:,2]) - gradE_anal = np.concatenate((Ex,Ey,Ez)) + gradE_anal = np.concatenate((Ex,Ey,Ez)) err = np.linalg.norm((gradE-gradE_anal), np.inf) if icount == 1: - print 'h | inf norm | error ratio' - print '---------------------------------------' + print 'h | inf norm | error ratio' + print '---------------------------------------' print '%6.4f | %8.2e |'% (h1[0,0], err) else: print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) - err_old = err - + err_old = err + diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/tests/test_tensorMesh.py index e1b74374..2a6cb0a1 100644 --- a/SimPEG/tests/test_tensorMesh.py +++ b/SimPEG/tests/test_tensorMesh.py @@ -10,8 +10,9 @@ class TestSequenceFunctions(unittest.TestCase): def setUp(self): a = np.array([1, 1, 1]) b = np.array([1, 2]) - x0 = np.array([3, 5]) - self.mesh2 = TensorMesh([a, b], x0) + c = np.array([1, 4]) + self.mesh2 = TensorMesh([a, b], np.array([3, 5])) + self.mesh3 = TensorMesh([a, b, c]) def test_vectorN_2D(self): testNx = np.array([3, 4, 5, 6]) @@ -29,6 +30,31 @@ class TestSequenceFunctions(unittest.TestCase): ytest = np.all(self.mesh2.vectorCCy == testNy) self.assertTrue(xtest and ytest) + def test_area_3D(self): + test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2]) + t1 = np.all(self.mesh3.area == test_area) + self.assertTrue(t1) + + def test_vol_3D(self): + test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8]) + t1 = np.all(self.mesh3.vol == test_vol) + self.assertTrue(t1) + + def test_vol_2D(self): + test_vol = np.array([1, 1, 1, 2, 2, 2]) + t1 = np.all(self.mesh2.vol == test_vol) + self.assertTrue(t1) + + def test_edge_3D(self): + test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]) + t1 = np.all(self.mesh3.edge == test_edge) + self.assertTrue(t1) + + def test_edge_2D(self): + test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2]) + t1 = np.all(self.mesh2.edge == test_edge) + self.assertTrue(t1) + if __name__ == '__main__': unittest.main() From 1c48365497c0c8d782a46e00e5c7e4aa411e5060 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 18 Jul 2013 14:03:41 -0700 Subject: [PATCH 04/11] Rename based on conventions discussed in meeting. --- SimPEG/DiffOperators.py | 78 +++++++++++++++++++-------------------- SimPEG/tests/test_curl.py | 2 +- SimPEG/tests/test_div.py | 2 +- SimPEG/tests/test_grad.py | 2 +- 4 files changed, 42 insertions(+), 42 deletions(-) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index f240a44f..f9d1f0aa 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -10,15 +10,15 @@ class DiffOperators(object): def __init__(self): raise Exception('You should use a Mesh class.') - def DIV(): - doc = "Construct the 3D divergence operator on Faces." + def faceDiv(): + doc = "Construct the 3D Divergence operator on Faces." def fget(self): - if(self._DIV is None): + if(self._faceDiv is None): # The number of cell centers in each direction - n = [x.size for x in self.h] - # Compute divergence operator on faces - dd = [ddx(x) for x in n] + n = [hk.size for hk in self.h] + # Compute faceDivergence operator on faces + dd = [ddx(k) for k in n] if(self.dim == 1): D = dd[0] elif(self.dim == 2): @@ -34,18 +34,18 @@ class DiffOperators(object): S = self.area # Compute cell volumes V = self.vol - self._DIV = sdiag(1/V)*D*sdiag(S) + self._faceDiv = sdiag(1/V)*D*sdiag(S) - return self._DIV + return self._faceDiv return locals() - _DIV = None - DIV = property(**DIV()) + _faceDiv = None + faceDiv = property(**faceDiv()) - def GRAD(): + def nodalGrad(): doc = "Construct the 3D nodal gradient operator." def fget(self): - if(self._GRAD is None): + if(self._nodalGrad is None): # The number of cell centers in each direction n1 = np.size(self.hx) n2 = np.size(self.hy) @@ -63,17 +63,17 @@ class DiffOperators(object): D3 = kron3(d3, speye(n2+1), speye(n1+1)) G = sparse.vstack((D1, D2, D3), format="csr") - self._GRAD = sdiag(1/L)*G - return self._GRAD + self._nodalGrad = sdiag(1/L)*G + return self._nodalGrad return locals() - _GRAD = None - GRAD = property(**GRAD()) + _nodalGrad = None + nodalGrad = property(**nodalGrad()) - def CURL(): + def edgeCurl(): doc = "Construct the 3D curl operator." def fget(self): - if(self._CURL is None): + if(self._edgeCurl is None): # The number of cell centers in each direction n1 = np.size(self.hx) n2 = np.size(self.hy) @@ -105,17 +105,17 @@ class DiffOperators(object): sparse.hstack((D31, O2, -D13)), sparse.hstack((-D21, D12, O3))), format="csr") - self._CURL = sdiag(1/S)*(C*sdiag(L)) - return self._CURL + self._edgeCurl = sdiag(1/S)*(C*sdiag(L)) + return self._edgeCurl return locals() - _CURL = None - CURL = property(**CURL()) + _edgeCurl = None + edgeCurl = property(**edgeCurl()) - def AVE_F(): - doc = "Construct the 3D averaging operator on cell faces." + def faceAve(): + doc = "Construct the 3D averaging operator on cell faces to cell centers." def fget(self): - if(self._AVE_F is None): + if(self._faceAve is None): # The number of cell centers in each direction n1 = np.size(self.hx) n2 = np.size(self.hy) @@ -125,19 +125,19 @@ class DiffOperators(object): av2 = av(n2) av3 = av(n3) - self._AVE_F = sparse.hstack(kron3(speye(n3), speye(n2), av1), - kron3(speye(n3), av2, speye(n3)), - kron3(av3, speye(n2), speye(n3)), format="csr") - return self._AVE_F + self._faceAve = sparse.hstack(kron3(speye(n3), speye(n2), av1), + kron3(speye(n3), av2, speye(n3)), + kron3(av3, speye(n2), speye(n3)), format="csr") + return self._faceAve return locals() - _AVE_F = None - AVE_F = property(**AVE_F()) + _faceAve = None + faceAve = property(**faceAve()) - def AVE_E(): + def edgeAve(): doc = "Construct the 3D averaging operator on cell edges." def fget(self): - if(self._AVE_E is None): + if(self._edgeAve is None): # The number of cell centers in each direction n1 = np.size(self.hx) n2 = np.size(self.hy) @@ -147,10 +147,10 @@ class DiffOperators(object): av2 = av(n2) av3 = av(n3) - self._AVE_E = sparse.hstack(kron3(av3, av2, speye(n1)), - kron3(av3, speye(n2), av1), - kron3(speye(n3), av2, av1), format="csr") - return self._AVE_E + self._edgeAve = sparse.hstack(kron3(av3, av2, speye(n1)), + kron3(av3, speye(n2), av1), + kron3(speye(n3), av2, av1), format="csr") + return self._edgeAve return locals() - _AVE_E = None - AVE_E = property(**AVE_E()) + _edgeAve = None + edgeAve = property(**edgeAve()) diff --git a/SimPEG/tests/test_curl.py b/SimPEG/tests/test_curl.py index e4d9812c..d8ee0018 100644 --- a/SimPEG/tests/test_curl.py +++ b/SimPEG/tests/test_curl.py @@ -20,7 +20,7 @@ for i in range(4): #n = M.plotGrid() # Generate DIV matrix - CURL = M.CURL + CURL = M.edgeCurl #Test function fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) diff --git a/SimPEG/tests/test_div.py b/SimPEG/tests/test_div.py index 35f4e408..dd469a3f 100644 --- a/SimPEG/tests/test_div.py +++ b/SimPEG/tests/test_div.py @@ -20,7 +20,7 @@ for i in range(4): #n = M.plotGrid() # Generate DIV matrix - DIV = M.DIV + DIV = M.faceDiv #Test function fun = lambda x: np.sin(x) diff --git a/SimPEG/tests/test_grad.py b/SimPEG/tests/test_grad.py index 245d9111..58401dc7 100644 --- a/SimPEG/tests/test_grad.py +++ b/SimPEG/tests/test_grad.py @@ -20,7 +20,7 @@ for i in range(4): #n = M.plotGrid() # Generate DIV matrix - GRAD = M.GRAD + GRAD = M.nodalGrad #Test function fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) From f626cedfb8d2dd1d9f2038a72da3b0c7b4ebd1f3 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 19 Jul 2013 11:24:14 -0700 Subject: [PATCH 05/11] Bug fixes to do with array sizes. Incorporated eldads outer product code --- SimPEG/DiffOperators.py | 75 +++++++++++++++++++++------------------ SimPEG/TensorMesh.py | 44 +++++++++++------------ SimPEG/sputils.py | 39 +++++++------------- SimPEG/tests/test_curl.py | 13 ++++--- SimPEG/tests/test_div.py | 16 ++++----- SimPEG/tests/test_grad.py | 13 ++++--- 6 files changed, 96 insertions(+), 104 deletions(-) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index f9d1f0aa..a755662f 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -1,6 +1,16 @@ import numpy as np -from scipy import sparse -from sputils import ddx, sdiag, speye, kron3, spzeros, av +from scipy import sparse as sp +from sputils import sdiag, speye, kron3, spzeros + + +def ddx(n): + """Define 1D derivatives""" + return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr") + + +def av(n): + """Define 1D averaging operator""" + return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") class DiffOperators(object): @@ -16,20 +26,20 @@ class DiffOperators(object): def fget(self): if(self._faceDiv is None): # The number of cell centers in each direction - n = [hk.size for hk in self.h] + n = self.n # Compute faceDivergence operator on faces dd = [ddx(k) for k in n] if(self.dim == 1): D = dd[0] elif(self.dim == 2): - D1 = sparse.kron(speye(n[1]), dd[0]) - D2 = sparse.kron(dd[1], speye(n[0])) - D = sparse.hstack((D1, D2), format="csr") + D1 = sp.kron(speye(n[1]), dd[0]) + D2 = sp.kron(dd[1], speye(n[0])) + D = sp.hstack((D1, D2), format="csr") elif(self.dim == 3): D1 = kron3(speye(n[2]), speye(n[1]), dd[0]) D2 = kron3(speye(n[2]), dd[1], speye(n[0])) D3 = kron3(dd[2], speye(n[1]), speye(n[0])) - D = sparse.hstack((D1, D2, D3), format="csr") + D = sp.hstack((D1, D2, D3), format="csr") # Compute areas of cell faces S = self.area # Compute cell volumes @@ -62,7 +72,7 @@ class DiffOperators(object): D2 = kron3(speye(n3+1), d2, speye(n1+1)) D3 = kron3(d3, speye(n2+1), speye(n1+1)) - G = sparse.vstack((D1, D2, D3), format="csr") + G = sp.vstack((D1, D2, D3), format="csr") self._nodalGrad = sdiag(1/L)*G return self._nodalGrad return locals() @@ -101,9 +111,9 @@ class DiffOperators(object): O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) - C = sparse.vstack((sparse.hstack((O1, -D32, D23)), - sparse.hstack((D31, O2, -D13)), - sparse.hstack((-D21, D12, O3))), format="csr") + C = sp.vstack((sp.hstack((O1, -D32, D23)), + sp.hstack((D31, O2, -D13)), + sp.hstack((-D21, D12, O3))), format="csr") self._edgeCurl = sdiag(1/S)*(C*sdiag(L)) return self._edgeCurl @@ -116,18 +126,16 @@ class DiffOperators(object): def fget(self): if(self._faceAve is None): - # The number of cell centers in each direction - n1 = np.size(self.hx) - n2 = np.size(self.hy) - n3 = np.size(self.hz) - - av1 = av(n1) - av2 = av(n2) - av3 = av(n3) - - self._faceAve = sparse.hstack(kron3(speye(n3), speye(n2), av1), - kron3(speye(n3), av2, speye(n3)), - kron3(av3, speye(n2), speye(n3)), format="csr") + n = self.n + if(self.dim == 1): + self._faceAve = av(n[0]) + elif(self.dim == 2): + self._faceAve = sp.hstack((sp.kron(speye(n[1]), av(n[0])), + sp.kron(av(n[1]), speye(n[0]))), format="csr") + elif(self.dim == 3): + self._faceAve = 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") return self._faceAve return locals() _faceAve = None @@ -139,17 +147,16 @@ class DiffOperators(object): def fget(self): if(self._edgeAve is None): # The number of cell centers in each direction - n1 = np.size(self.hx) - n2 = np.size(self.hy) - n3 = np.size(self.hz) - - av1 = av(n1) - av2 = av(n2) - av3 = av(n3) - - self._edgeAve = sparse.hstack(kron3(av3, av2, speye(n1)), - kron3(av3, speye(n2), av1), - kron3(speye(n3), av2, av1), format="csr") + n = self.n + if(self.dim == 1): + raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') + elif(self.dim == 2): + self._edgeAve = sp.hstack((sp.kron(av(n[1]), speye(n[0])), + sp.kron(speye(n[1]), av(n[0]))), format="csr") + elif(self.dim == 3): + self._edgeAve = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), + kron3(av(n[2]), speye(n[1]), av(n[0])), + kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") return self._edgeAve return locals() _edgeAve = None diff --git a/SimPEG/TensorMesh.py b/SimPEG/TensorMesh.py index 023362ae..8395c54d 100644 --- a/SimPEG/TensorMesh.py +++ b/SimPEG/TensorMesh.py @@ -22,15 +22,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators): """ def __init__(self, h, x0=None): - super(TensorMesh, self).__init__(np.array([len(x) for x in h]), x0) + 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): 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) # Ensure h contains 1D vectors - self._h = [x.ravel() for x in h] + self._h = [mkvc(x) for x in h] def h(): doc = "h is a list containing the cell widths of the tensor mesh in each dimension." @@ -193,17 +194,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators): def fget(self): if(self._vol is None): - vh = [mkvc(x, 2) for x in self.h] + vh = self.h # Compute cell volumes if(self.dim == 1): self._vol = mkvc(vh[0]) elif(self.dim == 2): # Cell sizes in each direction - self._vol = mkvc(vh[0]*vh[1].T) + self._vol = mkvc(np.outer(vh[0], vh[1])) elif(self.dim == 3): # Cell sizes in each direction - v12 = vh[0]*vh[1].T - self._vol = mkvc(mkvc(v12, 2)*vh[2].T) + self._vol = mkvc(np.outer(mkvc(np.outer(vh[0], vh[1])), vh[2])) return self._vol return locals() _vol = None @@ -215,20 +215,20 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators): def fget(self): if(self._area is None): # Ensure that we are working with column vectors - vh = [mkvc(x, 2) for x in self.h] + vh = self.h # The number of cell centers in each direction - n = [x.size for x in self.h] + n = self.n # Compute areas of cell faces if(self.dim == 1): - self._area = np.ones((n[0]+1, 1)) + self._area = np.ones(n[0]+1) elif(self.dim == 2): - area1 = np.ones((n[0]+1, 1))*vh[1].T - area2 = vh[0]*np.ones((n[1]+1, 1)).T + area1 = np.outer(np.ones(n[0]+1), vh[1]) + area2 = np.outer(vh[0], np.ones(n[1]+1)) self._area = np.r_[mkvc(area1), mkvc(area2)] elif(self.dim == 3): - area1 = np.ones((n[0]+1, 1))*mkvc(vh[1]*vh[2].T) - area2 = vh[0]*mkvc(np.ones((n[1]+1, 1))*vh[2].T) - area3 = vh[0]*mkvc(vh[1]*np.ones((n[2]+1, 1)).T) + area1 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], vh[2]))) + area2 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + area3 = np.outer(vh[0], mkvc(np.outer(vh[1], np.ones(n[2]+1)))) self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)] return self._area return locals() @@ -239,22 +239,22 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators): doc = "Construct edge legnths of the 3D model as 1d array." def fget(self): - if(self._area is None): + if(self._edge is None): # Ensure that we are working with column vectors - vh = [mkvc(x, 2) for x in self.h] + vh = self.h # The number of cell centers in each direction - n = [x.size for x in self.h] + n = self.n # Compute edge lengths if(self.dim == 1): self._edge = mkvc(vh[0]) elif(self.dim == 2): - l1 = vh[0]*np.ones((n[1]+1, 1)).T - l2 = np.ones((n[0]+1, 1))*vh[1].T + l1 = np.outer(vh[0], np.ones(n[1]+1)) + l2 = np.outer(np.ones(n[0]+1), vh[1]) self._edge = np.r_[mkvc(l1), mkvc(l2)] elif(self.dim == 3): - l1 = vh[0]*mkvc(np.ones((n[1]+1, 1))*np.ones((n[2]+1, 1)).T) - l2 = np.ones((n[0]+1, 1))*mkvc(vh[1]*np.ones((n[2]+1, 1)).T) - l3 = np.ones((n[0]+1, 1))*mkvc(np.ones((n[1]+1, 1))*vh[2].T) + l1 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) + l2 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + l3 = np.outer(np.ones(n[0]+1), mkvc(np.outer(np.ones(n[1]+1), vh[2]))) self._edge = np.r_[mkvc(l1), mkvc(l2), mkvc(l3)] return self._edge return locals() diff --git a/SimPEG/sputils.py b/SimPEG/sputils.py index f1e9a677..79bbdf04 100644 --- a/SimPEG/sputils.py +++ b/SimPEG/sputils.py @@ -1,43 +1,30 @@ import numpy as np -from scipy import sparse -from numpy import ones +from scipy import sparse as sp -def ddx(n): - """Define 1D derivatives""" - return sparse.spdiags((np.ones((n+1,1))*[-1,1]).T, [0,1], n, n+1, format="csr") - def sdiag(h): """Sparse diagonal matrix""" - return sparse.spdiags(h, 0, np.size(h), np.size(h), format="csr") + return sp.spdiags(h, 0, np.size(h), np.size(h), format="csr") + def speye(n): """Sparse identity""" - return sparse.identity(n, format="csr") + return sp.identity(n, format="csr") + def kron3(A, B, C): """Two kron prods""" - return sparse.kron(sparse.kron(A, B), C, format="csr") + return sp.kron(sp.kron(A, B), C, format="csr") + def spzeros(n1, n2): """spzeros""" - return sparse.coo_matrix((n1, n2)).tocsr() + return sp.coo_matrix((n1, n2)).tocsr() -def av(n): - """Define 1D averaging operator""" - return sparse.spdiags((0.5*np.ones((n+1,1))*[1,1]).T, [0,1], n, n+1, format="csr") - -def av(n): - """Define 1D average""" - return 0.5*(sparse.spdiags(ones(n+1), 0, n, n+1) + sparse.spdiags(ones(n+1), 1, n, n+1)) - -def spzeros(n1, n2): - """spzeros""" - return sparse.coo_matrix((n1, n2)) def appendBottom(A, B): """append on bottom""" - C = sparse.vstack((A, B)) + C = sp.vstack((A, B)) C = C.tocsr() return C @@ -51,7 +38,7 @@ def appendBottom3(A, B, C): def appendRight(A, B): """append on right""" - C = sparse.hstack((A, B)) + C = sp.hstack((A, B)) C = C.tocsr() return C @@ -65,9 +52,9 @@ def appendRight3(A, B, C): def blkDiag(A, B): """blockdigonal""" - O12 = sparse.coo_matrix((np.shape(A)[0], np.shape(B)[1])) - O21 = sparse.coo_matrix((np.shape(B)[0], np.shape(A)[1])) - C = sparse.vstack((sparse.hstack((A, O12)), sparse.hstack((O21, B)))) + O12 = sp.coo_matrix((np.shape(A)[0], np.shape(B)[1])) + O21 = sp.coo_matrix((np.shape(B)[0], np.shape(A)[1])) + C = sp.vstack((sp.hstack((A, O12)), sp.hstack((O21, B)))) C = C.tocsr() return C diff --git a/SimPEG/tests/test_curl.py b/SimPEG/tests/test_curl.py index d8ee0018..e90f8cfa 100644 --- a/SimPEG/tests/test_curl.py +++ b/SimPEG/tests/test_curl.py @@ -11,12 +11,11 @@ for i in range(4): icount=i+1 nc = 2**icount # Define the mesh - h1 = np.ones((1,nc))/nc - h2 = np.ones((1,nc))/nc - h3 = np.ones((1,nc))/nc + h1 = np.ones(nc)/nc + h2 = np.ones(nc)/nc + h3 = np.ones(nc)/nc h = [h1, h2, h3] - x0 = np.zeros((3, 1)) - M = TensorMesh(h, x0) + M = TensorMesh(h) #n = M.plotGrid() # Generate DIV matrix @@ -41,7 +40,7 @@ for i in range(4): if icount == 1: print 'h | inf norm | error ratio' print '---------------------------------------' - print '%6.4f | %8.2e |'% (h1[0,0], err) + print '%6.4f | %8.2e |'% (h1[0], err) else: - print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) + print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err) err_old = err \ No newline at end of file diff --git a/SimPEG/tests/test_div.py b/SimPEG/tests/test_div.py index dd469a3f..53ca506b 100644 --- a/SimPEG/tests/test_div.py +++ b/SimPEG/tests/test_div.py @@ -11,9 +11,9 @@ for i in range(4): icount=i+1 nc = 2**icount # Define the mesh - h1 = np.ones((1,nc))/nc - h2 = np.ones((1,nc))/nc - h3 = np.ones((1,nc))/nc + h1 = np.ones(nc)/nc + h2 = np.ones(nc)/nc + h3 = np.ones(nc)/nc h = [h1, h2, h3] x0 = np.zeros((3, 1)) M = TensorMesh(h, x0) @@ -24,9 +24,9 @@ for i in range(4): #Test function fun = lambda x: np.sin(x) - Fx = fun(M.gridFx[:,0]) - Fy = fun(M.gridFy[:,1]) - Fz = fun(M.gridFz[:,2]) + Fx = fun(M.gridFx[:, 0]) + Fy = fun(M.gridFy[:, 1]) + Fz = fun(M.gridFz[:, 2]) F = np.concatenate((Fx,Fy,Fz)) divF = DIV*F @@ -39,7 +39,7 @@ for i in range(4): if icount == 1: print 'h | inf norm | error ratio' print '---------------------------------------' - print '%6.4f | %8.2e |'% (h1[0,0], err) + print '%6.4f | %8.2e |'% (h1[0], err) else: - print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) + print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err) err_old = err diff --git a/SimPEG/tests/test_grad.py b/SimPEG/tests/test_grad.py index 58401dc7..fbc36aeb 100644 --- a/SimPEG/tests/test_grad.py +++ b/SimPEG/tests/test_grad.py @@ -11,12 +11,11 @@ for i in range(4): icount=i+1 nc = 2**icount # Define the mesh - h1 = np.ones((1,nc))/nc - h2 = np.ones((1,nc))/nc - h3 = np.ones((1,nc))/nc + h1 = np.ones(nc)/nc + h2 = np.ones(nc)/nc + h3 = np.ones(nc)/nc h = [h1, h2, h3] - x0 = np.zeros((3, 1)) - M = TensorMesh(h, x0) + M = TensorMesh(h) #n = M.plotGrid() # Generate DIV matrix @@ -38,8 +37,8 @@ for i in range(4): if icount == 1: print 'h | inf norm | error ratio' print '---------------------------------------' - print '%6.4f | %8.2e |'% (h1[0,0], err) + print '%6.4f | %8.2e |'% (h1[0], err) else: - print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) + print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err) err_old = err From ffb78fea10e47191f62804365800037b1dbf0a82 Mon Sep 17 00:00:00 2001 From: ehaber99 Date: Fri, 19 Jul 2013 14:43:15 -0700 Subject: [PATCH 06/11] added mass matrix --- SimPEG/DiffOperators.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index a755662f..0fb36c2b 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -161,3 +161,23 @@ class DiffOperators(object): return locals() _edgeAve = None edgeAve = property(**edgeAve()) + + + +def getEdgeMassMatrix(h,sigma): + # mass matix for products of edge functions w'*M(sigma)*e + + Av = getEdgeToCellAverge(h) + v = getVol(h) + sigma = mkvc(sigma) + + return sdiag(Av.T*(v*sigma)) + +def getFaceMassMatrix(h,sigma): + # mass matix for products of edge functions w'*M(sigma)*e + + Av = getFaceToCellAverge(h) + v = getVol(h) + sigma = mkvc(sigma) + + return sdiag(Av.T*(v*sigma)) \ No newline at end of file From 906d7d210f41dfa211daafbbde156abf35679910 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 19 Jul 2013 14:51:23 -0700 Subject: [PATCH 07/11] Cell Centered Grad, and Div in 3D --- SimPEG/DiffOperators.py | 130 ++++++++++++++++++++++++++++++++-------- 1 file changed, 104 insertions(+), 26 deletions(-) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index a755662f..064c4844 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -4,10 +4,40 @@ from sputils import sdiag, speye, kron3, spzeros def ddx(n): - """Define 1D derivatives""" + """Define 1D derivatives, inner, this means we go from n+1 to n+1""" return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr") +def checkBC(bc): + if(type(bc) is str): + bc = [bc, bc] + assert type(bc) is list, 'bc must be a list' + assert len(bc) == 2, 'bc must have two elements' + + for bc_i in bc: + assert type(bc_i) is str, "each bc must be a string" + assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'" + return bc + + +def ddxCellGrad(n, bc): + """Define 1D derivatives, outer, this means we go from n to n+1""" + bc = checkBC(bc) + + D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr") + # Set the first side + if(bc[0] == 'dirichlet'): + D[0, 0] = 2 + elif(bc[0] == 'neumann'): + D[0, 0] = 0 + # Set the second side + if(bc[1] == 'dirichlet'): + D[-1, -1] = -2 + elif(bc[1] == 'neumann'): + D[-1, -1] = 0 + return D + + def av(n): """Define 1D averaging operator""" return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") @@ -28,21 +58,19 @@ class DiffOperators(object): # The number of cell centers in each direction n = self.n # Compute faceDivergence operator on faces - dd = [ddx(k) for k in n] if(self.dim == 1): - D = dd[0] + D = ddx(n[0]) elif(self.dim == 2): - D1 = sp.kron(speye(n[1]), dd[0]) - D2 = sp.kron(dd[1], speye(n[0])) + D1 = sp.kron(speye(n[1]), ddx(n[0])) + D2 = sp.kron(ddx(n[1]), speye(n[0])) D = sp.hstack((D1, D2), format="csr") elif(self.dim == 3): - D1 = kron3(speye(n[2]), speye(n[1]), dd[0]) - D2 = kron3(speye(n[2]), dd[1], speye(n[0])) - D3 = kron3(dd[2], speye(n[1]), speye(n[0])) + D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0])) + D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0])) + D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0])) D = sp.hstack((D1, D2, D3), format="csr") - # Compute areas of cell faces + # Compute areas of cell faces & volumes S = self.area - # Compute cell volumes V = self.vol self._faceDiv = sdiag(1/V)*D*sdiag(S) @@ -57,28 +85,78 @@ class DiffOperators(object): def fget(self): if(self._nodalGrad is None): # The number of cell centers in each direction - n1 = np.size(self.hx) - n2 = np.size(self.hy) - n3 = np.size(self.hz) - + n = self.n + # Compute divergence operator on faces + if(self.dim == 1): + G = ddx(n[0]) + elif(self.dim == 2): + D1 = sp.kron(speye(n[1]+1), ddx(n[0])) + D2 = sp.kron(ddx(n[1]), speye(n[0]+1)) + G = sp.vstack((D1, D2), format="csr") + elif(self.dim == 3): + D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0])) + D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1)) + D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1)) + G = sp.vstack((D1, D2, D3), format="csr") # Compute lengths of cell edges L = self.edge - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - D1 = kron3(speye(n3+1), speye(n2+1), d1) - D2 = kron3(speye(n3+1), d2, speye(n1+1)) - D3 = kron3(d3, speye(n2+1), speye(n1+1)) - - G = sp.vstack((D1, D2, D3), format="csr") self._nodalGrad = sdiag(1/L)*G return self._nodalGrad return locals() _nodalGrad = None nodalGrad = property(**nodalGrad()) + def setCellGradBC(self, BC): + """ + e.g. + + BC = 'neumann' + BC = ['neumann', 'dirichlet', 'neumann'] + BC = [['neumann', 'dirichlet'], 'dirichlet', 'neumann'] + + """ + if(type(BC) is str): + BC = [BC for _ in self.n] # Repeat the str self.dim times + elif(type(BC) is list): + assert len(BC) == self.dim, 'BC list must be the size of your mesh' + else: + raise Exception("BC must be a str or a list.") + + 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 + return BC + _cellGradBC = '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) + n = self.n + if(self.dim == 1): + G = ddxCellGrad(n[0], BC[0]) + elif(self.dim == 2): + G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0])) + G2 = sp.kron(ddxCellGrad(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]), ddxCellGrad(n[0], BC[0])) + G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0])) + G3 = kron3(ddxCellGrad(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._cellGrad = sdiag(S)*G*sdiag(1/V) + return self._cellGrad + return locals() + _cellGrad = None + cellGrad = property(**cellGrad()) + def edgeCurl(): doc = "Construct the 3D curl operator." @@ -122,7 +200,7 @@ class DiffOperators(object): edgeCurl = property(**edgeCurl()) def faceAve(): - doc = "Construct the 3D averaging operator on cell faces to cell centers." + doc = "Construct the averaging operator on cell faces to cell centers." def fget(self): if(self._faceAve is None): @@ -142,7 +220,7 @@ class DiffOperators(object): faceAve = property(**faceAve()) def edgeAve(): - doc = "Construct the 3D averaging operator on cell edges." + doc = "Construct the averaging operator on cell edges to cell centers." def fget(self): if(self._edgeAve is None): From b1f41c14a50aeb979c3e460df40829a49a679a60 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 19 Jul 2013 15:51:53 -0700 Subject: [PATCH 08/11] Mass Matrices. Merged eldad's code. --- SimPEG/DiffOperators.py | 31 +++++++++++++------------------ 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index 7c194e76..bec2492e 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -1,6 +1,7 @@ import numpy as np from scipy import sparse as sp from sputils import sdiag, speye, kron3, spzeros +from utils import mkvc def ddx(n): @@ -240,22 +241,16 @@ class DiffOperators(object): _edgeAve = None edgeAve = property(**edgeAve()) + def getEdgeMass(self, materialProp=None): + """mass matix for products of edge functions w'*M(materialProp)*e""" + if(materialProp is None): + materialProp = np.ones(self.nC) + Av = self.edgeAve + return sdiag(Av.T * (self.vol * mkvc(materialProp))) - -def getEdgeMassMatrix(h,sigma): - # mass matix for products of edge functions w'*M(sigma)*e - - Av = getEdgeToCellAverge(h) - v = getVol(h) - sigma = mkvc(sigma) - - return sdiag(Av.T*(v*sigma)) - -def getFaceMassMatrix(h,sigma): - # mass matix for products of edge functions w'*M(sigma)*e - - Av = getFaceToCellAverge(h) - v = getVol(h) - sigma = mkvc(sigma) - - return sdiag(Av.T*(v*sigma)) + def getFaceMass(self, materialProp=None): + """mass matix for products of edge functions w'*M(materialProp)*e""" + if(materialProp is None): + materialProp = np.ones(self.nC) + Av = self.faceAve + return sdiag(Av.T*(self.vol*mkvc(materialProp))) From d0efdb509c094f093b8a5059c71f5ded80a66434 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 22 Jul 2013 12:26:17 -0700 Subject: [PATCH 09/11] Created an OrderTest class This creates a testable framework so that we can easily test order of convergence on things like operators by just writing the guts of the code. Merged Seogi's code into this framework. --- SimPEG/tests/OrderTest.py | 50 ++++++++++++++++ SimPEG/tests/test_curl.py | 46 --------------- SimPEG/tests/test_div.py | 45 --------------- SimPEG/tests/test_grad.py | 44 --------------- SimPEG/tests/test_operatorOrders.py | 88 +++++++++++++++++++++++++++++ 5 files changed, 138 insertions(+), 135 deletions(-) create mode 100644 SimPEG/tests/OrderTest.py delete mode 100644 SimPEG/tests/test_curl.py delete mode 100644 SimPEG/tests/test_div.py delete mode 100644 SimPEG/tests/test_grad.py create mode 100644 SimPEG/tests/test_operatorOrders.py diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py new file mode 100644 index 00000000..68837f7a --- /dev/null +++ b/SimPEG/tests/OrderTest.py @@ -0,0 +1,50 @@ +import sys +sys.path.append('../../') +from SimPEG import TensorMesh +import numpy as np +import unittest + + +class OrderTest(unittest.TestCase): + """Order test sets up the basics for testing order of decrease for a function on a mesh.""" + + name = "Order Test" + expectedOrder = 2 + meshSizes = [4, 8, 16, 32] + + def setupMesh(self, nc): + # Define the mesh + h1 = np.ones(nc)/nc + h2 = np.ones(nc)/nc + h3 = np.ones(nc)/nc + h = [h1, h2, h3] + self.M = TensorMesh(h) + + def getError(self): + """Overwrite this function with the guts of the test.""" + return 1. + + def orderTest(self): + order = [] + err_old = 0. + nc_old = 0. + for ii, nc in enumerate(self.meshSizes): + self.setupMesh(nc) + err = self.getError() + if ii == 0: + print '' + print 'Testing order of: ' + self.name + print ' h | inf norm | ratio | order' + print '------------------------------------------' + print '%4i | %8.2e |' % (nc, err) + else: + order.append(np.log(err/err_old)/np.log(float(nc_old)/float(nc))) + print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1]) + err_old = err + nc_old = nc + + self.assertTrue(len(np.where(np.array(order) > 0.9*self.expectedOrder)[0]) > np.floor(0.75*len(order))) + + +if __name__ == '__main__': + unittest.main() diff --git a/SimPEG/tests/test_curl.py b/SimPEG/tests/test_curl.py deleted file mode 100644 index e90f8cfa..00000000 --- a/SimPEG/tests/test_curl.py +++ /dev/null @@ -1,46 +0,0 @@ -import numpy as np - -import sys -sys.path.append('../') -from TensorMesh import TensorMesh - - -err=0. -print '>> Test Curl operator' -for i in range(4): - icount=i+1 - nc = 2**icount - # Define the mesh - h1 = np.ones(nc)/nc - h2 = np.ones(nc)/nc - h3 = np.ones(nc)/nc - h = [h1, h2, h3] - M = TensorMesh(h) - #n = M.plotGrid() - - # Generate DIV matrix - CURL = M.edgeCurl - #Test function - fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) - sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) - - Ex = fun(M.gridEx[:,1]) - Ey = fun(M.gridEy[:,2]) - Ez = fun(M.gridEz[:,0]) - E = np.concatenate((Ex,Ey,Ez)) - - Fx = sol(M.gridFx[:,2]) - Fy = sol(M.gridFy[:,0]) - Fz = sol(M.gridFz[:,1]) - curlE_anal = np.concatenate((Fx,Fy,Fz)) - - curlE = CURL*E - err = np.linalg.norm((curlE-curlE_anal), np.inf) - - if icount == 1: - print 'h | inf norm | error ratio' - print '---------------------------------------' - print '%6.4f | %8.2e |'% (h1[0], err) - else: - print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err) - err_old = err \ No newline at end of file diff --git a/SimPEG/tests/test_div.py b/SimPEG/tests/test_div.py deleted file mode 100644 index 53ca506b..00000000 --- a/SimPEG/tests/test_div.py +++ /dev/null @@ -1,45 +0,0 @@ -import numpy as np - -import sys -sys.path.append('../') -from TensorMesh import TensorMesh - - -err=0. -print '>> Test face Divergence operator' -for i in range(4): - icount=i+1 - nc = 2**icount - # Define the mesh - h1 = np.ones(nc)/nc - h2 = np.ones(nc)/nc - h3 = np.ones(nc)/nc - h = [h1, h2, h3] - x0 = np.zeros((3, 1)) - M = TensorMesh(h, x0) - #n = M.plotGrid() - - # Generate DIV matrix - DIV = M.faceDiv - - #Test function - fun = lambda x: np.sin(x) - Fx = fun(M.gridFx[:, 0]) - Fy = fun(M.gridFy[:, 1]) - Fz = fun(M.gridFz[:, 2]) - - F = np.concatenate((Fx,Fy,Fz)) - divF = DIV*F - sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) - divF_anal = sol(M.gridCC[:,0], M.gridCC[:,1], M.gridCC[:,2]) - - #err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2) - err = np.linalg.norm((divF-divF_anal), np.inf) - - if icount == 1: - print 'h | inf norm | error ratio' - print '---------------------------------------' - print '%6.4f | %8.2e |'% (h1[0], err) - else: - print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err) - err_old = err diff --git a/SimPEG/tests/test_grad.py b/SimPEG/tests/test_grad.py deleted file mode 100644 index fbc36aeb..00000000 --- a/SimPEG/tests/test_grad.py +++ /dev/null @@ -1,44 +0,0 @@ -import numpy as np - -import sys -sys.path.append('../') -from TensorMesh import TensorMesh - - -err=0. -print '>> Test nodal Gradient operator' -for i in range(4): - icount=i+1 - nc = 2**icount - # Define the mesh - h1 = np.ones(nc)/nc - h2 = np.ones(nc)/nc - h3 = np.ones(nc)/nc - h = [h1, h2, h3] - M = TensorMesh(h) - #n = M.plotGrid() - - # Generate DIV matrix - GRAD = M.nodalGrad - #Test function - fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) - sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) - - phi = fun(M.gridN[:,0], M.gridN[:,1], M.gridN[:,2]) - gradE = GRAD*phi - - Ex = sol(M.gridEx[:,0]) - Ey = sol(M.gridEy[:,1]) - Ez = sol(M.gridEz[:,2]) - - gradE_anal = np.concatenate((Ex,Ey,Ez)) - err = np.linalg.norm((gradE-gradE_anal), np.inf) - - if icount == 1: - print 'h | inf norm | error ratio' - print '---------------------------------------' - print '%6.4f | %8.2e |'% (h1[0], err) - else: - print '%6.4f | %8.2e | %6.4f' % (h1[0], err, err_old/err) - err_old = err - diff --git a/SimPEG/tests/test_operatorOrders.py b/SimPEG/tests/test_operatorOrders.py new file mode 100644 index 00000000..9bb9fa03 --- /dev/null +++ b/SimPEG/tests/test_operatorOrders.py @@ -0,0 +1,88 @@ +import numpy as np +from OrderTest import OrderTest +import unittest + + +class TestCurl(OrderTest): + + name = "Curl" + + def getError(self): + fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) + sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) + + Ex = fun(self.M.gridEx[:, 1]) + Ey = fun(self.M.gridEy[:, 2]) + Ez = fun(self.M.gridEz[:, 0]) + E = np.concatenate((Ex, Ey, Ez)) + + Fx = sol(self.M.gridFx[:, 2]) + Fy = sol(self.M.gridFy[:, 0]) + Fz = sol(self.M.gridFz[:, 1]) + curlE_anal = np.concatenate((Fx, Fy, Fz)) + + # Generate DIV matrix + CURL = self.M.edgeCurl + + curlE = CURL*E + err = np.linalg.norm((curlE-curlE_anal), np.inf) + return err + + def test_order(self): + self.orderTest() + + +class TestFaceDiv(OrderTest): + + name = "Face Divergence" + + def getError(self): + DIV = self.M.faceDiv + + #Test function + fun = lambda x: np.sin(x) + Fx = fun(self.M.gridFx[:, 0]) + Fy = fun(self.M.gridFy[:, 1]) + Fz = fun(self.M.gridFz[:, 2]) + + F = np.concatenate((Fx, Fy, Fz)) + divF = DIV*F + sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + divF_anal = sol(self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2]) + + err = np.linalg.norm((divF-divF_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestNodalGrad(OrderTest): + + name = "Nodal Gradient" + + def getError(self): + GRAD = self.M.nodalGrad + #Test function + fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) + + phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1], self.M.gridN[:, 2]) + gradE = GRAD*phi + + Ex = sol(self.M.gridEx[:, 0]) + Ey = sol(self.M.gridEy[:, 1]) + Ez = sol(self.M.gridEz[:, 2]) + + gradE_anal = np.concatenate((Ex, Ey, Ez)) + err = np.linalg.norm((gradE-gradE_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +if __name__ == '__main__': + unittest.main() From f86ebe12dbdcd83aa217cc2b1c61939c57363e7b Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 22 Jul 2013 13:40:18 -0700 Subject: [PATCH 10/11] Added Poisson Equation Tests. Forwards and Backwards. --- SimPEG/tests/OrderTest.py | 5 ++-- SimPEG/tests/test_operatorOrders.py | 37 ++++++++++++++++++++++++++--- 2 files changed, 37 insertions(+), 5 deletions(-) diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py index 68837f7a..652b224a 100644 --- a/SimPEG/tests/OrderTest.py +++ b/SimPEG/tests/OrderTest.py @@ -34,15 +34,16 @@ class OrderTest(unittest.TestCase): if ii == 0: print '' print 'Testing order of: ' + self.name + print '__________________________________________' print ' h | inf norm | ratio | order' - print '------------------------------------------' + print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~|~~~~~~~~~~' print '%4i | %8.2e |' % (nc, err) else: order.append(np.log(err/err_old)/np.log(float(nc_old)/float(nc))) print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1]) err_old = err nc_old = nc - + print '------------------------------------------' self.assertTrue(len(np.where(np.array(order) > 0.9*self.expectedOrder)[0]) > np.floor(0.75*len(order))) diff --git a/SimPEG/tests/test_operatorOrders.py b/SimPEG/tests/test_operatorOrders.py index 9bb9fa03..9fdebc10 100644 --- a/SimPEG/tests/test_operatorOrders.py +++ b/SimPEG/tests/test_operatorOrders.py @@ -1,10 +1,10 @@ import numpy as np from OrderTest import OrderTest import unittest +from scipy.sparse.linalg import dsolve class TestCurl(OrderTest): - name = "Curl" def getError(self): @@ -33,7 +33,6 @@ class TestCurl(OrderTest): class TestFaceDiv(OrderTest): - name = "Face Divergence" def getError(self): @@ -59,7 +58,6 @@ class TestFaceDiv(OrderTest): class TestNodalGrad(OrderTest): - name = "Nodal Gradient" def getError(self): @@ -84,5 +82,38 @@ class TestNodalGrad(OrderTest): self.orderTest() +class TestPoissonEqn(OrderTest): + name = "Poisson Equation" + meshSizes = [16, 20, 24] + + def getError(self): + # Create some functions to integrate + fun = lambda x: np.sin(2*np.pi*x[:, 0])*np.sin(2*np.pi*x[:, 1])*np.sin(2*np.pi*x[:, 2]) + sol = lambda x: -3.*((2*np.pi)**2)*fun(x) + + self.M.setCellGradBC('dirichlet') + + D = self.M.faceDiv + G = self.M.cellGrad + if self.forward: + sA = sol(self.M.gridCC) + sN = D*G*fun(self.M.gridCC) + err = np.linalg.norm((sA - sN), np.inf) + else: + fA = fun(self.M.gridCC) + fN = dsolve.spsolve(D*G, sol(self.M.gridCC)) + err = np.linalg.norm((fA - fN), np.inf) + return err + + def test_orderForward(self): + self.name = "Poisson Equation - Forward" + self.forward = True + self.orderTest() + + def test_orderBackward(self): + self.name = "Poisson Equation - Backward" + self.forward = False + self.orderTest() + if __name__ == '__main__': unittest.main() From 50f40c743dbedb140b65b558952087486b5984a5 Mon Sep 17 00:00:00 2001 From: Lars Ruthotto Date: Mon, 22 Jul 2013 15:12:34 -0700 Subject: [PATCH 11/11] added comments --- SimPEG/DiffOperators.py | 22 +-- SimPEG/getDiffOpps.py | 207 ---------------------------- SimPEG/sputils.py | 2 +- SimPEG/tests/OrderTest.py | 78 +++++++++-- SimPEG/tests/test_operatorOrders.py | 119 ---------------- SimPEG/tests/test_tensorMesh.py | 113 ++++++++++++++- 6 files changed, 196 insertions(+), 345 deletions(-) delete mode 100644 SimPEG/getDiffOpps.py delete mode 100644 SimPEG/tests/test_operatorOrders.py diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index bec2492e..c5daa96d 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -10,6 +10,7 @@ def ddx(n): def checkBC(bc): + """ Checks if boundary condition 'bc' is valid. """ if(type(bc) is str): bc = [bc, bc] assert type(bc) is list, 'bc must be a list' @@ -22,7 +23,7 @@ def checkBC(bc): def ddxCellGrad(n, bc): - """Define 1D derivatives, outer, this means we go from n to n+1""" + """Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1""" bc = checkBC(bc) D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr") @@ -40,7 +41,7 @@ def ddxCellGrad(n, bc): def av(n): - """Define 1D averaging operator""" + """Define 1D averaging operator from cell-centres to nodes.""" return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") @@ -49,10 +50,10 @@ class DiffOperators(object): Class creates the differential operators that you need! """ def __init__(self): - raise Exception('You should use a Mesh class.') + raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.') def faceDiv(): - doc = "Construct the 3D Divergence operator on Faces." + doc = "Construct divergence operator (face-stg to cell-centres)." def fget(self): if(self._faceDiv is None): @@ -81,7 +82,7 @@ class DiffOperators(object): faceDiv = property(**faceDiv()) def nodalGrad(): - doc = "Construct the 3D nodal gradient operator." + doc = "Construct gradient operator (nodes to edges)." def fget(self): if(self._nodalGrad is None): @@ -109,11 +110,14 @@ class DiffOperators(object): def setCellGradBC(self, BC): """ - e.g. + Function that sets the boundary conditions for cell-centred derivative operators. - BC = 'neumann' - BC = ['neumann', 'dirichlet', 'neumann'] - BC = [['neumann', 'dirichlet'], 'dirichlet', 'neumann'] + Examples: + + BC = 'neumann' # Neumann in all directions + BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else + BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain, + # Dirichlet else """ if(type(BC) is str): diff --git a/SimPEG/getDiffOpps.py b/SimPEG/getDiffOpps.py deleted file mode 100644 index bcfe2c8a..00000000 --- a/SimPEG/getDiffOpps.py +++ /dev/null @@ -1,207 +0,0 @@ -import numpy as np -from scipy import sparse -from utils import mkvc -from sputils import * -#from sputils import ddx, sdiag, speye, kron3, spzeros, appendBottom3, - -def getVol(h): - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # Compute cell volumes - V = mkvc(np.outer(mkvc(np.outer(h1,h2)),h3)) - - return V - -def getArea(h): - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - # Compute areas of cell faces - area1 = mkvc(np.outer(np.ones(n1+1),np.outer(h2,h3))) - area2 = mkvc(np.outer(h1,mkvc(np.outer(np.ones(n2+1),h3)))) - area3 = mkvc(np.outer(h1,mkvc(np.outer(h2,np.ones(n3+1))))) - area = np.hstack((np.hstack((area1, area2)), area3)) - - return area - -def getLength(h): - - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - # compute the length of each edge - Length1 = mkvc(np.outer(h1,mkvc(np.outer(np.ones(n2+1),np.ones(n3+1))))) - Length2 = mkvc(np.outer(np.ones(n1+1),mkvc(np.outer(h2,np.ones(n3+1))))) - Length3 = mkvc(np.outer(np.ones(n1+1),mkvc(np.outer(np.ones(n2+1),h3)))) - - Length = np.hstack((np.hstack((Length1, Length2)), Length3)) - - return Length - - -def getDivMatrix(h): - """Consturct the 3D divergence operator on Faces.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - area = getArea(h) - S = sdiag(area) - - # Compute cell volumes - V = getVol(h) - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - D1 = kron3(speye(n3), speye(n2), d1) - D2 = kron3(speye(n3), d2, speye(n1)) - D3 = kron3(d3, speye(n2), speye(n1)) - - D = sparse.hstack((sparse.hstack((D1, D2)), D3)) - return sdiag(1/V)*D*S - - -def getCurlMatrix(h): - """Edge CURL """ - - # Cell sizes in each direction - h1 = h[0]; h2 = h[1]; h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3) - - d1 = ddx(n1); d2 = ddx(n2); d3 = ddx(n3) - # derivatives on x-edge variables - D32 = kron3(d3, speye(n2), speye(n1+1)) - D23 = kron3(speye(n3), d2, speye(n1+1)) - D31 = kron3(d3, speye(n2+1), speye(n1)) - D13 = kron3(speye(n3), speye(n2+1), d1) - D21 = kron3(speye(n3+1), d2, speye(n1)) - D12 = kron3(speye(n3+1), speye(n2), d1) - - O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) - O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) - O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) - - CURL = appendBottom3( - appendRight3(O1, -D32, D23), - appendRight3(D31, O2, -D13), - appendRight3(-D21, D12, O3)) - - - area = getArea(h) - S = sdiag(1/area) - - # Compute edge length - lngth = getLength(h) - L = sdiag(lngth) - - return S*(CURL*L) - - -def getNodalGradient(h): - """Nodal Gradients""" - - # Cell sizes in each direction - h1 = h[0]; h2 = h[1]; h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3) - - - D1 = kron3(speye(n3+1), speye(n2+1), ddx(n1)) - D2 = kron3(speye(n3+1), ddx(n2), speye(n1+1)) - D3 = kron3(ddx(n3), speye(n2+1), speye(n1+1)) - - # topological gradient - GRAD = appendBottom3(D1, D2, D3) - - # scale for non-uniform mesh - # Compute edge length - lngth = getLength(h) - L = sdiag(1/lngth) - - return L*GRAD - - -def getEdgeToCellAverge(h): - - """Average from Edge to Cell center """ - - # Cell sizes in each direction - h1 = h[0]; h2 = h[1]; h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3) - - a1 = av(n1); a2 = av(n2); a3 = av(n3) - # derivatives on x-edge variables - A1 = kron3(a3, a2, speye(n1)) - A2 = kron3(a3, speye(n2), a1) - A3 = kron3(speye(n3), a2, a1) - - return appendRight3(A1, A2, A3) - - -def getFaceToCellAverge(h): - - """Average from Edge to Cell center """ - - # Cell sizes in each direction - h1 = h[0]; h2 = h[1]; h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1); n2 = np.size(h2); n3 = np.size(h3) - - a1 = av(n1); a2 = av(n2); a3 = av(n3) - # derivatives on x-edge variables - A1 = kron3(speye(n3), speye(n2), a1) - A2 = kron3(speye(n3), a2, speye(n1)) - A3 = kron3(a3, speye(n2), speye(n1)) - - return appendRight3(A1, A2, A3) - - -def getEdgeMassMatrix(h,sigma): - # mass matix for products of edge functions w'*M(sigma)*e - - Av = getEdgeToCellAverge(h) - v = getVol(h) - sigma = mkvc(sigma) - - return sdiag(Av.T*(v*sigma)) - -def getFaceMassMatrix(h,sigma): - # mass matix for products of edge functions w'*M(sigma)*e - - Av = getFaceToCellAverge(h) - v = getVol(h) - sigma = mkvc(sigma) - - return sdiag(Av.T*(v*sigma)) \ No newline at end of file diff --git a/SimPEG/sputils.py b/SimPEG/sputils.py index 79bbdf04..b078fff2 100644 --- a/SimPEG/sputils.py +++ b/SimPEG/sputils.py @@ -13,7 +13,7 @@ def speye(n): def kron3(A, B, C): - """Two kron prods""" + """Three kron prods""" return sp.kron(sp.kron(A, B), C, format="csr") diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py index 652b224a..7f9f97ca 100644 --- a/SimPEG/tests/OrderTest.py +++ b/SimPEG/tests/OrderTest.py @@ -6,14 +6,69 @@ import unittest class OrderTest(unittest.TestCase): - """Order test sets up the basics for testing order of decrease for a function on a mesh.""" + """ + + OrderTest is a base class for testing convergence orders with respect to mesh + sizes of integral/differential operators. + + Mathematical Problem: + + Given are an operator A and its discretization A[h]. For a given test function f + and h --> 0 we compare: + + error(h) = \| A[h](f) - A(f) \|_{\infty} + + Note that you can provide any norm. + + Test is passed when estimated rate order of convergence is at least 90% of the + estimated rate supplied by the user. + + Minimal example for a curl operator: + + class TestCURL(OrderTest): + name = "Curl" + + def getError(self): + # For given Mesh, generate A[h], f and A(f) and return norm of error. + + + fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) + sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) + + + Ex = fun(self.M.gridEx[:, 1]) + Ey = fun(self.M.gridEy[:, 2]) + Ez = fun(self.M.gridEz[:, 0]) + f = np.concatenate((Ex, Ey, Ez)) + + Fx = sol(self.M.gridFx[:, 2]) + Fy = sol(self.M.gridFy[:, 0]) + Fz = sol(self.M.gridFz[:, 1]) + Af = np.concatenate((Fx, Fy, Fz)) + + # Generate DIV matrix + Ah = self.M.edgeCurl + + curlE = Ah*E + err = np.linalg.norm((Ah*f -Af), np.inf) + return err + + def test_order(self): + # runs the test + self.orderTest() + + See also: test_operatorOrder.py + + """ name = "Order Test" expectedOrder = 2 meshSizes = [4, 8, 16, 32] def setupMesh(self, nc): - # Define the mesh + """ + For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc. + """ h1 = np.ones(nc)/nc h2 = np.ones(nc)/nc h3 = np.ones(nc)/nc @@ -21,10 +76,17 @@ class OrderTest(unittest.TestCase): self.M = TensorMesh(h) def getError(self): - """Overwrite this function with the guts of the test.""" + """For given h, generate A[h], f and A(f) and return norm of error.""" return 1. def orderTest(self): + """ + For number of cells specified in meshSizes setup mesh, call getError + and prints mesh size, error, ratio between current and previous error, + and estimated order of convergence. + + + """ order = [] err_old = 0. nc_old = 0. @@ -34,16 +96,16 @@ class OrderTest(unittest.TestCase): if ii == 0: print '' print 'Testing order of: ' + self.name - print '__________________________________________' - print ' h | inf norm | ratio | order' - print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~|~~~~~~~~~~' + print '_____________________________________________' + print ' h | error | e(i-1)/e(i) | order' + print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~' print '%4i | %8.2e |' % (nc, err) else: order.append(np.log(err/err_old)/np.log(float(nc_old)/float(nc))) - print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1]) + print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1]) err_old = err nc_old = nc - print '------------------------------------------' + print '---------------------------------------------' self.assertTrue(len(np.where(np.array(order) > 0.9*self.expectedOrder)[0]) > np.floor(0.75*len(order))) diff --git a/SimPEG/tests/test_operatorOrders.py b/SimPEG/tests/test_operatorOrders.py deleted file mode 100644 index 9fdebc10..00000000 --- a/SimPEG/tests/test_operatorOrders.py +++ /dev/null @@ -1,119 +0,0 @@ -import numpy as np -from OrderTest import OrderTest -import unittest -from scipy.sparse.linalg import dsolve - - -class TestCurl(OrderTest): - name = "Curl" - - def getError(self): - fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) - sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) - - Ex = fun(self.M.gridEx[:, 1]) - Ey = fun(self.M.gridEy[:, 2]) - Ez = fun(self.M.gridEz[:, 0]) - E = np.concatenate((Ex, Ey, Ez)) - - Fx = sol(self.M.gridFx[:, 2]) - Fy = sol(self.M.gridFy[:, 0]) - Fz = sol(self.M.gridFz[:, 1]) - curlE_anal = np.concatenate((Fx, Fy, Fz)) - - # Generate DIV matrix - CURL = self.M.edgeCurl - - curlE = CURL*E - err = np.linalg.norm((curlE-curlE_anal), np.inf) - return err - - def test_order(self): - self.orderTest() - - -class TestFaceDiv(OrderTest): - name = "Face Divergence" - - def getError(self): - DIV = self.M.faceDiv - - #Test function - fun = lambda x: np.sin(x) - Fx = fun(self.M.gridFx[:, 0]) - Fy = fun(self.M.gridFy[:, 1]) - Fz = fun(self.M.gridFz[:, 2]) - - F = np.concatenate((Fx, Fy, Fz)) - divF = DIV*F - sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) - divF_anal = sol(self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2]) - - err = np.linalg.norm((divF-divF_anal), np.inf) - - return err - - def test_order(self): - self.orderTest() - - -class TestNodalGrad(OrderTest): - name = "Nodal Gradient" - - def getError(self): - GRAD = self.M.nodalGrad - #Test function - fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) - sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) - - phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1], self.M.gridN[:, 2]) - gradE = GRAD*phi - - Ex = sol(self.M.gridEx[:, 0]) - Ey = sol(self.M.gridEy[:, 1]) - Ez = sol(self.M.gridEz[:, 2]) - - gradE_anal = np.concatenate((Ex, Ey, Ez)) - err = np.linalg.norm((gradE-gradE_anal), np.inf) - - return err - - def test_order(self): - self.orderTest() - - -class TestPoissonEqn(OrderTest): - name = "Poisson Equation" - meshSizes = [16, 20, 24] - - def getError(self): - # Create some functions to integrate - fun = lambda x: np.sin(2*np.pi*x[:, 0])*np.sin(2*np.pi*x[:, 1])*np.sin(2*np.pi*x[:, 2]) - sol = lambda x: -3.*((2*np.pi)**2)*fun(x) - - self.M.setCellGradBC('dirichlet') - - D = self.M.faceDiv - G = self.M.cellGrad - if self.forward: - sA = sol(self.M.gridCC) - sN = D*G*fun(self.M.gridCC) - err = np.linalg.norm((sA - sN), np.inf) - else: - fA = fun(self.M.gridCC) - fN = dsolve.spsolve(D*G, sol(self.M.gridCC)) - err = np.linalg.norm((fA - fN), np.inf) - return err - - def test_orderForward(self): - self.name = "Poisson Equation - Forward" - self.forward = True - self.orderTest() - - def test_orderBackward(self): - self.name = "Poisson Equation - Backward" - self.forward = False - self.orderTest() - -if __name__ == '__main__': - unittest.main() diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/tests/test_tensorMesh.py index 2a6cb0a1..c0b194db 100644 --- a/SimPEG/tests/test_tensorMesh.py +++ b/SimPEG/tests/test_tensorMesh.py @@ -3,9 +3,11 @@ import unittest import sys sys.path.append('../') from TensorMesh import TensorMesh +from OrderTest import OrderTest +from scipy.sparse.linalg import dsolve -class TestSequenceFunctions(unittest.TestCase): +class BasicTensorMeshTests(unittest.TestCase): def setUp(self): a = np.array([1, 1, 1]) @@ -55,6 +57,115 @@ class TestSequenceFunctions(unittest.TestCase): t1 = np.all(self.mesh2.edge == test_edge) self.assertTrue(t1) +class TestCurl(OrderTest): + name = "Curl" + def getError(self): + fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) + sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) + + Ex = fun(self.M.gridEx[:, 1]) + Ey = fun(self.M.gridEy[:, 2]) + Ez = fun(self.M.gridEz[:, 0]) + E = np.concatenate((Ex, Ey, Ez)) + + Fx = sol(self.M.gridFx[:, 2]) + Fy = sol(self.M.gridFy[:, 0]) + Fz = sol(self.M.gridFz[:, 1]) + curlE_anal = np.concatenate((Fx, Fy, Fz)) + + # Generate DIV matrix + CURL = self.M.edgeCurl + + curlE = CURL*E + err = np.linalg.norm((curlE-curlE_anal), np.inf) + return err + + def test_order(self): + self.orderTest() + + +class TestFaceDiv(OrderTest): + name = "Face Divergence" + + def getError(self): + DIV = self.M.faceDiv + + #Test function + fun = lambda x: np.sin(x) + Fx = fun(self.M.gridFx[:, 0]) + Fy = fun(self.M.gridFy[:, 1]) + Fz = fun(self.M.gridFz[:, 2]) + + F = np.concatenate((Fx, Fy, Fz)) + divF = DIV*F + sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + divF_anal = sol(self.M.gridCC[:, 0], self.M.gridCC[:, 1], self.M.gridCC[:, 2]) + + err = np.linalg.norm((divF-divF_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestNodalGrad(OrderTest): + name = "Nodal Gradient" + + def getError(self): + GRAD = self.M.nodalGrad + #Test function + fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) + + phi = fun(self.M.gridN[:, 0], self.M.gridN[:, 1], self.M.gridN[:, 2]) + gradE = GRAD*phi + + Ex = sol(self.M.gridEx[:, 0]) + Ey = sol(self.M.gridEy[:, 1]) + Ez = sol(self.M.gridEz[:, 2]) + + gradE_anal = np.concatenate((Ex, Ey, Ez)) + err = np.linalg.norm((gradE-gradE_anal), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestPoissonEqn(OrderTest): + name = "Poisson Equation" + meshSizes = [16, 20, 24] + + def getError(self): + # Create some functions to integrate + fun = lambda x: np.sin(2*np.pi*x[:, 0])*np.sin(2*np.pi*x[:, 1])*np.sin(2*np.pi*x[:, 2]) + sol = lambda x: -3.*((2*np.pi)**2)*fun(x) + + self.M.setCellGradBC('dirichlet') + + D = self.M.faceDiv + G = self.M.cellGrad + if self.forward: + sA = sol(self.M.gridCC) + sN = D*G*fun(self.M.gridCC) + err = np.linalg.norm((sA - sN), np.inf) + else: + fA = fun(self.M.gridCC) + fN = dsolve.spsolve(D*G, sol(self.M.gridCC)) + err = np.linalg.norm((fA - fN), np.inf) + return err + + def test_orderForward(self): + self.name = "Poisson Equation - Forward" + self.forward = True + self.orderTest() + + def test_orderBackward(self): + self.name = "Poisson Equation - Backward" + self.forward = False + self.orderTest() if __name__ == '__main__': unittest.main()