mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-28 21:47:21 +08:00
Testing boundary conditoins
- I am not sure what happend... but this is what I did for last push
This commit is contained in:
seogi
@@ -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()
|
||||
|
||||
Reference in New Issue
Block a user