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

Field projection derivatives working

Browse Source
This commit is contained in:
GudniRos
2015-11-02 17:26:52 -08:00
parent 960cb0a3e5
commit 0e81faf217
2 changed files with 24 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegMT/FieldsMT.py
+6 -4
View File
@@ -175,10 +175,10 @@ class FieldsMT_3D(FieldsMT):
return self._e_pyPrimary(e_pySolution,srcList) + self._e_pySecondary(e_pySolution,srcList)
def _e_pxDeriv_u(self, src, v, adjoint = False):
return v
return v[:len(v)/2]
def _e_pyDeriv_u(self, src, v, adjoint = False):
return v
return v[len(v)/2::]
def _e_pxDeriv_m(self, src, v, adjoint = False):
# assuming primary does not depend on the model
@@ -233,13 +233,15 @@ class FieldsMT_3D(FieldsMT):
# NOTE: v needs to be length 2*nE to account for both polarizations
def _b_pxSecondaryDeriv_u(self, src, v, adjoint = False):
C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,0]])
# C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,0]])
C = sp.hstack((self.mesh.edgeCurl,Utils.spzeros(self.mesh.nF,self.mesh.nE))) # This works for adjoint = None
if adjoint:
return - 1./(1j*omega(src.freq)) * (C.T * v)
return - 1./(1j*omega(src.freq)) * (C * v)
def _b_pySecondaryDeriv_u(self, src, v, adjoint = False):
C = sp.kron(self.mesh.edgeCurl,[[0,0],[0,1]])
# C = sp.kron(self.mesh.edgeCurl,[[0,0],[0,1]])
C = sp.hstack((Utils.spzeros(self.mesh.nF,self.mesh.nE),self.mesh.edgeCurl)) # This works for adjoint = None
if adjoint:
return - 1./(1j*omega(src.freq)) * (C.T * v)
return - 1./(1j*omega(src.freq)) * (C * v)
simpegMT/SurveyMT.py
+18 -20
View File
@@ -197,36 +197,34 @@ class RxMT(Survey.BaseRx):
hx_py = Pbx*f[src,'b_py']/mu_0
hy_py = Pby*f[src,'b_py']/mu_0
# Derivatives as lambda functions
spPe = Utils.spzeros(self.nD,mesh.nE)
spPb = Utils.spzeros(self.nD,mesh.nF)
# The size of the diratives should be nD,nU
ex_px_u = lambda vec: sp.hstack((Pex,spPe))*f._e_pxDeriv_u(src,vec)
ey_px_u = lambda vec: sp.hstack((Pey,spPe))*f._e_pxDeriv_u(src,vec)
ex_py_u = lambda vec: sp.hstack((spPe,Pex))*f._e_pyDeriv_u(src,vec)
ey_py_u = lambda vec: sp.hstack((spPe,Pey))*f._e_pyDeriv_u(src,vec)
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)
# NOTE: Think b_p?Deriv_u should return a 2*nF size matrix
hx_px_u = lambda vec: sp.hstack((Pbx,spPb))*f._b_pxDeriv_u(src,vec)/mu_0
hy_px_u = lambda vec: sp.hstack((Pby,spPb))*f._b_pxDeriv_u(src,vec)/mu_0
hx_py_u = lambda vec: sp.hstack((spPb,Pbx))*f._b_pyDeriv_u(src,vec)/mu_0
hy_py_u = lambda vec: sp.hstack((spPb,Pby))*f._b_pyDeriv_u(src,vec)/mu_0
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
# Update the input vector
# v = mkvc(v,2) # Make v into a column vector
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
# Define the components of the derivative
Hd = 1./(hx_px*hy_py - hx_py*hy_px)
Hd = sDiag(1./(sDiag(hx_px)*hy_py - sDiag(hx_py)*hy_px))
Hd_uV = hx_px_u(v)*hy_py + hx_px*hy_py_u(v) - hx_py*hy_px_u(v) - hx_py_u(v)*hy_px
# Calculate components
if 'zxx' in self.rxType:
Zij = ( ex_px*hy_py - ex_py*hy_px)*Hd
ZijN_uV = ex_px_u(v)*hy_py + ex_px*hy_py_u(v) - ex_py*hy_px_u(v) - ex_py_u(v)*hy_px
Zij = sDiag(Hd*( sDiag(ex_px)*hy_py - sDiag(ex_py)*hy_px ))
ZijN_uV = sDiag(hy_py)*ex_px_u(v) + sDiag(ex_px)*hy_py_u(v) - sDiag(ex_py)*hy_px_u(v) - sDiag(hy_px)*ex_py_u(v)
elif 'zxy' in self.rxType:
Zij = (-ex_px*hx_py + ex_py*hx_px)*Hd
ZijN_uV = -ex_px_u(v)*hx_py - ex_px*hx_py_u(v) + ex_py*hx_px_u(v) + ex_py_u(v)*hx_px
Zij = sDiag(Hd*(-sDiag(ex_px)*hx_py + sDiag(ex_py)*hx_px ))
ZijN_uV = -sDiag(hx_py)*ex_px_u(v) - sDiag(ex_px)*hx_py_u(v) + sDiag(ex_py)*hx_px_u(v) + sDiag(hx_px)*ex_py_u(v)
elif 'zyx' in self.rxType:
Zij = ( ey_px*hy_py - ey_py*hy_px)*Hd
ZijN_uV = ey_px_u(v)*hy_py + ey_px*hy_py_u(v) - ey_py*hy_px_u(v) - ey_py_u(v)*hy_px
Zij = sDiag(Hd*( sDiag(ey_px)*hy_py - sDiag(ey_py)*hy_px ))
ZijN_uV = sDiag(hy_py)*ey_px_u(v) + sDiag(ey_px)*hy_py_u(v) - sDiag(ey_py)*hy_px_u(v) - sDiag(hy_px)*ey_py_u(v)
elif 'zyy' in self.rxType:
Zij = (-ey_px*hx_py + ey_py*hx_px)*Hd
ZijN_uV = -ey_px_u(v)*hx_py - ey_px*hx_py_u(v) + ey_py*hx_px_u(v) + ey_py_u(v)*hx_px
Zij = sDiag(Hd*(-sDiag(ey_px)*hx_py + sDiag(ey_py)*hx_px ))
ZijN_uV = -sDiag(hx_py)*ey_px_u(v) - sDiag(ey_px)*hx_py_u(v) + sDiag(ey_py)*hx_px_u(v) + sDiag(hx_px)*ey_py_u(v)
# Calculate the complex derivative
PDeriv_complex = Hd * (ZijN_uV - Zij * Hd_uV )