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

use f where we are talking about fields

Browse Source
This commit is contained in:
Lindsey Heagy
2016-03-20 11:47:57 -07:00
parent 6a2caf5efc
commit 3d11431f2f
1 changed files with 56 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/FDEM/FDEM.py
+56 -56
View File
@@ -18,9 +18,9 @@ class BaseFDEMProblem(BaseEMProblem):
{\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
if using the E-B formulation (:code:`Problem_e`
or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
If we write Maxwell's equations in terms of
If we write Maxwell's equations in terms of
\\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
.. math ::
@@ -28,7 +28,7 @@ class BaseFDEMProblem(BaseEMProblem):
\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
\mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
"""
@@ -39,73 +39,73 @@ class BaseFDEMProblem(BaseEMProblem):
def fields(self, m=None):
"""
Solve the forward problem for the fields.
:param numpy.array m: inversion model (nP,)
:rtype numpy.array:
:return F: forward solution
:return f: forward solution
"""
self.curModel = m
F = self.fieldsPair(self.mesh, self.survey)
f = self.fieldsPair(self.mesh, self.survey)
for freq in self.survey.freqs:
A = self.getA(freq)
rhs = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
sol = Ainv * rhs
u = Ainv * rhs
Srcs = self.survey.getSrcByFreq(freq)
F[Srcs, self._solutionType] = sol
f[Srcs, self._solutionType] = u
Ainv.clean()
return F
return f
def Jvec(self, m, v, u=None):
def Jvec(self, m, v, f=None):
"""
Sensitivity times a vector.
:param numpy.array m: inversion model (nP,)
:param numpy.array v: vector which we take sensitivity product with (nP,)
:param SimPEG.EM.FDEM.Fields u: fields object
:param SimPEG.EM.FDEM.Fields u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
:return: Jv (ndata,)
"""
if u is None:
u = self.fields(m)
if f is None:
f = self.fields(m)
self.curModel = m
Jv = self.dataPair(self.survey)
for freq in self.survey.freqs:
A = self.getA(freq)
Ainv = self.Solver(A, **self.solverOpts)
A = self.getA(freq)
Ainv = self.Solver(A, **self.solverOpts) # create the concept of Ainv (actually a solve)
for src in self.survey.getSrcByFreq(freq):
u_src = u[src, self._solutionType]
dA_dm_v = self.getADeriv(freq, u_src, v)
dRHS_dm_v = self.getRHSDeriv(freq, src, v)
f_src = f[src, self._solutionType]
dA_dm_v = self.getADeriv(freq, f_src, v)
dRHS_dm_v = self.getRHSDeriv(freq, src, v)
du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
for rx in src.rxList:
df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None)
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
Jv[src, rx] = rx.evalDeriv(src, self.mesh, u, df_dm_v)
Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
Ainv.clean()
return Utils.mkvc(Jv)
def Jtvec(self, m, v, u=None):
def Jtvec(self, m, v, f=None):
"""
Sensitivity transpose times a vector
:param numpy.array m: inversion model (nP,)
:param numpy.array v: vector which we take adjoint product with (nP,)
:param SimPEG.EM.FDEM.Fields u: fields object
:param SimPEG.EM.FDEM.Fields u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
:return: Jv (ndata,)
"""
if u is None:
u = self.fields(m)
if f is None:
f = self.fields(m)
self.curModel = m
@@ -120,17 +120,17 @@ class BaseFDEMProblem(BaseEMProblem):
ATinv = self.Solver(AT, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
u_src = u[src, self._solutionType]
f_src = f[src, self._solutionType]
for rx in src.rxList:
PTv = rx.evalDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
ATinvdf_duT = ATinv * df_duT
dA_dmT = self.getADeriv(freq, u_src, ATinvdf_duT, adjoint=True)
dA_dmT = self.getADeriv(freq, f_src, ATinvdf_duT, adjoint=True)
dRHS_dmT = self.getRHSDeriv(freq, src, ATinvdf_duT, adjoint=True)
du_dmT = -dA_dmT + dRHS_dmT
@@ -144,7 +144,7 @@ class BaseFDEMProblem(BaseEMProblem):
Jtv += - np.array(df_dmT, dtype=complex).real
else:
raise Exception('Must be real or imag')
ATinv.clean()
return Utils.mkvc(Jtv)
@@ -154,7 +154,7 @@ class BaseFDEMProblem(BaseEMProblem):
Evaluates the sources for a given frequency and puts them in matrix form
:param float freq: Frequency
:rtype: (numpy.ndarray, numpy.ndarray)
:rtype: (numpy.ndarray, numpy.ndarray)
:return: S_m, S_e (nE or nF, nSrc)
"""
Srcs = self.survey.getSrcByFreq(freq)
@@ -207,7 +207,7 @@ class Problem_e(BaseFDEMProblem):
def getA(self, freq):
"""
System matrix
.. math ::
\mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}}
@@ -230,12 +230,12 @@ class Problem_e(BaseFDEMProblem):
.. math ::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = i \omega \\frac{d \mathbf{M^e_{\sigma}}\mathbf{v} }{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
dsig_dm = self.curModel.sigmaDeriv
@@ -248,13 +248,13 @@ class Problem_e(BaseFDEMProblem):
def getRHS(self, freq):
"""
Right hand side for the system
Right hand side for the system
.. math ::
\mathbf{RHS} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
@@ -266,7 +266,7 @@ class Problem_e(BaseFDEMProblem):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
@@ -346,12 +346,12 @@ class Problem_b(BaseFDEMProblem):
.. math ::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \\frac{\mathbf{M^e_{\sigma}} \mathbf{v}}{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
MfMui = self.MfMui
@@ -373,13 +373,13 @@ class Problem_b(BaseFDEMProblem):
def getRHS(self, freq):
"""
Right hand side for the system
Right hand side for the system
.. math ::
\mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
@@ -497,12 +497,12 @@ class Problem_j(BaseFDEMProblem):
\\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \mathbf{M^e_{mu^{-1}}} \mathbf{C^{\\top}} \\frac{d \mathbf{M^f_{\sigma^{-1}}}\mathbf{v} }{d \mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nF,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
MeMuI = self.MeMuI
@@ -522,7 +522,7 @@ class Problem_j(BaseFDEMProblem):
def getRHS(self, freq):
"""
Right hand side for the system
Right hand side for the system
.. math ::
@@ -546,7 +546,7 @@ class Problem_j(BaseFDEMProblem):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
@@ -626,12 +626,12 @@ class Problem_h(BaseFDEMProblem):
.. math::
\\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C}^{\\top}\\frac{d \mathbf{M^f_{\\rho}}\mathbf{v} }{d\mathbf{m}}
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param float freq: frequency
:param numpy.ndarray u: solution vector (nE,)
:param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
:return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
"""
MeMu = self.MeMu
@@ -644,14 +644,14 @@ class Problem_h(BaseFDEMProblem):
def getRHS(self, freq):
"""
Right hand side for the system
Right hand side for the system
.. math ::
\mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
@@ -663,7 +663,7 @@ class Problem_h(BaseFDEMProblem):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source