wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-13 17:45:30 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5 from simpeg/MapUpdates

Browse Source 
Map updates (there are still serious(?) problems with TDEM), but not with these updates.
This commit is contained in:
Rowan Cockett
2014-04-15 16:43:58 -07:00
parent 3f7f95de8c 2d90799415
commit f4e210a61d
10 changed files with 153 additions and 68 deletions
Download Patch File Download Diff File Expand all files Collapse all files
README.md
+1 -1
View File
@@ -15,7 +15,7 @@ Tests:
[https://travis-ci.org/simpeg/simpegem](https://travis-ci.org/simpeg/simpegem)
Build Status:
[![Build Status](https://travis-ci.org/simpeg/simpegem.png)](https://travis-ci.org/simpeg/simpegem)
[![Build Status](https://travis-ci.org/simpeg/simpegem.svg?branch=master)](https://travis-ci.org/simpeg/simpegem)
Bugs & Issues:
[https://github.com/simpeg/simpegem/issues](https://github.com/simpeg/simpegem/issues)
docs/index.rst
+1 -1
View File
@@ -52,7 +52,7 @@ Testing simpegEM
================
* Master Branch
.. image:: https://travis-ci.org/simpeg/simpegem.png?branch=master
.. image:: https://travis-ci.org/simpeg/simpegem.svg?branch=master
:target: https://travis-ci.org/simpeg/simpegem
:alt: Master Branch
:align: center
simpegEM/FDEM/FDEM.py
+3 -3
View File
@@ -66,13 +66,13 @@ class BaseProblemFDEM(Problem.BaseProblem):
@property
def curTModel(self):
if getattr(self, '_curTModel', None) is None:
self._curTModel = self.model.transform(self.curModel)
self._curTModel = self.mapping.transform(self.curModel)
return self._curTModel
@property
def curTModelDeriv(self):
if getattr(self, '_curTModelDeriv', None) is None:
self._curTModelDeriv = self.model.transformDeriv(self.curModel)
self._curTModelDeriv = self.mapping.transformDeriv(self.curModel)
return self._curTModelDeriv
def fields(self, m):
@@ -132,7 +132,7 @@ class BaseProblemFDEM(Problem.BaseProblem):
if not isinstance(v, self.dataPair):
v = self.dataPair(self.survey, v)
Jtv = np.zeros(self.model.nP)
Jtv = np.zeros(self.mapping.nP)
for freq in self.survey.freqs:
AT = self.getA(freq).T
simpegEM/TDEM/BaseTDEM.py
+17 -10
View File
@@ -22,13 +22,19 @@ class MixinInitialFieldCalc(object):
return F
def _getInitialFields_VMD_MVP(self):
if self.mesh._meshType is 'CYL1D':
MVP = Sources.MagneticDipoleVectorPotential(np.r_[0,0,self.survey.txLoc], np.c_[np.zeros(self.mesh.nN), self.mesh.gridN], 'x')
if self.mesh._meshType is 'CYL':
if self.mesh.isSymmetric:
MVP = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEy, 'y')
# MVP = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, np.c_[np.zeros(self.mesh.nN), self.mesh.gridN], 'x')
else:
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
elif self.mesh._meshType is 'TENSOR':
MVPx = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEx, 'x')
MVPy = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEy, 'y')
MVPz = Sources.MagneticDipoleVectorPotential(self.survey.txLoc, self.mesh.gridEz, 'z')
MVP = np.concatenate((MVPx, MVPy, MVPz))
else:
raise Exception('Unknown mesh for VMD')
# Initialize field object
F = FieldsTDEM(self.mesh, 1, self.times.size, store=self.storeTheseFields)
@@ -60,7 +66,7 @@ class MixinTimeStuff(object):
nsteps = property(**nsteps())
def times():
doc = "Modelling times"
doc = "Modeling times"
def fget(self):
t = np.r_[1:self.nsteps[0]+1]*self.dt[0]
for i in range(1,self.dt.size):
@@ -86,8 +92,8 @@ class MixinTimeStuff(object):
class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
"""docstring for ProblemTDEM1D"""
def __init__(self, model, **kwargs):
BaseProblem.__init__(self, model, **kwargs)
def __init__(self, mesh, mapping=None, **kwargs):
BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs)
####################################################
@@ -116,10 +122,10 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
def MeSigmaI(self): return self._MeSigmaI
def makeMassMatrices(self, m):
m = self.model.transform(m)
self._MeSigma = self.mesh.getMass(m, loc='e')
sig = self.mapping.transform(m)
self._MeSigma = self.mesh.getEdgeInnerProduct(sig)
self._MeSigmaI = sdiag(1/self.MeSigma.diagonal())
self._MfMui = self.mesh.getMass(1/mu_0, loc='f')
self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0)
def calcFields(self, sol, solType, tInd):
@@ -134,7 +140,8 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
return {'b':b, 'e':e}
solveOpts = {'factorize':True,'backend':'scipy'}
Solver = Solver
solveOpts = {}
def fields(self, m):
self.makeMassMatrices(m)
@@ -153,7 +160,7 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
dtFact = dt
A = self.getA(tInd)
# print 'Factoring... (dt = ' + str(dt) + ')'
Asolve = Solver(A, options=self.solveOpts)
Asolve = self.Solver(A, **self.solveOpts)
# print 'Done'
rhs = RHS(tInd, F)
sol = Asolve.solve(rhs)
simpegEM/TDEM/SurveyTDEM.py
+1 -1
View File
@@ -47,7 +47,7 @@ class SurveyTDEM1D(BaseSurvey):
def Qrx(self):
if self._Qrx is None:
if self.rxType == 'bz':
locType = 'fz'
locType = 'Fz'
self._Qrx = self.prob.mesh.getInterpolationMat(self.rxLoc, locType=locType)
return self._Qrx
_Qrx = None
simpegEM/TDEM/TDEM_b.py
+13 -10
View File
@@ -16,8 +16,8 @@ class ProblemTDEM_b(ProblemBaseTDEM):
with \\\(\\b\\\) defined on cell faces and \\\(\e\\\) defined on edges.
"""
def __init__(self, model, **kwargs):
ProblemBaseTDEM.__init__(self, model, **kwargs)
def __init__(self, mesh, mapping=None, **kwargs):
ProblemBaseTDEM.__init__(self, mesh, mapping=mapping, **kwargs)
solType = 'b'
@@ -35,11 +35,11 @@ class ProblemTDEM_b(ProblemBaseTDEM):
"""
dt = self.getDt(tInd)
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + (1/dt)*self.MfMui
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + (1.0/dt)*self.MfMui
def getRHS(self, tInd, F):
dt = self.getDt(tInd)
return (1/dt)*self.MfMui*F.get_b(tInd-1)
return (1.0/dt)*self.MfMui*F.get_b(tInd-1)
####################################################
@@ -74,7 +74,8 @@ class ProblemTDEM_b(ProblemBaseTDEM):
if u is None:
u = self.fields(m)
p = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
c = self.mesh.getEdgeMassDeriv()*self.model.transformDeriv(m)*vec
curModel = self.mapping.transform(m)
c = self.mesh.getEdgeInnerProductDeriv(curModel)*self.mapping.transformDeriv(m)*vec
for i in range(self.times.size):
ei = u.get_e(i)
pVal = np.empty_like(ei)
@@ -91,7 +92,9 @@ class ProblemTDEM_b(ProblemBaseTDEM):
tmp = np.zeros((self.mesh.nE,self.survey.nTx))
for i in range(self.nTimes):
tmp += v.get_e(i)*u.get_e(i)
p = -mkvc(self.model.transformDeriv(m).T*self.mesh.getEdgeMassDeriv().T*tmp)
curModel = self.mapping.transform(m)
p = -mkvc(self.mapping.transformDeriv(m).T*self.mesh.getEdgeInnerProductDeriv(curModel).T*tmp)
return p
def solveAh(self, m, p):
@@ -100,7 +103,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
if tInd == 0:
return rhs
dt = self.getDt(tInd)
return rhs + 1./dt*self.MfMui*u.get_b(tInd-1)
return rhs + 1.0/dt*self.MfMui*u.get_b(tInd-1)
def AhCalcFields(sol, solType, tInd):
b = sol
@@ -117,7 +120,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
if tInd == self.nTimes-1:
return rhs
dt = self.getDt(tInd+1)
return rhs + 1./dt*self.MfMui*u.get_b(tInd+1)
return rhs + 1.0/dt*self.MfMui*u.get_b(tInd+1)
def AhtCalcFields(sol, solType, tInd):
b = sol
@@ -167,14 +170,14 @@ class ProblemTDEM_b(ProblemBaseTDEM):
self.makeMassMatrices(sigma)
dt = self.getDt(0)
b = 1/dt*self.MfMui*vec.get_b(0) + self.MfMui*self.mesh.edgeCurl*vec.get_e(0)
b = 1.0/dt*self.MfMui*vec.get_b(0) + self.MfMui*self.mesh.edgeCurl*vec.get_e(0)
e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(0) - self.MeSigma*vec.get_e(0)
f = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
f.set_b(b, 0)
f.set_e(e, 0)
for i in range(1,self.nTimes):
dt = self.getDt(i)
b = 1/dt*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1/dt*self.MfMui*vec.get_b(i-1)
b = 1.0/dt*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1.0/dt*self.MfMui*vec.get_b(i-1)
e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(i) - self.MeSigma*vec.get_e(i)
f.set_b(b, i)
f.set_e(e, i)
simpegEM/Tests/test_FDEM.py
+10 -4
View File
@@ -15,7 +15,7 @@ def getProblem(fdemType, comp):
hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs)))
mesh = Mesh.TensorMesh([hx,hy,hz],[-hx.sum()/2., -hy.sum()/2., -hz.sum()/2.])
model = Model.LogModel(mesh)
mapping = Maps.ExpMap(mesh)
x = np.linspace(-30,30,6)
XYZ = Utils.ndgrid(x,x,np.r_[0])
@@ -25,13 +25,19 @@ def getProblem(fdemType, comp):
survey = EM.FDEM.SurveyFDEM([Tx0])
if fdemType == 'e':
prb = EM.FDEM.ProblemFDEM_e(model)
prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
elif fdemType == 'b':
prb = EM.FDEM.ProblemFDEM_b(model)
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
else:
raise NotImplementedError()
prb.pair(survey)
try:
from mumpsSCI import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
pass
return prb
def adjointTest(fdemType, comp):
@@ -42,7 +48,7 @@ def adjointTest(fdemType, comp):
survey = prb.survey
v = np.random.rand(survey.nD)
w = np.random.rand(prb.model.nP)
w = np.random.rand(prb.mapping.nP)
u = prb.fields(m)
vJw = v.dot(prb.Jvec(m, w, u=u))
simpegEM/Tests/test_FDEM_analytics.py
+2 -2
View File
@@ -17,7 +17,7 @@ class FDEM_analyticTests(unittest.TestCase):
hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs)))
mesh = Mesh.TensorMesh([hx,hy,hz], x0=[-hx.sum()/2.,-hy.sum()/2.,-hz.sum()/2.,])
model = Model.LogModel(mesh)
mapping = Maps.ExpMap(mesh)
x = np.linspace(-10,10,5)
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
@@ -26,7 +26,7 @@ class FDEM_analyticTests(unittest.TestCase):
survey = EM.FDEM.SurveyFDEM([Tx0])
prb = EM.FDEM.ProblemFDEM_b(model)
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
prb.pair(survey)
prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.splu, checkAccuracy=False)
simpegEM/Tests/test_TDEM_b_DerivAdjoint.py
+17 -16
View File
@@ -12,23 +12,23 @@ class TDEM_bDerivTests(unittest.TestCase):
npad = 20
hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs)))
hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2)
mesh = Mesh.CylMesh([hx,1,hy], [0,0,-hy.sum()/2])
active = mesh.vectorCCz<0.
model = Model.ActiveModel(mesh, active, -8, nC=mesh.nCz)
model = Model.ComboModel(mesh,
[Model.LogModel, Model.Vertical1DModel, model])
activeMap = Maps.ActiveCells(mesh, active, -8, nC=mesh.nCz)
mapping = Maps.ComboMap(mesh,
[Maps.ExpMap, Maps.Vertical1DMap, activeMap])
opts = {'txLoc':0.,
'txType':'VMD_MVP',
'rxLoc':np.r_[150., 0.],
'txType': 'VMD_MVP',
'rxLoc':np.r_[40., 0., 0.],
'rxType':'bz',
'timeCh':np.logspace(-4,-2,20),
}
self.dat = EM.TDEM.SurveyTDEM1D(**opts)
self.prb = EM.TDEM.ProblemTDEM_b(model)
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10])
self.sigma = np.ones(mesh.nCz)*1e-8
@@ -50,7 +50,8 @@ class TDEM_bDerivTests(unittest.TestCase):
Ahu = prb.AhVec(sigma, u)
V1 = Ahu.get_b(0)
V2 = 1/prb.getDt(0)*prb.MfMui*u.get_b(-1)
V2 = 1./prb.getDt(0)*prb.MfMui*u.get_b(-1)
print np.linalg.norm(V1-V2), np.linalg.norm(V2), np.linalg.norm(V1-V2)/np.linalg.norm(V2)
self.assertTrue(np.linalg.norm(V1-V2)/np.linalg.norm(V2) < 1.e-6)
V1 = Ahu.get_e(0)
@@ -137,10 +138,10 @@ class TDEM_bDerivTests(unittest.TestCase):
"""
# Random model and perturbation
sigma = np.random.rand(self.prb.model.nP)
sigma = np.random.rand(self.prb.mapping.nP)
f = self.prb.fields(sigma)
dm = 1000*np.random.rand(self.prb.model.nP)
dm = 1000*np.random.rand(self.prb.mapping.nP)
h = 0.01
derChk = lambda m: [self.prb.AhVec(m, f).fieldVec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).fieldVec()]
@@ -155,7 +156,7 @@ class TDEM_bDerivTests(unittest.TestCase):
mesh = self.mesh
sigma = self.sigma
dm = 10*np.random.rand(prb.model.nP)
dm = 10*np.random.rand(prb.mapping.nP)
f = prb.fields(sigma)
derChk = lambda m: [self.prb.fields(m).fieldVec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).fieldVec()]
@@ -172,7 +173,7 @@ class TDEM_bDerivTests(unittest.TestCase):
sigma = self.sigma
# d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
d_sig = 10*np.random.rand(prb.model.nP)
d_sig = 10*np.random.rand(prb.mapping.nP)
derChk = lambda m: [prb.survey.dpred(m), lambda mx: -prb.Jvec(sigma, mx)]
@@ -221,7 +222,7 @@ class TDEM_bDerivTests(unittest.TestCase):
def test_solveAhtVsAhtVec(self):
prb = self.prb
mesh = self.mesh
sigma = np.random.rand(prb.model.nP)
sigma = np.random.rand(prb.mapping.nP)
f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nTimes, 'b')
for i in range(f1.nTimes):
@@ -258,8 +259,8 @@ class TDEM_bDerivTests(unittest.TestCase):
mesh = self.mesh
prb = self.prb
m = np.random.rand(prb.model.nP)
sigma = np.random.rand(prb.model.nP)
m = np.random.rand(prb.mapping.nP)
sigma = np.random.rand(prb.mapping.nP)
u = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b')
for i in range(u.nTimes):
@@ -280,7 +281,7 @@ class TDEM_bDerivTests(unittest.TestCase):
prb = self.prb
sigma = self.sigma
m = np.random.rand(prb.model.nP)
m = np.random.rand(prb.mapping.nP)
d = np.random.rand(prb.nTimes)
V1 = d.dot(prb.Jvec(sigma, m))
simpegEM/Tests/test_forward_EMproblem.py
+88 -20
View File
@@ -5,54 +5,122 @@ from scipy.constants import mu_0
from simpegEM.Utils.Ana import hzAnalyticDipoleT
import matplotlib.pyplot as plt
import simpegem1d as EM1D
class TDEM_bTests(unittest.TestCase):
def setUp(self):
cs = 10.
ncx = 15
ncy = 10
npad = 20
hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs)))
hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2)
# cs = 20.
# ncx = 15
# ncy = 10
# npad = 15
# hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs)))
# hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
# mesh = Mesh.CylMesh([hx,1,hz], [0,0,-hz.sum()/2])
cs, nc, npad = 20., 15, 10
hx = Utils.meshTensors(((npad,cs), (nc,cs), (npad,cs)))
hy = Utils.meshTensors(((npad,cs), (nc,cs), (npad,cs)))
hz = Utils.meshTensors(((npad,cs), (nc,cs), (npad,cs)))
mesh = Mesh.TensorMesh([hx,hy,hz], [-hx.sum()/2.,-hy.sum()/2.,-hz.sum()/2.])
active = mesh.vectorCCz<0.
model = Model.ActiveModel(mesh, active, -8, nC=mesh.nCz)
model = Model.ComboModel(mesh,
[Model.LogModel, Model.Vertical1DModel, model])
actMap = Maps.ActiveCells(mesh, active, -8, nC=mesh.nCz)
mapping = Maps.ComboMap(mesh,
[Maps.ExpMap, Maps.Vertical1DMap, actMap])
opts = {'txLoc':0.,
opts = {'txLoc':np.array([0., 0., 0.]),
'txType':'VMD_MVP',
'rxLoc':np.r_[30., 0.],
'rxLoc':np.array([[10., 0., 0.]]),
'rxType':'bz',
'timeCh':np.logspace(-4,-2.5, 21),
'timeCh':np.logspace(-5,-4, 21),
}
sig_half = 1e-3
self.sig_half = sig_half
self.dat = EM.TDEM.SurveyTDEM1D(**opts)
self.prb = EM.TDEM.ProblemTDEM_b(model)
self.prb.setTimes([1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4], [40, 40, 40, 40, 40, 40])
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
try:
from mumpsSCI import MumpsSolver
self.prb.Solver = MumpsSolver
except ImportError:
pass
self.prb.setTimes([1e-6], [100])
# self.prb.setTimes([1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4], [40, 40, 40, 40, 40, 40])
self.sigma = np.ones(mesh.nCz)*1e-8
self.sigma[mesh.vectorCCz<0] = 1e-3
self.sigma[active] = sig_half
self.sigma = np.log(self.sigma[active])
self.showIt = False
self.showIt = True
self.prb.pair(self.dat)
TDsurvey = EM1D.BaseEM1D.EM1DSurveyTD()
TDsurvey.rxLoc = np.array([0., 0., 30.])
TDsurvey.txLoc = np.array([0., 0., 80.])
TDsurvey.fieldtype = 'secondary'
TDsurvey.waveType = 'stepoff'
TDsurvey.time = self.prb.times #np.logspace(-5, -2, 64)
TDsurvey.setFrequency(TDsurvey.time)
nearthick = np.logspace(-1, 1, 5)
deepthick = np.logspace(1, 2, 10)
hx = np.r_[nearthick, deepthick]
mesh1D = Mesh.TensorMesh([hx], [0.])
depth = -mesh1D.gridN
LocSigZ = -mesh1D.gridCC
nlay = depth.size
topo = np.r_[0., 0., 0.]
TDsurvey.depth = depth
TDsurvey.topo = topo
TDsurvey.LocSigZ = LocSigZ
TDsurvey.HalfSwitch = True
TDsurvey.Setup1Dsystem()
chi_half = 0.
expmap = EM1D.BaseEM1D.BaseEM1DMap(mesh1D)
mappingReal = Maps.ComboMap(mesh1D, [expmap])
m_1D = np.log(np.ones(nlay)*sig_half)
TDsurvey.rxType = 'Bz'
WT0, WT1, YBASE = EM1D.DigFilter.LoadWeights()
options = {'WT0': WT0, 'WT1': WT1, 'YBASE': YBASE}
prob = EM1D.EM1D.EM1D(mesh1D, mapping=mappingReal, **options)
prob.pair(TDsurvey)
prob.chi = np.zeros(TDsurvey.nlay)
survey = TDsurvey
options = options
prob.CondType = 'Real'
prob.survey.txType = 'VMD'
prob.survey.offset = 1e-5
m_1D = np.log(np.ones(prob.survey.nlay)*sig_half)
Bz = survey.dpred(m_1D)
self.Bzanal = Bz
def test_analitic_b(self):
bz_calc = self.dat.dpred(self.sigma)
bz_ana = mu_0*hzAnalyticDipoleT(self.dat.rxLoc[0], self.prb.times, np.exp(self.sigma[0]))
# bz_ana = self.Bzanal
bz_ana = mu_0*hzAnalyticDipoleT(self.dat.rxLoc[0,0], self.prb.times, self.sig_half)
ind = self.prb.times > 1e-5
diff = np.linalg.norm(bz_calc[ind].flatten() - bz_ana[ind].flatten())/np.linalg.norm(bz_ana[ind].flatten())
if self.showIt == True:
plt.loglog(self.prb.times[bz_calc>0], bz_calc[bz_calc>0], 'b', self.prb.times[bz_calc<0], -bz_calc[bz_calc<0], 'b--')
plt.loglog(self.prb.times[bz_calc>0], bz_calc[bz_calc>0], 'r', self.prb.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--')
plt.loglog(self.prb.times, abs(bz_ana), 'b*')
plt.xlim(1e-5, 1e-2)
plt.show()
print diff
print 'Difference: ', diff
self.assertTrue(diff < 0.10)