From 50f40c743dbedb140b65b558952087486b5984a5 Mon Sep 17 00:00:00 2001 From: Lars Ruthotto Date: Mon, 22 Jul 2013 15:12:34 -0700 Subject: [PATCH] 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()