mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-27 19:48:52 +08:00
sketching out TDEM problem
This commit is contained in:
Lindsey Heagy
@@ -0,0 +1,68 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import SimPEG
|
||||
from SimPEG import Utils
|
||||
from SimPEG.EM.Utils import omega
|
||||
from SimPEG.Utils import Zero, Identity
|
||||
|
||||
class Fields(SimPEG.Problem.Fields):
|
||||
"""
|
||||
|
||||
Fancy Field Storage for a TDEM survey. Only one field type is stored for
|
||||
each problem, the rest are computed. The fields obejct acts like an array and is indexed by
|
||||
|
||||
.. code-block:: python
|
||||
|
||||
f = problem.fields(m)
|
||||
e = f[srcList,'e']
|
||||
b = f[srcList,'b']
|
||||
|
||||
If accessing all sources for a given field, use the :code:`:`
|
||||
|
||||
.. code-block:: python
|
||||
|
||||
f = problem.fields(m)
|
||||
e = f[:,'e']
|
||||
b = f[:,'b']
|
||||
|
||||
The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies)
|
||||
"""
|
||||
|
||||
knownFields = {}
|
||||
dtype = float
|
||||
|
||||
|
||||
|
||||
|
||||
class Fields_b(Fields):
|
||||
"""Fancy Field Storage for a TDEM survey."""
|
||||
knownFields = {'bSolution': 'F'}
|
||||
aliasFields = {
|
||||
'b': ['bSolution', 'F', '_b'],
|
||||
'e': ['bSolution', 'E', '_e'],
|
||||
}
|
||||
|
||||
def startup(self):
|
||||
self.MeSigmaI = self.survey.prob.MeSigmaI
|
||||
self.edgeCurl = self.survey.prob.mesh.edgeCurl
|
||||
self.MfMui = self.survey.prob.MfMui
|
||||
|
||||
def _b(self, bSolution, srcList, tInd):
|
||||
return bSolution
|
||||
|
||||
def _bDeriv_u(self, src, du_dm_v, adjoint = False):
|
||||
return Identity()*v
|
||||
|
||||
def _bDeriv_m(self, src, v, adjoint = False):
|
||||
return Zero()
|
||||
|
||||
def _e(self, bSolution, srcList, tInd):
|
||||
return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) )
|
||||
|
||||
def _eDeriv_u(self, src, du_dm_v, adjoint = False):
|
||||
raise NotImplementedError
|
||||
|
||||
def _eDeriv_m(self, src, v, adjoint = False):
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
@@ -1,10 +1,7 @@
|
||||
from SimPEG import Utils, Survey, np
|
||||
from SimPEG.Survey import BaseSurvey
|
||||
from SimPEG.EM.Utils import *
|
||||
from BaseTDEM import FieldsTDEM
|
||||
|
||||
|
||||
class RxTDEM(Survey.BaseTimeRx):
|
||||
class Rx(Survey.BaseTimeRx):
|
||||
|
||||
knownRxTypes = {
|
||||
'ex':['e', 'Ex', 'N'],
|
||||
@@ -62,138 +59,4 @@ class RxTDEM(Survey.BaseTimeRx):
|
||||
if not adjoint:
|
||||
return P * Utils.mkvc(v[src, self.projField, :])
|
||||
elif adjoint:
|
||||
return P.T * v[src, self]
|
||||
|
||||
|
||||
class SrcTDEM(Survey.BaseSrc):
|
||||
rxPair = RxTDEM
|
||||
radius = None
|
||||
|
||||
def getInitialFields(self, mesh):
|
||||
F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
|
||||
return F0
|
||||
|
||||
def getJs(self, mesh, time):
|
||||
return None
|
||||
|
||||
|
||||
class SrcTDEM_VMD_MVP(SrcTDEM):
|
||||
|
||||
def __init__(self,rxList,loc,waveformType="STEPOFF"):
|
||||
self.loc = loc
|
||||
self.waveformType = waveformType
|
||||
SrcTDEM.__init__(self,rxList)
|
||||
|
||||
def getInitialFields(self, mesh):
|
||||
"""Vertical magnetic dipole, magnetic vector potential"""
|
||||
if self.waveformType == "STEPOFF":
|
||||
print ">> Step waveform: Non-zero initial condition"
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
return {"b": mesh.edgeCurl*MVP}
|
||||
elif self.waveformType == "GENERAL":
|
||||
print ">> General waveform: Zero initial condition"
|
||||
return {"b": np.zeros(mesh.nF)}
|
||||
else:
|
||||
raise NotImplementedError("Only use STEPOFF or GENERAL")
|
||||
|
||||
def getMeS(self, mesh, MfMui):
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
|
||||
|
||||
|
||||
class SrcTDEM_CircularLoop_MVP(SrcTDEM):
|
||||
def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
|
||||
self.loc = loc
|
||||
self.radius = radius
|
||||
self.waveformType = waveformType
|
||||
SrcTDEM.__init__(self,rxList)
|
||||
|
||||
def getInitialFields(self, mesh):
|
||||
"""Circular Loop, magnetic vector potential"""
|
||||
if self.waveformType == "STEPOFF":
|
||||
print ">> Step waveform: Non-zero initial condition"
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
|
||||
else:
|
||||
raise Exception('Unknown mesh for CircularLoop')
|
||||
return {"b": mesh.edgeCurl*MVP}
|
||||
elif self.waveformType == "GENERAL":
|
||||
print ">> General waveform: Zero initial condition"
|
||||
return {"b": np.zeros(mesh.nF)}
|
||||
else:
|
||||
raise NotImplementedError("Only use STEPOFF or GENERAL")
|
||||
|
||||
def getMeS(self, mesh, MfMui):
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
|
||||
else:
|
||||
raise Exception('Unknown mesh for CircularLoop')
|
||||
return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
|
||||
|
||||
|
||||
class SurveyTDEM(Survey.BaseSurvey):
|
||||
"""
|
||||
docstring for SurveyTDEM
|
||||
"""
|
||||
srcPair = SrcTDEM
|
||||
|
||||
def __init__(self, srcList, **kwargs):
|
||||
# Sort these by frequency
|
||||
self.srcList = srcList
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
|
||||
def projectFields(self, u):
|
||||
data = Survey.Data(self)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.projectFields(src, self.mesh, self.prob.timeMesh, u)
|
||||
return data
|
||||
|
||||
def projectFieldsDeriv(self, u, v=None, adjoint=False):
|
||||
assert v is not None, 'v to multiply must be provided.'
|
||||
|
||||
if not adjoint:
|
||||
data = Survey.Data(self)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.projectFieldsDeriv(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.projectFieldsDeriv(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
|
||||
|
||||
|
||||
return P.T * v[src, self]
|
||||
@@ -0,0 +1,180 @@
|
||||
from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
|
||||
from SimPEG.EM.Base import BaseEMProblem
|
||||
from SimPEG.EM.TDEM.SurveyTDEM import Survey as SurveyTDEM
|
||||
from SimPEG.EM.TDEM.FieldsTDEM import *
|
||||
from scipy.constants import mu_0
|
||||
import time
|
||||
|
||||
class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
|
||||
"""
|
||||
We start with the first order form of Maxwell's equations
|
||||
"""
|
||||
surveyPair = SurveyTDEM
|
||||
fieldsPair = Fields
|
||||
|
||||
def __init__(self, mesh, mapping=None, **kwargs):
|
||||
Problem.BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
|
||||
|
||||
|
||||
# _FieldsForward_pair = FieldsTDEM #: used for the forward calculation only
|
||||
|
||||
def fields(self, m):
|
||||
"""
|
||||
Solve the forward problem for the fields.
|
||||
|
||||
:param numpy.array m: inversion model (nP,)
|
||||
:rtype numpy.array:
|
||||
:return F: fields
|
||||
"""
|
||||
|
||||
tic = time.time()
|
||||
self.curModel = m
|
||||
F = self.fieldsPair(self.mesh, self.survey)
|
||||
|
||||
# for
|
||||
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
return None
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
return None
|
||||
|
||||
def getSourceTerm(self, tInd):
|
||||
return None
|
||||
# return S_m, S_e
|
||||
|
||||
|
||||
|
||||
##########################################################################################
|
||||
################################ E-B Formulation #########################################
|
||||
##########################################################################################
|
||||
|
||||
class Problem_b(BaseTDEMProblem):
|
||||
"""
|
||||
Starting from the quasi-static E-B formulation of Maxwell's equations (semi-discretized)
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{C} \mathbf{e} + \\frac{\partial \mathbf{b}}{\partial t} = \mathbf{s_m} \\\\
|
||||
\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}
|
||||
|
||||
where :math:`\mathbf{s_e}` is an integrated quantity, we eliminate :math:`\mathbf{e}` using
|
||||
|
||||
.. math::
|
||||
\mathbf{e} = \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}
|
||||
|
||||
to obtain a second order semi-discretized system in :math:`\mathbf{b}`
|
||||
|
||||
.. math::
|
||||
\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} + \\frac{\partial \mathbf{b}}{\partial t} = \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e} + \mathbf{s_m}
|
||||
|
||||
and moving everything except the time derivative to the rhs gives
|
||||
|
||||
.. math::
|
||||
\\frac{\partial \mathbf{b}}{\partial t} = -\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} + \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e} + \mathbf{s_m}
|
||||
|
||||
For the time discretization, we use backward euler. To solve for the :math:`n+1`th time step, we have
|
||||
|
||||
.. math::
|
||||
\\frac{\mathbf{b}^{n+1} - \mathbf{b}^{n}}{\mathbf{dt}} = -\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b}^{n+1} + \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}^{n+1} + \mathbf{s_m}^{n+1}
|
||||
|
||||
re-arranging to put :math:`\mathbf{b}^{n+1}` on the left hand side gives
|
||||
|
||||
.. math::
|
||||
(\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}) \mathbf{b}^{n+1} = \mathbf{b}^{n} + \mathbf{dt}(\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}^{n+1} + \mathbf{s_m}^{n+1})
|
||||
|
||||
:param Mesh mesh: mesh
|
||||
:param Mapping mapping: mapping
|
||||
"""
|
||||
|
||||
_fieldType = 'b'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = Fields_b
|
||||
surveyPair = SurveyTDEM
|
||||
|
||||
def __init__(self, mesh, mapping=None, **kwargs):
|
||||
BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs)
|
||||
|
||||
def getA(self, tInd):
|
||||
"""
|
||||
System matrix at a given time index
|
||||
|
||||
.. math::
|
||||
(\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f})
|
||||
|
||||
"""
|
||||
|
||||
dt = self.timeSteps[tInd]
|
||||
C = self.mesh.edgeCurl
|
||||
MeSigmaI = self.MeSigmaI
|
||||
MfMui = self.MfMui
|
||||
I = Utils.speye(self.mesh.nF)
|
||||
|
||||
A = I + dt * ( C * ( MeSigmaI * (C.T * MfMui ) ) )
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return MeMui.T * A
|
||||
return A
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
dt = self.timeSteps[tInd]
|
||||
C = self.mesh.edgeCurl
|
||||
MeSigmaI = self.MeSigmaIDeriv
|
||||
MfMui = self.MfMui
|
||||
I = Utils.speye(self.mesh.nF)
|
||||
|
||||
if adjoint:
|
||||
if self._makeASymmetric is True:
|
||||
v = MfMui * v
|
||||
return dt * MfMui.T * ( C * ( MeSigmaIDeriv.T * ( C.T * v ) ) )
|
||||
|
||||
ADeriv = dt * ( C * ( MeSigmaIDeriv * (C.T * ( MfMui * v ) ) ) )
|
||||
if self._makeASymmetric is True:
|
||||
return MeMui.T * ADeriv
|
||||
return ADeriv
|
||||
|
||||
|
||||
def getRHS(self, tInd):
|
||||
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 ?
|
||||
|
||||
B_n = np.c_[[F[src,'b',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)
|
||||
|
||||
def getRHSDeriv(self, tInd, src, v, adjoint=False):
|
||||
|
||||
raise NotImplementedError
|
||||
|
||||
dt = self.timeSteps[tInd]
|
||||
C = self.mesh.edgeCurl
|
||||
MeSigmaI = self.MeSigmaI
|
||||
MeSigmaIDeriv = self.MeSigmaIDeriv
|
||||
MfMui = self.MfMui
|
||||
|
||||
S_m, S_e = src.eval(tInd+1, self) # I think this is tInd+1 ?
|
||||
S_m, S_e = src.evalDeriv(tInd+1, self, adjoint=adjoint) # I think this is tInd+1 ?
|
||||
|
||||
B_n = np.c_[[F[src,'b',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)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -1,3 +1,3 @@
|
||||
from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
|
||||
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
|
||||
from TDEM_b import ProblemTDEM_b
|
||||
from TDEM import BaseTDEMProblem, Problem_b
|
||||
from FieldsTDEM import Fields, Fields_b
|
||||
from SurveyTDEM import Survey, Src, Rx
|
||||
@@ -0,0 +1,199 @@
|
||||
from SimPEG import Utils, Survey, np
|
||||
from SimPEG.Survey import BaseSurvey
|
||||
from SimPEG.EM.Utils import *
|
||||
from BaseTDEM import FieldsTDEM
|
||||
|
||||
|
||||
class RxTDEM(Survey.BaseTimeRx):
|
||||
|
||||
knownRxTypes = {
|
||||
'ex':['e', 'Ex', 'N'],
|
||||
'ey':['e', 'Ey', 'N'],
|
||||
'ez':['e', 'Ez', 'N'],
|
||||
|
||||
'bx':['b', 'Fx', 'N'],
|
||||
'by':['b', 'Fy', 'N'],
|
||||
'bz':['b', 'Fz', 'N'],
|
||||
|
||||
'dbxdt':['b', 'Fx', 'CC'],
|
||||
'dbydt':['b', 'Fy', 'CC'],
|
||||
'dbzdt':['b', 'Fz', 'CC'],
|
||||
}
|
||||
|
||||
def __init__(self, locs, times, rxType):
|
||||
Survey.BaseTimeRx.__init__(self, locs, times, rxType)
|
||||
|
||||
@property
|
||||
def projField(self):
|
||||
"""Field Type projection (e.g. e b ...)"""
|
||||
return self.knownRxTypes[self.rxType][0]
|
||||
|
||||
@property
|
||||
def projGLoc(self):
|
||||
"""Grid Location projection (e.g. Ex Fy ...)"""
|
||||
return self.knownRxTypes[self.rxType][1]
|
||||
|
||||
@property
|
||||
def projTLoc(self):
|
||||
"""Time Location projection (e.g. CC N)"""
|
||||
return self.knownRxTypes[self.rxType][2]
|
||||
|
||||
def getTimeP(self, timeMesh):
|
||||
"""
|
||||
Returns the time projection matrix.
|
||||
|
||||
.. note::
|
||||
|
||||
This is not stored in memory, but is created on demand.
|
||||
"""
|
||||
if self.rxType in ['dbxdt','dbydt','dbzdt']:
|
||||
return timeMesh.getInterpolationMat(self.times, self.projTLoc)*timeMesh.faceDiv
|
||||
else:
|
||||
return timeMesh.getInterpolationMat(self.times, self.projTLoc)
|
||||
|
||||
def projectFields(self, src, mesh, timeMesh, u):
|
||||
P = self.getP(mesh, timeMesh)
|
||||
u_part = Utils.mkvc(u[src, self.projField, :])
|
||||
return P*u_part
|
||||
|
||||
def projectFieldsDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
|
||||
P = self.getP(mesh, timeMesh)
|
||||
|
||||
if not adjoint:
|
||||
return P * Utils.mkvc(v[src, self.projField, :])
|
||||
elif adjoint:
|
||||
return P.T * v[src, self]
|
||||
|
||||
|
||||
class SrcTDEM(Survey.BaseSrc):
|
||||
rxPair = RxTDEM
|
||||
radius = None
|
||||
|
||||
def getInitialFields(self, mesh):
|
||||
F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
|
||||
return F0
|
||||
|
||||
def getJs(self, mesh, time):
|
||||
return None
|
||||
|
||||
|
||||
class SrcTDEM_VMD_MVP(SrcTDEM):
|
||||
|
||||
def __init__(self,rxList,loc,waveformType="STEPOFF"):
|
||||
self.loc = loc
|
||||
self.waveformType = waveformType
|
||||
SrcTDEM.__init__(self,rxList)
|
||||
|
||||
def getInitialFields(self, mesh):
|
||||
"""Vertical magnetic dipole, magnetic vector potential"""
|
||||
if self.waveformType == "STEPOFF":
|
||||
print ">> Step waveform: Non-zero initial condition"
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
return {"b": mesh.edgeCurl*MVP}
|
||||
elif self.waveformType == "GENERAL":
|
||||
print ">> General waveform: Zero initial condition"
|
||||
return {"b": np.zeros(mesh.nF)}
|
||||
else:
|
||||
raise NotImplementedError("Only use STEPOFF or GENERAL")
|
||||
|
||||
def getMeS(self, mesh, MfMui):
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
|
||||
|
||||
|
||||
class SrcTDEM_CircularLoop_MVP(SrcTDEM):
|
||||
def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
|
||||
self.loc = loc
|
||||
self.radius = radius
|
||||
self.waveformType = waveformType
|
||||
SrcTDEM.__init__(self,rxList)
|
||||
|
||||
def getInitialFields(self, mesh):
|
||||
"""Circular Loop, magnetic vector potential"""
|
||||
if self.waveformType == "STEPOFF":
|
||||
print ">> Step waveform: Non-zero initial condition"
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
|
||||
else:
|
||||
raise Exception('Unknown mesh for CircularLoop')
|
||||
return {"b": mesh.edgeCurl*MVP}
|
||||
elif self.waveformType == "GENERAL":
|
||||
print ">> General waveform: Zero initial condition"
|
||||
return {"b": np.zeros(mesh.nF)}
|
||||
else:
|
||||
raise NotImplementedError("Only use STEPOFF or GENERAL")
|
||||
|
||||
def getMeS(self, mesh, MfMui):
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
|
||||
else:
|
||||
raise Exception('Unknown mesh for CircularLoop')
|
||||
return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
|
||||
|
||||
|
||||
class SurveyTDEM(Survey.BaseSurvey):
|
||||
"""
|
||||
docstring for SurveyTDEM
|
||||
"""
|
||||
srcPair = SrcTDEM
|
||||
|
||||
def __init__(self, srcList, **kwargs):
|
||||
# Sort these by frequency
|
||||
self.srcList = srcList
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
|
||||
def projectFields(self, u):
|
||||
data = Survey.Data(self)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.projectFields(src, self.mesh, self.prob.timeMesh, u)
|
||||
return data
|
||||
|
||||
def projectFieldsDeriv(self, u, v=None, adjoint=False):
|
||||
assert v is not None, 'v to multiply must be provided.'
|
||||
|
||||
if not adjoint:
|
||||
data = Survey.Data(self)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.projectFieldsDeriv(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.projectFieldsDeriv(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
|
||||
|
||||
|
||||
@@ -0,0 +1,3 @@
|
||||
from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
|
||||
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
|
||||
from TDEM_b import ProblemTDEM_b
|
||||
+12
-12
@@ -347,10 +347,10 @@ and
|
||||
|
||||
|
||||
|
||||
TDEM - B formulation
|
||||
====================
|
||||
TDEM Problem
|
||||
============
|
||||
|
||||
.. automodule:: SimPEG.EM.TDEM.TDEM_b
|
||||
.. automodule:: SimPEG.EM.TDEM.TDEM
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -359,7 +359,7 @@ TDEM - B formulation
|
||||
Field Storage
|
||||
=============
|
||||
|
||||
.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.FieldsTDEM
|
||||
.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.Fields
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -369,19 +369,19 @@ Field Storage
|
||||
TDEM Survey Classes
|
||||
===================
|
||||
|
||||
.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.SurveyTDEM
|
||||
.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.Survey
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
:inherited-members:
|
||||
|
||||
|
||||
Base Classes
|
||||
============
|
||||
.. Base Classes
|
||||
.. ============
|
||||
|
||||
.. automodule:: SimPEG.EM.TDEM.BaseTDEM
|
||||
:show-inheritance:
|
||||
:members:
|
||||
:undoc-members:
|
||||
:inherited-members:
|
||||
.. .. automodule:: SimPEG.EM.TDEM.BaseTDEM
|
||||
.. :show-inheritance:
|
||||
.. :members:
|
||||
.. :undoc-members:
|
||||
.. :inherited-members:
|
||||
|
||||
|
||||
Reference in New Issue
Block a user