mirror of
https://github.com/wassname/simpeg.git
synced 2026-08-12 12:30:37 +08:00
Merge branch 'master' of https://github.com/simpeg/simpegem into em/dev
Conflicts: .coveragerc .gitignore .travis.yml docs/api_Utils.rst docs/conf.py docs/index.rst requirements.txt setup.py
This commit is contained in:
Rowan Cockett
@@ -0,0 +1,11 @@
|
||||
if __name__ == '__main__':
|
||||
import os
|
||||
import glob
|
||||
import unittest
|
||||
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,10 @@
|
||||
import unittest, os
|
||||
from SimPEG.EM import Examples
|
||||
|
||||
class EM_ExamplesRunning(unittest.TestCase):
|
||||
|
||||
def test_CylInversion(self):
|
||||
Examples.CylInversion.run(plotIt=False)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,478 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
|
||||
testDerivs = True
|
||||
testCrossCheck = True
|
||||
testAdjoint = True
|
||||
testEB = True
|
||||
testHJ = True
|
||||
|
||||
verbose = False
|
||||
|
||||
TOL = 1e-5
|
||||
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
|
||||
CONDUCTIVITY = 1e1
|
||||
MU = mu_0
|
||||
freq = 1e-1
|
||||
addrandoms = True
|
||||
|
||||
SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
|
||||
|
||||
|
||||
def getProblem(fdemType, comp):
|
||||
cs = 5.
|
||||
ncx, ncy, ncz = 6, 6, 6
|
||||
npad = 3
|
||||
hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
|
||||
mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
|
||||
|
||||
mapping = Maps.ExpMap(mesh)
|
||||
|
||||
x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0.])
|
||||
Rx0 = EM.FDEM.RxFDEM(XYZ, comp)
|
||||
|
||||
if SrcType is 'MagDipole':
|
||||
Src = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.])
|
||||
elif SrcType is 'MagDipole_Bfield':
|
||||
Src = EM.FDEM.SrcFDEM_MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.])
|
||||
elif SrcType is 'CircularLoop':
|
||||
Src2 = EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.])
|
||||
|
||||
if verbose:
|
||||
print ' Fetching %s problem' % (fdemType)
|
||||
|
||||
if fdemType == 'e':
|
||||
if SrcType is 'RawVec':
|
||||
S_m = np.zeros(mesh.nF)
|
||||
S_e = np.zeros(mesh.nE)
|
||||
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
|
||||
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
|
||||
Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Src])
|
||||
prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
|
||||
|
||||
elif fdemType == 'b':
|
||||
if SrcType is 'RawVec':
|
||||
S_m = np.zeros(mesh.nF)
|
||||
S_e = np.zeros(mesh.nE)
|
||||
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
|
||||
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
|
||||
Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Src])
|
||||
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
|
||||
|
||||
elif fdemType == 'j':
|
||||
if SrcType is 'RawVec':
|
||||
S_m = np.zeros(mesh.nE)
|
||||
S_e = np.zeros(mesh.nF)
|
||||
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
|
||||
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
|
||||
Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Src])
|
||||
prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping)
|
||||
|
||||
elif fdemType == 'h':
|
||||
if SrcType is 'RawVec':
|
||||
S_m = np.zeros(mesh.nE)
|
||||
S_e = np.zeros(mesh.nF)
|
||||
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
|
||||
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
|
||||
Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Src])
|
||||
prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
|
||||
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
prb.pair(survey)
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
pass
|
||||
|
||||
return prb
|
||||
|
||||
def adjointTest(fdemType, comp):
|
||||
prb = getProblem(fdemType, comp)
|
||||
print 'Adjoint %s formulation - %s' % (fdemType, comp)
|
||||
|
||||
m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
mu = np.ones(prb.mesh.nC)*MU
|
||||
|
||||
if addrandoms is True:
|
||||
m = m + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
|
||||
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
|
||||
|
||||
survey = prb.survey
|
||||
# prb.PropMap.PropModel.mu = mu
|
||||
# prb.PropMap.PropModel.mui = 1./mu
|
||||
u = prb.fields(m)
|
||||
|
||||
v = np.random.rand(survey.nD)
|
||||
w = np.random.rand(prb.mesh.nC)
|
||||
|
||||
vJw = v.dot(prb.Jvec(m, w, u))
|
||||
wJtv = w.dot(prb.Jtvec(m, v, u))
|
||||
tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
|
||||
print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
|
||||
return np.abs(vJw - wJtv) < tol
|
||||
|
||||
|
||||
def derivTest(fdemType, comp):
|
||||
|
||||
prb = getProblem(fdemType, comp)
|
||||
print '%s formulation - %s' % (fdemType, comp)
|
||||
x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
mu = np.log(np.ones(prb.mesh.nC)*MU)
|
||||
|
||||
if addrandoms is True:
|
||||
x0 = x0 + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
|
||||
mu = mu + np.random.randn(prb.mapping.nP)*MU*1e-1
|
||||
|
||||
# prb.PropMap.PropModel.mu = mu
|
||||
# prb.PropMap.PropModel.mui = 1./mu
|
||||
|
||||
survey = prb.survey
|
||||
def fun(x):
|
||||
return survey.dpred(x), lambda x: prb.Jvec(x0, x)
|
||||
return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
|
||||
|
||||
|
||||
def crossCheckTest(fdemType, comp):
|
||||
|
||||
l2norm = lambda r: np.sqrt(r.dot(r))
|
||||
|
||||
prb1 = getProblem(fdemType, comp)
|
||||
mesh = prb1.mesh
|
||||
print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp)
|
||||
m = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
|
||||
mu = np.log(np.ones(mesh.nC)*MU)
|
||||
|
||||
if addrandoms is True:
|
||||
m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
|
||||
mu = mu + np.random.randn(mesh.nC)*MU*1e-1
|
||||
|
||||
# prb1.PropMap.PropModel.mu = mu
|
||||
# prb1.PropMap.PropModel.mui = 1./mu
|
||||
survey1 = prb1.survey
|
||||
d1 = survey1.dpred(m)
|
||||
|
||||
if verbose:
|
||||
print ' Problem 1 solved'
|
||||
|
||||
if fdemType == 'e':
|
||||
prb2 = getProblem('b', comp)
|
||||
elif fdemType == 'b':
|
||||
prb2 = getProblem('e', comp)
|
||||
elif fdemType == 'j':
|
||||
prb2 = getProblem('h', comp)
|
||||
elif fdemType == 'h':
|
||||
prb2 = getProblem('j', comp)
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
# prb2.mu = mu
|
||||
survey2 = prb2.survey
|
||||
d2 = survey2.dpred(m)
|
||||
|
||||
if verbose:
|
||||
print ' Problem 2 solved'
|
||||
|
||||
r = d2-d1
|
||||
l2r = l2norm(r)
|
||||
|
||||
tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR])
|
||||
print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol
|
||||
return l2r < tol
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
class FDEM_DerivTests(unittest.TestCase):
|
||||
|
||||
if testDerivs:
|
||||
if testEB:
|
||||
def test_Jvec_exr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'exr'))
|
||||
def test_Jvec_eyr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'eyr'))
|
||||
def test_Jvec_ezr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'ezr'))
|
||||
def test_Jvec_exi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'exi'))
|
||||
def test_Jvec_eyi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'eyi'))
|
||||
def test_Jvec_ezi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'ezi'))
|
||||
|
||||
def test_Jvec_bxr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bxr'))
|
||||
def test_Jvec_byr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'byr'))
|
||||
def test_Jvec_bzr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bzr'))
|
||||
def test_Jvec_bxi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bxi'))
|
||||
def test_Jvec_byi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'byi'))
|
||||
def test_Jvec_bzi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bzi'))
|
||||
|
||||
def test_Jvec_exr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'exr'))
|
||||
def test_Jvec_eyr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'eyr'))
|
||||
def test_Jvec_ezr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'ezr'))
|
||||
def test_Jvec_exi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'exi'))
|
||||
def test_Jvec_eyi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'eyi'))
|
||||
def test_Jvec_ezi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'ezi'))
|
||||
|
||||
def test_Jvec_bxr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bxr'))
|
||||
def test_Jvec_byr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'byr'))
|
||||
def test_Jvec_bzr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bzr'))
|
||||
def test_Jvec_bxi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bxi'))
|
||||
def test_Jvec_byi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'byi'))
|
||||
def test_Jvec_bzi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bzi'))
|
||||
|
||||
if testHJ:
|
||||
def test_Jvec_jxr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jxr'))
|
||||
def test_Jvec_jyr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jyr'))
|
||||
def test_Jvec_jzr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jzr'))
|
||||
def test_Jvec_jxi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jxi'))
|
||||
def test_Jvec_jyi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jyi'))
|
||||
def test_Jvec_jzi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jzi'))
|
||||
|
||||
def test_Jvec_hxr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hxr'))
|
||||
def test_Jvec_hyr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hyr'))
|
||||
def test_Jvec_hzr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hzr'))
|
||||
def test_Jvec_hxi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hxi'))
|
||||
def test_Jvec_hyi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hyi'))
|
||||
def test_Jvec_hzi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hzi'))
|
||||
|
||||
def test_Jvec_hxr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hxr'))
|
||||
def test_Jvec_hyr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hyr'))
|
||||
def test_Jvec_hzr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hzr'))
|
||||
def test_Jvec_hxi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hxi'))
|
||||
def test_Jvec_hyi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hyi'))
|
||||
def test_Jvec_hzi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hzi'))
|
||||
|
||||
def test_Jvec_hxr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jxr'))
|
||||
def test_Jvec_hyr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jyr'))
|
||||
def test_Jvec_hzr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jzr'))
|
||||
def test_Jvec_hxi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jxi'))
|
||||
def test_Jvec_hyi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jyi'))
|
||||
def test_Jvec_hzi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jzi'))
|
||||
|
||||
|
||||
if testAdjoint:
|
||||
if testEB:
|
||||
def test_Jtvec_adjointTest_exr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'exr'))
|
||||
def test_Jtvec_adjointTest_eyr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'eyr'))
|
||||
def test_Jtvec_adjointTest_ezr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'ezr'))
|
||||
def test_Jtvec_adjointTest_exi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'exi'))
|
||||
def test_Jtvec_adjointTest_eyi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'eyi'))
|
||||
def test_Jtvec_adjointTest_ezi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'ezi'))
|
||||
|
||||
def test_Jtvec_adjointTest_bxr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bxr'))
|
||||
def test_Jtvec_adjointTest_byr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'byr'))
|
||||
def test_Jtvec_adjointTest_bzr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bzr'))
|
||||
def test_Jtvec_adjointTest_bxi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bxi'))
|
||||
def test_Jtvec_adjointTest_byi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'byi'))
|
||||
def test_Jtvec_adjointTest_bzi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_exr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'exr'))
|
||||
def test_Jtvec_adjointTest_eyr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'eyr'))
|
||||
def test_Jtvec_adjointTest_ezr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'ezr'))
|
||||
def test_Jtvec_adjointTest_exi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'exi'))
|
||||
def test_Jtvec_adjointTest_eyi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'eyi'))
|
||||
def test_Jtvec_adjointTest_ezi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'ezi'))
|
||||
def test_Jtvec_adjointTest_bxr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bxr'))
|
||||
def test_Jtvec_adjointTest_byr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'byr'))
|
||||
def test_Jtvec_adjointTest_bzr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bzr'))
|
||||
def test_Jtvec_adjointTest_bxi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bxi'))
|
||||
def test_Jtvec_adjointTest_byi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'byi'))
|
||||
def test_Jtvec_adjointTest_bzi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bzi'))
|
||||
|
||||
|
||||
if testHJ:
|
||||
def test_Jtvec_adjointTest_jxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxr'))
|
||||
def test_Jtvec_adjointTest_jyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyr'))
|
||||
def test_Jtvec_adjointTest_jzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzr'))
|
||||
def test_Jtvec_adjointTest_jxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxi'))
|
||||
def test_Jtvec_adjointTest_jyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyi'))
|
||||
def test_Jtvec_adjointTest_jzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzi'))
|
||||
|
||||
|
||||
if testCrossCheck:
|
||||
if testEB:
|
||||
def test_EB_CrossCheck_exr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exr'))
|
||||
def test_EB_CrossCheck_eyr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyr'))
|
||||
def test_EB_CrossCheck_ezr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezr'))
|
||||
def test_EB_CrossCheck_exi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exi'))
|
||||
def test_EB_CrossCheck_eyi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyi'))
|
||||
def test_EB_CrossCheck_ezi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezi'))
|
||||
|
||||
def test_EB_CrossCheck_bxr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxr'))
|
||||
def test_EB_CrossCheck_byr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byr'))
|
||||
def test_EB_CrossCheck_bzr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzr'))
|
||||
def test_EB_CrossCheck_bxi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxi'))
|
||||
def test_EB_CrossCheck_byi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byi'))
|
||||
def test_EB_CrossCheck_bzi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzi'))
|
||||
|
||||
if testHJ:
|
||||
def test_HJ_CrossCheck_jxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxr'))
|
||||
def test_HJ_CrossCheck_jyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyr'))
|
||||
def test_HJ_CrossCheck_jzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzr'))
|
||||
def test_HJ_CrossCheck_jxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxi'))
|
||||
def test_HJ_CrossCheck_jyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyi'))
|
||||
def test_HJ_CrossCheck_jzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzi'))
|
||||
|
||||
def test_HJ_CrossCheck_hxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxr'))
|
||||
def test_HJ_CrossCheck_hyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyr'))
|
||||
def test_HJ_CrossCheck_hzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzr'))
|
||||
def test_HJ_CrossCheck_hxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxi'))
|
||||
def test_HJ_CrossCheck_hyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyi'))
|
||||
def test_HJ_CrossCheck_hzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzi'))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,62 @@
|
||||
from SimPEG import Tests, Utils, np
|
||||
import SimPEG.EM.Analytics.FDEMcasing as Casing
|
||||
import unittest
|
||||
from scipy.constants import mu_0
|
||||
|
||||
|
||||
n = 50.
|
||||
freq = 1.
|
||||
a = 5e-2
|
||||
b = a + 1e-2
|
||||
sigma = np.r_[10., 5.5e6, 1e-1]
|
||||
mu = mu_0*np.r_[1.,100.,1.]
|
||||
srcloc = np.r_[0., 0., 0.]
|
||||
xobs = np.random.rand(n)+10.
|
||||
yobs = np.zeros(n)
|
||||
zobs = np.random.randn(n)
|
||||
plotit = False
|
||||
|
||||
def CasingMagDipoleDeriv_r(x):
|
||||
obsloc = np.vstack([x, yobs, zobs]).T
|
||||
|
||||
f = Casing._getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu)
|
||||
g = Utils.sdiag(Casing._getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu))
|
||||
|
||||
return f,g
|
||||
|
||||
def CasingMagDipoleDeriv_z(z):
|
||||
obsloc = np.vstack([xobs, yobs, z]).T
|
||||
|
||||
f = Casing._getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu)
|
||||
g = Utils.sdiag(Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu))
|
||||
|
||||
return f,g
|
||||
|
||||
def CasingMagDipole2Deriv_z_r(x):
|
||||
obsloc = np.vstack([x, yobs, zobs]).T
|
||||
|
||||
f = Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu)
|
||||
g = Utils.sdiag(Casing._getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu))
|
||||
|
||||
return f,g
|
||||
|
||||
def CasingMagDipole2Deriv_z_z(z):
|
||||
obsloc = np.vstack([xobs, yobs, z]).T
|
||||
|
||||
f = Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu)
|
||||
g = Utils.sdiag(Casing._getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu))
|
||||
|
||||
return f,g
|
||||
|
||||
|
||||
|
||||
class Casing_DerivTest(unittest.TestCase):
|
||||
def test_derivs(self):
|
||||
Tests.checkDerivative(CasingMagDipoleDeriv_r,np.ones(n)*10+np.random.randn(n),plotIt=False)
|
||||
Tests.checkDerivative(CasingMagDipoleDeriv_z,np.random.randn(n),plotIt=False)
|
||||
Tests.checkDerivative(CasingMagDipole2Deriv_z_r,np.ones(n)*10+np.random.randn(n),plotIt=False)
|
||||
Tests.checkDerivative(CasingMagDipole2Deriv_z_z,np.random.randn(n),plotIt=False)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,243 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
from scipy.constants import mu_0
|
||||
|
||||
plotIt = False
|
||||
tol_EBdipole = 1e-2
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pylab
|
||||
|
||||
|
||||
class FDEM_analyticTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 10.
|
||||
ncx, ncy, ncz = 10, 10, 10
|
||||
npad = 4
|
||||
freq = 1e2
|
||||
|
||||
hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
|
||||
mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC')
|
||||
|
||||
mapping = Maps.ExpMap(mesh)
|
||||
|
||||
x = np.linspace(-10,10,5)
|
||||
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
|
||||
rxList = EM.FDEM.RxFDEM(XYZ, 'exi')
|
||||
Src0 = EM.FDEM.SrcFDEM_MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Src0])
|
||||
|
||||
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
|
||||
prb.pair(survey)
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
prb.Solver = SolverLU
|
||||
|
||||
sig = 1e-1
|
||||
sigma = np.ones(mesh.nC)*sig
|
||||
sigma[mesh.gridCC[:,2] > 0] = 1e-8
|
||||
m = np.log(sigma)
|
||||
|
||||
self.prb = prb
|
||||
self.mesh = mesh
|
||||
self.m = m
|
||||
self.Src0 = Src0
|
||||
self.sig = sig
|
||||
|
||||
def test_Transect(self):
|
||||
print 'Testing Transect for analytic'
|
||||
|
||||
u = self.prb.fields(self.m)
|
||||
|
||||
bfz = self.mesh.r(u[self.Src0, 'b'],'F','Fz','M')
|
||||
x = np.linspace(-55,55,12)
|
||||
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
|
||||
|
||||
P = self.mesh.getInterpolationMat(XYZ, 'Fz')
|
||||
|
||||
an = EM.Analytics.FDEM.hzAnalyticDipoleF(x, self.Src0.freq, self.sig)
|
||||
|
||||
diff = np.log10(np.abs(P*np.imag(u[self.Src0, 'b']) - mu_0*np.imag(an)))
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
plt.plot(x,np.log10(np.abs(P*np.imag(u[self.Src0, 'b']))))
|
||||
plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r')
|
||||
plt.plot(x,diff,'g')
|
||||
plt.show()
|
||||
|
||||
# We want the difference to be an orderMag less
|
||||
# than the analytic solution. Note that right at
|
||||
# the source, both the analytic and the numerical
|
||||
# solution will be poor. Use plotIt up top to see that...
|
||||
orderMag = 1.6
|
||||
passed = np.abs(np.mean(diff - np.log10(np.abs(mu_0*np.imag(an))))) > orderMag
|
||||
self.assertTrue(passed)
|
||||
|
||||
|
||||
def test_CylMeshEBDipoles(self):
|
||||
print 'Testing CylMesh Electric and Magnetic Dipoles in a wholespace- Analytic: J-formulation'
|
||||
sigmaback = 1.
|
||||
mur = 2.
|
||||
freq = 1.
|
||||
skdpth = 500./np.sqrt(sigmaback*freq)
|
||||
|
||||
csx, ncx, npadx = 5, 50, 25
|
||||
csz, ncz, npadz = 5, 50, 25
|
||||
hx = Utils.meshTensor([(csx,ncx), (csx,npadx,1.3)])
|
||||
hz = Utils.meshTensor([(csz,npadz,-1.3), (csz,ncz), (csz,npadz,1.3)])
|
||||
mesh = Mesh.CylMesh([hx,1,hz], [0.,0.,-hz.sum()/2]) # define the cylindrical mesh
|
||||
|
||||
if plotIt:
|
||||
mesh.plotGrid()
|
||||
|
||||
# make sure mesh is big enough
|
||||
self.assertTrue(mesh.hz.sum() > skdpth*2.)
|
||||
self.assertTrue(mesh.hx.sum() > skdpth*2.)
|
||||
|
||||
SigmaBack = sigmaback*np.ones((mesh.nC))
|
||||
MuBack = mur*mu_0*np.ones((mesh.nC))
|
||||
|
||||
# set up source
|
||||
# test electric dipole
|
||||
src_loc = np.r_[0.,0.,0.]
|
||||
s_ind = Utils.closestPoints(mesh,src_loc,'Fz') + mesh.nFx
|
||||
|
||||
de = np.zeros(mesh.nF,dtype=complex)
|
||||
de[s_ind] = 1./csz
|
||||
de_p = [EM.FDEM.SrcFDEM_RawVec_e([],freq,de/mesh.area)]
|
||||
|
||||
dm_p = [EM.FDEM.SrcFDEM_MagDipole([],freq,src_loc)]
|
||||
|
||||
|
||||
# Pair the problem and survey
|
||||
surveye = EM.FDEM.SurveyFDEM(de_p)
|
||||
surveym = EM.FDEM.SurveyFDEM(dm_p)
|
||||
|
||||
mapping = [('sigma', Maps.IdentityMap(mesh)),('mu', Maps.IdentityMap(mesh))]
|
||||
|
||||
prbe = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
|
||||
prbm = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
|
||||
|
||||
prbe.pair(surveye) # pair problem and survey
|
||||
prbm.pair(surveym)
|
||||
|
||||
# solve
|
||||
fieldsBackE = prbe.fields(np.r_[SigmaBack, MuBack]) # Done
|
||||
fieldsBackM = prbm.fields(np.r_[SigmaBack, MuBack]) # Done
|
||||
|
||||
|
||||
rlim = [20.,500.]
|
||||
lookAtTx = de_p
|
||||
r = mesh.vectorCCx[np.argmin(np.abs(mesh.vectorCCx-rlim[0])):np.argmin(np.abs(mesh.vectorCCx-rlim[1]))]
|
||||
z = 100.
|
||||
|
||||
# where we choose to measure
|
||||
XYZ = Utils.ndgrid(r, np.r_[0.], np.r_[z])
|
||||
|
||||
Pf = mesh.getInterpolationMat(XYZ, 'CC')
|
||||
Zero = sp.csr_matrix(Pf.shape)
|
||||
Pfx,Pfz = sp.hstack([Pf,Zero]),sp.hstack([Zero,Pf])
|
||||
|
||||
jn = fieldsBackE[de_p,'j']
|
||||
bn = fieldsBackM[dm_p,'b']
|
||||
|
||||
Rho = Utils.sdiag(1./SigmaBack)
|
||||
Rho = sp.block_diag([Rho,Rho])
|
||||
|
||||
en = Rho*mesh.aveF2CCV*jn
|
||||
bn = mesh.aveF2CCV*bn
|
||||
|
||||
ex,ez = Pfx*en, Pfz*en
|
||||
bx,bz = Pfx*bn, Pfz*bn
|
||||
|
||||
# get analytic solution
|
||||
exa, eya, eza = EM.Analytics.FDEM.ElectricDipoleWholeSpace(XYZ, src_loc, sigmaback, freq,orientation='Z',mu= mur*mu_0)
|
||||
exa, eya, eza = Utils.mkvc(exa,2), Utils.mkvc(eya,2), Utils.mkvc(eza,2)
|
||||
|
||||
bxa, bya, bza = EM.Analytics.FDEM.MagneticDipoleWholeSpace(XYZ, src_loc, sigmaback, freq,orientation='Z',mu= mur*mu_0)
|
||||
bxa, bya, bza = Utils.mkvc(bxa,2), Utils.mkvc(bya,2), Utils.mkvc(bza,2)
|
||||
|
||||
print ' comp, anayltic, numeric, num - ana, (num - ana)/ana'
|
||||
print ' ex:', np.linalg.norm(exa), np.linalg.norm(ex), np.linalg.norm(exa-ex), np.linalg.norm(exa-ex)/np.linalg.norm(exa)
|
||||
print ' ez:', np.linalg.norm(eza), np.linalg.norm(ez), np.linalg.norm(eza-ez), np.linalg.norm(eza-ez)/np.linalg.norm(eza)
|
||||
|
||||
print ' bx:', np.linalg.norm(bxa), np.linalg.norm(bx), np.linalg.norm(bxa-bx), np.linalg.norm(bxa-bx)/np.linalg.norm(bxa)
|
||||
print ' bz:', np.linalg.norm(bza), np.linalg.norm(bz), np.linalg.norm(bza-bz), np.linalg.norm(bza-bz)/np.linalg.norm(bza)
|
||||
|
||||
if plotIt:
|
||||
# Edipole
|
||||
plt.subplot(221)
|
||||
plt.plot(r,ex.real,'o',r,exa.real,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Ex Real')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(222)
|
||||
plt.plot(r,ex.imag,'o',r,exa.imag,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Ex Imag')
|
||||
plt.legend(['Num','Ana'],bbox_to_anchor=(1.5,0.5))
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(223)
|
||||
plt.plot(r,ez.real,'o',r,eza.real,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Ez Real')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(224)
|
||||
plt.plot(r,ez.imag,'o',r,eza.imag,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Ez Imag')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.tight_layout()
|
||||
|
||||
# Bdipole
|
||||
plt.subplot(221)
|
||||
plt.plot(r,bx.real,'o',r,bxa.real,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Bx Real')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(222)
|
||||
plt.plot(r,bx.imag,'o',r,bxa.imag,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Bx Imag')
|
||||
plt.legend(['Num','Ana'],bbox_to_anchor=(1.5,0.5))
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(223)
|
||||
plt.plot(r,bz.real,'o',r,bza.real,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Bz Real')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.subplot(224)
|
||||
plt.plot(r,bz.imag,'o',r,bza.imag,linewidth=2)
|
||||
plt.grid(which='both')
|
||||
plt.title('Bz Imag')
|
||||
plt.xlabel('r (m)')
|
||||
|
||||
plt.tight_layout()
|
||||
|
||||
self.assertTrue(np.linalg.norm(exa-ex)/np.linalg.norm(exa) < tol_EBdipole)
|
||||
self.assertTrue(np.linalg.norm(eza-ez)/np.linalg.norm(eza) < tol_EBdipole)
|
||||
|
||||
self.assertTrue(np.linalg.norm(bxa-bx)/np.linalg.norm(bxa) < tol_EBdipole)
|
||||
self.assertTrue(np.linalg.norm(bza-bz)/np.linalg.norm(bza) < tol_EBdipole)
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,314 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
|
||||
plotIt = False
|
||||
tol = 1e-6
|
||||
|
||||
class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
ncy = 6
|
||||
npad = 20
|
||||
hx = [(cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
mesh = Mesh.CylMesh([hx,1,hy], '00C')
|
||||
|
||||
active = mesh.vectorCCz<0.
|
||||
activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
|
||||
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap
|
||||
|
||||
rxOffset = 40.
|
||||
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
|
||||
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM([src])
|
||||
|
||||
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
# self.prb.timeSteps = [1e-5]
|
||||
self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
|
||||
# self.prb.timeSteps = [(1e-05, 100)]
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
self.prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
self.prb.Solver = SolverLU
|
||||
|
||||
self.sigma = np.ones(mesh.nCz)*1e-8
|
||||
self.sigma[mesh.vectorCCz<0] = 1e-1
|
||||
self.sigma = np.log(self.sigma[active])
|
||||
|
||||
self.prb.pair(survey)
|
||||
self.mesh = mesh
|
||||
|
||||
def test_AhVec(self):
|
||||
"""
|
||||
Test that fields and AhVec produce consistent results
|
||||
"""
|
||||
|
||||
prb = self.prb
|
||||
sigma = self.sigma
|
||||
|
||||
u = prb.fields(sigma)
|
||||
Ahu = prb._AhVec(sigma, u)
|
||||
|
||||
V1 = Ahu[:,'b',1]
|
||||
V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0]
|
||||
self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6)
|
||||
|
||||
V1 = Ahu[:,'e',1]
|
||||
return np.linalg.norm(V1) < 1.e-6
|
||||
|
||||
for i in range(2,prb.nT):
|
||||
|
||||
dt = prb.timeSteps[i]
|
||||
|
||||
V1 = Ahu[:,'b',i]
|
||||
V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1]
|
||||
# print np.linalg.norm(V1), np.linalg.norm(V2)
|
||||
self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
|
||||
|
||||
V1 = Ahu[:,'e',i]
|
||||
V2 = prb.MeSigma*u[:,'e',i]
|
||||
# print np.linalg.norm(V1), np.linalg.norm(V2)
|
||||
return np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6
|
||||
|
||||
def test_AhVecVSMat_OneTS(self):
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [1e-05]
|
||||
sigma = self.sigma
|
||||
prb.curModel = sigma
|
||||
|
||||
dt = prb.timeSteps[0]
|
||||
a11 = 1/dt*prb.MfMui*sp.identity(prb.mesh.nF)
|
||||
a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
a22 = -prb.MeSigma
|
||||
A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
|
||||
f = prb.fields(sigma)
|
||||
u1 = A*f.tovec()
|
||||
u2 = prb._AhVec(sigma,f).tovec()
|
||||
|
||||
self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
|
||||
|
||||
def test_solveAhVSMat_OneTS(self):
|
||||
prb = self.prb
|
||||
|
||||
prb.timeSteps = [1e-05]
|
||||
|
||||
sigma = self.sigma
|
||||
prb.curModel = sigma
|
||||
|
||||
dt = prb.timeSteps[0]
|
||||
a11 = 1.0/dt*prb.MfMui*sp.identity(prb.mesh.nF)
|
||||
a12 = prb.MfMui*prb.mesh.edgeCurl
|
||||
a21 = prb.mesh.edgeCurl.T*prb.MfMui
|
||||
a22 = -prb.MeSigma
|
||||
A = sp.bmat([[a11,a12],[a21,a22]])
|
||||
|
||||
f = prb.fields(sigma)
|
||||
f[:,:,0] = {'b':0}
|
||||
f[:,'b',1] = 0
|
||||
|
||||
self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec()))
|
||||
|
||||
u1 = prb.solveAh(sigma,f).tovec().flatten()
|
||||
u2 = sp.linalg.spsolve(A.tocsr(),f.tovec())
|
||||
|
||||
self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
|
||||
|
||||
def test_solveAhVsAhVec(self):
|
||||
|
||||
prb = self.prb
|
||||
mesh = self.prb.mesh
|
||||
sigma = self.sigma
|
||||
self.prb.curModel = sigma
|
||||
|
||||
f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
f[:,'b',:] = 0.0
|
||||
for i in range(prb.nT):
|
||||
f[:,'e', i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
Ahf = prb._AhVec(sigma, f)
|
||||
f_test = prb.solveAh(sigma, Ahf)
|
||||
|
||||
u1 = f.tovec()
|
||||
u2 = f_test.tovec()
|
||||
self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
|
||||
|
||||
def test_DerivG(self):
|
||||
"""
|
||||
Test the derivative of c with respect to sigma
|
||||
"""
|
||||
|
||||
# Random model and perturbation
|
||||
sigma = np.random.rand(self.prb.mapping.nP)
|
||||
|
||||
f = self.prb.fields(sigma)
|
||||
dm = 1000*np.random.rand(self.prb.mapping.nP)
|
||||
h = 0.01
|
||||
|
||||
derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
|
||||
print '\ntest_DerivG'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
return passed
|
||||
|
||||
def test_Deriv_dUdM(self):
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
|
||||
dm = 10*np.random.rand(prb.mapping.nP)
|
||||
f = prb.fields(sigma)
|
||||
|
||||
derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
|
||||
print '\n'
|
||||
print 'test_Deriv_dUdM'
|
||||
Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
|
||||
def test_Deriv_J(self):
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
|
||||
# d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
|
||||
d_sig = 10*np.random.rand(prb.mapping.nP)
|
||||
|
||||
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
|
||||
print '\n'
|
||||
print 'test_Deriv_J'
|
||||
Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
|
||||
|
||||
def test_projectAdjoint(self):
|
||||
prb = self.prb
|
||||
survey = prb.survey
|
||||
mesh = self.mesh
|
||||
|
||||
# Generate random fields and data
|
||||
f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(prb.nT):
|
||||
f[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
f[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
d_vec = np.random.rand(survey.nD)
|
||||
d = Survey.Data(survey,v=d_vec)
|
||||
|
||||
# Check that d.T*Q*f = f.T*Q.T*d
|
||||
V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec())
|
||||
V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec())
|
||||
|
||||
self.assertTrue((V1-V2)/np.abs(V1) < tol)
|
||||
|
||||
def test_adjointAhVsAht(self):
|
||||
prb = self.prb
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
|
||||
f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
f1[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
f1[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
f2[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
f2[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec())
|
||||
V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec())
|
||||
self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
|
||||
|
||||
# def test_solveAhtVsAhtVec(self):
|
||||
# prb = self.prb
|
||||
# mesh = self.mesh
|
||||
# sigma = np.random.rand(prb.mapping.nP)
|
||||
|
||||
# f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey)
|
||||
# for i in range(1,prb.nT+1):
|
||||
# f1[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f1[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# f2 = prb.solveAht(sigma, f1)
|
||||
# f3 = prb._AhtVec(sigma, f2)
|
||||
|
||||
# if True:
|
||||
# import matplotlib.pyplot as plt
|
||||
# plt.plot(f3.tovec(),'b')
|
||||
# plt.plot(f1.tovec(),'r')
|
||||
# plt.show()
|
||||
# V1 = np.linalg.norm(f3.tovec()-f1.tovec())
|
||||
# V2 = np.linalg.norm(f1.tovec())
|
||||
# print 'AhtVsAhtVec', V1, V2, f1.tovec()
|
||||
# print 'I am gunna fail this one: boo. :('
|
||||
# self.assertLess(V1/V2, 1e-6)
|
||||
|
||||
# def test_adjointsolveAhVssolveAht(self):
|
||||
# prb = self.prb
|
||||
# mesh = self.mesh
|
||||
# sigma = self.sigma
|
||||
|
||||
# f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(1,prb.nT+1):
|
||||
# f1[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f1[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
# for i in range(1,prb.nT+1):
|
||||
# f2[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
# f2[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
# V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
|
||||
# V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
|
||||
# print V1, V2
|
||||
# self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
|
||||
def test_adjointGvecVsGtvec(self):
|
||||
mesh = self.mesh
|
||||
prb = self.prb
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
sigma = np.random.rand(prb.mapping.nP)
|
||||
|
||||
u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
u[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
u[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
v[:,'b',i] = np.random.rand(mesh.nF, 1)
|
||||
v[:,'e',i] = np.random.rand(mesh.nE, 1)
|
||||
|
||||
V1 = m.dot(prb.Gtvec(sigma, v, u))
|
||||
V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
|
||||
self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
|
||||
|
||||
def test_adjointJvecVsJtvec(self):
|
||||
mesh = self.mesh
|
||||
prb = self.prb
|
||||
sigma = self.sigma
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
d = np.random.rand(prb.survey.nD)
|
||||
|
||||
V1 = d.dot(prb.Jvec(sigma, m))
|
||||
V2 = m.dot(prb.Jtvec(sigma, d))
|
||||
passed = np.abs(V1-V2)/np.abs(V1) < tol
|
||||
print 'AdjointTest', V1, V2, passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,153 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
|
||||
plotIt = False
|
||||
|
||||
class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
ncy = 6
|
||||
npad = 20
|
||||
hx = [(cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
mesh = Mesh.CylMesh([hx,1,hy], '00C')
|
||||
|
||||
active = mesh.vectorCCz<0.
|
||||
activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
|
||||
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap
|
||||
|
||||
rxOffset = 40.
|
||||
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
|
||||
src = EM.TDEM.SrcTDEM_VMD_MVP( [rx], loc=np.array([0., 0., 0.]))
|
||||
rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz')
|
||||
src2 = EM.TDEM.SrcTDEM_VMD_MVP( [rx2], loc=np.array([0., 0., 0.]))
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM([src,src2])
|
||||
|
||||
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
# self.prb.timeSteps = [1e-5]
|
||||
self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
|
||||
# self.prb.timeSteps = [(1e-05, 100)]
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
self.prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
self.prb.Solver = SolverLU
|
||||
|
||||
self.sigma = np.ones(mesh.nCz)*1e-8
|
||||
self.sigma[mesh.vectorCCz<0] = 1e-1
|
||||
self.sigma = np.log(self.sigma[active])
|
||||
|
||||
self.prb.pair(survey)
|
||||
self.mesh = mesh
|
||||
|
||||
def test_DerivG(self):
|
||||
"""
|
||||
Test the derivative of c with respect to sigma
|
||||
"""
|
||||
|
||||
# Random model and perturbation
|
||||
sigma = np.random.rand(self.prb.mapping.nP)
|
||||
|
||||
f = self.prb.fields(sigma)
|
||||
dm = 1000*np.random.rand(self.prb.mapping.nP)
|
||||
h = 0.01
|
||||
|
||||
derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
|
||||
print '\ntest_DerivG'
|
||||
Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
|
||||
def test_Deriv_dUdM(self):
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
|
||||
dm = 10*np.random.rand(prb.mapping.nP)
|
||||
f = prb.fields(sigma)
|
||||
|
||||
derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
|
||||
print '\n'
|
||||
print 'test_Deriv_dUdM'
|
||||
Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
|
||||
def test_Deriv_J(self):
|
||||
|
||||
prb = self.prb
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
mesh = self.mesh
|
||||
sigma = self.sigma
|
||||
|
||||
# d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
|
||||
d_sig = 10*np.random.rand(prb.mapping.nP)
|
||||
|
||||
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
|
||||
print '\n'
|
||||
print 'test_Deriv_J'
|
||||
Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
|
||||
|
||||
def test_projectAdjoint(self):
|
||||
prb = self.prb
|
||||
survey = prb.survey
|
||||
nSrc = survey.nSrc
|
||||
mesh = self.mesh
|
||||
|
||||
# Generate random fields and data
|
||||
f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(prb.nT):
|
||||
f[:,'b',i] = np.random.rand(mesh.nF, nSrc)
|
||||
f[:,'e',i] = np.random.rand(mesh.nE, nSrc)
|
||||
d_vec = np.random.rand(survey.nD)
|
||||
d = Survey.Data(survey,v=d_vec)
|
||||
|
||||
# Check that d.T*Q*f = f.T*Q.T*d
|
||||
V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec())
|
||||
V2 = np.sum((f.tovec())*(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec()))
|
||||
|
||||
self.assertTrue((V1-V2)/np.abs(V1) < 1e-6)
|
||||
|
||||
def test_adjointGvecVsGtvec(self):
|
||||
mesh = self.mesh
|
||||
prb = self.prb
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
sigma = np.random.rand(prb.mapping.nP)
|
||||
|
||||
u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
u[:,'b',i] = np.random.rand(mesh.nF, 2)
|
||||
u[:,'e',i] = np.random.rand(mesh.nE, 2)
|
||||
|
||||
v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(1,prb.nT+1):
|
||||
v[:,'b',i] = np.random.rand(mesh.nF, 2)
|
||||
v[:,'e',i] = np.random.rand(mesh.nE, 2)
|
||||
|
||||
V1 = m.dot(prb.Gtvec(sigma, v, u))
|
||||
V2 = np.sum(v.tovec()*prb.Gvec(sigma, m, u).tovec())
|
||||
self.assertTrue(np.abs(V1-V2)/np.abs(V1) <1e-6)
|
||||
|
||||
def test_adjointJvecVsJtvec(self):
|
||||
mesh = self.mesh
|
||||
prb = self.prb
|
||||
sigma = self.sigma
|
||||
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
d = np.random.rand(prb.survey.nD)
|
||||
|
||||
V1 = d.dot(prb.Jvec(sigma, m))
|
||||
V2 = m.dot(prb.Jtvec(sigma, d))
|
||||
print 'AdjointTest', V1, V2
|
||||
self.assertTrue(np.abs(V1-V2)/np.abs(V1) < 1e-6)
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,94 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
|
||||
plotIt = False
|
||||
|
||||
def getProb(meshType='CYL',rxTypes='bx,bz',nSrc=1):
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
ncy = 6
|
||||
npad = 20
|
||||
hx = [(cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
mesh = Mesh.CylMesh([hx,1,hy], '00C')
|
||||
|
||||
active = mesh.vectorCCz<0.
|
||||
activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
|
||||
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap
|
||||
|
||||
rxOffset = 40.
|
||||
|
||||
srcs = []
|
||||
for ii in range(nSrc):
|
||||
rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')]
|
||||
srcs += [EM.TDEM.SrcTDEM_VMD_MVP(rxs,np.array([0., 0., 0.]))]
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM(srcs)
|
||||
|
||||
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
# prb.timeSteps = [1e-5]
|
||||
prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
|
||||
# prb.timeSteps = [(1e-05, 100)]
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
prb.Solver = SolverLU
|
||||
|
||||
sigma = np.ones(mesh.nCz)*1e-8
|
||||
sigma[mesh.vectorCCz<0] = 1e-1
|
||||
sigma = np.log(sigma[active])
|
||||
|
||||
prb.pair(survey)
|
||||
return prb, mesh, sigma
|
||||
|
||||
def dotestJvec(prb, mesh, sigma):
|
||||
prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
|
||||
# d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
|
||||
d_sig = 10*np.random.rand(prb.mapping.nP)
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
|
||||
return Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20)
|
||||
|
||||
def dotestAdjoint(prb, mesh, sigma):
|
||||
m = np.random.rand(prb.mapping.nP)
|
||||
d = np.random.rand(prb.survey.nD)
|
||||
|
||||
V1 = d.dot(prb.Jvec(sigma, m))
|
||||
V2 = m.dot(prb.Jtvec(sigma, d))
|
||||
print 'AdjointTest', V1, V2
|
||||
return np.abs(V1-V2)/np.abs(V1), 1e-6
|
||||
|
||||
class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
def test_Jvec_bx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx')))
|
||||
def test_Adjoint_bx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx')))
|
||||
|
||||
def test_Jvec_bxbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz')))
|
||||
def test_Adjoint_bxbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz')))
|
||||
|
||||
def test_Jvec_bxbz_2src(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nSrc=2)))
|
||||
def test_Adjoint_bxbz_2src(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nSrc=2)))
|
||||
|
||||
def test_Jvec_bxbzbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz,bz')))
|
||||
def test_Adjoint_bxbzbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz,bz')))
|
||||
|
||||
def test_Jvec_dbxdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbxdt')))
|
||||
def test_Adjoint_dbxdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbxdt')))
|
||||
|
||||
def test_Jvec_dbzdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbzdt')))
|
||||
def test_Adjoint_dbzdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbzdt')))
|
||||
|
||||
def test_Jvec_dbxdtbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbxdt,bz')))
|
||||
def test_Adjoint_dbxdtbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbxdt,bz')))
|
||||
|
||||
def test_Jvec_ey(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey')))
|
||||
def test_Adjoint_ey(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey')))
|
||||
|
||||
def test_Jvec_eybzdbxdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey,bz,dbxdt')))
|
||||
def test_Adjoint_eybzdbxdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey,bz,dbxdt')))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,89 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
from scipy.constants import mu_0
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
except ImportError, e:
|
||||
MumpsSolver = SolverLU
|
||||
|
||||
|
||||
def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,1e-3], showIt=False):
|
||||
if meshType == 'CYL':
|
||||
cs, ncx, ncz, npad = 5., 30, 10, 15
|
||||
hx = [(cs,ncx), (cs,npad,1.3)]
|
||||
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
|
||||
mesh = Mesh.CylMesh([hx,1,hz], '00C')
|
||||
elif meshType == 'TENSOR':
|
||||
cs, nc, npad = 20., 13, 5
|
||||
hx = [(cs,npad,-1.3), (cs,nc), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,nc), (cs,npad,1.3)]
|
||||
hz = [(cs,npad,-1.3), (cs,nc), (cs,npad,1.3)]
|
||||
mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC')
|
||||
|
||||
active = mesh.vectorCCz<0.
|
||||
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
|
||||
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
|
||||
|
||||
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
|
||||
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
|
||||
# src = EM.TDEM.SrcTDEM([rx], loc=np.array([0., 0., 0.]))
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM([src])
|
||||
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
prb.Solver = MumpsSolver
|
||||
|
||||
prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40), (0.0001, 40), (0.0005, 40)]
|
||||
|
||||
sigma = np.ones(mesh.nCz)*1e-8
|
||||
sigma[active] = sig_half
|
||||
sigma = np.log(sigma[active])
|
||||
prb.pair(survey)
|
||||
|
||||
bz_ana = mu_0*EM.Analytics.hzAnalyticDipoleT(rx.locs[0][0]+1e-3, rx.times, sig_half)
|
||||
|
||||
bz_calc = survey.dpred(sigma)
|
||||
|
||||
ind = np.logical_and(rx.times > bounds[0],rx.times < bounds[1])
|
||||
log10diff = np.linalg.norm(np.log10(np.abs(bz_calc[ind])) - np.log10(np.abs(bz_ana[ind])))/np.linalg.norm(np.log10(np.abs(bz_ana[ind])))
|
||||
print 'Difference: ', log10diff
|
||||
|
||||
if showIt == True:
|
||||
plt.loglog(rx.times[bz_calc>0], bz_calc[bz_calc>0], 'r', rx.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--')
|
||||
plt.loglog(rx.times, abs(bz_ana), 'b*')
|
||||
plt.title('sig_half = %e'%sig_half)
|
||||
plt.show()
|
||||
|
||||
return log10diff
|
||||
|
||||
|
||||
class TDEM_bTests(unittest.TestCase):
|
||||
|
||||
def test_analytic_p2_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+2) < 0.01)
|
||||
def test_analytic_p1_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+1) < 0.01)
|
||||
def test_analytic_p0_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+0) < 0.01)
|
||||
def test_analytic_m1_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-1) < 0.01)
|
||||
def test_analytic_m2_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-2) < 0.01)
|
||||
def test_analytic_m3_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-3) < 0.02)
|
||||
|
||||
def test_analytic_p0_CYL_1m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e+0) < 0.01)
|
||||
def test_analytic_m1_CYL_1m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-1) < 0.01)
|
||||
def test_analytic_m2_CYL_1m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-2) < 0.01)
|
||||
def test_analytic_m3_CYL_1m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-3) < 0.02)
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
Reference in New Issue
Block a user