From 72fe771071e088b3e9fe8d4b42009351ea604f2c Mon Sep 17 00:00:00 2001 From: seogi Date: Wed, 19 Feb 2014 15:24:37 -0800 Subject: [PATCH] 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()