mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-31 12:40:46 +08:00
Merge branch 'dev' into maps/feat-parametrizedBlock
This commit is contained in:
Lindsey Heagy
@@ -1,6 +1,7 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG import Tests
|
||||
|
||||
|
||||
class MyPropMap(Maps.PropMap):
|
||||
@@ -187,6 +188,34 @@ class TestPropMaps(unittest.TestCase):
|
||||
|
||||
MyReciprocalPropMap([('sigma', iMap), ('mu', iMap)]) # This should be fine
|
||||
|
||||
def test_linked_derivs_sigma(self):
|
||||
mesh = Mesh.TensorMesh([4,5], x0='CC')
|
||||
|
||||
mapping = Maps.ExpMap(mesh)
|
||||
propmap = MyReciprocalPropMap([('rho', mapping)])
|
||||
|
||||
x0 = np.random.rand(mesh.nC)
|
||||
m = propmap(x0)
|
||||
|
||||
# test Sigma
|
||||
testme = lambda v: [1./(m.rhoMap*v), m.sigmaDeriv]
|
||||
print 'Testing Rho from Sigma'
|
||||
Tests.checkDerivative(testme, x0, dx=0.01*x0, num=5, plotIt=False)
|
||||
|
||||
def test_linked_derivs_rho(self):
|
||||
mesh = Mesh.TensorMesh([4,5], x0='CC')
|
||||
|
||||
mapping = Maps.ExpMap(mesh)
|
||||
propmap = MyReciprocalPropMap([('sigma', mapping)])
|
||||
|
||||
x0 = np.random.rand(mesh.nC)
|
||||
m = propmap(x0)
|
||||
|
||||
# test Sigma
|
||||
testme = lambda v: [1./(m.sigmaMap*v), m.rhoDeriv]
|
||||
print 'Testing Rho from Sigma'
|
||||
Tests.checkDerivative(testme, x0, dx=0.01*x0, num=5, plotIt=False)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
|
||||
@@ -0,0 +1,12 @@
|
||||
import os
|
||||
import glob
|
||||
import unittest
|
||||
|
||||
if __name__ == '__main__':
|
||||
test_file_strings = glob.glob('test_*.py')
|
||||
module_strings = [str[0:len(str)-3] for str in test_file_strings]
|
||||
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
|
||||
in module_strings]
|
||||
testSuite = unittest.TestSuite(suites)
|
||||
|
||||
unittest.TextTestRunner(verbosity=2).run(testSuite)
|
||||
@@ -0,0 +1,69 @@
|
||||
import unittest
|
||||
from SimPEG import Mesh, Utils, EM, Maps, np
|
||||
import SimPEG.EM.Static.DC as DC
|
||||
|
||||
class DCProblemAnalyticTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 12.5
|
||||
hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)]
|
||||
hy = [(cs,7, -1.3),(cs,20)]
|
||||
mesh = Mesh.TensorMesh([hx, hy],x0="CN")
|
||||
sighalf = 1e-2
|
||||
sigma = np.ones(mesh.nC)*sighalf
|
||||
x = np.linspace(-135, 250., 20)
|
||||
M = Utils.ndgrid(x-12.5, np.r_[0.])
|
||||
N = Utils.ndgrid(x+12.5, np.r_[0.])
|
||||
A0loc = np.r_[-150, 0.]
|
||||
A1loc = np.r_[-130, 0.]
|
||||
rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
|
||||
data_anal = EM.Analytics.DCAnalyticHalf(np.r_[A0loc, 0.], rxloc, sighalf, earth_type="halfspace")
|
||||
|
||||
rx = DC.Rx.Dipole_ky(M, N)
|
||||
src0 = DC.Src.Pole([rx], A0loc)
|
||||
survey = DC.Survey_ky([src0])
|
||||
|
||||
self.survey = survey
|
||||
self.mesh = mesh
|
||||
self.sigma = sigma
|
||||
self.data_anal = data_anal
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
self.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
self.Solver = SolverLU
|
||||
|
||||
def test_Problem3D_N(self):
|
||||
|
||||
problem = DC.Problem2D_N(self.mesh)
|
||||
problem.Solver = self.Solver
|
||||
problem.pair(self.survey)
|
||||
data = self.survey.dpred(self.sigma)
|
||||
err= np.linalg.norm((data-self.data_anal)/self.data_anal)**2 / self.data_anal.size
|
||||
if err < 0.05:
|
||||
passed = True
|
||||
print ">> DC analytic test for Problem3D_N is passed"
|
||||
else:
|
||||
passed = False
|
||||
print ">> DC analytic test for Problem3D_N is failed"
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_Problem3D_CC(self):
|
||||
problem = DC.Problem2D_CC(self.mesh)
|
||||
problem.Solver = self.Solver
|
||||
problem.pair(self.survey)
|
||||
data = self.survey.dpred(self.sigma)
|
||||
err= np.linalg.norm((data-self.data_anal)/self.data_anal)**2 / self.data_anal.size
|
||||
if err < 0.05:
|
||||
passed = True
|
||||
print ">> DC analytic test for Problem3D_CC is passed"
|
||||
else:
|
||||
passed = False
|
||||
print ">> DC analytic test for Problem3D_CC is failed"
|
||||
self.assertTrue(passed)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -0,0 +1,127 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import SimPEG.EM.Static.DC as DC
|
||||
|
||||
|
||||
class DCProblem_2DTestsCC(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 12.5
|
||||
hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)]
|
||||
hy = [(cs,7, -1.3),(cs,20)]
|
||||
mesh = Mesh.TensorMesh([hx, hy],x0="CN")
|
||||
x = np.linspace(-135, 250., 20)
|
||||
M = Utils.ndgrid(x-12.5, np.r_[0.])
|
||||
N = Utils.ndgrid(x+12.5, np.r_[0.])
|
||||
A0loc = np.r_[-150, 0.]
|
||||
A1loc = np.r_[-130, 0.]
|
||||
rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
|
||||
rx = DC.Rx.Dipole_ky(M, N)
|
||||
src0 = DC.Src.Pole([rx], A0loc)
|
||||
src1 = DC.Src.Pole([rx], A1loc)
|
||||
survey = DC.Survey_ky([src0, src1])
|
||||
problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
|
||||
problem.pair(survey)
|
||||
|
||||
mSynth = np.ones(mesh.nC)*1.
|
||||
survey.makeSyntheticData(mSynth)
|
||||
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-10
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
class DCProblemTestsN(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 12.5
|
||||
hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)]
|
||||
hy = [(cs,7, -1.3),(cs,20)]
|
||||
mesh = Mesh.TensorMesh([hx, hy],x0="CN")
|
||||
x = np.linspace(-135, 250., 20)
|
||||
M = Utils.ndgrid(x-12.5, np.r_[0.])
|
||||
N = Utils.ndgrid(x+12.5, np.r_[0.])
|
||||
A0loc = np.r_[-150, 0.]
|
||||
A1loc = np.r_[-130, 0.]
|
||||
rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
|
||||
rx = DC.Rx.Dipole_ky(M, N)
|
||||
src0 = DC.Src.Pole([rx], A0loc)
|
||||
src1 = DC.Src.Pole([rx], A1loc)
|
||||
survey = DC.Survey_ky([src0, src1])
|
||||
problem = DC.Problem2D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
|
||||
problem.pair(survey)
|
||||
|
||||
mSynth = np.ones(mesh.nC)*1.
|
||||
survey.makeSyntheticData(mSynth)
|
||||
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-8
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,71 @@
|
||||
import unittest
|
||||
from SimPEG import Mesh, Utils, EM, Maps, np
|
||||
import SimPEG.EM.Static.DC as DC
|
||||
|
||||
class DCProblemAnalyticTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 25.
|
||||
hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
|
||||
hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
|
||||
hz = [(cs,7, -1.3),(cs,20)]
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
|
||||
sigma = np.ones(mesh.nC)*1e-2
|
||||
|
||||
x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
|
||||
y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
|
||||
Aloc = np.r_[-200., 0., 0.]
|
||||
Bloc = np.r_[200., 0., 0.]
|
||||
M = Utils.ndgrid(x-25.,y, np.r_[0.])
|
||||
N = Utils.ndgrid(x+25.,y, np.r_[0.])
|
||||
phiA = EM.Analytics.DCAnalyticHalf(Aloc, [M,N], 1e-2, earth_type="halfspace")
|
||||
phiB = EM.Analytics.DCAnalyticHalf(Bloc, [M,N], 1e-2, earth_type="halfspace")
|
||||
data_anal = phiA-phiB
|
||||
|
||||
rx = DC.Rx.Dipole(M, N)
|
||||
src = DC.Src.Dipole([rx], Aloc, Bloc)
|
||||
survey = DC.Survey([src])
|
||||
|
||||
self.survey = survey
|
||||
self.mesh = mesh
|
||||
self.sigma = sigma
|
||||
self.data_anal = data_anal
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
self.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
self.Solver = SolverLU
|
||||
|
||||
def test_Problem3D_N(self):
|
||||
problem = DC.Problem3D_N(self.mesh)
|
||||
problem.Solver = self.Solver
|
||||
problem.pair(self.survey)
|
||||
data = self.survey.dpred(self.sigma)
|
||||
err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal)
|
||||
if err < 0.2:
|
||||
passed = True
|
||||
print ">> DC analytic test for Problem3D_N is passed"
|
||||
else:
|
||||
passed = False
|
||||
print ">> DC analytic test for Problem3D_N is failed"
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_Problem3D_CC(self):
|
||||
problem = DC.Problem3D_CC(self.mesh)
|
||||
problem.Solver = self.Solver
|
||||
problem.pair(self.survey)
|
||||
data = self.survey.dpred(self.sigma)
|
||||
err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal)
|
||||
if err < 0.2:
|
||||
passed = True
|
||||
print ">> DC analytic test for Problem3D_CC is passed"
|
||||
else:
|
||||
passed = False
|
||||
print ">> DC analytic test for Problem3D_CC is failed"
|
||||
self.assertTrue(passed)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -0,0 +1,127 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import SimPEG.EM.Static.DC as DC
|
||||
|
||||
|
||||
class DCProblemTestsCC(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
aSpacing=2.5
|
||||
nElecs=5
|
||||
|
||||
surveySize = nElecs*aSpacing - aSpacing
|
||||
cs = surveySize/nElecs/4
|
||||
|
||||
mesh = Mesh.TensorMesh([
|
||||
[(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
|
||||
[(cs,3, -1.3),(cs,3,1.3)],
|
||||
# [(cs,5, -1.3),(cs,10)]
|
||||
],'CN')
|
||||
|
||||
srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
|
||||
survey = DC.Survey(srcList)
|
||||
problem = DC.Problem3D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
|
||||
problem.pair(survey)
|
||||
|
||||
mSynth = np.ones(mesh.nC)
|
||||
survey.makeSyntheticData(mSynth)
|
||||
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-10
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
class DCProblemTestsN(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
aSpacing=2.5
|
||||
nElecs=10
|
||||
|
||||
surveySize = nElecs*aSpacing - aSpacing
|
||||
cs = surveySize/nElecs/4
|
||||
|
||||
mesh = Mesh.TensorMesh([
|
||||
[(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
|
||||
[(cs,3, -1.3),(cs,3,1.3)],
|
||||
# [(cs,5, -1.3),(cs,10)]
|
||||
],'CN')
|
||||
|
||||
srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
|
||||
survey = DC.Survey(srcList)
|
||||
problem = DC.Problem3D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
|
||||
problem.pair(survey)
|
||||
|
||||
mSynth = np.ones(mesh.nC)
|
||||
survey.makeSyntheticData(mSynth)
|
||||
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-8
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
|
||||
self.assertTrue(passed)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,96 @@
|
||||
import unittest
|
||||
from SimPEG import Mesh, Utils, EM, Maps, np
|
||||
import SimPEG.EM.Static.DC as DC
|
||||
import SimPEG.EM.Static.IP as IP
|
||||
|
||||
class IPProblemAnalyticTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 12.5
|
||||
npad=2
|
||||
hx = [(cs,npad, -1.3),(cs,21),(cs,npad, 1.3)]
|
||||
hy = [(cs,npad, -1.3),(cs,21),(cs,npad, 1.3)]
|
||||
hz = [(cs,npad, -1.3),(cs,20)]
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
|
||||
|
||||
x = mesh.vectorCCx[(mesh.vectorCCx>-80.)&(mesh.vectorCCx<80.)]
|
||||
y = mesh.vectorCCx[(mesh.vectorCCy>-80.)&(mesh.vectorCCy<80.)]
|
||||
Aloc = np.r_[-100., 0., 0.]
|
||||
Bloc = np.r_[100., 0., 0.]
|
||||
M = Utils.ndgrid(x-12.5,y, np.r_[0.])
|
||||
N = Utils.ndgrid(x+12.5,y, np.r_[0.])
|
||||
radius = 50.
|
||||
xc = np.r_[0., 0., -100]
|
||||
blkind = Utils.ModelBuilder.getIndicesSphere(xc, radius, mesh.gridCC)
|
||||
sigmaInf = np.ones(mesh.nC)*1e-2
|
||||
eta = np.zeros(mesh.nC)
|
||||
eta[blkind] = 0.1
|
||||
sigma0 = sigmaInf*(1.-eta)
|
||||
|
||||
rx = DC.Rx.Dipole(M, N)
|
||||
src = DC.Src.Dipole([rx], Aloc, Bloc)
|
||||
surveyDC = DC.Survey([src])
|
||||
|
||||
self.surveyDC = surveyDC
|
||||
self.mesh = mesh
|
||||
self.sigmaInf = sigmaInf
|
||||
self.sigma0 = sigma0
|
||||
self.src = src
|
||||
self.eta = eta
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
self.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
self.Solver = SolverLU
|
||||
|
||||
def test_Problem3D_N(self):
|
||||
|
||||
problemDC = DC.Problem3D_N(self.mesh)
|
||||
problemDC.Solver = self.Solver
|
||||
problemDC.pair(self.surveyDC)
|
||||
data0 = self.surveyDC.dpred(self.sigma0)
|
||||
finf = problemDC.fields(self.sigmaInf)
|
||||
datainf = self.surveyDC.dpred(self.sigmaInf, f=finf)
|
||||
problemIP = IP.Problem3D_N(self.mesh, sigma=self.sigmaInf, Ainv=problemDC.Ainv, f=finf)
|
||||
problemIP.Solver = self.Solver
|
||||
surveyIP = IP.Survey([self.src])
|
||||
problemIP.pair(surveyIP)
|
||||
data_full = data0 - datainf
|
||||
data = surveyIP.dpred(self.eta)
|
||||
err= np.linalg.norm((data-data_full)/data_full)**2 / data_full.size
|
||||
if err < 0.05:
|
||||
passed = True
|
||||
print ">> IP forward test for Problem3D_N is passed"
|
||||
else:
|
||||
passed = False
|
||||
print ">> IP forward test for Problem3D_N is failed"
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_Problem3D_CC(self):
|
||||
|
||||
problemDC = DC.Problem3D_CC(self.mesh)
|
||||
problemDC.Solver = self.Solver
|
||||
problemDC.pair(self.surveyDC)
|
||||
data0 = self.surveyDC.dpred(self.sigma0)
|
||||
finf = problemDC.fields(self.sigmaInf)
|
||||
datainf = self.surveyDC.dpred(self.sigmaInf, f=finf)
|
||||
problemIP = IP.Problem3D_CC(self.mesh, rho=1./self.sigmaInf, Ainv=problemDC.Ainv, f=finf)
|
||||
problemIP.Solver = self.Solver
|
||||
surveyIP = IP.Survey([self.src])
|
||||
problemIP.pair(surveyIP)
|
||||
data_full = data0 - datainf
|
||||
data = surveyIP.dpred(self.eta)
|
||||
err= np.linalg.norm((data-data_full)/data_full)**2 / data_full.size
|
||||
if err < 0.05:
|
||||
passed = True
|
||||
print ">> IP forward test for Problem3D_CC is passed"
|
||||
else:
|
||||
passed = False
|
||||
print ">> IP forward test for Problem3D_CC is failed"
|
||||
self.assertTrue(passed)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -0,0 +1,126 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import SimPEG.EM.Static.DC as DC
|
||||
import SimPEG.EM.Static.IP as IP
|
||||
|
||||
|
||||
class IPProblemTestsCC(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
aSpacing=2.5
|
||||
nElecs=5
|
||||
|
||||
surveySize = nElecs*aSpacing - aSpacing
|
||||
cs = surveySize/nElecs/4
|
||||
|
||||
mesh = Mesh.TensorMesh([
|
||||
[(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
|
||||
[(cs,3, -1.3),(cs,3,1.3)],
|
||||
# [(cs,5, -1.3),(cs,10)]
|
||||
],'CN')
|
||||
|
||||
srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
|
||||
survey = IP.Survey(srcList)
|
||||
sigma = np.ones(mesh.nC)
|
||||
problem = IP.Problem3D_CC(mesh, rho=1./sigma)
|
||||
problem.pair(survey)
|
||||
mSynth = np.ones(mesh.nC)*0.1
|
||||
survey.makeSyntheticData(mSynth)
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-10
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
class IPProblemTestsN(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
aSpacing=2.5
|
||||
nElecs=5
|
||||
|
||||
surveySize = nElecs*aSpacing - aSpacing
|
||||
cs = surveySize/nElecs/4
|
||||
|
||||
mesh = Mesh.TensorMesh([
|
||||
[(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
|
||||
[(cs,3, -1.3),(cs,3,1.3)],
|
||||
# [(cs,5, -1.3),(cs,10)]
|
||||
],'CN')
|
||||
|
||||
srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
|
||||
survey = IP.Survey(srcList)
|
||||
sigma = np.ones(mesh.nC)
|
||||
problem = IP.Problem3D_N(mesh, sigma=sigma)
|
||||
problem.pair(survey)
|
||||
mSynth = np.ones(mesh.nC)*0.1
|
||||
survey.makeSyntheticData(mSynth)
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-8
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,232 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import SimPEG
|
||||
from SimPEG import Mesh, Utils, EM, Maps, np, Survey
|
||||
from SimPEG.EM.Static import SIP, DC, IP
|
||||
from pymatsolver import MumpsSolver
|
||||
|
||||
|
||||
class IPProblemTestsCC(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 25.
|
||||
hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hz = [(cs,0, -1.3),(cs,20)]
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
|
||||
blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
|
||||
blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
|
||||
sigma = np.ones(mesh.nC)*1e-2
|
||||
eta = np.zeros(mesh.nC)
|
||||
tau = np.ones_like(sigma)*1.
|
||||
eta[blkind0] = 0.1
|
||||
eta[blkind1] = 0.1
|
||||
tau[blkind0] = 0.1
|
||||
tau[blkind1] = 0.01
|
||||
|
||||
x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
|
||||
y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
|
||||
Aloc = np.r_[-200., 0., 0.]
|
||||
Bloc = np.r_[200., 0., 0.]
|
||||
M = Utils.ndgrid(x-25.,y, np.r_[0.])
|
||||
N = Utils.ndgrid(x+25.,y, np.r_[0.])
|
||||
|
||||
times = np.arange(10)*1e-3 + 1e-3
|
||||
rx = SIP.Rx.Dipole(M, N, times)
|
||||
src = SIP.Src.Dipole([rx], Aloc, Bloc)
|
||||
survey = SIP.Survey([src])
|
||||
colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))]
|
||||
problem = SIP.Problem3D_CC(mesh, rho=1./sigma, mapping=colemap)
|
||||
problem.Solver = MumpsSolver
|
||||
problem.pair(survey)
|
||||
mSynth = np.r_[eta, 1./tau]
|
||||
survey.makeSyntheticData(mSynth)
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC*2)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-10
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
class IPProblemTestsN(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 25.
|
||||
hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hz = [(cs,0, -1.3),(cs,20)]
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
|
||||
blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
|
||||
blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
|
||||
sigma = np.ones(mesh.nC)*1e-2
|
||||
eta = np.zeros(mesh.nC)
|
||||
tau = np.ones_like(sigma)*1.
|
||||
eta[blkind0] = 0.1
|
||||
eta[blkind1] = 0.1
|
||||
tau[blkind0] = 0.1
|
||||
tau[blkind1] = 0.01
|
||||
|
||||
x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
|
||||
y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
|
||||
Aloc = np.r_[-200., 0., 0.]
|
||||
Bloc = np.r_[200., 0., 0.]
|
||||
M = Utils.ndgrid(x-25.,y, np.r_[0.])
|
||||
N = Utils.ndgrid(x+25.,y, np.r_[0.])
|
||||
|
||||
times = np.arange(10)*1e-3 + 1e-3
|
||||
rx = SIP.Rx.Dipole(M, N, times)
|
||||
src = SIP.Src.Dipole([rx], Aloc, Bloc)
|
||||
survey = SIP.Survey([src])
|
||||
colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))]
|
||||
problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap)
|
||||
problem.Solver = MumpsSolver
|
||||
problem.pair(survey)
|
||||
mSynth = np.r_[eta, 1./tau]
|
||||
survey.makeSyntheticData(mSynth)
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC*2)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-8
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
class IPProblemTestsN_air(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 25.
|
||||
hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hz = [(cs,0, -1.3),(cs,20),(cs,0, 1.3)]
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCC")
|
||||
blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
|
||||
blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
|
||||
sigma = np.ones(mesh.nC)*1e-2
|
||||
airind = mesh.gridCC[:,2]>0.
|
||||
sigma[airind] = 1e-8
|
||||
eta = np.zeros(mesh.nC)
|
||||
tau = np.ones_like(sigma)*1.
|
||||
eta[blkind0] = 0.1
|
||||
eta[blkind1] = 0.1
|
||||
tau[blkind0] = 0.1
|
||||
tau[blkind1] = 0.01
|
||||
|
||||
actmapeta = Maps.InjectActiveCells(mesh, ~airind, 0.)
|
||||
actmaptau = Maps.InjectActiveCells(mesh, ~airind, 1.)
|
||||
|
||||
x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
|
||||
y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
|
||||
Aloc = np.r_[-200., 0., 0.]
|
||||
Bloc = np.r_[200., 0., 0.]
|
||||
M = Utils.ndgrid(x-25.,y, np.r_[0.])
|
||||
N = Utils.ndgrid(x+25.,y, np.r_[0.])
|
||||
|
||||
times = np.arange(10)*1e-3 + 1e-3
|
||||
rx = SIP.Rx.Dipole(M, N, times)
|
||||
src = SIP.Src.Dipole([rx], Aloc, Bloc)
|
||||
survey = SIP.Survey([src])
|
||||
colemap = [("eta", Maps.IdentityMap(mesh)*actmapeta), ("taui", Maps.IdentityMap(mesh)*actmaptau)]
|
||||
problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap)
|
||||
problem.Solver = MumpsSolver
|
||||
problem.pair(survey)
|
||||
mSynth = np.r_[eta[~airind], 1./tau[~airind]]
|
||||
survey.makeSyntheticData(mSynth)
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
regmap = Maps.IdentityMap(nP=int(mSynth[~airind].size*2))
|
||||
reg = SIP.MultiRegularization(mesh, mapping=regmap, nModels=2, indActive=~airind)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-8
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,411 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import unittest
|
||||
import matplotlib.pyplot as plt
|
||||
from SimPEG import *
|
||||
|
||||
MESHTYPES = ['uniformTensorMesh']
|
||||
|
||||
def getxBCyBC_CC(mesh, alpha, beta, gamma):
|
||||
# def getxBCyBC(mesh, alpha, beta, gamma):
|
||||
"""
|
||||
This is a subfunction generating mixed-boundary condition:
|
||||
|
||||
.. math::
|
||||
|
||||
\nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q
|
||||
|
||||
\rho \vec{j} = -\nabla \phi \phi
|
||||
|
||||
\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 (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]
|
||||
h_xm, h_xp = mesh.hx[0], mesh.hx[-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)
|
||||
|
||||
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")
|
||||
|
||||
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]
|
||||
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]
|
||||
|
||||
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)
|
||||
|
||||
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 = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
|
||||
|
||||
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[4], beta[4], gamma[4]
|
||||
alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5]
|
||||
|
||||
# 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]
|
||||
|
||||
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)
|
||||
|
||||
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_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)
|
||||
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
|
||||
# Mrhoj = D.T V phi + P_BC*Utils.sdiag(y_BC)*M phi - P_BC*x_BC
|
||||
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_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)
|
||||
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[:,2]), np.ones_like(gBFzp[:,2])
|
||||
beta_zm, beta_zp = np.ones_like(gBFzm[:,2]), 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)
|
||||
|
||||
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_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)
|
||||
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()
|
||||
Reference in New Issue
Block a user