mirror of
https://github.com/wassname/simpeg.git
synced 2026-08-11 11:26:01 +08:00
forward is running, but not passing
This commit is contained in:
Lindsey Heagy
@@ -5,7 +5,7 @@ from SimPEG import Utils
|
||||
from SimPEG.EM.Utils import omega
|
||||
from SimPEG.Utils import Zero, Identity
|
||||
|
||||
class Fields(SimPEG.Problem.Fields):
|
||||
class Fields(SimPEG.Problem.TimeFields):
|
||||
"""
|
||||
|
||||
Fancy Field Storage for a TDEM survey. Only one field type is stored for
|
||||
@@ -32,8 +32,6 @@ class Fields(SimPEG.Problem.Fields):
|
||||
dtype = float
|
||||
|
||||
|
||||
|
||||
|
||||
class Fields_b(Fields):
|
||||
"""Fancy Field Storage for a TDEM survey."""
|
||||
knownFields = {'bSolution': 'F'}
|
||||
@@ -57,7 +55,11 @@ class Fields_b(Fields):
|
||||
return Zero()
|
||||
|
||||
def _e(self, bSolution, srcList, tInd):
|
||||
return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) )
|
||||
e = self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) )
|
||||
for i, src in enumerate(srcList):
|
||||
_, S_e = src.eval(self.prob, tInd)
|
||||
e[:,i,tInd] = e[:,i,tInd] - self.MeSigmaI * S_e
|
||||
return e
|
||||
|
||||
def _eDeriv_u(self, src, du_dm_v, adjoint = False):
|
||||
raise NotImplementedError
|
||||
|
||||
+106
-13
@@ -1,13 +1,14 @@
|
||||
from SimPEG import Survey, np
|
||||
from SimPEG.Survey import BaseSurvey
|
||||
import SimPEG
|
||||
from SimPEG import np, Utils
|
||||
from SimPEG.Utils import Zero, Identity
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils import *
|
||||
|
||||
####################################################
|
||||
# Receivers
|
||||
####################################################
|
||||
|
||||
class Rx(Survey.BaseTimeRx):
|
||||
class Rx(SimPEG.Survey.BaseTimeRx):
|
||||
|
||||
knownRxTypes = {
|
||||
'ex':['e', 'Ex', 'N'],
|
||||
@@ -24,7 +25,7 @@ class Rx(Survey.BaseTimeRx):
|
||||
}
|
||||
|
||||
def __init__(self, locs, times, rxType):
|
||||
Survey.BaseTimeRx.__init__(self, locs, times, rxType)
|
||||
SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType)
|
||||
|
||||
@property
|
||||
def projField(self):
|
||||
@@ -77,6 +78,10 @@ class BaseWaveform(object):
|
||||
self.offTime = offTime
|
||||
self.hasInitialFields = hasInitialFields
|
||||
|
||||
def _assertMatchesPair(self, pair):
|
||||
assert (isinstance(self, pair)
|
||||
), "Waveform object must be an instance of a %s BaseWaveform class."%(pair.__name__)
|
||||
|
||||
def eval(self, time):
|
||||
raise NotImplementedError
|
||||
|
||||
@@ -91,11 +96,11 @@ class StepOffWaveform(BaseWaveform):
|
||||
|
||||
|
||||
|
||||
class BaseSrc(Survey.BaseSrc):
|
||||
class BaseSrc(SimPEG.Survey.BaseSrc):
|
||||
|
||||
rxPair = Rx
|
||||
integrate = True
|
||||
waveformPair = StepOffWaveform
|
||||
waveformPair = BaseWaveform
|
||||
|
||||
@property
|
||||
def waveform(self):
|
||||
@@ -103,17 +108,17 @@ class BaseSrc(Survey.BaseSrc):
|
||||
return getattr(self, '_waveform', None)
|
||||
@waveform.setter
|
||||
def waveform(self, val):
|
||||
if self.waveform is None:
|
||||
val._assertMatchesPair(self.waveformPair)
|
||||
self._waveform = val
|
||||
# if self.waveform is None:
|
||||
val._assertMatchesPair(self.waveformPair)
|
||||
self._waveform = val
|
||||
|
||||
|
||||
def __init__(self, rxList, waveform = None):
|
||||
self.waveform = waveform
|
||||
Survey.BaseSrc.__init__(self, rxList)
|
||||
SimPEG.Survey.BaseSrc.__init__(self, rxList)
|
||||
|
||||
|
||||
def bInitial(self, mesh):
|
||||
def bInitial(self):
|
||||
return Zero()
|
||||
|
||||
def eval(self, prob, time):
|
||||
@@ -135,7 +140,7 @@ class BaseSrc(Survey.BaseSrc):
|
||||
|
||||
|
||||
class MagDipole(BaseSrc):
|
||||
def __init__(self, rxList, waveform, loc, orientation='Z', moment=1., mu=mu_0, **kwargs):
|
||||
def __init__(self, rxList, waveform=None, loc=None, orientation='Z', moment=1., mu=mu_0):
|
||||
|
||||
self.loc = loc
|
||||
self.orientation = orientation
|
||||
@@ -143,10 +148,98 @@ class MagDipole(BaseSrc):
|
||||
self.moment = moment
|
||||
self.mu = mu
|
||||
self.integrate = False
|
||||
BaseSrc.__init__(self, rxList)
|
||||
BaseSrc.__init__(self, rxList, waveform)
|
||||
|
||||
def _bfromVectorPotential(self, prob):
|
||||
if prob._eqLocs is 'FE':
|
||||
gridX = prob.mesh.gridEx
|
||||
gridY = prob.mesh.gridEy
|
||||
gridZ = prob.mesh.gridEz
|
||||
C = prob.mesh.edgeCurl
|
||||
|
||||
elif prob._eqLocs is 'EF':
|
||||
gridX = prob.mesh.gridFx
|
||||
gridY = prob.mesh.gridFy
|
||||
gridZ = prob.mesh.gridFz
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
|
||||
if prob.mesh._meshType is 'CYL':
|
||||
if not prob.mesh.isSymmetric:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
|
||||
|
||||
else:
|
||||
srcfct = MagneticDipoleVectorPotential
|
||||
ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
|
||||
ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
|
||||
az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
|
||||
a = np.concatenate((ax, ay, az))
|
||||
|
||||
return C*a
|
||||
|
||||
|
||||
def bInitial(self, prob):
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
if self.waveform.hasInitialFields is False:
|
||||
return Zero()
|
||||
|
||||
return self._bfromVectorPotential(prob)
|
||||
|
||||
def S_m(self, prob, time):
|
||||
if self.waveform.hasInitialFields is False:
|
||||
raise NotImplementedError
|
||||
return Zero()
|
||||
|
||||
def S_e(self, prob, time):
|
||||
if self.waveform.hasInitialFields is False:
|
||||
raise NotImplementedError
|
||||
return Zero()
|
||||
|
||||
####################################################
|
||||
# Survey
|
||||
####################################################
|
||||
|
||||
class Survey(SimPEG.Survey.BaseSurvey):
|
||||
"""
|
||||
Time domain electromagnetic survey
|
||||
"""
|
||||
|
||||
srcPair = BaseSrc
|
||||
rxPair = Rx
|
||||
|
||||
def __init__(self, srcList, **kwargs):
|
||||
# Sort these by frequency
|
||||
self.srcList = srcList
|
||||
SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
|
||||
def eval(self, u):
|
||||
data = SimPEG.Survey.Data(self)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.eval(src, self.mesh, self.prob.timeMesh, u)
|
||||
return data
|
||||
|
||||
def evalDeriv(self, u, v=None, adjoint=False):
|
||||
assert v is not None, 'v to multiply must be provided.'
|
||||
|
||||
if not adjoint:
|
||||
data = SimPEG.Survey.Data(self)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
|
||||
return data
|
||||
else:
|
||||
f = FieldsTDEM(self.mesh, self)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
|
||||
Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
|
||||
if rx.projField not in f: # first time we are projecting
|
||||
f[src, rx.projField, :] = Ptv
|
||||
else: # there are already fields, so let's add to them!
|
||||
f[src, rx.projField, :] += Ptv
|
||||
return f
|
||||
|
||||
|
||||
|
||||
+53
-11
@@ -29,20 +29,59 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
|
||||
tic = time.time()
|
||||
self.curModel = m
|
||||
|
||||
F = self.fieldsPair(self.mesh, self.survey)
|
||||
|
||||
# for
|
||||
# set initial fields
|
||||
for i, src in enumerate(self.survey.srcList):
|
||||
F[src,'bSolution',0] = src.bInitial(self)
|
||||
|
||||
# timestep to solve forward
|
||||
Ainv = None
|
||||
for tInd, dt in enumerate(self.timeSteps):
|
||||
print dt, self.timeSteps[tInd]
|
||||
if Ainv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same
|
||||
Ainv.clean()
|
||||
|
||||
if Ainv is None:
|
||||
A = self.getA(tInd)
|
||||
if self.verbose: print 'Factoring... (dt = %e)'%dt
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
if self.verbose: print 'Done'
|
||||
|
||||
rhs = self.getRHS(tInd, F)
|
||||
if self.verbose: print ' Solving... (tInd = %d)'%tInd
|
||||
sol = Ainv * rhs
|
||||
if self.verbose: print ' Done...'
|
||||
if sol.ndim == 1:
|
||||
sol.shape = (sol.size,1)
|
||||
F[:,self._fieldType+'Solution',tInd+1] = sol
|
||||
Ainv.clean()
|
||||
return F
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
return None
|
||||
raise NotImplementedError
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
return None
|
||||
raise NotImplementedError
|
||||
|
||||
def getSourceTerm(self, tInd):
|
||||
return None
|
||||
# return S_m, S_e
|
||||
|
||||
Srcs = self.survey.srcList
|
||||
|
||||
if self._eqLocs is 'FE':
|
||||
S_m = np.zeros((self.mesh.nF,len(Srcs)))
|
||||
S_e = np.zeros((self.mesh.nE,len(Srcs)))
|
||||
elif self._eqLocs is 'EF':
|
||||
S_m = np.zeros((self.mesh.nE,len(Srcs)))
|
||||
S_e = np.zeros((self.mesh.nF,len(Srcs)))
|
||||
|
||||
for i, src in enumerate(Srcs):
|
||||
smi, sei = src.eval(self, self.times[tInd])
|
||||
S_m[:,i] = S_m[:,i] + smi
|
||||
S_e[:,i] = S_e[:,i] + sei
|
||||
|
||||
return S_m, S_e
|
||||
|
||||
|
||||
|
||||
@@ -114,7 +153,7 @@ class Problem_b(BaseTDEMProblem):
|
||||
A = I + dt * ( C * ( MeSigmaI * (C.T * MfMui ) ) )
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return MeMui.T * A
|
||||
return MfMui.T * A
|
||||
return A
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
@@ -135,21 +174,24 @@ class Problem_b(BaseTDEMProblem):
|
||||
return ADeriv
|
||||
|
||||
|
||||
def getRHS(self, tInd):
|
||||
def getRHS(self, tInd, F):
|
||||
dt = self.timeSteps[tInd]
|
||||
C = self.mesh.edgeCurl
|
||||
MeSigmaI = self.MeSigmaI
|
||||
MfMui = self.MfMui
|
||||
|
||||
S_m, S_e = self.getSourceTerm(tInd+1) # I think this is tInd+1 ?
|
||||
S_m, S_e = self.getSourceTerm(tInd+1)
|
||||
|
||||
B_n = np.c_[[F[src,'b',tInd] for src in self.survey.srcList]].T
|
||||
B_n = np.c_[[F[src,'bSolution',tInd] for src in self.survey.srcList]].T
|
||||
if B_n.shape[0] is not 1:
|
||||
raise NotImplementedError('getRHS not implemented for this shape of B_n')
|
||||
|
||||
return B_n + dt * (C * (MeSigmaIDeriv * S_e) + S_m)
|
||||
rhs = B_n[:,:,0].T + dt * (C * (MeSigmaI * S_e) + S_m)
|
||||
if self._makeASymmetric:
|
||||
return MfMui.T * rhs
|
||||
return rhs
|
||||
|
||||
def getRHSDeriv(self, tInd, src, v, adjoint=False):
|
||||
def getRHSDeriv(self, tInd, F, src, v, adjoint=False):
|
||||
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
@@ -10,7 +10,7 @@ except ImportError, e:
|
||||
MumpsSolver = SolverLU
|
||||
|
||||
|
||||
def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,1e-3], showIt=False):
|
||||
def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,1e-3], showIt=True):
|
||||
if meshType == 'CYL':
|
||||
cs, ncx, ncz, npad = 5., 30, 10, 15
|
||||
hx = [(cs,ncx), (cs,npad,1.3)]
|
||||
@@ -28,11 +28,11 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
|
||||
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
|
||||
|
||||
rx = EM.TDEM.Rx(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.SurveyTDEM.MagDipole([rx], waveform= EM.TDEM.SurveyTDEM.StepOffWaveform(), 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)
|
||||
survey = EM.TDEM.Survey([src])
|
||||
prb = EM.TDEM.Problem_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)]
|
||||
@@ -48,6 +48,8 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
|
||||
|
||||
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 ' |bz_ana| = ',np.linalg.norm(bz_ana), ' |bz_num| = ', np.linalg.norm(bz_calc), ' |bz_ana - bz_num| =', np.linalg.norm(bz_ana-bz_calc)
|
||||
print 'Difference: ', log10diff
|
||||
|
||||
if showIt == True:
|
||||
@@ -61,35 +63,35 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
|
||||
|
||||
class TDEM_SimpleSrcTests(unittest.TestCase):
|
||||
def test_source(self):
|
||||
waveform = EM.TDEM.SurveyTDEM.StepOffWaveform([])
|
||||
waveform = EM.TDEM.SurveyTDEM.StepOffWaveform()
|
||||
assert waveform.eval(0.) == 0.
|
||||
|
||||
|
||||
|
||||
|
||||
# class TDEM_bTests(unittest.TestCase):
|
||||
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_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)
|
||||
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)
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user