From 002554209b0aeecea262f3aedc67cfa661138545 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 21 Jul 2016 22:28:20 -0700 Subject: [PATCH] listing and remove wildcard inputs from tes_Mixed_boundaryPoisson --- tests/mesh/test_Mixed_boundaryPoisson.py | 222 ++++++++++++++--------- 1 file changed, 137 insertions(+), 85 deletions(-) diff --git a/tests/mesh/test_Mixed_boundaryPoisson.py b/tests/mesh/test_Mixed_boundaryPoisson.py index 3aa1dbcd..3ba2e692 100644 --- a/tests/mesh/test_Mixed_boundaryPoisson.py +++ b/tests/mesh/test_Mixed_boundaryPoisson.py @@ -2,12 +2,12 @@ import numpy as np import scipy.sparse as sp import unittest import matplotlib.pyplot as plt -from SimPEG import * +from SimPEG import Mesh, Tests, Utils, Solver MESHTYPES = ['uniformTensorMesh'] + def getxBCyBC_CC(mesh, alpha, beta, gamma): -# def getxBCyBC(mesh, alpha, beta, gamma): """ This is a subfunction generating mixed-boundary condition: @@ -17,17 +17,19 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): \rho \vec{j} = -\nabla \phi \phi - \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega + \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r + = \partial \Omega xBC = f_1(\alpha, \beta, \gamma) yBC = f(\alpha, \beta, \gamma) Computes xBC and yBC for cell-centered discretizations """ - if mesh.dim == 1: #1D + + 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 + fCCxm, fCCxp = mesh.cellBoundaryInd nBC = fCCxm.sum()+fCCxp.sum() h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] @@ -50,11 +52,11 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): xBC = np.r_[xBC_xm, xBC_xp] yBC = np.r_[yBC_xm, yBC_xp] - elif mesh.dim == 2: #2D + 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") - fxm,fxp,fym,fyp = mesh.faceBoundaryInd + fxm, fxp, fym, fyp = mesh.faceBoundaryInd nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] @@ -65,8 +67,10 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): # 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_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) - h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) + h_xm = mesh.hx[0]*np.ones_like(alpha_xm) + h_xp = mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym = mesh.hy[0]*np.ones_like(alpha_ym) + h_yp = mesh.hy[-1]*np.ones_like(alpha_yp) 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) @@ -87,8 +91,10 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): 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]]) + 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] @@ -98,11 +104,11 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): xBC = np.r_[xBC_x, xBC_y] yBC = np.r_[yBC_x, yBC_y] - elif mesh.dim == 3: #3D + 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 + fxm, fxp, fym, fyp, fzm, fzp = mesh.faceBoundaryInd nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] @@ -116,9 +122,12 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] - h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) - h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) - h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp) + h_xm = mesh.hx[0]*np.ones_like(alpha_xm) + h_xp = mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym = mesh.hy[0]*np.ones_like(alpha_ym) + h_yp = mesh.hy[-1]*np.ones_like(alpha_yp) + h_zm = mesh.hz[0]*np.ones_like(alpha_zm) + h_zp = mesh.hz[-1]*np.ones_like(alpha_zp) 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) @@ -148,9 +157,12 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): 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]]) + 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] @@ -165,6 +177,7 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma): return xBC, yBC + class Test1D_InhomogeneousMixed(Tests.OrderTest): name = "1D - Mixed" meshTypes = MESHTYPES @@ -173,11 +186,14 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest): 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) + # Test function + def phi_fun(x): return np.cos(np.pi*x) + + def j_fun(x): return np.pi*np.sin(np.pi*x) + + def phi_deriv(x): return -j_fun(x) + + def q_fun(x): return (np.pi**2)*np.cos(np.pi*x) xc_ana = phi_fun(self.M.gridCC) q_ana = q_fun(self.M.gridCC) @@ -189,17 +205,16 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest): # Setup Mixed B.C (alpha, beta, gamma) alpha_xm, alpha_xp = 1., 1. - beta_xm, beta_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]]) + 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_CC(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) @@ -214,7 +229,7 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest): A = Div*MfrhoI*G if self.myTest == 'xc': - #TODO: fix the null space + # TODO: fix the null space Ainv = Solver(A) xc = Ainv*rhs err = np.linalg.norm((xc-xc_ana), np.inf) @@ -222,13 +237,13 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest): 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 @@ -237,40 +252,59 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest): 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) + # Test function + def phi_fun(x): + return np.cos(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) + + def j_funX(x): + return +np.pi*np.sin(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) + + def j_funY(x): + return +np.pi*np.cos(np.pi*x[:, 0])*np.sin(np.pi*x[:, 1]) + + def phideriv_funX(x): + return -j_funX(x) + + def phideriv_funY(x): + return -j_funY(x) + + def q_fun(x): + return +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] + 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,:] + 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]) + alpha_xm = np.ones_like(gBFxm[:, 0]) + alpha_xp = np.ones_like(gBFxp[:, 0]) + beta_xm = np.ones_like(gBFxm[:, 0]) + beta_xp = np.ones_like(gBFxp[:, 0]) + alpha_ym = np.ones_like(gBFym[:, 1]) + alpha_yp = np.ones_like(gBFyp[:, 1]) + beta_ym = np.ones_like(gBFym[:, 1]) + beta_yp = 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) + phiderivX_bc_xm = phideriv_funX(gBFxm) + phiderivX_bc_xp = phideriv_funX(gBFxp) + phiderivY_bc_ym = phideriv_funY(gBFym) + phiderivY_bc_yp = phideriv_funY(gBFyp) + + def gamma_fun(alpha, beta, phi, phi_deriv): + return alpha*phi + beta*phi_deriv - 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) @@ -282,7 +316,6 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest): x_BC, y_BC = getxBCyBC_CC(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) @@ -303,13 +336,13 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest): 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 @@ -318,51 +351,74 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest): 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]) + # Test function + def phi_fun(x): + return (np.cos(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) * + np.cos(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) + def j_funX(x): + return (np.pi*np.sin(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) * + np.cos(np.pi*x[:, 2])) - q_fun = lambda x: 3*(np.pi**2)*phi_fun(x) + def j_funY(x): + return (np.pi*np.cos(np.pi*x[:, 0])*np.sin(np.pi*x[:, 1]) * + np.cos(np.pi*x[:, 2])) + + def j_funZ(x): + return (np.pi*np.cos(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) * + np.sin(np.pi*x[:, 2])) + + def phideriv_funX(x): return -j_funX(x) + + def phideriv_funY(x): return -j_funY(x) + + def phideriv_funZ(x): return -j_funZ(x) + + def q_fun(x): return 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] + 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,:] + 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[:,2]), np.ones_like(gBFzp[:,2]) - beta_zm, beta_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) - + alpha_xm = np.ones_like(gBFxm[:, 0]) + alpha_xp = np.ones_like(gBFxp[:, 0]) + beta_xm = np.ones_like(gBFxm[:, 0]) + beta_xp = np.ones_like(gBFxp[:, 0]) + alpha_ym = np.ones_like(gBFym[:, 1]) + alpha_yp = np.ones_like(gBFyp[:, 1]) + beta_ym = np.ones_like(gBFym[:, 1]) + beta_yp = np.ones_like(gBFyp[:, 1]) + alpha_zm = np.ones_like(gBFzm[:, 2]) + alpha_zp = np.ones_like(gBFzp[:, 2]) + beta_zm = np.ones_like(gBFzm[:, 2]) + beta_zp = np.ones_like(gBFzp[:, 2]) 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_funZ(gBFzm), phideriv_funZ(gBFzp) + phiderivX_bc_xm = phideriv_funX(gBFxm) + phiderivX_bc_xp = phideriv_funX(gBFxp) + phiderivY_bc_ym = phideriv_funY(gBFym) + phiderivY_bc_yp = phideriv_funY(gBFyp) + phiderivY_bc_zm = phideriv_funZ(gBFzm) + phiderivY_bc_zp = phideriv_funZ(gBFzp) + + def gamma_fun(alpha, beta, phi, phi_deriv): + return alpha*phi + beta*phi_deriv - 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) @@ -376,7 +432,6 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest): x_BC, y_BC = getxBCyBC_CC(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) @@ -390,7 +445,7 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest): A = Div*MfrhoI*G if self.myTest == 'xc': - #TODO: fix the null space + # TODO: fix the null space Ainv = Solver(A) xc = Ainv*rhs err = np.linalg.norm((xc-xc_ana), np.inf) @@ -398,14 +453,11 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest): 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()