mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-05 22:18:42 +08:00
Added derivative support for the 3D problem.
This commit is contained in:
GudniRos
+195
-19
@@ -45,7 +45,7 @@ class FieldsMT_1D(FieldsMT):
|
||||
return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
|
||||
|
||||
def _eDeriv_u(self, src, v, adjoint = False):
|
||||
return None
|
||||
return v
|
||||
|
||||
def _eDeriv_m(self, src, v, adjoint = False):
|
||||
# assuming primary does not depend on the model
|
||||
@@ -124,23 +124,199 @@ class FieldsMT_3D(FieldsMT):
|
||||
"""
|
||||
Fields storage for the 3D MT solution.
|
||||
"""
|
||||
knownFields = {'e_px':'E','e_py':'E','b_px':'F','b_py':'F'}
|
||||
aliasFields = { }
|
||||
# 'e_1d' : ['e_1dSolution','F','_e'],
|
||||
# 'e_1dPrimary' : ['e_1dSolution','F','_ePrimary'],
|
||||
# 'e_1dSecondary' : ['e_1dSolution','F','_eSecondary'],
|
||||
# 'b_1d' : ['e_1dSolution','E','_b'],
|
||||
# 'b_1dPrimary' : ['e_1dSolution','E','_bPrimary'],
|
||||
# 'b_1dSecondary' : ['e_1dSolution','E','_bSecondary']
|
||||
# }
|
||||
# Define the known the alias fields
|
||||
# Assume that the solution of e on the E.
|
||||
## NOTE: Need to make this more general, to allow for other solutions formats.
|
||||
knownFields = {'e_pxSolution':'E','e_pySoluiton':'E'}
|
||||
aliasFields = {
|
||||
'e_px' : ['e_pxSolution','E','_e_px'],
|
||||
'e_pxPrimary' : ['e_pxSolution','E','_e_pxPrimary'],
|
||||
'e_pxSecondary' : ['e_pxSolution','E','_e_pxSecondary'],
|
||||
'b_px' : ['e_pxSolution','F','_b_px'],
|
||||
'b_pxPrimary' : ['e_pxSolution','F','_b_pxPrimary'],
|
||||
'b_pxSecondary' : ['e_pxSolution','F','_b_pxSecondary']
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsMT.__init__(self,mesh,survey,**kwargs)
|
||||
|
||||
# knownFields = {'e_pxSolution':'E','e_pySoluiton':'E'}
|
||||
# aliasFields = {
|
||||
# 'e_px' : ['e_pxSolution','E','_epx'],
|
||||
# 'e_pxPrimary' : ['e_pxSolution','E','_epxPrimary'],
|
||||
# 'e_pxSecondary' : ['e_pxSolution','E','_epxSecondary'],
|
||||
# 'b_px' : ['e_pxSolution','F','_bpx'],
|
||||
# 'b_pxPrimary' : ['e_pxSolution','F','_bpxPrimary'],
|
||||
# 'b_pxSecondary' : ['e_pxSolution','F','_bpxSecondary']
|
||||
# }
|
||||
def _e_pxPrimary(self, e_pxSolution, srcList):
|
||||
e_pxPrimary = np.zeros_like(e_pxSolution)
|
||||
for i, src in enumerate(srcList):
|
||||
ep = src.ePrimary(self.survey.prob)
|
||||
if ep is not None:
|
||||
e_pxPrimary[:,i] = ep[:,0]
|
||||
return e_pxPrimary
|
||||
|
||||
def _e_pyPrimary(self, e_pySolution, srcList):
|
||||
e_pyPrimary = np.zeros_like(e_pySolution)
|
||||
for i, src in enumerate(srcList):
|
||||
ep = src.ePrimary(self.survey.prob)
|
||||
if ep is not None:
|
||||
e_pyPrimary[:,i] = ep[:,1]
|
||||
return e_pyPrimary
|
||||
|
||||
def _e_pxSecondary(self, e_pxSolution, srcList):
|
||||
return e_pxSolution
|
||||
|
||||
def _e_pySecondary(self, e_pySolution, srcList):
|
||||
return e_pySolution
|
||||
|
||||
def _e_px(self, e_pxSolution, srcList):
|
||||
return self._e_pxPrimary(e_pxSolution,srcList) + self._e_pxSecondary(e_pxSolution,srcList)
|
||||
|
||||
def _e_py(self, e_pySolution, srcList):
|
||||
return self._e_pyPrimary(e_pySolution,srcList) + self._e_pySecondary(e_pySolution,srcList)
|
||||
|
||||
def _e_pxDeriv_u(self, src, v, adjoint = False):
|
||||
return v
|
||||
|
||||
def _e_pyDeriv_u(self, src, v, adjoint = False):
|
||||
return v
|
||||
|
||||
def _e_pxDeriv_m(self, src, v, adjoint = False):
|
||||
# assuming primary does not depend on the model
|
||||
return None
|
||||
def _e_pyDeriv_m(self, src, v, adjoint = False):
|
||||
# assuming primary does not depend on the model
|
||||
return None
|
||||
|
||||
def _b_pxPrimary(self, e_pxSolution, srcList):
|
||||
b_pxPrimary = np.zeros([self.survey.mesh.nE,e_pxSolution.shape[1]], dtype = complex)
|
||||
for i, src in enumerate(srcList):
|
||||
bp = src.bPrimary(self.survey.prob)
|
||||
if bp is not None:
|
||||
b_pxPrimary[:,i] += bp[:,-1]
|
||||
return b_pxPrimary
|
||||
|
||||
def _b_pyPrimary(self, e_pySolution, srcList):
|
||||
b_pyPrimary = np.zeros([self.survey.mesh.nE,e_pySolution.shape[1]], dtype = complex)
|
||||
for i, src in enumerate(srcList):
|
||||
bp = src.bPrimary(self.survey.prob)
|
||||
if bp is not None:
|
||||
b_pyPrimary[:,i] += bp[:,-1]
|
||||
return b_pyPrimary
|
||||
|
||||
def _b_pxSecondary(self, e_pxSolution, srcList):
|
||||
C = self.mesh.edgeCurl
|
||||
b = (C * e_pxSolution)
|
||||
for i, src in enumerate(srcList):
|
||||
b[:,i] *= - 1./(1j*omega(src.freq))
|
||||
# There is no magnetic source in the MT problem
|
||||
# S_m, _ = src.eval(self.survey.prob)
|
||||
# if S_m is not None:
|
||||
# b[:,i] += 1./(1j*omega(src.freq)) * S_m
|
||||
return b
|
||||
|
||||
def _b_pySecondary(self, e_pySolution, srcList):
|
||||
C = self.mesh.edgeCurl
|
||||
b = (C * e_pySolution)
|
||||
for i, src in enumerate(srcList):
|
||||
b[:,i] *= - 1./(1j*omega(src.freq))
|
||||
# There is no magnetic source in the MT problem
|
||||
# S_m, _ = src.eval(self.survey.prob)
|
||||
# if S_m is not None:
|
||||
# b[:,i] += 1./(1j*omega(src.freq)) * S_m
|
||||
return b
|
||||
|
||||
def _b_px(self, eSolution, srcList):
|
||||
return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
|
||||
|
||||
def _b_pxSecondaryDeriv_u(self, src, v, adjoint = False):
|
||||
C = self.mesh.edgeCurl
|
||||
if adjoint:
|
||||
return - 1./(1j*omega(src.freq)) * (C.T * v)
|
||||
return - 1./(1j*omega(src.freq)) * (C * v)
|
||||
|
||||
def _b_py(self, eSolution, srcList):
|
||||
return self._b_pyPrimary(eSolution, srcList) + self._b_pySecondary(eSolution, srcList)
|
||||
|
||||
def _b_pxSecondaryDeriv_u(self, src, v, adjoint = False):
|
||||
C = self.mesh.edgeCurl
|
||||
if adjoint:
|
||||
return - 1./(1j*omega(src.freq)) * (C.T * v)
|
||||
return - 1./(1j*omega(src.freq)) * (C * v)
|
||||
|
||||
def _b_pySecondaryDeriv_m(self, src, v, adjoint = False):
|
||||
# Doesn't depend on m
|
||||
# _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
|
||||
# S_eDeriv = S_eDeriv(v)
|
||||
# if S_eDeriv is not None:
|
||||
# return 1./(1j * omega(src.freq)) * S_eDeriv
|
||||
return None
|
||||
|
||||
def _b_pxSecondaryDeriv_m(self, src, v, adjoint = False):
|
||||
# Doesn't depend on m
|
||||
# _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
|
||||
# S_eDeriv = S_eDeriv(v)
|
||||
# if S_eDeriv is not None:
|
||||
# return 1./(1j * omega(src.freq)) * S_eDeriv
|
||||
return None
|
||||
|
||||
def _b_pxDeriv_u(self, src, v, adjoint=False):
|
||||
# Primary does not depend on u
|
||||
return self._b_pxSecondaryDeriv_u(src, v, adjoint)
|
||||
|
||||
def _b_pyDeriv_u(self, src, v, adjoint=False):
|
||||
# Primary does not depend on u
|
||||
return self._b_pySecondaryDeriv_u(src, v, adjoint)
|
||||
|
||||
def _b_pxDeriv_m(self, src, v, adjoint=False):
|
||||
# Assuming the primary does not depend on the model
|
||||
return self._b_pxSecondaryDeriv_m(src, v, adjoint)
|
||||
|
||||
def _b_pyDeriv_m(self, src, v, adjoint=False):
|
||||
# Assuming the primary does not depend on the model
|
||||
return self._b_pySecondaryDeriv_m(src, v, adjoint)
|
||||
|
||||
def _f_pxDeriv_u(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the fields object wrt u.
|
||||
|
||||
:param MTsrc src: MT source
|
||||
:param numpy.ndarray v: random vector of f_sol.size
|
||||
This function stacks the fields derivatives appropriately
|
||||
|
||||
return a vector of size (nreEle+nrbEle)
|
||||
"""
|
||||
|
||||
de_du = v #Utils.spdiag(np.ones((self.nF,)))
|
||||
db_du = self._b_pxDeriv_u(src, v, adjoint)
|
||||
# Return the stack
|
||||
# This doesn't work...
|
||||
return np.vstack((de_du,db_du))
|
||||
|
||||
def _f_pyDeriv_u(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the fields object wrt u.
|
||||
|
||||
:param MTsrc src: MT source
|
||||
:param numpy.ndarray v: random vector of f_sol.size
|
||||
This function stacks the fields derivatives appropriately
|
||||
|
||||
return a vector of size (nreEle+nrbEle)
|
||||
"""
|
||||
|
||||
de_du = v #Utils.spdiag(np.ones((self.nF,)))
|
||||
db_du = self._b_pyDeriv_u(src, v, adjoint)
|
||||
# Return the stack
|
||||
# This doesn't work...
|
||||
return np.vstack((de_du,db_du))
|
||||
|
||||
def _f_pxDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the fields object wrt m.
|
||||
|
||||
This function stacks the fields derivatives appropriately
|
||||
"""
|
||||
# The fields have no dependance to the model.
|
||||
return None
|
||||
|
||||
def _f_pyDeriv_m(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the fields object wrt m.
|
||||
|
||||
This function stacks the fields derivatives appropriately
|
||||
"""
|
||||
# The fields have no dependance to the model.
|
||||
return None
|
||||
@@ -104,16 +104,13 @@ class eForm_ps(BaseMTProblem):
|
||||
|
||||
# Store the fields
|
||||
Src = self.survey.getSrcByFreq(freq)[0]
|
||||
# Calculate total e
|
||||
|
||||
e = Src.ePrimary(self) + e_s
|
||||
# Store the fieldss
|
||||
F[Src, 'e_px'] = e[:,0]
|
||||
F[Src, 'e_py'] = e[:,1]
|
||||
F[Src, 'e_pxSolution'] = e_s[:,0]
|
||||
F[Src, 'e_pySolution'] = e_s[:,1]
|
||||
# Note curl e = -iwb so b = -curl/iw
|
||||
b = -( self.mesh.edgeCurl * e )/( 1j*omega(freq) )
|
||||
F[Src, 'b_px'] = b[:,0]
|
||||
F[Src, 'b_py'] = b[:,1]
|
||||
# b = -( self.mesh.edgeCurl * e )/( 1j*omega(freq) )
|
||||
# F[Src, 'b_px'] = b[:,0]
|
||||
# F[Src, 'b_py'] = b[:,1]
|
||||
if self.verbose:
|
||||
print 'Ran for {:f} seconds'.format(time.time()-startTime)
|
||||
sys.stdout.flush()
|
||||
|
||||
+101
-2
@@ -175,7 +175,57 @@ class RxMT(Survey.BaseRx):
|
||||
elif self.projType is 'Z2D':
|
||||
raise NotImplementedError('Has not be implement for 2D impedance tensor')
|
||||
elif self.projType is 'Z3D':
|
||||
raise NotImplementedError('Has not be implement for full 3D impedance tensor')
|
||||
if self.locs.ndim == 3:
|
||||
eFLocs = self.locs[:,:,0]
|
||||
bFLocs = self.locs[:,:,1]
|
||||
else:
|
||||
eFLocs = self.locs
|
||||
bFLocs = self.locs
|
||||
# Get the projection
|
||||
Pex = mesh.getInterpolationMat(eFLocs,'Ex')
|
||||
Pey = mesh.getInterpolationMat(eFLocs,'Ey')
|
||||
Pbx = mesh.getInterpolationMat(bFLocs,'Fx')
|
||||
Pby = mesh.getInterpolationMat(bFLocs,'Fy')
|
||||
# Get the fields at location
|
||||
# px: x-polaration and py: y-polaration.
|
||||
ex_px = Pex*f[src,'e_px']
|
||||
ey_px = Pey*f[src,'e_px']
|
||||
ex_py = Pex*f[src,'e_py']
|
||||
ey_py = Pey*f[src,'e_py']
|
||||
hx_px = Pbx*f[src,'b_px']/mu_0
|
||||
hy_px = Pby*f[src,'b_px']/mu_0
|
||||
hx_py = Pbx*f[src,'b_py']/mu_0
|
||||
hy_py = Pby*f[src,'b_py']/mu_0
|
||||
# Derivatives as lambda functions
|
||||
ex_px_u = lambda vec: Pex*f._e_pxDeriv_u(src,vec)
|
||||
ey_px_u = lambda vec: Pey*f._e_pxDeriv_u(src,vec)
|
||||
ex_py_u = lambda vec: Pex*f._e_pyDeriv_u(src,vec)
|
||||
ey_py_u = lambda vec: Pey*f._e_pyDeriv_u(src,vec)
|
||||
hx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec)/mu_0
|
||||
hy_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec)/mu_0
|
||||
hx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec)/mu_0
|
||||
hy_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec)/mu_0
|
||||
|
||||
# Define the components of the derivative
|
||||
Hd = (hx_px*hy_py - hx_py*hy_px)
|
||||
Hd_uV = hx_px_u(hy_py)*v - hx_py*hy_px_u(v) + hx_px*hy_py_u(v) - hx_py_u(hy_px)
|
||||
|
||||
if 'zxx' in self.rxType:
|
||||
Zij = ( ex_px*hy_py - ex_py*hy_px)/Hd
|
||||
ZijN_uV = ex_px_u(hy_py)*v - ex_py*hy_px_u(v) + ex_px*hy_py_u(v) - ex_py_u(hy_px)*v
|
||||
elif 'zxy' in self.rxType:
|
||||
Zij = (-ex_px*hx_py + ex_py*hx_px)/Hd
|
||||
ZijN_uV = -ex_px_u(hx_py)*v + ex_py*hx_px_u(v) - ex_px*hx_py_u(v) + ex_py_u(hx_px)*v
|
||||
elif 'zyx' in self.rxType:
|
||||
Zij = ( ey_px*hy_py - ey_py*hy_px)/Hd
|
||||
ZijN_uV = ey_px_u(hy_py)*v - ey_py*hy_px_u(v) + ey_px*hy_py_u(v) - ey_py_u(hy_px)*v
|
||||
elif 'zyy' in self.rxType:
|
||||
Zij = (-ey_px*hx_py + ey_py*hx_px)/Hd
|
||||
ZijN_uV = -ey_px_u(hx_py)*v + ey_py*hx_px_u(v) - ey_px*hx_py_u(v) +ey_py_u(hx_px)*v
|
||||
|
||||
# Calculate the complex derivative
|
||||
PDeriv_complex = ZijN_uV/Hd + Zij * (Hd_uV/Hd)
|
||||
|
||||
# Extract the real number for the real/imag components.
|
||||
Pv = np.array(getattr(PDeriv_complex, real_or_imag))
|
||||
elif adjoint:
|
||||
@@ -192,7 +242,56 @@ class RxMT(Survey.BaseRx):
|
||||
elif self.projType is 'Z2D':
|
||||
raise NotImplementedError('Has not be implement for 2D impedance tensor')
|
||||
elif self.projType is 'Z3D':
|
||||
raise NotImplementedError('Has not be implement for full 3D impedance tensor')
|
||||
if self.locs.ndim == 3:
|
||||
eFLocs = self.locs[:,:,0]
|
||||
bFLocs = self.locs[:,:,1]
|
||||
else:
|
||||
eFLocs = self.locs
|
||||
bFLocs = self.locs
|
||||
# Get the projection
|
||||
Pex = mesh.getInterpolationMat(eFLocs,'Ex')
|
||||
Pey = mesh.getInterpolationMat(eFLocs,'Ey')
|
||||
Pbx = mesh.getInterpolationMat(bFLocs,'Fx')
|
||||
Pby = mesh.getInterpolationMat(bFLocs,'Fy')
|
||||
# Get the fields at location
|
||||
# px: x-polaration and py: y-polaration.
|
||||
ex_px = Pex*f[src,'e_px']
|
||||
ey_px = Pey*f[src,'e_px']
|
||||
ex_py = Pex*f[src,'e_py']
|
||||
ey_py = Pey*f[src,'e_py']
|
||||
hx_px = Pbx*f[src,'b_px']/mu_0
|
||||
hy_px = Pby*f[src,'b_px']/mu_0
|
||||
hx_py = Pbx*f[src,'b_py']/mu_0
|
||||
hy_py = Pby*f[src,'b_py']/mu_0
|
||||
# Derivatives as lambda functions
|
||||
ex_px_u = lambda vec: Pex*f._e_pxDeriv_u(src,vec)
|
||||
ey_px_u = lambda vec: Pey*f._e_pxDeriv_u(src,vec)
|
||||
ex_py_u = lambda vec: Pex*f._e_pyDeriv_u(src,vec)
|
||||
ey_py_u = lambda vec: Pey*f._e_pyDeriv_u(src,vec)
|
||||
hx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec)/mu_0
|
||||
hy_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec)/mu_0
|
||||
hx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec)/mu_0
|
||||
hy_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec)/mu_0
|
||||
|
||||
# Define the components of the derivative
|
||||
Hd = (hx_px*hy_py - hx_py*hy_px)
|
||||
Hd_uV = hx_px_u(hy_py)*v - hx_py*hy_px_u(v) + hx_px*hy_py_u(v) - hx_py_u(hy_px)
|
||||
# Need to fix this to reflect the adjoint
|
||||
if 'zxx' in self.rxType:
|
||||
Zij = ( ex_px*hy_py - ex_py*hy_px)/Hd
|
||||
ZijN_uV = ex_px_u(hy_py)*v - ex_py*hy_px_u(v) + ex_px*hy_py_u(v) - ex_py_u(hy_px)*v
|
||||
elif 'zxy' in self.rxType:
|
||||
Zij = (-ex_px*hx_py + ex_py*hx_px)/Hd
|
||||
ZijN_uV = -ex_px_u(hx_py)*v + ex_py*hx_px_u(v) - ex_px*hx_py_u(v) + ex_py_u(hx_px)*v
|
||||
elif 'zyx' in self.rxType:
|
||||
Zij = ( ey_px*hy_py - ey_py*hy_px)/Hd
|
||||
ZijN_uV = ey_px_u(hy_py)*v - ey_py*hy_px_u(v) + ey_px*hy_py_u(v) - ey_py_u(hy_px)*v
|
||||
elif 'zyy' in self.rxType:
|
||||
Zij = (-ey_px*hx_py + ey_py*hx_px)/Hd
|
||||
ZijN_uV = -ey_px_u(hx_py)*v + ey_py*hx_px_u(v) - ey_px*hx_py_u(v) +ey_py_u(hx_px)*v
|
||||
|
||||
# Calculate the complex derivative
|
||||
PDeriv_real = ZijN_uV/Hd + (Hd_uV/Hd)*Zij^T
|
||||
# Extract the data
|
||||
if real_or_imag == 'imag':
|
||||
Pv = 1j*PDeriv_real
|
||||
|
||||
Reference in New Issue
Block a user