mirror of
https://github.com/wassname/simpeg.git
synced 2026-08-15 12:54:52 +08:00
458 lines
16 KiB
Python
458 lines
16 KiB
Python
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)
|
|
|
|
|
|
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)
|
|
|
|
# set initial fields
|
|
F[:,self._fieldType+'Solution',0] = self.getInitialFields()
|
|
|
|
# timestep to solve forward
|
|
Ainv = None
|
|
for tInd, dt in enumerate(self.timeSteps):
|
|
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()
|
|
Ainv = None
|
|
|
|
if Ainv is None:
|
|
A = self.getAdiag(tInd)
|
|
if self.verbose: print 'Factoring... (dt = %e)'%dt
|
|
Ainv = self.Solver(A, **self.solverOpts)
|
|
if self.verbose: print 'Done'
|
|
|
|
rhs = self.getRHS(tInd+1) # this is on the nodes of the time mesh
|
|
Asubdiag = self.getAsubdiag(tInd)
|
|
|
|
if self.verbose: print ' Solving... (tInd = %d)'%tInd+1
|
|
sol = Ainv * (rhs - Asubdiag * F[:,self._fieldType+'Solution',tInd]) # taking a step
|
|
|
|
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):
|
|
|
|
if u is None:
|
|
u = self.fields(m)
|
|
|
|
ftype = self._fieldType + 'Solution' # the thing we solved for
|
|
self.curModel = m
|
|
|
|
Jv = self.dataPair(self.survey)
|
|
|
|
# mat to store previous time-step's solution deriv times a vector for each source
|
|
# size: nu x nSrc
|
|
dun_dm_v = self.getInitialFieldsDeriv(v) # can over-write this at each timestep
|
|
|
|
#
|
|
df_dm_v = Fields_Derivs(self.mesh, self.survey) # store the field derivs we need to project to calc full deriv
|
|
|
|
Adiaginv = None
|
|
|
|
for tInd, dt in zip(range(self.nT+1), self.timeSteps):
|
|
if Adiaginv 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
|
|
Adiaginv.clean()
|
|
Adiaginv = None
|
|
|
|
if Adiaginv is None:
|
|
A = self.getAdiag(tInd)
|
|
Adiaginv = self.Solver(A, **self.solverOpts)
|
|
|
|
Asubdiag = self.getAsubdiag(tInd)
|
|
|
|
for i, src in enumerate(self.survey.srcList):
|
|
|
|
# here, we are lagging by a timestep, so filling in as we go
|
|
for projField in set([rx.projField for rx in src.rxList]):
|
|
df_dmFun = getattr(u, '_%sDeriv'%projField, None)
|
|
# df_dm_v is dense, but we only need the times at (rx.P.T * ones > 0)
|
|
# This should be called rx.footprint
|
|
df_dm_v[src, '%sDeriv'%projField , tInd] = df_dmFun(tInd, src, dun_dm_v[:,i], v)
|
|
|
|
un_src = u[src,ftype,tInd+1]
|
|
|
|
dA_dm_v = self.getAdiagDeriv(tInd, un_src, v) # cell centered on time mesh
|
|
dRHS_dm_v = self.getRHSDeriv(tInd+1, src, v) # on nodes of time mesh
|
|
# dAsubdiag_dm_v = 0
|
|
|
|
JRHS = dRHS_dm_v - dA_dm_v # - dAsubdiag_dm_v (which is zero)
|
|
|
|
# step in time and overwrite
|
|
if tInd != len(self.timeSteps+1):
|
|
dun_dm_v[:,i] = Adiaginv * (JRHS - Asubdiag * dun_dm_v[:,i])
|
|
|
|
for src in self.survey.srcList:
|
|
for rx in src.rxList:
|
|
Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, df_dm_v)
|
|
|
|
Adiaginv.clean()
|
|
return Utils.mkvc(Jv)
|
|
|
|
|
|
def Jtvec(self, m, v, u=None):
|
|
|
|
if u is None:
|
|
u = self.fields(m)
|
|
|
|
self.curModel = m
|
|
ftype = self._fieldType + 'Solution' # the thing we solved for
|
|
|
|
# Ensure v is a data object.
|
|
if not isinstance(v, self.dataPair):
|
|
v = self.dataPair(self.survey, v)
|
|
|
|
# TODO: make this general
|
|
# if self._fieldType is 'b':
|
|
# dun_dmT_v = np.zeros((len(m), self.survey.nSrc))
|
|
|
|
|
|
PT_v = Fields_Derivs(self.mesh, self.survey) #PT_v is a fields object
|
|
|
|
# TODO: This will only work for b formulation right now b/c of the mesh.nF
|
|
df_duT_v = np.zeros((self.mesh.nF,self.nT+1))
|
|
ATinv_df_duT_v = np.zeros((self.mesh.nF,self.nT))
|
|
|
|
for src in self.survey.srcList:
|
|
# for rx in src.rxList:
|
|
for projField in set([rx.projField for rx in src.rxList]):
|
|
PT_v[src,'%sDeriv'%projField, :] = rx.evalDeriv(src, self.mesh, self.timeMesh, v, adjoint = True) # All the fields for a given src, reciever.
|
|
|
|
df_duTFun = getattr(u, '_%sDeriv'%projField, None)
|
|
|
|
# TODO: don't need to recompute df_dmT_v every time... only need it once
|
|
df_duT_v_cur, df_dmT_v = df_duTFun(None, src, None, PT_v[src,'%sDeriv'%projField,:], adjoint=True) # this seems odd
|
|
|
|
df_duT_v = df_duT_v + df_duT_v_cur
|
|
|
|
AdiagTinv = None
|
|
JTv = df_dmT_v
|
|
|
|
for tInd in reversed(range(self.nT)): #enumerate(reversed(list(self.timeSteps))):
|
|
if AdiagTinv is not None and (tInd < self.nT and self.timeSteps[tInd] != self.timeSteps[tInd + 1]):# keep factors if dt is the same as previous step b/c A will be the same
|
|
AdiagTinv.clean()
|
|
AdiagTinv = None
|
|
|
|
if AdiagTinv is None:
|
|
Adiag = self.getAdiag(tInd)
|
|
AdiagTinv = self.Solver(Adiag.T, **self.solverOpts)
|
|
|
|
Asubdiag = self.getAsubdiag(tInd)
|
|
|
|
# solve against df_duT_v
|
|
|
|
if tInd < self.nT:
|
|
# print Utils.mkvc(AdiagTinv * df_duT_v[:,tInd],2).shape, ATinv_df_duT_v[:,tInd].shape
|
|
ATinv_df_duT_v[:,tInd] = AdiagTinv * df_duT_v[:,tInd]
|
|
else:
|
|
ATinv_df_duT_v[:,tInd] = AdiagTinv * (df_duT_v[:,tInd+1] - Asubdiag.T * df_duT_v[:,tInd])
|
|
|
|
un_src = u[src,ftype,tInd]
|
|
dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[:,tInd], adjoint=True) # cell centered on time mesh
|
|
|
|
dRHST_dm_v = self.getRHSDeriv(tInd, src, ATinv_df_duT_v[:,tInd], adjoint=True) # on nodes of time mesh
|
|
# dAsubdiag_dm_v = 0
|
|
|
|
JTv = JTv + (-dAT_dm_v + dRHST_dm_v)
|
|
|
|
return JTv
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
# for i, src in enumerate(self.survey.srcList):
|
|
|
|
# un_src = u[src,ftype,tInd+1] # fields for this source at tInd
|
|
|
|
# for rx in src.rxList:
|
|
|
|
# df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
|
|
# df_duT_v, df_dmT_v = df_duTFun(tInd, src, None, PT_v[src,'%sDeriv'%rx.projField,tInd], adjoint=True)
|
|
|
|
# ATinv_df_duT_v = AdiagTinv * df_duT_v
|
|
|
|
# dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v, adjoint=True)
|
|
# dRHST_dm_v = self.getRHSDeriv(tInd, src, ATinv_df_duT_v)
|
|
# # dAsubdiagT_dm_v = 0
|
|
|
|
# print dAT_dm_v.shape, Asubdiag.shape, ATinv_df_duT_v.shape, dun_dmT_v.shape,
|
|
# dun_dmT_v[:,i] = (dRHST_dm_v - dAT_dm_v - Asubdiag.T*dun_dmT_v[:,i])
|
|
|
|
# # rhsT_v = self.getJRHS(tInd, src, u_src, ATinv_df_duT_v, dun_dmT_v[:,i], adjoint = True)
|
|
|
|
# JTv = JTv + dun_dmT_v
|
|
|
|
# return Utils.mkvc(JTv)
|
|
|
|
|
|
|
|
# def getJRHS(self, tInd, src, u, v, adjoint = False):
|
|
|
|
# dA_dm = self.getADeriv(tInd, u, v, adjoint)
|
|
# dRHS_dm = self.getRHSDeriv(tInd, src, v, adjoint)
|
|
|
|
# b = - dA_dm + dRHS_dm
|
|
|
|
# return b
|
|
|
|
|
|
def getSourceTerm(self, tInd):
|
|
|
|
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
|
|
|
|
def getInitialFields(self):
|
|
|
|
Srcs = self.survey.srcList
|
|
|
|
if self._fieldType is 'b' or self._fieldType is 'j':
|
|
ifields = np.zeros((self.mesh.nF, len(Srcs)))
|
|
elif self._fieldType is 'e' or self._fieldType is 'h':
|
|
ifields = np.zeros((self.mesh.nE, len(Srcs)))
|
|
|
|
for i,src in enumerate(Srcs):
|
|
ifields[:,i] = ifields[:,i] + getattr(src, '%sInitial'%self._fieldType, None)(self)
|
|
|
|
return ifields
|
|
|
|
def getInitialFieldsDeriv(self, v, adjoint=False):
|
|
|
|
Srcs = self.survey.srcList
|
|
|
|
if self._fieldType is 'b' or self._fieldType is 'j':
|
|
ifieldsDeriv = np.zeros((self.mesh.nF, len(Srcs)))
|
|
elif self._fieldType is 'e' or self._fieldType is 'h':
|
|
ifieldsDeriv = np.zeros((self.mesh.nE, len(Srcs)))
|
|
|
|
for i,src in enumerate(Srcs):
|
|
ifieldsDeriv[:,i] = ifieldsDeriv[:,i] + getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)
|
|
|
|
return ifieldsDeriv
|
|
|
|
|
|
##########################################################################################
|
|
################################ 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 getAdiag(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 = 1./dt * I + ( C * ( MeSigmaI * (C.T * MfMui ) ) )
|
|
|
|
if self._makeASymmetric is True:
|
|
return MfMui.T * A
|
|
return A
|
|
|
|
def getAdiagDeriv(self, tInd, u, v, adjoint=False):
|
|
C = self.mesh.edgeCurl
|
|
MeSigmaIDeriv = lambda x: self.MeSigmaIDeriv(x)
|
|
MfMui = self.MfMui
|
|
|
|
if adjoint:
|
|
if self._makeASymmetric is True:
|
|
v = MfMui * v
|
|
return MeSigmaIDeriv(C.T * ( MfMui * u )).T * ( C.T * v )
|
|
|
|
ADeriv = ( C * ( MeSigmaIDeriv(C.T * ( MfMui * u )) * v ) )
|
|
if self._makeASymmetric is True:
|
|
return MfMui.T * ADeriv
|
|
return ADeriv
|
|
|
|
|
|
def getAsubdiag(self, tInd):
|
|
|
|
dt = self.timeSteps[tInd]
|
|
MfMui = self.MfMui
|
|
Asubdiag = - 1./dt * sp.eye(self.mesh.nF)
|
|
|
|
if self._makeASymmetric is True:
|
|
return MfMui.T * Asubdiag
|
|
|
|
return Asubdiag
|
|
|
|
|
|
|
|
def getRHS(self, tInd):
|
|
C = self.mesh.edgeCurl
|
|
MeSigmaI = self.MeSigmaI
|
|
MfMui = self.MfMui
|
|
|
|
S_m, S_e = self.getSourceTerm(tInd)
|
|
|
|
# B_n = np.c_[[F[src,'bSolution',tInd] for src in self.survey.srcList]]
|
|
# if B_n.shape[0] is not 1:
|
|
# raise NotImplementedError('getRHS not implemented for this shape of B_n')
|
|
|
|
rhs = (C * (MeSigmaI * S_e) + S_m) # + 1./dt * B_n[:,:,0].T
|
|
if self._makeASymmetric is True:
|
|
return MfMui.T * rhs
|
|
return rhs
|
|
|
|
def getRHSDeriv(self, tInd, src, v, adjoint=False):
|
|
|
|
C = self.mesh.edgeCurl
|
|
MeSigmaI = self.MeSigmaI
|
|
MeSigmaIDeriv = lambda u: self.MeSigmaIDeriv(u)
|
|
MfMui = self.MfMui
|
|
|
|
_, S_e = src.eval(tInd, self)
|
|
S_mDeriv, S_eDeriv = src.evalDeriv(self.times[tInd], self, adjoint=adjoint) # I think this is tInd+1 ?
|
|
|
|
if adjoint:
|
|
if self._makeASymmetric is True:
|
|
v = self.MfMui * v
|
|
if isinstance(S_e, Utils.Zero):
|
|
MeSigmaIDerivT_v = Utils.Zero()
|
|
else:
|
|
MeSigmaIDerivT_v = MeSigmaIDeriv(S_e).T * v
|
|
RHSDeriv = MeSigmaIDerivT_v + S_eDeriv( MeSigmaI.T * ( C.T * v ) ) + S_mDeriv(v) #+ dbn_dm_v / dt #this will be given the transposed version
|
|
return RHSDeriv
|
|
|
|
if isinstance(S_e, Utils.Zero):
|
|
MeSigmaIDeriv_v = Utils.Zero()
|
|
else:
|
|
MeSigmaIDeriv_v = MeSigmaIDeriv(S_e) * v
|
|
|
|
RHSDeriv = (C * (MeSigmaIDeriv_v + MeSigmaI * S_eDeriv(v) + S_mDeriv(v))) #+ dbn_dm_v / dt
|
|
|
|
if self._makeASymmetric is True:
|
|
return self.MfMui.T * RHSDeriv
|
|
return RHSDeriv
|
|
|
|
|
|
@Utils.timeIt
|
|
def getJdiags(self, tInd, adjoint = False):
|
|
# The matrix that we are computing has the form:
|
|
#
|
|
# - - - - - -
|
|
# | Adiag | | uderiv1 | | b1 |
|
|
# | Asub Adiag | | uderiv2 | | b2 |
|
|
# | Asub Adiag | | uderiv3 | = | b3 |
|
|
# | ... ... | | ... | | .. |
|
|
# | Asub Adiag | | uderivn | | bn |
|
|
# - - - - - -
|
|
|
|
if adjoint:
|
|
raise NotImplementedError
|
|
|
|
dt = self.timeSteps[tInd]
|
|
|
|
Adiag = self.getA(tInd)
|
|
Asub = - 1./dt * Utils.speye(self.mesh.nF)
|
|
|
|
if self._makeASymmetric:
|
|
Asub = self.MfMui.T * Asub
|
|
|
|
return Adiag, Asub
|
|
|
|
|
|
|
|
|
|
|