From fbdde2ac245a5b002c67064065472f3d28bb557e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 17 Feb 2014 12:00:50 -0800 Subject: [PATCH 01/15] boundaryCondition initial work. --- SimPEG/Mesh/DiffOperators.py | 98 ++++++++++++++++ SimPEG/Mesh/TensorMesh.py | 50 ++++++++ SimPEG/Mesh/TensorView.py | 6 +- SimPEG/Tests/TestUtils.py | 12 +- SimPEG/Tests/test_boundaryPoisson.py | 166 +++++++++++++++++++++++++++ 5 files changed, 325 insertions(+), 7 deletions(-) create mode 100644 SimPEG/Tests/test_boundaryPoisson.py diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 1712964e..8ca498d7 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -460,6 +460,104 @@ class DiffOperators(object): _edgeCurl = None edgeCurl = property(**edgeCurl()) + def getBCProjWF(self, BC, discretization='CC'): + """ + + The weak form boundary condition projection matrices. + + 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): + BC = [BC for _ in self.vnC] # 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) + + + def projDirichlet(n, bc): + bc = checkBC(bc) + ij = ([0,n], [0,1]) + vals = [0,0] + if(bc[0] == 'dirichlet'): + vals[0] = -1 + if(bc[1] == 'dirichlet'): + vals[1] = 1 + return sp.csr_matrix((vals, ij), shape=(n+1,2)) + + def projNeumannIn(n, bc): + bc = checkBC(bc) + P = sp.identity(n+1).tocsr() + if(bc[0] == 'neumann'): + P = P[1:,:] + if(bc[0] == 'neumann'): + P = P[:-1,:] + return P + + def projNeumannOut(n, bc): + bc = checkBC(bc) + ij = ([0, 1],[0, n]) + vals = [0,0] + if(bc[0] == 'neumann'): + vals[0] = 1 + if(bc[1] == 'neumann'): + vals[1] = 1 + return sp.csr_matrix((vals, ij), shape=(2,n+1)) + + n = self.vnC + indF = self.faceBoundaryInd + if(self.dim == 1): + Pbc = projDirichlet(n[0], BC[0]) + indF = indF[0] | indF[1] + Pbc = Pbc*sdiag(self.area[indF]) + + Pin = projNeumannIn(n[0], BC[0]) + + Pout = projNeumannOut(n[0], BC[0]) + elif(self.dim == 2): + Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0])) + Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0])) + Pbc = sp.block_diag((Pbc1, Pbc2), format="csr") + indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])] + Pbc = Pbc*sdiag(self.area[indF]) + + P1 = sp.kron(speye(n[1]), projNeumannIn(n[0], BC[0])) + P2 = sp.kron(projNeumannIn(n[1], BC[1]), speye(n[0])) + Pin = sp.block_diag((P1, P2), format="csr") + + P1 = sp.kron(speye(n[1]), projNeumannOut(n[0], BC[0])) + P2 = sp.kron(projNeumannOut(n[1], BC[1]), speye(n[0])) + Pout = sp.block_diag((P1, P2), format="csr") + elif(self.dim == 3): + Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0])) + Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0])) + Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0])) + Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr") + indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])] + Pbc = Pbc*sdiag(self.area[indF]) + + P1 = kron3(speye(n[2]), speye(n[1]), projNeumannIn(n[0], BC[0])) + P2 = kron3(speye(n[2]), projNeumannIn(n[1], BC[1]), speye(n[0])) + P3 = kron3(projNeumannIn(n[2], BC[2]), speye(n[1]), speye(n[0])) + Pin = sp.block_diag((P1, P2, P3), format="csr") + + P1 = kron3(speye(n[2]), speye(n[1]), projNeumannOut(n[0], BC[0])) + P2 = kron3(speye(n[2]), projNeumannOut(n[1], BC[1]), speye(n[0])) + P3 = kron3(projNeumannOut(n[2], BC[2]), speye(n[1]), speye(n[0])) + Pout = sp.block_diag((P1, P2, P3), format="csr") + + return Pbc, Pin, Pout + + # --------------- Averaging --------------------- @property diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 12253bb5..8f83e793 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -419,6 +419,56 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) return Q + + @property + def faceBoundaryInd(self): + """ + Find indices of boundary faces in each direction + """ + if self.dim==1: + indxd = (self.gridFx==min(self.gridFx)) + indxu = (self.gridFx==max(self.gridFx)) + return indxd, indxu + elif self.dim==2: + indxd = (self.gridFx[:,0]==min(self.gridFx[:,0])) + indxu = (self.gridFx[:,0]==max(self.gridFx[:,0])) + indyd = (self.gridFy[:,1]==min(self.gridFy[:,1])) + indyu = (self.gridFy[:,1]==max(self.gridFy[:,1])) + return indxd, indxu, indyd, indyu + elif self.dim==3: + indxd = (self.gridFx[:,0]==min(self.gridFx[:,0])) + indxu = (self.gridFx[:,0]==max(self.gridFx[:,0])) + indyd = (self.gridFy[:,1]==min(self.gridFy[:,1])) + indyu = (self.gridFy[:,1]==max(self.gridFy[:,1])) + indzd = (self.gridFz[:,2]==min(self.gridFz[:,2])) + indzu = (self.gridFz[:,2]==max(self.gridFz[:,2])) + return indxd, indxu, indyd, indyu, indzd, indzu + + @property + def cellBoundaryInd(self): + """ + Find indices of boundary faces in each direction + """ + if self.dim==1: + indxd = (self.gridCC==min(self.gridCC)) + indxu = (self.gridCC==max(self.gridCC)) + return indxd, indxu + elif self.dim==2: + indxd = (self.gridCC[:,0]==min(self.gridCC[:,0])) + indxu = (self.gridCC[:,0]==max(self.gridCC[:,0])) + indyd = (self.gridCC[:,1]==min(self.gridCC[:,1])) + indyu = (self.gridCC[:,1]==max(self.gridCC[:,1])) + return indxd, indxu, indyd, indyu + elif self.dim==3: + indxd = (self.gridCC[:,0]==min(self.gridCC[:,0])) + indxu = (self.gridCC[:,0]==max(self.gridCC[:,0])) + indyd = (self.gridCC[:,1]==min(self.gridCC[:,1])) + indyu = (self.gridCC[:,1]==max(self.gridCC[:,1])) + indzd = (self.gridCC[:,2]==min(self.gridCC[:,2])) + indzu = (self.gridCC[:,2]==max(self.gridCC[:,2])) + return indxd, indxu, indyd, indyu, indzd, indzu + + if __name__ == '__main__': print('Welcome to tensor mesh!') diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index e4dd9243..597bf989 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -74,7 +74,7 @@ class TensorView(object): if I.size != np.prod(self.vnEz): ex, ey, I = self.r(I,'E','E','M') elif imageType[0] == 'E': plotAll = len(imageType) == 1 - options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} + options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":False} fig = plt.figure(figNum) # Determine the subplot number: 131, 121 numPlots = 130 if plotAll else len(imageType)/2*10+100 @@ -92,10 +92,11 @@ class TensorView(object): ax_z = plt.subplot(numPlots+pltNum) self.plotImage(ez, imageType='Ez', ax=ax_z, **options) pltNum +=1 + if showIt: plt.show() return elif imageType[0] == 'F': plotAll = len(imageType) == 1 - options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt} + options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":False} fig = plt.figure(figNum) # Determine the subplot number: 131, 121 numPlots = 130 if plotAll else len(imageType)/2*10+100 @@ -113,6 +114,7 @@ class TensorView(object): ax_z = plt.subplot(numPlots+pltNum) self.plotImage(fxyz[2], imageType='Fz', ax=ax_z, **options) pltNum +=1 + if showIt: plt.show() return else: raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'") diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index d84408d4..b46fbf10 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -93,9 +93,9 @@ class OrderTest(unittest.TestCase): h3 = np.ones(nc)/nc h = [h1, h2, h3] elif 'random' in self._meshType: - h1 = np.random.rand(nc) - h2 = np.random.rand(nc) - h3 = np.random.rand(nc) + h1 = np.random.rand(nc)*nc*0.5 + nc*0.5 + h2 = np.random.rand(nc)*nc*0.5 + nc*0.5 + h3 = np.random.rand(nc)*nc*0.5 + nc*0.5 h = [hi/np.sum(hi) for hi in [h1, h2, h3]] # normalize else: raise Exception('Unexpected meshType') @@ -111,10 +111,12 @@ class OrderTest(unittest.TestCase): kwrd = 'rotate' else: raise Exception('Unexpected meshType') - if self.meshDimension == 2: + if self.meshDimension == 1: + raise Exception('Lom not supported for 1D') + elif self.meshDimension == 2: X, Y = Utils.exampleLomGird([nc, nc], kwrd) self.M = LogicallyOrthogonalMesh([X, Y]) - if self.meshDimension == 3: + elif self.meshDimension == 3: X, Y, Z = Utils.exampleLomGird([nc, nc, nc], kwrd) self.M = LogicallyOrthogonalMesh([X, Y, Z]) return 1./nc diff --git a/SimPEG/Tests/test_boundaryPoisson.py b/SimPEG/Tests/test_boundaryPoisson.py new file mode 100644 index 00000000..fd6ede5d --- /dev/null +++ b/SimPEG/Tests/test_boundaryPoisson.py @@ -0,0 +1,166 @@ +import numpy as np +import scipy.sparse as sp +import unittest +from TestUtils import OrderTest +import matplotlib.pyplot as plt +from SimPEG import Utils, Solver + +MESHTYPES = ['uniformTensorMesh'] + +class Test1D_InhomogeneousDirichlet(OrderTest): + name = "1D - Dirichlet" + meshTypes = MESHTYPES + meshDimension = 1 + expectedOrders = 2 + meshSizes = [4, 8, 16, 32, 64, 128] + + def getError(self): + #Test function + phi = lambda x: np.cos(np.pi*x) + j_fun = lambda x: -np.pi*np.sin(np.pi*x) + q_fun = lambda x: -(np.pi**2)*np.cos(np.pi*x) + + xc_anal = phi(self.M.gridCC) + q_anal = q_fun(self.M.gridCC) + j_anal = j_fun(self.M.gridFx) + + #TODO: Check where our boundary conditions are CCx or Nx + # vec = self.M.vectorNx + vec = self.M.vectorCCx + bc = phi(vec[[0,-1]]) + + P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'dirichlet']]) + + Mc = self.M.getFaceInnerProduct() + McI = Utils.sdInv(self.M.getFaceInnerProduct()) + G = -self.M.faceDiv.T * Utils.sdiag(self.M.vol) + D = self.M.faceDiv + j = McI*(G*xc_anal + P*bc) + q = D*j + + # Rearrange if we know q to solve for x + A = D*McI*G + rhs = q_anal - D*McI*P*bc + + + if self.myTest == 'j': + err = np.linalg.norm((j-j_anal), np.inf) + elif self.myTest == 'q': + err = np.linalg.norm((q-q_anal), np.inf) + elif self.myTest == 'xc': + xc = Solver(A).solve(rhs) + err = np.linalg.norm((xc-xc_anal), np.inf) + elif self.myTest == 'xcJ': + xc = Solver(A).solve(rhs) + j = McI*(G*xc + P*bc) + err = np.linalg.norm((j-j_anal), np.inf) + + return err + + def test_orderJ(self): + self.name = "1D - InhomogeneousDirichlet_Forward j" + self.myTest = 'j' + self.orderTest() + + def test_orderQ(self): + self.name = "1D - InhomogeneousDirichlet_Forward q" + self.myTest = 'q' + self.orderTest() + + def test_orderX(self): + self.name = "1D - InhomogeneousDirichlet_Inverse" + self.myTest = 'xc' + self.orderTest() + + def test_orderXJ(self): + self.name = "1D - InhomogeneousDirichlet_Inverse J" + self.myTest = 'xcJ' + self.orderTest() + + +class Test2D_InhomogeneousDirichlet(OrderTest): + name = "2D - Dirichlet" + meshTypes = MESHTYPES + meshDimension = 2 + expectedOrders = 2 + meshSizes = [4, 8, 16, 32] + + def getError(self): + #Test function + phi = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1]) + j_funX = lambda x: -np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1]) + j_funY = lambda x: -np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1]) + q_fun = lambda x: -2*(np.pi**2)*phi(x) + + xc_anal = phi(self.M.gridCC) + q_anal = q_fun(self.M.gridCC) + jX_anal = j_funX(self.M.gridFx) + jY_anal = j_funY(self.M.gridFy) + j_anal = np.r_[jX_anal,jY_anal] + + #TODO: Check where our boundary conditions are CCx or Nx + # fxm,fxp,fym,fyp = self.M.faceBoundaryInd + # gBFx = self.M.gridFx[(fxm|fxp),:] + # gBFy = self.M.gridFy[(fym|fyp),:] + fxm,fxp,fym,fyp = self.M.cellBoundaryInd + gBFx = self.M.gridCC[(fxm|fxp),:] + gBFy = self.M.gridCC[(fym|fyp),:] + + bc = phi(np.r_[gBFx,gBFy]) + + # P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2)) + + P, Pin, Pout = self.M.getBCProjWF('dirichlet') + + Mc = self.M.getFaceInnerProduct() + McI = Utils.sdInv(self.M.getFaceInnerProduct()) + G = -self.M.faceDiv.T * Utils.sdiag(self.M.vol) + D = self.M.faceDiv + j = McI*(G*xc_anal + P*bc) + q = D*j + + # self.M.plotImage(j, 'FxFy', showIt=True) + + # Rearrange if we know q to solve for x + A = D*McI*G + rhs = q_anal - D*McI*P*bc + + if self.myTest == 'j': + err = np.linalg.norm((j-j_anal), np.inf) + elif self.myTest == 'q': + err = np.linalg.norm((q-q_anal), np.inf) + elif self.myTest == 'xc': + xc = Solver(A).solve(rhs) + err = np.linalg.norm((xc-xc_anal), np.inf) + elif self.myTest == 'xcJ': + xc = Solver(A).solve(rhs) + j = McI*(G*xc + P*bc) + err = np.linalg.norm((j-j_anal), np.inf) + + return err + + def test_orderJ(self): + self.name = "2D - InhomogeneousDirichlet_Forward j" + self.myTest = 'j' + self.orderTest() + + def test_orderQ(self): + self.name = "2D - InhomogeneousDirichlet_Forward q" + self.myTest = 'q' + self.orderTest() + + def test_orderX(self): + self.name = "2D - InhomogeneousDirichlet_Inverse" + self.myTest = 'xc' + self.orderTest() + + def test_orderXJ(self): + self.name = "2D - InhomogeneousDirichlet_Inverse J" + self.myTest = 'xcJ' + self.orderTest() + + + + +if __name__ == '__main__': + unittest.main() From d8c41e664ba6308082c741cd2c68ed274db6a891 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 17 Feb 2014 12:26:53 -0800 Subject: [PATCH 02/15] make general. --- SimPEG/Tests/test_boundaryPoisson.py | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/SimPEG/Tests/test_boundaryPoisson.py b/SimPEG/Tests/test_boundaryPoisson.py index fd6ede5d..fd8984a1 100644 --- a/SimPEG/Tests/test_boundaryPoisson.py +++ b/SimPEG/Tests/test_boundaryPoisson.py @@ -27,20 +27,24 @@ class Test1D_InhomogeneousDirichlet(OrderTest): #TODO: Check where our boundary conditions are CCx or Nx # vec = self.M.vectorNx vec = self.M.vectorCCx - bc = phi(vec[[0,-1]]) + + phi_bc = phi(vec[[0,-1]]) + j_bc = j_fun(vec[[0,-1]]) P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'dirichlet']]) Mc = self.M.getFaceInnerProduct() McI = Utils.sdInv(self.M.getFaceInnerProduct()) - G = -self.M.faceDiv.T * Utils.sdiag(self.M.vol) + G = -Pin.T*Pin*self.M.faceDiv.T * Utils.sdiag(self.M.vol) D = self.M.faceDiv - j = McI*(G*xc_anal + P*bc) - q = D*j + j = McI*(G*xc_anal + P*phi_bc) + q = D*Pin.T*Pin*j + D*Pout.T*j_bc # Rearrange if we know q to solve for x - A = D*McI*G - rhs = q_anal - D*McI*P*bc + A = D*Pin.T*Pin*McI*G + rhs = q_anal - D*Pin.T*Pin*McI*P*phi_bc - D*Pout.T*j_bc + # A = D*McI*G + # rhs = q_anal - D*McI*P*phi_bc if self.myTest == 'j': @@ -48,11 +52,15 @@ class Test1D_InhomogeneousDirichlet(OrderTest): elif self.myTest == 'q': err = np.linalg.norm((q-q_anal), np.inf) elif self.myTest == 'xc': + #TODO: fix the null space xc = Solver(A).solve(rhs) + print np.linalg.norm(Utils.mkvc(A*xc) - rhs) err = np.linalg.norm((xc-xc_anal), np.inf) elif self.myTest == 'xcJ': + #TODO: fix the null space xc = Solver(A).solve(rhs) - j = McI*(G*xc + P*bc) + print np.linalg.norm(Utils.mkvc(A*xc) - rhs) + j = McI*(G*xc + P*phi_bc) err = np.linalg.norm((j-j_anal), np.inf) return err From 6f52e9b50d59bb3c6946c0ca159f8ceca1e09f4c Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 17 Feb 2014 12:34:15 -0800 Subject: [PATCH 03/15] make symmetric and solve with CG --- SimPEG/Tests/test_boundaryPoisson.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/SimPEG/Tests/test_boundaryPoisson.py b/SimPEG/Tests/test_boundaryPoisson.py index fd8984a1..5f3fc9ac 100644 --- a/SimPEG/Tests/test_boundaryPoisson.py +++ b/SimPEG/Tests/test_boundaryPoisson.py @@ -35,14 +35,15 @@ class Test1D_InhomogeneousDirichlet(OrderTest): Mc = self.M.getFaceInnerProduct() McI = Utils.sdInv(self.M.getFaceInnerProduct()) - G = -Pin.T*Pin*self.M.faceDiv.T * Utils.sdiag(self.M.vol) + V = Utils.sdiag(self.M.vol) + G = -Pin.T*Pin*self.M.faceDiv.T * V D = self.M.faceDiv j = McI*(G*xc_anal + P*phi_bc) - q = D*Pin.T*Pin*j + D*Pout.T*j_bc + q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc # Rearrange if we know q to solve for x - A = D*Pin.T*Pin*McI*G - rhs = q_anal - D*Pin.T*Pin*McI*P*phi_bc - D*Pout.T*j_bc + A = V*D*Pin.T*Pin*McI*G + rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc # A = D*McI*G # rhs = q_anal - D*McI*P*phi_bc @@ -53,8 +54,9 @@ class Test1D_InhomogeneousDirichlet(OrderTest): err = np.linalg.norm((q-q_anal), np.inf) elif self.myTest == 'xc': #TODO: fix the null space - xc = Solver(A).solve(rhs) - print np.linalg.norm(Utils.mkvc(A*xc) - rhs) + solver = Solver(A, doDirect=False, options={'M':'J','iterSolver':'CG','backend':'scipy','maxIter':1000}) + xc = solver.solve(rhs) + print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs) err = np.linalg.norm((xc-xc_anal), np.inf) elif self.myTest == 'xcJ': #TODO: fix the null space From 5a971baaf427053e4ecca9b953472d3108008fd7 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 17 Feb 2014 12:58:02 -0800 Subject: [PATCH 04/15] bug fix for comparison with q --- SimPEG/Tests/test_boundaryPoisson.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Tests/test_boundaryPoisson.py b/SimPEG/Tests/test_boundaryPoisson.py index 5f3fc9ac..42a7bc59 100644 --- a/SimPEG/Tests/test_boundaryPoisson.py +++ b/SimPEG/Tests/test_boundaryPoisson.py @@ -51,7 +51,7 @@ class Test1D_InhomogeneousDirichlet(OrderTest): if self.myTest == 'j': err = np.linalg.norm((j-j_anal), np.inf) elif self.myTest == 'q': - err = np.linalg.norm((q-q_anal), np.inf) + err = np.linalg.norm((q-V*q_anal), np.inf) elif self.myTest == 'xc': #TODO: fix the null space solver = Solver(A, doDirect=False, options={'M':'J','iterSolver':'CG','backend':'scipy','maxIter':1000}) From 9c8b64b9f1c2fd3b5b02a0a44770604d16ae5f36 Mon Sep 17 00:00:00 2001 From: seogi Date: Wed, 19 Feb 2014 14:25:55 -0800 Subject: [PATCH 05/15] Testing BC for cell-centered system - inhomogenous neumann BC - Mixed BC were tested Fixing bug in the function "projNeumannIn" --- SimPEG/Mesh/DiffOperators.py | 2 +- SimPEG/Tests/test_boundaryPoisson.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 8ca498d7..70a7e657 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -499,7 +499,7 @@ class DiffOperators(object): P = sp.identity(n+1).tocsr() if(bc[0] == 'neumann'): P = P[1:,:] - if(bc[0] == 'neumann'): + if(bc[1] == 'neumann'): P = P[:-1,:] return P diff --git a/SimPEG/Tests/test_boundaryPoisson.py b/SimPEG/Tests/test_boundaryPoisson.py index 42a7bc59..a14ab2a9 100644 --- a/SimPEG/Tests/test_boundaryPoisson.py +++ b/SimPEG/Tests/test_boundaryPoisson.py @@ -171,6 +171,5 @@ class Test2D_InhomogeneousDirichlet(OrderTest): - if __name__ == '__main__': unittest.main() From 72fe771071e088b3e9fe8d4b42009351ea604f2c Mon Sep 17 00:00:00 2001 From: seogi Date: Wed, 19 Feb 2014 15:24:37 -0800 Subject: [PATCH 06/15] Testing boundary conditoins - I am not sure what happend... but this is what I did for last push --- SimPEG/Tests/test_boundaryPoisson.py | 326 +++++++++++++++++++++++++++ 1 file changed, 326 insertions(+) diff --git a/SimPEG/Tests/test_boundaryPoisson.py b/SimPEG/Tests/test_boundaryPoisson.py index a14ab2a9..5e028dee 100644 --- a/SimPEG/Tests/test_boundaryPoisson.py +++ b/SimPEG/Tests/test_boundaryPoisson.py @@ -169,7 +169,333 @@ class Test2D_InhomogeneousDirichlet(OrderTest): self.myTest = 'xcJ' self.orderTest() +class Test1D_InhomogeneousNeumann(OrderTest): + name = "1D - Neumann" + meshTypes = MESHTYPES + meshDimension = 1 + expectedOrders = 2 + meshSizes = [4, 8, 16, 32, 64, 128] + def getError(self): + #Test function + phi = lambda x: np.sin(np.pi*x) + j_fun = lambda x: np.pi*np.cos(np.pi*x) + q_fun = lambda x: -(np.pi**2)*np.sin(np.pi*x) + + xc_anal = phi(self.M.gridCC) + q_anal = q_fun(self.M.gridCC) + j_anal = j_fun(self.M.gridFx) + + #TODO: Check where our boundary conditions are CCx or Nx + vecN = self.M.vectorNx + vecC = self.M.vectorCCx + + phi_bc = phi(vecC[[0,-1]]) + j_bc = j_fun(vecN[[0,-1]]) + + P, Pin, Pout = self.M.getBCProjWF([['neumann', 'neumann']]) + + Mc = self.M.getFaceInnerProduct() + McI = Utils.sdInv(self.M.getFaceInnerProduct()) + V = Utils.sdiag(self.M.vol) + G = -Pin.T*Pin*self.M.faceDiv.T * V + D = self.M.faceDiv + j = McI*(G*xc_anal + P*phi_bc) + q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc + + # Rearrange if we know q to solve for x + A = V*D*Pin.T*Pin*McI*G + rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc + # A = D*McI*G + # rhs = q_anal - D*McI*P*phi_bc + + + if self.myTest == 'j': + err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf) + elif self.myTest == 'q': + err = np.linalg.norm((q-V*q_anal), np.inf) + elif self.myTest == 'xc': + #TODO: fix the null space + xc, info = sp.linalg.minres(A, rhs, tol = 1e-6) + err = np.linalg.norm((xc-xc_anal), np.inf) + if info > 0: + print 'Solve does not work well' + print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs) + elif self.myTest == 'xcJ': + #TODO: fix the null space + xc, info = sp.linalg.minres(A, rhs, tol = 1e-6) + j = McI*(G*xc + P*phi_bc) + err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf) + if info > 0: + print 'Solve does not work well' + print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs) + return err + + def test_orderJ(self): + self.name = "1D - InhomogeneousNeumann_Forward j" + self.myTest = 'j' + self.orderTest() + + def test_orderQ(self): + self.name = "1D - InhomogeneousNeumann_Forward q" + self.myTest = 'q' + self.orderTest() + + def test_orderXJ(self): + self.name = "1D - InhomogeneousNeumann_Inverse J" + self.myTest = 'xcJ' + self.orderTest() + +class Test2D_InhomogeneousNeumann(OrderTest): + name = "2D - Neumann" + meshTypes = MESHTYPES + meshDimension = 2 + expectedOrders = 2 + meshSizes = [4, 8, 16, 32] + # meshSizes = [4] + + def getError(self): + #Test function + phi = lambda x: np.sin(np.pi*x[:,0])*np.sin(np.pi*x[:,1]) + j_funX = lambda x: np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1]) + j_funY = lambda x: np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1]) + q_fun = lambda x: -2*(np.pi**2)*phi(x) + + xc_anal = phi(self.M.gridCC) + q_anal = q_fun(self.M.gridCC) + jX_anal = j_funX(self.M.gridFx) + jY_anal = j_funY(self.M.gridFy) + j_anal = np.r_[jX_anal,jY_anal] + + #TODO: Check where our boundary conditions are CCx or Nx + + cxm,cxp,cym,cyp = self.M.cellBoundaryInd + fxm,fxp,fym,fyp = self.M.faceBoundaryInd + + gBFx = self.M.gridFx[(fxm|fxp),:] + gBFy = self.M.gridFy[(fym|fyp),:] + + gBCx = self.M.gridCC[(cxm|cxp),:] + gBCy = self.M.gridCC[(cym|cyp),:] + + phi_bc = phi(np.r_[gBFx,gBFy]) + j_bc = np.r_[j_funX(gBFx), j_funY(gBFy)] + + # P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2)) + + P, Pin, Pout = self.M.getBCProjWF('neumann') + + Mc = self.M.getFaceInnerProduct() + McI = Utils.sdInv(self.M.getFaceInnerProduct()) + V = Utils.sdiag(self.M.vol) + G = -Pin.T*Pin*self.M.faceDiv.T * V + D = self.M.faceDiv + j = McI*(G*xc_anal + P*phi_bc) + q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc + + # Rearrange if we know q to solve for x + A = V*D*Pin.T*Pin*McI*G + rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc + + if self.myTest == 'j': + err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf) + elif self.myTest == 'q': + err = np.linalg.norm((q-V*q_anal), np.inf) + elif self.myTest == 'xc': + #TODO: fix the null space + xc, info = sp.linalg.minres(A, rhs, tol = 1e-6) + err = np.linalg.norm((xc-xc_anal), np.inf) + if info > 0: + print 'Solve does not work well' + print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs) + elif self.myTest == 'xcJ': + #TODO: fix the null space + xc, info = sp.linalg.minres(A, rhs, tol = 1e-6) + j = McI*(G*xc + P*phi_bc) + err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf) + if info > 0: + print 'Solve does not work well' + print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs) + return err + + def test_orderJ(self): + self.name = "2D - InhomogeneousNeumann_Forward j" + self.myTest = 'j' + self.orderTest() + + def test_orderQ(self): + self.name = "2D - InhomogeneousNeumann_Forward q" + self.myTest = 'q' + self.orderTest() + + def test_orderXJ(self): + self.name = "2D - InhomogeneousNeumann_Inverse J" + self.myTest = 'xcJ' + self.orderTest() + +class Test1D_InhomogeneousMixed(OrderTest): + name = "1D - Mixed" + meshTypes = MESHTYPES + meshDimension = 1 + expectedOrders = 2 + meshSizes = [4, 8, 16, 32, 64, 128] + + def getError(self): + #Test function + phi = lambda x: np.cos(0.5*np.pi*x) + j_fun = lambda x: -0.5*np.pi*np.sin(0.5*np.pi*x) + q_fun = lambda x: -0.25*(np.pi**2)*np.cos(0.5*np.pi*x) + + xc_anal = phi(self.M.gridCC) + q_anal = q_fun(self.M.gridCC) + j_anal = j_fun(self.M.gridFx) + + #TODO: Check where our boundary conditions are CCx or Nx + vecN = self.M.vectorNx + vecC = self.M.vectorCCx + + phi_bc = phi(vecC[[0,-1]]) + j_bc = j_fun(vecN[[0,-1]]) + + P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'neumann']]) + + Mc = self.M.getFaceInnerProduct() + McI = Utils.sdInv(self.M.getFaceInnerProduct()) + V = Utils.sdiag(self.M.vol) + G = -Pin.T*Pin*self.M.faceDiv.T * V + D = self.M.faceDiv + j = McI*(G*xc_anal + P*phi_bc) + q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc + + # Rearrange if we know q to solve for x + A = V*D*Pin.T*Pin*McI*G + rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc + # A = D*McI*G + # rhs = q_anal - D*McI*P*phi_bc + + + if self.myTest == 'j': + err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf) + elif self.myTest == 'q': + err = np.linalg.norm((q-V*q_anal), np.inf) + elif self.myTest == 'xc': + #TODO: fix the null space + xc, info = sp.linalg.minres(A, rhs, tol = 1e-6) + err = np.linalg.norm((xc-xc_anal), np.inf) + if info > 0: + print 'Solve does not work well' + print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs) + elif self.myTest == 'xcJ': + #TODO: fix the null space + xc, info = sp.linalg.minres(A, rhs, tol = 1e-6) + j = McI*(G*xc + P*phi_bc) + err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf) + if info > 0: + print 'Solve does not work well' + print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs) + return err + + def test_orderJ(self): + self.name = "1D - InhomogeneousMixed_Forward j" + self.myTest = 'j' + self.orderTest() + + def test_orderQ(self): + self.name = "1D - InhomogeneousMixed_Forward q" + self.myTest = 'q' + self.orderTest() + + def test_orderXJ(self): + self.name = "1D - InhomogeneousMixed_Inverse J" + self.myTest = 'xcJ' + self.orderTest() + +class Test2D_InhomogeneousMixed(OrderTest): + name = "2D - Mixed" + meshTypes = MESHTYPES + meshDimension = 2 + expectedOrders = 2 + meshSizes = [2, 4, 8, 16] + # meshSizes = [4] + + def getError(self): + #Test function + phi = lambda x: np.cos(0.5*np.pi*x[:,0])*np.cos(0.5*np.pi*x[:,1]) + j_funX = lambda x: -0.5*np.pi*np.sin(0.5*np.pi*x[:,0])*np.cos(0.5*np.pi*x[:,1]) + j_funY = lambda x: -0.5*np.pi*np.cos(0.5*np.pi*x[:,0])*np.sin(0.5*np.pi*x[:,1]) + q_fun = lambda x: -2*((0.5*np.pi)**2)*phi(x) + + xc_anal = phi(self.M.gridCC) + q_anal = q_fun(self.M.gridCC) + jX_anal = j_funX(self.M.gridFx) + jY_anal = j_funY(self.M.gridFy) + j_anal = np.r_[jX_anal,jY_anal] + + #TODO: Check where our boundary conditions are CCx or Nx + + cxm,cxp,cym,cyp = self.M.cellBoundaryInd + fxm,fxp,fym,fyp = self.M.faceBoundaryInd + + gBFx = self.M.gridFx[(fxm|fxp),:] + gBFy = self.M.gridFy[(fym|fyp),:] + + gBCx = self.M.gridCC[(cxm|cxp),:] + gBCy = self.M.gridCC[(cym|cyp),:] + + phi_bc = phi(np.r_[gBCx,gBCy]) + j_bc = np.r_[j_funX(gBFx), j_funY(gBFy)] + + # P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2)) + + P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'neumann'], ['dirichlet', 'neumann']]) + + Mc = self.M.getFaceInnerProduct() + McI = Utils.sdInv(self.M.getFaceInnerProduct()) + V = Utils.sdiag(self.M.vol) + G = -Pin.T*Pin*self.M.faceDiv.T * V + D = self.M.faceDiv + j = McI*(G*xc_anal + P*phi_bc) + q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc + + # Rearrange if we know q to solve for x + A = V*D*Pin.T*Pin*McI*G + rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc + + if self.myTest == 'j': + err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf) + elif self.myTest == 'q': + err = np.linalg.norm((q-V*q_anal), np.inf) + elif self.myTest == 'xc': + #TODO: fix the null space + xc, info = sp.linalg.minres(A, rhs, tol = 1e-6) + err = np.linalg.norm((xc-xc_anal), np.inf) + if info > 0: + print 'Solve does not work well' + print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs) + elif self.myTest == 'xcJ': + #TODO: fix the null space + xc, info = sp.linalg.minres(A, rhs, tol = 1e-6) + j = McI*(G*xc + P*phi_bc) + err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf) + if info > 0: + print 'Solve does not work well' + print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs) + return err + + def test_orderJ(self): + self.name = "2D - InhomogeneousMixed_Forward j" + self.myTest = 'j' + self.orderTest() + + def test_orderQ(self): + self.name = "2D - InhomogeneousMixed_Forward q" + self.myTest = 'q' + self.orderTest() + + def test_orderXJ(self): + self.name = "2D - InhomogeneousMixed_Inverse J" + self.myTest = 'xcJ' + self.orderTest() if __name__ == '__main__': unittest.main() From aad54e64ddd70953f21911bbcb3bff9cdff9bbe5 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 22 Feb 2014 13:03:05 -0800 Subject: [PATCH 07/15] Updates to plotSlice to work on irregular grid --- SimPEG/Mesh/TensorView.py | 55 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 51 insertions(+), 4 deletions(-) diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index ef0957eb..e3967c4c 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -241,8 +241,27 @@ class TensorView(object): def plotSlice(self, v, vType='CC', normal='Z', ind=None, grid=False, view='real', ax=None, clim=None, showIt=False, + pcolorOpts={}, + streamOpts={'color':'k'}, gridOpts={'color':'k'} ): + + """ + Plots a slice of a 3D mesh. + + .. plot:: + + from SimPEG import * + mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2))) + M = Mesh.TensorMesh(mT) + q = np.zeros(M.vnC) + q[[4,4],[4,4],[2,6]]=[-1,1] + q = Utils.mkvc(q) + A = M.faceDiv*M.cellGrad + b = Solver(A).solve(q) + M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8}) + + """ viewOpts = ['real','imag','abs','vec'] normalOpts = ['X', 'Y', 'Z'] vTypeOpts = ['CC','F','E'] @@ -302,13 +321,31 @@ class TensorView(object): v = doSlice(v) if clim is None: clim = [v.min(),v.max()] - out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1]),) + out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1], **pcolorOpts),) elif view in ['vec']: U, V = doSlice(v) if clim is None: - clim = [v.min(),v.max()] - out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, 0.5*(U+V).T, vmin=clim[0], vmax=clim[1]),) - out += (plt.streamplot(tM.vectorCCx, tM.vectorCCy, U.T, V.T),) + uv = np.r_[mkvc(U), mkvc(V)] + uv = np.sqrt(uv**2) + clim = [uv.min(),uv.max()] + + # Matplotlib seems to not support irregular + # spaced vectors at the moment. So we will + # Interpolate down to a regular mesh at the + # smallest mesh size in this 2D slice. + nxi = int(tM.hx.sum()/tM.hx.min()) + nyi = int(tM.hy.sum()/tM.hy.min()) + tMi = self.__class__([np.ones(nxi)*tM.hx.sum()/nxi, + np.ones(nyi)*tM.hy.sum()/nyi]) + P = tM.getInterpolationMat(tMi.gridCC,'CC',zerosOutside=True) + Ui = P*mkvc(U) + Vi = P*mkvc(V) + Ui = tMi.r(Ui, 'CC', 'CC', 'M') + Vi = tMi.r(Vi, 'CC', 'CC', 'M') + # End Interpolation + + out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, np.sqrt(U**2+V**2).T, vmin=clim[0], vmax=clim[1], **pcolorOpts),) + out += (plt.streamplot(tMi.vectorCCx, tMi.vectorCCy, Ui.T, Vi.T, **streamOpts),) if grid: xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten() @@ -512,3 +549,13 @@ class TensorView(object): return animate(fig, animateFrame, frames=len(frames)) +if __name__ == '__main__': + from SimPEG import * + mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2))) + M = Mesh.TensorMesh(mT) + q = np.zeros(M.vnC) + q[[4,4],[4,4],[2,6]]=[-1,1] + q = Utils.mkvc(q) + A = M.faceDiv*M.cellGrad + b = Solver(A).solve(q) + M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8}) From d41e6dd13df263b5f87de261a9d6b43027e9922d Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 22 Feb 2014 13:26:03 -0800 Subject: [PATCH 08/15] updates to plotSlice --- SimPEG/Mesh/TensorView.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index e3967c4c..d345f0c3 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -264,7 +264,7 @@ class TensorView(object): """ viewOpts = ['real','imag','abs','vec'] normalOpts = ['X', 'Y', 'Z'] - vTypeOpts = ['CC','F','E'] + vTypeOpts = ['CC', 'CCv','F','E'] # Some user error checking assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts) @@ -296,11 +296,15 @@ class TensorView(object): def doSlice(v): if vType == 'CC': return getIndSlice(self.r(v,'CC','CC','M')) - # Now just deal with 'F' and 'E' - aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC') - v = getattr(self,aveOp)*v # average to cell centers (might be a vector) - if view == 'vec': + elif vType == 'CCv': v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M') + assert view == 'vec', 'Other types for CCv not yet supported' + else: + # Now just deal with 'F' and 'E' + aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC') + v = getattr(self,aveOp)*v # average to cell centers (might be a vector) + v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M') + if view == 'vec': outSlice = [] if 'X' not in normal: outSlice.append(getIndSlice(v[0])) if 'Y' not in normal: outSlice.append(getIndSlice(v[1])) @@ -357,6 +361,8 @@ class TensorView(object): ax.set_xlabel('y' if normal == 'X' else 'x') ax.set_ylabel('y' if normal == 'Z' else 'z') ax.set_title('Slice %d' % ind) + ax.set_xlim(*tM.vectorNx[[0,-1]]) + ax.set_ylim(*tM.vectorNy[[0,-1]]) if showIt: plt.show() return out From 33431c4e5f5e165dd773fd12a8c9e4714baba776 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 22 Feb 2014 13:33:10 -0800 Subject: [PATCH 09/15] minor bug fix. --- SimPEG/Mesh/TensorView.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index d345f0c3..5757b00b 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -349,7 +349,7 @@ class TensorView(object): # End Interpolation out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, np.sqrt(U**2+V**2).T, vmin=clim[0], vmax=clim[1], **pcolorOpts),) - out += (plt.streamplot(tMi.vectorCCx, tMi.vectorCCy, Ui.T, Vi.T, **streamOpts),) + out += (ax.streamplot(tMi.vectorCCx, tMi.vectorCCy, Ui.T, Vi.T, **streamOpts),) if grid: xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten() From 48b00400b3a59a7d238ffdf19432817c72f595dc Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 24 Feb 2014 12:24:35 -0800 Subject: [PATCH 10/15] Lowered tolerance on utility tests --- SimPEG/Tests/test_utils.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index 4a68b139..429c797c 100644 --- a/SimPEG/Tests/test_utils.py +++ b/SimPEG/Tests/test_utils.py @@ -3,6 +3,7 @@ from SimPEG.Utils import * from SimPEG import Mesh, np, sp from SimPEG.Tests import checkDerivative +TOL = 1e-8 class TestCheckDerivative(unittest.TestCase): @@ -104,7 +105,7 @@ class TestSequenceFunctions(unittest.TestCase): sp.hstack((sdiag(a[2]), sdiag(a[3]))))) Z2 = B*A - sp.eye(10, 10) - self.assertTrue(np.linalg.norm(Z2.todense().ravel(), 2) < 1e-12) + self.assertTrue(np.linalg.norm(Z2.todense().ravel(), 2) < TOL) a = [np.random.rand(5, 1) for i in range(9)] B = inv3X3BlockDiagonal(*a) @@ -115,7 +116,7 @@ class TestSequenceFunctions(unittest.TestCase): Z3 = B*A - sp.eye(15, 15) - self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-12) + self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < TOL) def test_invPropertyTensor2D(self): @@ -134,9 +135,9 @@ class TestSequenceFunctions(unittest.TestCase): B2 = invPropertyTensor(M, prop, returnMatrix=True) Z = B1*A - sp.identity(M.nC*2) - self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12) + self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL) Z = B2*A - sp.identity(M.nC*2) - self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12) + self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL) def test_invPropertyTensor3D(self): @@ -158,9 +159,9 @@ class TestSequenceFunctions(unittest.TestCase): B2 = invPropertyTensor(M, prop, returnMatrix=True) Z = B1*A - sp.identity(M.nC*3) - self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12) + self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL) Z = B2*A - sp.identity(M.nC*3) - self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12) + self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL) if __name__ == '__main__': From 90fa9c372794a7c36a85a1bb50a611710819636b Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 24 Feb 2014 14:38:19 -0800 Subject: [PATCH 11/15] Added check on modelPair in Problem. --- SimPEG/Problem.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index c6f381f5..f93c0c7a 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -1,5 +1,5 @@ import Utils, Data, numpy as np, scipy.sparse as sp - +import Model class BaseProblem(object): """ @@ -39,10 +39,12 @@ class BaseProblem(object): counter = None #: A SimPEG.Utils.Counter object dataPair = Data.BaseData + modelPair = Model.BaseModel def __init__(self, mesh, model, *args, **kwargs): Utils.setKwargs(self, **kwargs) self.mesh = mesh + assert isinstance(d, self.modelPair), "Model object must be an instance of a %s class."%(self.modelPair.__name__) self.model = model @property From aef96902e9ecae7748728a1d2b5ea289647265ad Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 25 Feb 2014 09:50:09 -0800 Subject: [PATCH 12/15] Minor updates. --- SimPEG/Mesh/DiffOperators.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 70a7e657..03536f88 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -291,8 +291,8 @@ class DiffOperators(object): """ if(type(BC) is str): - BC = [BC for _ in self.vnC] # Repeat the str self.dim times - elif(type(BC) is list): + BC = [BC]*self.dim + if(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.") @@ -473,6 +473,10 @@ class DiffOperators(object): # Dirichlet else """ + + if discretization is not 'CC': + raise NotImplementedError('Boundary conditions only implemented for CC discretization.') + if(type(BC) is str): BC = [BC for _ in self.vnC] # Repeat the str self.dim times elif(type(BC) is list): From 30fd9a3c3d23236861c63d02e27c69771bf45f7c Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 25 Feb 2014 09:51:39 -0800 Subject: [PATCH 13/15] Docs reorganization. --- docs/api_BaseMesh.rst | 8 ----- docs/api_Cyl1DMesh.rst | 8 ----- docs/api_DiffOperators.rst | 8 ----- docs/api_Forward.rst | 6 ++-- docs/api_InnerProducts.rst | 8 ----- docs/api_Inverse.rst | 8 ++--- docs/api_LogicallyOrthogonalMesh.rst | 10 ------ docs/api_Mesh.rst | 53 ++++++++++++++++++++++++++++ docs/api_Parameters.rst | 2 +- docs/api_TensorMesh.rst | 10 ------ docs/api_Tests.rst | 2 +- docs/api_Utils.rst | 17 +++------ docs/index.rst | 33 ++++++++--------- 13 files changed, 82 insertions(+), 91 deletions(-) delete mode 100644 docs/api_BaseMesh.rst delete mode 100644 docs/api_Cyl1DMesh.rst delete mode 100644 docs/api_DiffOperators.rst delete mode 100644 docs/api_InnerProducts.rst delete mode 100644 docs/api_LogicallyOrthogonalMesh.rst create mode 100644 docs/api_Mesh.rst delete mode 100644 docs/api_TensorMesh.rst diff --git a/docs/api_BaseMesh.rst b/docs/api_BaseMesh.rst deleted file mode 100644 index fcf0e8c1..00000000 --- a/docs/api_BaseMesh.rst +++ /dev/null @@ -1,8 +0,0 @@ -.. _api_BaseMesh: - -Base Mesh -********* - -.. automodule:: SimPEG.Mesh.BaseMesh - :members: - :undoc-members: diff --git a/docs/api_Cyl1DMesh.rst b/docs/api_Cyl1DMesh.rst deleted file mode 100644 index 573e6c5e..00000000 --- a/docs/api_Cyl1DMesh.rst +++ /dev/null @@ -1,8 +0,0 @@ -.. _api_Cyl1DMesh: - -Cylindrical 1D Mesh -******************* - -.. automodule:: SimPEG.Mesh.Cyl1DMesh - :members: - :undoc-members: diff --git a/docs/api_DiffOperators.rst b/docs/api_DiffOperators.rst deleted file mode 100644 index a3a2ca98..00000000 --- a/docs/api_DiffOperators.rst +++ /dev/null @@ -1,8 +0,0 @@ -.. _api_DiffOperators: - -Differential Operators -********************** - -.. automodule:: SimPEG.Mesh.DiffOperators - :members: - :undoc-members: diff --git a/docs/api_Forward.rst b/docs/api_Forward.rst index 91009294..067a902f 100644 --- a/docs/api_Forward.rst +++ b/docs/api_Forward.rst @@ -2,7 +2,7 @@ Model -***** +===== .. automodule:: SimPEG.Model :show-inheritance: @@ -11,7 +11,7 @@ Model :inherited-members: Data -**** +==== .. automodule:: SimPEG.Data :show-inheritance: @@ -20,7 +20,7 @@ Data :inherited-members: Problem -******* +======= .. automodule:: SimPEG.Problem :show-inheritance: diff --git a/docs/api_InnerProducts.rst b/docs/api_InnerProducts.rst deleted file mode 100644 index e2f54c8e..00000000 --- a/docs/api_InnerProducts.rst +++ /dev/null @@ -1,8 +0,0 @@ -.. _api_InnerProducts: - -Inner Products -************** - -.. automodule:: SimPEG.Mesh.InnerProducts - :members: - :undoc-members: diff --git a/docs/api_Inverse.rst b/docs/api_Inverse.rst index 21f4386e..ab711f81 100644 --- a/docs/api_Inverse.rst +++ b/docs/api_Inverse.rst @@ -2,7 +2,7 @@ Regularization -************** +============== .. automodule:: SimPEG.Regularization :show-inheritance: @@ -11,7 +11,7 @@ Regularization Objective Function -****************** +================== .. automodule:: SimPEG.ObjFunction :members: @@ -19,7 +19,7 @@ Objective Function Optimize -******** +======== .. automodule:: SimPEG.Optimization :show-inheritance: @@ -28,7 +28,7 @@ Optimize :undoc-members: Inversion -********* +========= .. automodule:: SimPEG.Inversion :show-inheritance: diff --git a/docs/api_LogicallyOrthogonalMesh.rst b/docs/api_LogicallyOrthogonalMesh.rst deleted file mode 100644 index 449ed6f8..00000000 --- a/docs/api_LogicallyOrthogonalMesh.rst +++ /dev/null @@ -1,10 +0,0 @@ -.. _api_LogicallyOrthogonalMesh: - -Logically Orthogonal Mesh -************************* - -.. automodule:: SimPEG.Mesh.LogicallyOrthogonalMesh - :show-inheritance: - :members: - :undoc-members: - :inherited-members: diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst new file mode 100644 index 00000000..a964b24d --- /dev/null +++ b/docs/api_Mesh.rst @@ -0,0 +1,53 @@ +.. _api_Mesh: + + +Tensor Mesh +=========== + +.. automodule:: SimPEG.Mesh.TensorMesh + :show-inheritance: + :members: + :undoc-members: + :inherited-members: + + +Cylindrical 1D Mesh +=================== + +.. automodule:: SimPEG.Mesh.Cyl1DMesh + :members: + :undoc-members: + + +Logically Orthogonal Mesh +========================= + +.. automodule:: SimPEG.Mesh.LogicallyOrthogonalMesh + :show-inheritance: + :members: + :undoc-members: + :inherited-members: + + + +Base Mesh +========= + +.. automodule:: SimPEG.Mesh.BaseMesh + :members: + :undoc-members: + + +Inner Products +============== + +.. automodule:: SimPEG.Mesh.InnerProducts + :members: + :undoc-members: + +Differential Operators +====================== + +.. automodule:: SimPEG.Mesh.DiffOperators + :members: + :undoc-members: diff --git a/docs/api_Parameters.rst b/docs/api_Parameters.rst index e90b58a5..1a333b88 100644 --- a/docs/api_Parameters.rst +++ b/docs/api_Parameters.rst @@ -2,7 +2,7 @@ Parameters -********** +========== .. automodule:: SimPEG.Parameters :show-inheritance: diff --git a/docs/api_TensorMesh.rst b/docs/api_TensorMesh.rst deleted file mode 100644 index 9cf80c41..00000000 --- a/docs/api_TensorMesh.rst +++ /dev/null @@ -1,10 +0,0 @@ -.. _api_TensorMesh: - -Tensor Mesh -*********** - -.. automodule:: SimPEG.Mesh.TensorMesh - :show-inheritance: - :members: - :undoc-members: - :inherited-members: diff --git a/docs/api_Tests.rst b/docs/api_Tests.rst index 0ab45c51..614d8747 100644 --- a/docs/api_Tests.rst +++ b/docs/api_Tests.rst @@ -1,7 +1,7 @@ .. _api_Tests: Testing SimPEG -************** +============== .. automodule:: SimPEG.Tests.TestUtils :members: diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst index 9c40cabe..641e6d60 100644 --- a/docs/api_Utils.rst +++ b/docs/api_Utils.rst @@ -17,42 +17,35 @@ Utilities :undoc-members: Matrix Utilities -**************** +================ .. automodule:: SimPEG.Utils.matutils :members: :undoc-members: -Sparse Utilities -**************** - -.. automodule:: SimPEG.Utils.sputils - :members: - :undoc-members: - LOM Utilities -************* +============= .. automodule:: SimPEG.Utils.lomutils :members: :undoc-members: Mesh Utilities -************** +============== .. automodule:: SimPEG.Utils.meshutils :members: :undoc-members: Model Builder Utilities -*********************** +======================= .. automodule:: SimPEG.Utils.ModelBuilder :members: :undoc-members: Interpolation Utilities -*********************** +======================= .. automodule:: SimPEG.Utils.interputils :members: diff --git a/docs/index.rst b/docs/index.rst index 9ed98291..84e0f925 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -11,10 +11,12 @@ The vision is to create a package for finite volume simulation and parameter est applications to geophysical imaging and subsurface flow. To enable these goals, this package has the following features: -* is modular with respect to discretization, physics, optimization, and regularization -* is built with the (large-scale) inverse problem in mind -* provides a framework for geophysical and hydrogeologic problems -* supports 1D, 2D and 3D problems + + * is modular with respect to discretization, physics, optimization, and regularization + * is built with the (large-scale) inverse problem in mind + * provides a framework for geophysical and hydrogeologic problems + * supports 1D, 2D and 3D problems + .. image:: simpeg-framework.png :width: 400 px @@ -22,20 +24,15 @@ these goals, this package has the following features: :align: center Meshing & Operators -=================== +******************* .. toctree:: :maxdepth: 2 - api_BaseMesh - api_TensorMesh - api_LogicallyOrthogonalMesh - api_Cyl1DMesh - api_DiffOperators - api_InnerProducts + api_Mesh Forward Problems -================ +**************** .. toctree:: :maxdepth: 2 @@ -43,7 +40,7 @@ Forward Problems api_Forward Inversion -========= +********* .. toctree:: :maxdepth: 2 @@ -52,7 +49,7 @@ Inversion api_Parameters Testing SimPEG -============== +************** .. toctree:: :maxdepth: 2 @@ -60,20 +57,20 @@ Testing SimPEG api_Tests * Master Branch - .. image:: https://travis-ci.org/simpeg/simpeg.png?branch=master + .. image:: https://travis-ci.org/simpeg/simpeg.png?branch*master :target: https://travis-ci.org/simpeg/simpeg :alt: Master Branch :align: center * Develop Branch - .. image:: https://travis-ci.org/simpeg/simpeg.png?branch=develop + .. image:: https://travis-ci.org/simpeg/simpeg.png?branch*develop :target: https://travis-ci.org/simpeg/simpeg :alt: Develop Branch :align: center Utility Codes -============= +************* .. toctree:: :maxdepth: 2 @@ -82,7 +79,7 @@ Utility Codes Project Index & Search -====================== +********************** * :ref:`genindex` * :ref:`modindex` From e0ea38ff855d2675e1adec983dc3aed7198c0f13 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 25 Feb 2014 12:14:59 -0800 Subject: [PATCH 14/15] Added nonLinearModel Object --- SimPEG/Model.py | 62 +++++++++++++++++++++++++++++++++++++++++++++++ SimPEG/Problem.py | 2 +- 2 files changed, 63 insertions(+), 1 deletion(-) diff --git a/SimPEG/Model.py b/SimPEG/Model.py index f6d4d6d4..ce2349d0 100644 --- a/SimPEG/Model.py +++ b/SimPEG/Model.py @@ -64,6 +64,68 @@ class BaseModel(object): m = self.example() return checkDerivative(lambda m : [self.transform(m), self.transformDeriv(m)], m, plotIt=False) +class BaseNonLinearModel(object): + """ + SimPEG BaseNonLinearModel + + """ + + __metaclass__ = Utils.SimPEGMetaClass + + counter = None #: A SimPEG.Utils.Counter object + mesh = None #: A SimPEG Mesh + + def __init__(self, mesh): + self.mesh = mesh + + def transform(self, u, m): + """ + :param numpy.array u: fields + :param numpy.array m: model + :rtype: numpy.array + :return: transformed model + + The *transform* changes the model into the physical property. + + """ + return m + + def transformDerivU(self, u, m): + """ + :param numpy.array u: fields + :param numpy.array m: model + :rtype: scipy.csr_matrix + :return: derivative of transformed model + + The *transform* changes the model into the physical property. + The *transformDerivU* provides the derivative of the *transform* with respect to the fields. + """ + raise NotImplementedError('The transformDerivU is not implemented.') + + + def transformDerivM(self, u, m): + """ + :param numpy.array u: fields + :param numpy.array m: model + :rtype: scipy.csr_matrix + :return: derivative of transformed model + + The *transform* changes the model into the physical property. + The *transformDerivU* provides the derivative of the *transform* with respect to the model. + """ + raise NotImplementedError('The transformDerivM is not implemented.') + + @property + def nP(self): + """Number of parameters in the model.""" + return self.mesh.nC + + def example(self): + raise NotImplementedError('The example is not implemented.') + + def test(self, m=None): + raise NotImplementedError('The test is not implemented.') + class LogModel(BaseModel): """SimPEG LogModel""" diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index c6f381f5..7ea8141c 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -40,7 +40,7 @@ class BaseProblem(object): dataPair = Data.BaseData - def __init__(self, mesh, model, *args, **kwargs): + def __init__(self, mesh, model, **kwargs): Utils.setKwargs(self, **kwargs) self.mesh = mesh self.model = model From 425d472f31c2a039c4b78d8ad9ac3bed25c8264f Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 26 Feb 2014 09:16:41 -0800 Subject: [PATCH 15/15] updates to problem and data. --- SimPEG/Data.py | 20 +++++++------------- SimPEG/Problem.py | 2 +- 2 files changed, 8 insertions(+), 14 deletions(-) diff --git a/SimPEG/Data.py b/SimPEG/Data.py index 0290be36..12696406 100644 --- a/SimPEG/Data.py +++ b/SimPEG/Data.py @@ -52,6 +52,7 @@ class BaseData(object): instead of recalculating the fields (which may be expensive!). .. math:: + d_\\text{pred} = P(u(m)) Where P is a projection of the fields onto the data space. @@ -67,29 +68,22 @@ class BaseData(object): .. math:: + d_\\text{pred} = \mathbf{P} u(m) """ return u - @Utils.count - def projectFieldsAdjoint(self, d): + def projectFieldsDeriv(self, u): """ - This function is the adjoint of the projection. - **projectFieldsAdjoint** is used in the - calculation of the sensitivities. + This function projects the fields onto the data space. + .. math:: - u = \mathbf{P}^\\top d - :param numpy.array d: data - :param numpy.array u: fields (ish) - :rtype: fields like object - :return: data + \\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P} """ - return d - - #TODO: def projectFieldDeriv(self, u): Does this need to be made??! + return sp.identity(u.size) @Utils.count def residual(self, m, u=None): diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index f93c0c7a..350ff93a 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -44,7 +44,7 @@ class BaseProblem(object): def __init__(self, mesh, model, *args, **kwargs): Utils.setKwargs(self, **kwargs) self.mesh = mesh - assert isinstance(d, self.modelPair), "Model object must be an instance of a %s class."%(self.modelPair.__name__) + assert isinstance(model, self.modelPair), "Model object must be an instance of a %s class."%(self.modelPair.__name__) self.model = model @property