From 2c09be9fc12a8046feadcc240f25b74dd2d439ae Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 21 Apr 2016 14:44:37 -0700 Subject: [PATCH] Working Mixed boundary conditions and testing ... --- tests/mesh/test_Mixed_boundaryPoisson.py | 395 +++++++++++++++++++++++ 1 file changed, 395 insertions(+) create mode 100644 tests/mesh/test_Mixed_boundaryPoisson.py diff --git a/tests/mesh/test_Mixed_boundaryPoisson.py b/tests/mesh/test_Mixed_boundaryPoisson.py new file mode 100644 index 00000000..7a15f36d --- /dev/null +++ b/tests/mesh/test_Mixed_boundaryPoisson.py @@ -0,0 +1,395 @@ +import numpy as np +import scipy.sparse as sp +import unittest +import matplotlib.pyplot as plt +from SimPEG import * + +MESHTYPES = ['uniformTensorMesh'] + +def getxBCyBC(mesh, alpha, beta, gamma): +# def getxBCyBC(mesh, alpha, beta, gamma): + """ + + """ + if mesh.dim == 1: #1D + if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2): + raise Exception("Lenght of list, alpha should be 2") + fCCxm,fCCxp = mesh.cellBoundaryInd + nBC = fCCxm.sum()+fCCxp.sum() + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + + xBC = np.r_[xBC_xm, xBC_xp] + yBC = np.r_[yBC_xm, yBC_xp] + + elif mesh.dim == 2: #2D + if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4): + raise Exception("Lenght of list, alpha should be 4") + + fCCxm,fCCxp,fCCym,fCCyp = mesh.cellBoundaryInd + fxm,fxp,fym,fyp = mesh.faceBoundaryInd + nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] + alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] + + h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym) + b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym) + a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp) + b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + xBC_ym = 0.5*a_ym + xBC_yp = 0.5*a_yp/b_yp + yBC_ym = 0.5*(1.-b_ym) + yBC_yp = 0.5*(1.-1./b_yp) + + sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) + sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) + + xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx] + xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] + yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx] + yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy] + + xBC = np.r_[xBC_x, xBC_y] + yBC = np.r_[yBC_x, yBC_y] + + elif mesh.dim == 3: #3D + if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6): + raise Exception("Lenght of list, alpha should be 6") + fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd + fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd + nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum() + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp] + h_zm, h_zp = mesh.gridCC[fCCzm], mesh.gridCC[fCCzp] + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] + alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] + alpha_zm, beta_zm, gamma_zm = alpha[2], beta[2], gamma[2] + alpha_zp, beta_zp, gamma_zp = alpha[3], beta[3], gamma[3] + + h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym) + b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym) + a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp) + b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp) + + a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm) + b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm) + a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp) + b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + xBC_ym = 0.5*a_ym + xBC_yp = 0.5*a_yp/b_yp + yBC_ym = 0.5*(1.-b_ym) + yBC_yp = 0.5*(1.-1./b_yp) + xBC_zm = 0.5*a_zm + xBC_zp = 0.5*a_zp/b_zp + yBC_zm = 0.5*(1.-b_zm) + yBC_zp = 0.5*(1.-1./b_zp) + + sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) + sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) + sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]]) + + xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx] + xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] + xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz] + + yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx] + yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy] + yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz] + + xBC = np.r_[xBC_x, xBC_y, xBC_z] + yBC = np.r_[yBC_x, yBC_y, yBC_z] + + return xBC, yBC + +class Test1D_InhomogeneousMixed(Tests.OrderTest): + name = "1D - Mixed" + meshTypes = MESHTYPES + meshDimension = 1 + expectedOrders = 2 + meshSizes = [4, 8, 16, 32] + + def getError(self): + #Test function + phi_fun = lambda x: np.cos(np.pi*x) + j_fun = lambda x: np.pi*np.sin(np.pi*x) + phi_deriv = lambda x: -j_fun(x) + q_fun = lambda x: (np.pi**2)*np.cos(np.pi*x) + + xc_ana = phi_fun(self.M.gridCC) + q_ana = q_fun(self.M.gridCC) + j_ana = j_fun(self.M.gridFx) + + # Get boundary locations + vecN = self.M.vectorNx + vecC = self.M.vectorCCx + + # Setup Mixed B.C (alpha, beta, gamma) + alpha_xm, alpha_xp = 1., 1. + beta_xm, beta_xp = 1., 1. + alpha = np.r_[alpha_xm, alpha_xp] + beta = np.r_[beta_xm, beta_xp] + vecN = self.M.vectorNx + vecC = self.M.vectorCCx + phi_bc = phi_fun(vecN[[0,-1]]) + phi_deriv_bc = phi_deriv(vecN[[0,-1]]) + gamma = alpha*phi_bc + beta*phi_deriv_bc + x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + + + sigma = np.ones(self.M.nC) + Mfrho = self.M.getFaceInnerProduct(1./sigma) + MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) + V = Utils.sdiag(self.M.vol) + Div = V*self.M.faceDiv + P_BC, B = self.M.getBCProjWF_simple() + q = q_fun(self.M.gridCC) + M = B*self.M.aveCC2F + G = Div.T - P_BC*Utils.sdiag(y_BC)*M + rhs = V*q + Div*MfrhoI*P_BC*x_BC + A = Div*MfrhoI*G + + if self.myTest == 'xc': + #TODO: fix the null space + Ainv = Solver(A) + xc = Ainv*rhs + err = np.linalg.norm((xc-xc_ana), np.inf) + else: + NotImplementedError + return err + + + def test_order(self): + print "==== Testing Mixed boudary conduction for CC-problem ====" + self.name = "1D" + self.myTest = 'xc' + self.orderTest() + +class Test2D_InhomogeneousMixed(Tests.OrderTest): + name = "2D - Mixed" + meshTypes = MESHTYPES + meshDimension = 2 + expectedOrders = 2 + meshSizes = [4, 8, 16, 32] + + def getError(self): + #Test function + phi_fun = 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]) + phideriv_funX = lambda x: -j_funX(x) + phideriv_funY = lambda x: -j_funY(x) + q_fun = lambda x: +2*(np.pi**2)*phi_fun(x) + + xc_ana = phi_fun(self.M.gridCC) + q_ana = q_fun(self.M.gridCC) + jX_ana = j_funX(self.M.gridFx) + jY_ana = j_funY(self.M.gridFy) + j_ana = np.r_[jX_ana,jY_ana] + + # Get boundary locations + fxm,fxp,fym,fyp = self.M.faceBoundaryInd + gBFxm = self.M.gridFx[fxm,:] + gBFxp = self.M.gridFx[fxp,:] + gBFym = self.M.gridFy[fym,:] + gBFyp = self.M.gridFy[fyp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + + phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp) + phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp) + + phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp) + phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp) + + gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv + gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm) + gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp) + gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym) + gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp) + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] + + x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + + + sigma = np.ones(self.M.nC) + Mfrho = self.M.getFaceInnerProduct(1./sigma) + MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) + V = Utils.sdiag(self.M.vol) + Div = V*self.M.faceDiv + P_BC, B = self.M.getBCProjWF_simple() + q = q_fun(self.M.gridCC) + M = B*self.M.aveCC2F + G = Div.T - P_BC*Utils.sdiag(y_BC)*M + rhs = V*q + Div*MfrhoI*P_BC*x_BC + A = Div*MfrhoI*G + + if self.myTest == 'xc': + Ainv = Solver(A) + xc = Ainv*rhs + err = np.linalg.norm((xc-xc_ana), np.inf) + else: + NotImplementedError + return err + + + def test_order(self): + print "==== Testing Mixed boudary conduction for CC-problem ====" + self.name = "2D" + self.myTest = 'xc' + self.orderTest() + +class Test3D_InhomogeneousMixed(Tests.OrderTest): + name = "3D - Mixed" + meshTypes = MESHTYPES + meshDimension = 3 + expectedOrders = 2 + meshSizes = [4, 8, 16] + + def getError(self): + #Test function + phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2]) + j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2]) + j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])*np.cos(np.pi*x[:,2]) + j_funZ = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.sin(np.pi*x[:,2]) + + phideriv_funX = lambda x: -j_funX(x) + phideriv_funY = lambda x: -j_funY(x) + phideriv_funZ = lambda x: -j_funZ(x) + + q_fun = lambda x: 3*(np.pi**2)*phi_fun(x) + + xc_ana = phi_fun(self.M.gridCC) + q_ana = q_fun(self.M.gridCC) + jX_ana = j_funX(self.M.gridFx) + jY_ana = j_funY(self.M.gridFy) + j_ana = np.r_[jX_ana,jY_ana,jY_ana] + + # Get boundary locations + fxm,fxp,fym,fyp,fzm,fzp = self.M.faceBoundaryInd + gBFxm = self.M.gridFx[fxm,:] + gBFxp = self.M.gridFx[fxp,:] + gBFym = self.M.gridFy[fym,:] + gBFyp = self.M.gridFy[fyp,:] + gBFzm = self.M.gridFz[fzm,:] + gBFzp = self.M.gridFz[fzp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + alpha_zm, alpha_zp = np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1]) + beta_zm, beta_zp = np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1]) + + + phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp) + phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp) + phi_bc_zm, phi_bc_zp = phi_fun(gBFzm), phi_fun(gBFzp) + + phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp) + phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp) + phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funY(gBFzm), phideriv_funY(gBFzp) + + gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv + gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm) + gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp) + gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym) + gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp) + gamma_zm = gamma_fun(alpha_zm, beta_zm, phi_bc_zm, phiderivY_bc_zm) + gamma_zp = gamma_fun(alpha_zp, beta_zp, phi_bc_zp, phiderivY_bc_zp) + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] + + x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma) + + + sigma = np.ones(self.M.nC) + Mfrho = self.M.getFaceInnerProduct(1./sigma) + MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) + V = Utils.sdiag(self.M.vol) + Div = V*self.M.faceDiv + P_BC, B = self.M.getBCProjWF_simple() + q = q_fun(self.M.gridCC) + M = B*self.M.aveCC2F + G = Div.T - P_BC*Utils.sdiag(y_BC)*M + rhs = V*q + Div*MfrhoI*P_BC*x_BC + A = Div*MfrhoI*G + + if self.myTest == 'xc': + #TODO: fix the null space + Ainv = Solver(A) + xc = Ainv*rhs + err = np.linalg.norm((xc-xc_ana), np.inf) + else: + NotImplementedError + return err + + + def test_order(self): + print "==== Testing Mixed boudary conduction for CC-problem ====" + self.name = "3D" + self.myTest = 'xc' + self.orderTest() + + + +if __name__ == '__main__': + unittest.main()