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

Fixed dimension bugs in code.

Browse Source
This commit is contained in:
GudniRos
2015-11-04 23:57:11 -08:00
parent 9f69a33512
commit d7234cea9e
4 changed files with 1451 additions and 1426 deletions
Download Patch File Download Diff File Expand all files Collapse all files
notebooks/MT Script-3D_layerTest.ipynb
+1396 -1413
View File
File diff suppressed because one or more lines are too long
simpegMT/BaseMT.py
+2
View File
@@ -123,6 +123,8 @@ class BaseMTProblem(BaseFDEMProblem):
# Get the adjoint projectFieldsDeriv
# PTv needs to be nE,
PTv = rx.projectFieldsDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
# Need to reshape PTv
PTv = np.hstack((mkvc(PTv[:len(PTv)/2],2),mkvc(PTv[len(PTv)/2::],2)))
# Get the
dA_duIT = ATinv * PTv
dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True)
simpegMT/SurveyMT.py
+11 -11
View File
@@ -211,7 +211,7 @@ class RxMT(Survey.BaseRx):
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
# Define the components of the derivative
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
Hd_uV = sDiag(hy_py)*hx_px_u(v) + sDiag(hx_px)*hy_py_u(v) - sDiag(hx_py)*hy_px_u(v) - sDiag(hy_px)*hx_py_u(v)
# Calculate components
if 'zxx' in self.rxType:
Zij = sDiag(Hd*( sDiag(ex_px)*hy_py - sDiag(ex_py)*hy_px ))
@@ -228,7 +228,7 @@ class RxMT(Survey.BaseRx):
# Calculate the complex derivative
PDeriv_complex = Hd * (ZijN_uV - Zij * Hd_uV )
# ero
# Extract the real number for the real/imag components.
Pv = np.array(getattr(PDeriv_complex, real_or_imag))
elif adjoint:
@@ -258,14 +258,14 @@ class RxMT(Survey.BaseRx):
Pby = mesh.getInterpolationMat(bFLocs,'Fy')
# Get the fields at location
# px: x-polaration and py: y-polaration.
aex_px = mkvc(f[src,'e_px'],2).T*Pex.T
aey_px = mkvc(f[src,'e_px'],2).T*Pey.T
aex_py = mkvc(f[src,'e_py'],2).T*Pex.T
aey_py = mkvc(f[src,'e_py'],2).T*Pey.T
ahx_px = mkvc(f[src,'b_px'],2).T/mu_0*Pbx.T
ahy_px = mkvc(f[src,'b_px'],2).T/mu_0*Pby.T
ahx_py = mkvc(f[src,'b_py'],2).T/mu_0*Pbx.T
ahy_py = mkvc(f[src,'b_py'],2).T/mu_0*Pby.T
aex_px = mkvc(mkvc(f[src,'e_px'],2).T*Pex.T)
aey_px = mkvc(mkvc(f[src,'e_px'],2).T*Pey.T)
aex_py = mkvc(mkvc(f[src,'e_py'],2).T*Pex.T)
aey_py = mkvc(mkvc(f[src,'e_py'],2).T*Pey.T)
ahx_px = mkvc(mkvc(f[src,'b_px'],2).T/mu_0*Pbx.T)
ahy_px = mkvc(mkvc(f[src,'b_px'],2).T/mu_0*Pby.T)
ahx_py = mkvc(mkvc(f[src,'b_py'],2).T/mu_0*Pbx.T)
ahy_py = mkvc(mkvc(f[src,'b_py'],2).T/mu_0*Pby.T)
# Derivatives as lambda functions
aex_px_u = lambda vec: f._e_pxDeriv_u(src,Pex.T*vec,adjoint=True)
aey_px_u = lambda vec: f._e_pxDeriv_u(src,Pey.T*vec,adjoint=True)
@@ -281,7 +281,7 @@ class RxMT(Survey.BaseRx):
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2))
# Define the components of the derivative
aHd = sDiag(1./(sDiag(ahy_py)*ahx_px - sDiag(ahy_px)*ahx_py))
aHd = sDiag(1./(sDiag(ahx_px)*ahy_py - sDiag(ahx_py)*ahy_px))
aHd_uV = lambda x: ahx_px_u(sDiag(ahy_py)*x) + ahx_px_u(sDiag(ahy_py)*x) - ahy_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(ahy_px)*x)
# Need to fix this to reflect the adjoint
if 'zxx' in self.rxType:
simpegMT/Tests/test_Problem3D_againstAnalytic.py
+42 -2
View File
@@ -32,7 +32,7 @@ def getInputs():
## Setup the the survey object
# Receiver locations
rx_x, rx_y = np.meshgrid(np.arange(-3000,3001,500),np.arange(-1000,1001,500))
rx_loc = np.array([[0, 0, 0]]) #np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
rx_loc = np.array([[0, 0, 0]]) #,[-100,-100,0],[100,100,0]]) #np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
return M, freqs, rx_loc, elev
@@ -87,7 +87,7 @@ def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False,expMap=True
# Make a receiver list
rxList = []
if comp == 'All':
for rxType in ['zxyr','zxyi','zyxr','zyxi']:
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi',]:
rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,rxType))
else:
rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,comp))
@@ -121,6 +121,46 @@ def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False,expMap=True
return (survey, problem)
def setupSimpegMTfwd_eForm_ps_multiRx(inputSetup,comp='All',singleFreq=False,expMap=True):
M,freqs,sig,sigBG,rx_loc = inputSetup
# Add to the receiver list
rx_loc = np.vstack((rx_loc,np.array([[-100,100,0]])))# ,[-100,-100,0],[100,-100,0],[100,100,0]])))
# Make a receiver list
rxList = []
if comp == 'All':
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi',]:
rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,rxType))
else:
rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,comp))
# Source list
srcList =[]
if singleFreq:
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,singleFreq))
else:
for freq in freqs:
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
# Survey MT
survey = simpegmt.SurveyMT.SurveyMT(srcList)
## Setup the problem object
sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
if expMap:
problem = simpegmt.ProblemMT3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) )
problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
problem.curModel = np.log(sig)
else:
problem = simpegmt.ProblemMT3D.eForm_ps(M,sigmaPrimary= sigma1d)
problem.curModel = sig
problem.pair(survey)
problem.verbose = False
try:
from pymatsolver import MumpsSolver
problem.Solver = MumpsSolver
except:
pass
return (survey, problem)
def getAppResPhs(MTdata):
# Make impedance
def appResPhs(freq,z):