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/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 1712964e..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.") @@ -460,6 +460,108 @@ 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 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): + 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[1] == '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 21662d8f..8a80ffea 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -442,6 +442,56 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): Q[indZeros, :] = 0 return Q.tocsr() + + @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 5757b00b..4c04cda9 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/Problem.py b/SimPEG/Problem.py index 7ea8141c..414d0c33 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, **kwargs): Utils.setKwargs(self, **kwargs) self.mesh = mesh + assert isinstance(model, self.modelPair), "Model object must be an instance of a %s class."%(self.modelPair.__name__) self.model = model @property 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..5e028dee --- /dev/null +++ b/SimPEG/Tests/test_boundaryPoisson.py @@ -0,0 +1,501 @@ +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 + + 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()) + 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((j-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 + 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 + xc = Solver(A).solve(rhs) + 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 + + 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() + +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() 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__':