mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-28 21:12:41 +08:00
Reorignized to have u = [u_px,u_py].
Everything runs but tests are not passing.
This commit is contained in:
GudniRos
+19
-14
@@ -2,8 +2,7 @@ from simpegEM.FDEM import BaseFDEMProblem
|
||||
from SurveyMT import SurveyMT
|
||||
from DataMT import DataMT
|
||||
from FieldsMT import FieldsMT
|
||||
from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc
|
||||
import numpy as np
|
||||
from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
|
||||
|
||||
class BaseMTProblem(BaseFDEMProblem):
|
||||
|
||||
@@ -75,24 +74,27 @@ class BaseMTProblem(BaseFDEMProblem):
|
||||
# Loop all the frequenies
|
||||
for freq in self.survey.freqs:
|
||||
dA_du = self.getA(freq) #
|
||||
|
||||
dA_duI = self.Solver(dA_du, **self.solverOpts)
|
||||
|
||||
for src in self.survey.getSrcByFreq(freq):
|
||||
# We need fDeriv_m = df/du*du/dm + df/dm
|
||||
# Construct du/dm, it requires a solve
|
||||
ftype = self._fieldType + 'Solution'
|
||||
u_src = u[src, ftype]
|
||||
dA_dm = self.getADeriv_m(freq, u_src, v)
|
||||
dRHS_dm = self.getRHSDeriv_m(freq, v)
|
||||
# NOTE: need to account for the 2 polarizations in the derivatives.
|
||||
u_src = u[src,:]
|
||||
# dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
|
||||
dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
|
||||
dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
|
||||
if dRHS_dm is None:
|
||||
du_dm = dA_duI * ( - dA_dm )
|
||||
du_dm = dA_duI * ( -dA_dm )
|
||||
else:
|
||||
du_dm = dA_duI * ( - dA_dm + dRHS_dm )
|
||||
du_dm = dA_duI * ( -dA_dm + dRHS_dm )
|
||||
# Calculate the projection derivatives
|
||||
for rx in src.rxList:
|
||||
# Get the projection derivative
|
||||
PDeriv_u = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m
|
||||
Jv[src, rx] = PDeriv_u(du_dm)
|
||||
# v should be of size nE,2 (each column for 2 polarizations)
|
||||
PDeriv_u = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, we don't have have PDeriv wrt m
|
||||
Jv[src, rx] = PDeriv_u(mkvc(du_dm))
|
||||
# Return the vectorized sensitivities
|
||||
return mkvc(Jv)
|
||||
|
||||
@@ -110,25 +112,28 @@ class BaseMTProblem(BaseFDEMProblem):
|
||||
|
||||
for freq in self.survey.freqs:
|
||||
AT = self.getA(freq).T
|
||||
|
||||
ATinv = self.Solver(AT, **self.solverOpts)
|
||||
|
||||
for src in self.survey.getSrcByFreq(freq):
|
||||
ftype = self._fieldType + 'Solution'
|
||||
u_src = u[src, ftype]
|
||||
u_src = u[src, :]
|
||||
|
||||
for rx in src.rxList:
|
||||
# Get the adjoint projectFieldsDeriv
|
||||
PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
|
||||
# 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
|
||||
# Get the
|
||||
dA_duIT = ATinv * PTv
|
||||
dA_dmT = self.getADeriv_m(freq, u_src, dA_duIT, adjoint=True)
|
||||
dRHS_dmT = self.getRHSDeriv_m(freq, dA_duIT, adjoint=True)
|
||||
dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True)
|
||||
dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
|
||||
# Make du_dmT
|
||||
if dRHS_dmT is None:
|
||||
du_dmT = -dA_dmT
|
||||
else:
|
||||
du_dmT = -dA_dmT + dRHS_dmT
|
||||
# Select the correct component
|
||||
# du_dmT needs to be of size nC,
|
||||
real_or_imag = rx.projComp
|
||||
if real_or_imag == 'real':
|
||||
Jtv += du_dmT.real
|
||||
|
||||
@@ -228,23 +228,18 @@ class FieldsMT_3D(FieldsMT):
|
||||
def _b_px(self, eSolution, srcList):
|
||||
return self._b_pxPrimary(eSolution, srcList) + self._b_pxSecondary(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)
|
||||
|
||||
# NOTE: v needs to be length 2*nE to account for both polarizations
|
||||
def _b_pxSecondaryDeriv_u(self, src, v, adjoint = False):
|
||||
C = self.mesh.edgeCurl
|
||||
C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,1]])
|
||||
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 = self.mesh.edgeCurl
|
||||
C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,1]])
|
||||
if adjoint:
|
||||
return - 1./(1j*omega(src.freq)) * (C.T * v)
|
||||
return - 1./(1j*omega(src.freq)) * (C * v)
|
||||
|
||||
@@ -92,7 +92,7 @@ class eForm_psField(BaseMTProblem):
|
||||
"""
|
||||
|
||||
Src = self.survey.getSrcByFreq(freq)[0]
|
||||
S_eDeriv = Src.S_eDeriv(self, v, adjoint)
|
||||
S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
|
||||
return -1j * omega(freq) * S_eDeriv
|
||||
|
||||
def fields(self, m):
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
from SimPEG import Survey, Problem, Utils, Models, np, sp, SolverLU as SimpegSolver
|
||||
from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver
|
||||
from simpegEM.Utils.EMUtils import omega
|
||||
from scipy.constants import mu_0
|
||||
from simpegMT.BaseMT import BaseMTProblem
|
||||
@@ -20,7 +20,7 @@ class eForm_ps(BaseMTProblem):
|
||||
"""
|
||||
|
||||
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
|
||||
_fieldType = 'e_px'
|
||||
_fieldType = 'e'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = FieldsMT_3D
|
||||
_sigmaPrimary = None
|
||||
@@ -52,14 +52,32 @@ class eForm_ps(BaseMTProblem):
|
||||
|
||||
def getADeriv_m(self, freq, u, v, adjoint=False):
|
||||
|
||||
dMe_dsig = self.MeSigmaDeriv( u )
|
||||
# Nee to account for both the polarizations
|
||||
# dMe_dsig = (self.MeSigmaDeriv( u['e_pxSolution'] ) + self.MeSigmaDeriv( u['e_pySolution'] ))
|
||||
# dMe_dsig = (self.MeSigmaDeriv( u['e_pxSolution'] + u['e_pySolution'] ))
|
||||
|
||||
# # dMe_dsig = self.MeSigmaDeriv( u )
|
||||
# if adjoint:
|
||||
# return 1j * omega(freq) * ( dMe_dsig.T * v ) # As in simpegEM
|
||||
|
||||
# return 1j * omega(freq) * ( dMe_dsig * v ) # As in simpegEM
|
||||
|
||||
# This considers both polarizations and returns a nE,2 matrix for each polarization
|
||||
if adjoint:
|
||||
return 1j * omega(freq) * ( dMe_dsig.T * v ) # As in simpegEM
|
||||
# return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
|
||||
dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v
|
||||
else:
|
||||
# Need a nE,2 matrix to be returned
|
||||
dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
|
||||
return 1j * omega(freq) * dMe_dsigV
|
||||
# Stacking them
|
||||
|
||||
# if adjoint:
|
||||
# dMe_dsig = sp.vstack((self.MeSigmaDeriv( u['e_pxSolution'] ), self.MeSigmaDeriv( u['e_pxSolution'] ) )).T
|
||||
# # self.MeSigmaDeriv(u['e_pySolution'] ).T*v,2) ))
|
||||
# else:
|
||||
# dMe_dsig = sp.vstack((self.MeSigmaDeriv( u['e_pxSolution'] ), self.MeSigmaDeriv( u['e_pxSolution'] ) ))
|
||||
# return 1j * omega(freq) * dMe_dsig*v
|
||||
|
||||
return 1j * omega(freq) * ( dMe_dsig * v ) # As in simpegEM
|
||||
# return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
@@ -80,7 +98,7 @@ class eForm_ps(BaseMTProblem):
|
||||
"""
|
||||
|
||||
Src = self.survey.getSrcByFreq(freq)[0]
|
||||
S_eDeriv = Src.S_eDeriv(self, v, adjoint)
|
||||
S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
|
||||
return -1j * omega(freq) * S_eDeriv
|
||||
|
||||
def fields(self, m):
|
||||
|
||||
+27
-21
@@ -173,7 +173,7 @@ class RxMT(Survey.BaseRx):
|
||||
dP_db = mkvc( Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pbx*f._bDeriv_u(src,v)/mu_0),2)
|
||||
PDeriv_complex = np.sum(np.hstack((dP_de,dP_db)),1)
|
||||
elif self.projType is 'Z2D':
|
||||
raise NotImplementedError('Has not be implement for 2D impedance tensor')
|
||||
raise NotImplementedError('Has not been implement for 2D impedance tensor')
|
||||
elif self.projType is 'Z3D':
|
||||
if self.locs.ndim == 3:
|
||||
eFLocs = self.locs[:,:,0]
|
||||
@@ -197,14 +197,15 @@ class RxMT(Survey.BaseRx):
|
||||
hx_py = Pbx*mkvc(f[src,'b_py']/mu_0,2)
|
||||
hy_py = Pby*mkvc(f[src,'b_py']/mu_0,2)
|
||||
# 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
|
||||
ex_px_u = lambda vec: sp.hstack((Pex,Pex))*f._e_pxDeriv_u(src,vec)
|
||||
ey_px_u = lambda vec: sp.hstack((Pey,Pey))*f._e_pxDeriv_u(src,vec)
|
||||
ex_py_u = lambda vec: sp.hstack((Pex,Pex))*f._e_pyDeriv_u(src,vec)
|
||||
ey_py_u = lambda vec: sp.hstack((Pey,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,Pbx))*f._b_pxDeriv_u(src,vec)/mu_0
|
||||
hy_px_u = lambda vec: sp.hstack((Pby,Pby))*f._b_pxDeriv_u(src,vec)/mu_0
|
||||
hx_py_u = lambda vec: sp.hstack((Pbx,Pbx))*f._b_pyDeriv_u(src,vec)/mu_0
|
||||
hy_py_u = lambda vec: sp.hstack((Pby,Pby))*f._b_pyDeriv_u(src,vec)/mu_0
|
||||
|
||||
# Update the input vector
|
||||
v = mkvc(v,2) # Make v into a column vector
|
||||
@@ -226,7 +227,7 @@ class RxMT(Survey.BaseRx):
|
||||
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)
|
||||
PDeriv_complex = ZijN_uV*Hd.toarray() - Zij * (Hd_uV*Hd.toarray())
|
||||
|
||||
# Extract the real number for the real/imag components.
|
||||
Pv = np.array(getattr(PDeriv_complex, real_or_imag))
|
||||
@@ -266,14 +267,14 @@ class RxMT(Survey.BaseRx):
|
||||
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
|
||||
# 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)
|
||||
aex_py_u = lambda vec: f._e_pyDeriv_u(src,Pex.T*vec,adjoint=True)
|
||||
aey_py_u = lambda vec: f._e_pyDeriv_u(src,Pey.T*vec,adjoint=True)
|
||||
ahx_px_u = lambda vec: f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
|
||||
ahy_px_u = lambda vec: f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
|
||||
ahx_py_u = lambda vec: f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
|
||||
ahy_py_u = lambda vec: f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
|
||||
aex_px_u = lambda vec: f._e_pxDeriv_u(src,sp.hstack((Pex,Pex)).T*vec,adjoint=True)
|
||||
aey_px_u = lambda vec: f._e_pxDeriv_u(src,sp.hstack((Pey,Pey)).T*vec,adjoint=True)
|
||||
aex_py_u = lambda vec: f._e_pyDeriv_u(src,sp.hstack((Pex,Pex)).T*vec,adjoint=True)
|
||||
aey_py_u = lambda vec: f._e_pyDeriv_u(src,sp.hstack((Pey,Pey)).T*vec,adjoint=True)
|
||||
ahx_px_u = lambda vec: f._b_pxDeriv_u(src,sp.hstack((Pbx,Pbx)).T*vec,adjoint=True)/mu_0
|
||||
ahy_px_u = lambda vec: f._b_pxDeriv_u(src,sp.hstack((Pby,Pby)).T*vec,adjoint=True)/mu_0
|
||||
ahx_py_u = lambda vec: f._b_pyDeriv_u(src,sp.hstack((Pbx,Pbx)).T*vec,adjoint=True)/mu_0
|
||||
ahy_py_u = lambda vec: f._b_pyDeriv_u(src,sp.hstack((Pby,Pby)).T*vec,adjoint=True)/mu_0
|
||||
|
||||
# Update the input vector
|
||||
v = mkvc(v,2) # Make v into a column vector
|
||||
@@ -376,7 +377,7 @@ class srcMT_polxy_1Dprimary(srcMT):
|
||||
C = problem.mesh.nodalGrad
|
||||
elif problem.mesh.dim == 3:
|
||||
C = problem.mesh.edgeCurl
|
||||
bBG_bp = (- C * self.ePrimary(problem) )/( 1j*omega(self.freq) )
|
||||
bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
|
||||
return bBG_bp
|
||||
|
||||
def S_e(self,problem):
|
||||
@@ -399,7 +400,7 @@ class srcMT_polxy_1Dprimary(srcMT):
|
||||
Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p)
|
||||
return (Mesigma - Mesigma_p) * e_p
|
||||
|
||||
def S_eDeriv(self, problem, v, adjoint = False):
|
||||
def S_eDeriv_m(self, problem, v, adjoint = False):
|
||||
'''
|
||||
Get the derivative of S_e wrt to sigma (m)
|
||||
'''
|
||||
@@ -413,7 +414,12 @@ class srcMT_polxy_1Dprimary(srcMT):
|
||||
if problem.mesh.dim == 3:
|
||||
# Need to take the derivative of both u_px and u_py
|
||||
ePri = self.ePrimary(problem)
|
||||
MsigmaDeriv = problem.MeSigmaDeriv(ePri[:,0]) + problem.MeSigmaDeriv(ePri[:,1])
|
||||
# MsigmaDeriv = problem.MeSigmaDeriv(ePri[:,0]) + problem.MeSigmaDeriv(ePri[:,1])
|
||||
# MsigmaDeriv = problem.MeSigmaDeriv(np.sum(ePri,axis=1))
|
||||
if adjoint:
|
||||
return sp.hstack(( problem.MeSigmaDeriv(ePri[:,0]).T, problem.MeSigmaDeriv(ePri[:,1]).T ))*v
|
||||
else:
|
||||
return np.hstack(( mkvc(problem.MeSigmaDeriv(ePri[:,0]) * v,2), mkvc(problem.MeSigmaDeriv(ePri[:,1])*v,2) ))
|
||||
if adjoint:
|
||||
#
|
||||
return MsigmaDeriv.T * v
|
||||
|
||||
@@ -23,7 +23,8 @@ def getInputs():
|
||||
"""
|
||||
# Make a mesh
|
||||
# M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360])
|
||||
M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,10,-1.3),(1000.,2),(1000,10,1.3)]], x0=['C','C','C'])# Setup the model
|
||||
# M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model
|
||||
M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(500,8,-1.3),(500.,8),(500,8,1.3)]], x0=['C','C','C'])# Setup the model
|
||||
# Set the frequencies
|
||||
freqs = np.logspace(1,-3,5)
|
||||
elev = 0
|
||||
@@ -36,6 +37,17 @@ def getInputs():
|
||||
return M, freqs, rx_loc, elev
|
||||
|
||||
|
||||
def random(conds):
|
||||
''' Returns a halfspace model based on the inputs'''
|
||||
M, freqs, rx_loc, elev = getInputs()
|
||||
|
||||
# Backround
|
||||
sigBG = np.ones(M.nC)*conds
|
||||
# Add randomness to the model (10% of the value).
|
||||
sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 )
|
||||
|
||||
return (M, freqs, sig, sigBG, rx_loc)
|
||||
|
||||
def halfSpace(conds):
|
||||
''' Returns a halfspace model based on the inputs'''
|
||||
M, freqs, rx_loc, elev = getInputs()
|
||||
@@ -70,7 +82,7 @@ def twoLayer(conds):
|
||||
|
||||
return (M, freqs, sig, sigBG, rx_loc)
|
||||
|
||||
def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False):
|
||||
def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False,expMap=True):
|
||||
M,freqs,sig,sigBG,rx_loc = inputSetup
|
||||
# Make a receiver list
|
||||
rxList = []
|
||||
@@ -81,9 +93,9 @@ def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False):
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,comp))
|
||||
# Source list
|
||||
srcList =[]
|
||||
sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
|
||||
|
||||
if singleFreq:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freqs[2]))
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,singleFreq))
|
||||
else:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
|
||||
@@ -91,16 +103,21 @@ def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False):
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
|
||||
## Setup the problem object
|
||||
problem = simpegmt.ProblemMT3D.eForm_ps(M,sigmaPrimary=sigma1d)
|
||||
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
|
||||
problem.pair(survey)
|
||||
problem.curMod = sig
|
||||
problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
|
||||
|
||||
return (survey, problem)
|
||||
|
||||
@@ -113,34 +130,49 @@ def getAppResPhs(MTdata):
|
||||
recData = MTdata.toRecArray('Complex')
|
||||
return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx'])
|
||||
|
||||
def adjointTest(inputSetup):
|
||||
|
||||
survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup)
|
||||
print 'Adjoint test of eForm primary/secondary\n'
|
||||
|
||||
m = problem.curMod
|
||||
|
||||
# if addrandoms is True:
|
||||
# m = m + np.random.randn(problem.mesh.nC)*CONDUCTIVITY*1e-1
|
||||
def test_JvecAdjoint(inputSetup,comp='All',freq=False):
|
||||
(M, freqs, sig, sigBG, rx_loc) = inputSetup
|
||||
survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=freq)
|
||||
print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs))
|
||||
|
||||
m = sig
|
||||
u = problem.fields(m)
|
||||
|
||||
v = np.random.rand(survey.nD)
|
||||
v = np.random.rand(survey.nD,)
|
||||
# print problem.PropMap.PropModel.nP
|
||||
w = np.random.rand(problem.mesh.nC)
|
||||
w = np.random.rand(problem.mesh.nC,)
|
||||
|
||||
vJw = v.dot(problem.Jvec(m, w, u))
|
||||
wJtv = w.dot(problem.Jtvec(m, v, u))
|
||||
vJw = v.ravel().dot(problem.Jvec(m, w, u))
|
||||
wJtv = w.ravel().dot(problem.Jtvec(m, v, u))
|
||||
tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
|
||||
print ' vJw wJtv vJw - wJtv tol abs(vJw - wJtv) < tol'
|
||||
print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
|
||||
return np.abs(vJw - wJtv) < tol
|
||||
|
||||
# Test the Jvec derivative
|
||||
def test_DerivJvec(inputSetup,comp='All',freq=False,expMap=True):
|
||||
(M, freqs, sig, sigBG, rx_loc) = inputSetup
|
||||
survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp=comp,singleFreq=freq,expMap=expMap)
|
||||
print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs)
|
||||
# problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
|
||||
# problem.sigmaPrimary = np.log(sigBG)
|
||||
x0 = np.log(sigBG)
|
||||
# cond = sig[0]
|
||||
# x0 = np.log(np.ones(problem.mesh.nC)*cond)
|
||||
# problem.sigmaPrimary = x0
|
||||
# if True:
|
||||
# x0 = x0 + np.random.randn(problem.mesh.nC)*cond*1e-1
|
||||
survey = problem.survey
|
||||
def fun(x):
|
||||
return survey.dpred(x), lambda x: problem.Jvec(x0, x)
|
||||
return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
|
||||
|
||||
def derivProjfields(inputSetup):
|
||||
|
||||
survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup)
|
||||
def test_DerivProjfields(inputSetup,comp='All',freq=False):
|
||||
|
||||
survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp,freq)
|
||||
print 'Derivative test of data projection for eFormulation primary/secondary\n\n'
|
||||
# problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
|
||||
|
||||
# Define a src and rx
|
||||
src = survey.srcList[-1]
|
||||
@@ -158,39 +190,58 @@ def derivProjfields(inputSetup):
|
||||
return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)
|
||||
|
||||
|
||||
def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True):
|
||||
def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False):
|
||||
if appR:
|
||||
label = 'resistivity'
|
||||
else:
|
||||
label = 'phase'
|
||||
# Make the survey and the problem
|
||||
survey, problem = setupSimpegMTfwd_eForm_ps(halfSpace(sigmaHalf))
|
||||
survey, problem = setupSimpegMTfwd_eForm_ps(halfSpace(sigmaHalf),expMap=expMap)
|
||||
print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf)
|
||||
|
||||
data = problem.dataPair(survey,survey.dpred(problem.curMod))
|
||||
data = problem.dataPair(survey,survey.dpred(problem.curModel))
|
||||
# Calculate the app phs
|
||||
app_rpxy, app_rpyx = np.array(getAppResPhs(data))
|
||||
if appR:
|
||||
return np.all(np.abs(app_rpxy[0,:] - np.ones(survey.nFreq)/sigmaHalf) * sigmaHalf < .35)
|
||||
return np.all(np.abs(app_rpxy[0,:] - np.ones(survey.nFreq)/sigmaHalf) * sigmaHalf < .4)
|
||||
else:
|
||||
return np.all(np.abs(app_rpxy[1,:] + np.ones(survey.nFreq)*135) / 135 < .35)
|
||||
return np.all(np.abs(app_rpxy[1,:] + np.ones(survey.nFreq)*135) / 135 < .4)
|
||||
|
||||
class TestAnalytics(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
pass
|
||||
# def test_appRes2en1(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(2e-1))
|
||||
def test_appRes1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2))
|
||||
def test_appPhs1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False))
|
||||
# # Test apparent resistivity and phase
|
||||
# def test_appRes1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2))
|
||||
# def test_appPhs1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False))
|
||||
|
||||
def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3))
|
||||
def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False))
|
||||
# def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3))
|
||||
# def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False))
|
||||
|
||||
# Do a derivative test
|
||||
def test_deriv1(self):self.assertTrue(derivProjfields(halfSpace(1e-3)))
|
||||
def test_derivProj1(self):self.assertTrue(test_DerivProjfields(halfSpace(1e-2)))
|
||||
|
||||
# Do a derivative test of Jvec
|
||||
# def test_derivJvec_zxxr(self):self.assertTrue(test_DerivJvec(random(1e-2),'zxxr',1))
|
||||
# def test_derivJvec_zxxi(self):self.assertTrue(test_DerivJvec(random(1e-2),'zxxi',1))
|
||||
def test_derivJvec_zxyr(self):self.assertTrue(test_DerivJvec(random(1e-2),'zxyr',.1))
|
||||
# def test_derivJvec_zxyi(self):self.assertTrue(test_DerivJvec(random(1e-2),'zxyi',1))
|
||||
def test_derivJvec_zyxr(self):self.assertTrue(test_DerivJvec(random(1e-2),'zyxr',1))
|
||||
# def test_derivJvec_zyxi(self):self.assertTrue(test_DerivJvec(random(1e-2),'zyxi',1))
|
||||
# def test_derivJvec_zyyr(self):self.assertTrue(test_DerivJvec(random(1e-2),'zyyr',1))
|
||||
# def test_derivJvec_zyyi(self):self.assertTrue(test_DerivJvec(random(1e-2),'zyyi',1))
|
||||
# def test_derivJvec_All(self):self.assertTrue(test_DerivJvec(random(1e-2),'All',1))
|
||||
|
||||
# Test the adjoint of Jvec and Jtvec
|
||||
def test_adjointDeriv1(self):self.assertTrue(adjointTest(halfSpace(1e-3)))
|
||||
# def test_JvecAdjoint_zxxr(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zxxr',1))
|
||||
# def test_JvecAdjoint_zxxi(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zxxi',1))
|
||||
def test_JvecAdjoint_zxyr(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zxyr',.1))
|
||||
# def test_JvecAdjoint_zxyi(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zxyi',1))
|
||||
def test_JvecAdjoint_zyxr(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zyxr',1))
|
||||
# def test_JvecAdjoint_zyxi(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zyxi',1))
|
||||
# def test_JvecAdjoint_zyyr(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zyyr',1))
|
||||
# def test_JvecAdjoint_zyyi(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zyyi',1))
|
||||
# def test_JvecAdjoint_All(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'All',1))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -3,6 +3,8 @@ import numpy as np, matplotlib.pyplot as plt, sys
|
||||
import SimPEG as simpeg
|
||||
import simpegMT as simpegmt
|
||||
import numpy.lib.recfunctions as recFunc
|
||||
from scipy.constants import mu_0
|
||||
|
||||
def getAppRes(MTdata):
|
||||
# Make impedance
|
||||
zList = []
|
||||
@@ -23,6 +25,10 @@ def appResPhs(freq,z):
|
||||
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
|
||||
return app_res, app_phs
|
||||
|
||||
def skindepth(rho,freq):
|
||||
''' Function to calculate the skindepth of EM waves'''
|
||||
return np.sqrt( (rho*((1/(freq * mu_0 * np.pi )))))
|
||||
|
||||
def rec2ndarr(x,dt=float):
|
||||
return x.view((dt, len(x.dtype.names)))
|
||||
|
||||
|
||||
Reference in New Issue
Block a user