wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-27 18:44:32 +08:00
Code Issues Packages Projects Releases Wiki Activity

listing and remove wildcard inputs from tes_Mixed_boundaryPoisson

Browse Source
This commit is contained in:
Lindsey Heagy
2016-07-21 22:28:20 -07:00
parent 76f140f89b
commit 002554209b
1 changed files with 137 additions and 85 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/mesh/test_Mixed_boundaryPoisson.py
+137 -85
View File
@@ -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()