wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-27 18:25:42 +08:00
Code Issues Packages Projects Releases Wiki Activity

SimPEG.MT Mergify.

Browse Source 
Merge branch 'move2Simpeg' of https://github.com/simpeg/simpegmt into mt/dev

Conflicts:
	.gitignore
	.travis.yml
	LICENSE
	docs/conf.py
	docs/index.rst
	requirements.txt
	setup.py
This commit is contained in:
Rowan Cockett
2016-01-14 15:36:58 -08:00
parent e15913cf84 d77b393d42
commit 558e8879fe
28 changed files with 3595 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.travis.yml
+1
View File
@@ -19,6 +19,7 @@ env:
- TEST_DIR=tests/em/fdem/inverse/derivs
- TEST_DIR=tests/em/tdem
- TEST_DIR=tests/flow
- TEST_DIR=tests/mt
- TEST_DIR=tests/examples
- TEST_DIR=tests/em/fdem/inverse/adjoint
- TEST_DIR=tests/em/fdem/forward
SimPEG/MT/BaseMT.py
+116
View File
@@ -0,0 +1,116 @@
from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem
from SurveyMT import Survey, Data
from FieldsMT import BaseMTFields
class BaseMTProblem(BaseFDEMProblem):
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
Utils.setKwargs(self, **kwargs)
# Set the default pairs of the problem
surveyPair = Survey
dataPair = Data
fieldsPair = BaseMTFields
# Set the solver
Solver = SimpegSolver
solverOpts = {}
verbose = False
# Notes:
# Use the forward and devs from BaseFDEMProblem
# Might need to add more stuff here.
def Jvec(self, m, v, u=None):
"""
Function to calculate the data sensitivities dD/dm times a vector.
:param numpy.ndarray (nC, 1) - conductive model
:param numpy.ndarray (nC, 1) - random vector
:param MTfields object (optional) - MT fields object, if not given it is calculated
:rtype: MTdata object
:return: Data sensitivities wrt m
"""
# Calculate the fields
if u is None:
u = self.fields(m)
# Set current model
self.curModel = m
# Initiate the Jv object
Jv = self.dataPair(self.survey)
# 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
# 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 )
else:
du_dm = dA_duI * ( -dA_dm + dRHS_dm )
# Calculate the projection derivatives
for rx in src.rxList:
# Get the projection derivative
# v should be of size 2*nE (for 2 polarizations)
PDeriv_u = lambda t: rx.projectFieldsDeriv(src, self.mesh, u, t) # 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)
def Jtvec(self, m, v, u=None):
if u is None:
u = self.fields(m)
self.curModel = m
# Ensure v is a data object.
if not isinstance(v, self.dataPair):
v = self.dataPair(self.survey, v)
Jtv = np.zeros(m.size)
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, :]
for rx in src.rxList:
# 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
# Get the
dA_duIT = ATinv * PTv
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
elif real_or_imag == 'imag':
Jtv += -du_dmT.real
else:
raise Exception('Must be real or imag')
return Jtv
SimPEG/MT/Examples/simple3DfowardProblem.py
+41
View File
@@ -0,0 +1,41 @@
# Test script to use SimPEG.MT platform to forward model synthetic data.
# Import
import SimPEG as simpeg
from SimPEG import MT
import numpy as np
# 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])
# Setup the model
conds = [1e-2,1]
sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,[-1000,-1000,-400],[1000,1000,-200],conds)
sig[M.gridCC[:,2]>0] = 1e-8
sig[M.gridCC[:,2]<-600] = 1e-1
sigBG = np.zeros(M.nC) + conds[0]
sigBG[M.gridCC[:,2]>0] = 1e-8
## Setup the the survey object
# Receiver locations
rx_x, rx_y = np.meshgrid(np.arange(-500,501,50),np.arange(-500,501,50))
rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),np.zeros((np.prod(rx_x.shape),1))))
# Make a receiver list
rxList = []
for loc in rx_loc:
# NOTE: loc has to be a (1,3) np.ndarray otherwise errors accure
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
rxList.append(MT.RxMT(simpeg.mkvc(loc,2).T,rxType))
# Source list
srcList =[]
for freq in np.logspace(3,-3,7):
srcList.append(MT.SurveyMT.srcMT(freq,rxList))
# Survey MT
survey = MT.Survey(srcList)
## Setup the problem object
problem = MT.ProblemMT.MTProblem(M)
problem.pair(survey)
fields = problem.fields(sig,sigBG)
mtData = survey.projectFields(fields)
SimPEG/MT/FieldsMT.py
+348
View File
@@ -0,0 +1,348 @@
from SimPEG import Survey, Utils, Problem, np, sp, mkvc
from scipy.constants import mu_0
import sys
from numpy.lib import recfunctions as recFunc
from SimPEG.EM.Utils import omega
##############
### Fields ###
##############
class BaseMTFields(Problem.Fields):
"""Field Storage for a MT survey."""
knownFields = {}
dtype = complex
class Fields1D_e(BaseMTFields):
"""
Fields storage for the 1D MT solution.
"""
knownFields = {'e_1dSolution':'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']
}
def __init__(self,mesh,survey,**kwargs):
BaseMTFields.__init__(self,mesh,survey,**kwargs)
def _ePrimary(self, eSolution, srcList):
ePrimary = np.zeros_like(eSolution)
for i, src in enumerate(srcList):
ep = src.ePrimary(self.survey.prob)
if ep is not None:
ePrimary[:,i] = ep[:,-1]
return ePrimary
def _eSecondary(self, eSolution, srcList):
return eSolution
def _e(self, eSolution, srcList):
return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
def _eDeriv_u(self, src, v, adjoint = False):
return v
def _eDeriv_m(self, src, v, adjoint = False):
# assuming primary does not depend on the model
return None
def _bPrimary(self, eSolution, srcList):
bPrimary = np.zeros([self.survey.mesh.nE,eSolution.shape[1]], dtype = complex)
for i, src in enumerate(srcList):
bp = src.bPrimary(self.survey.prob)
if bp is not None:
bPrimary[:,i] += bp[:,-1]
return bPrimary
def _bSecondary(self, eSolution, srcList):
C = self.mesh.nodalGrad
b = (C * eSolution)
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(self, eSolution, srcList):
return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
def _bSecondaryDeriv_u(self, src, v, adjoint = False):
C = self.mesh.nodalGrad
if adjoint:
return - 1./(1j*omega(src.freq)) * (C.T * v)
return - 1./(1j*omega(src.freq)) * (C * v)
def _bSecondaryDeriv_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 _bDeriv_u(self, src, v, adjoint=False):
# Primary does not depend on u
return self._bSecondaryDeriv_u(src, v, adjoint)
def _bDeriv_m(self, src, v, adjoint=False):
# Assuming the primary does not depend on the model
return self._bSecondaryDeriv_m(src, v, adjoint)
def _fDeriv_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._bDeriv_u(src, v, adjoint)
# Return the stack
# This doesn't work...
return np.vstack((de_du,db_du))
def _fDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the fields object wrt m.
This function stacks the fields derivatives appropriately
"""
return None
class Fields3D_e(BaseMTFields):
"""
Fields storage for the 3D MT solution.
"""
# 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_pySolution':'E'}
aliasFields = {
'e_px' : ['e_pxSolution','E','_e_px'],
'e_pxPrimary' : ['e_pxSolution','E','_e_pxPrimary'],
'e_pxSecondary' : ['e_pxSolution','E','_e_pxSecondary'],
'e_py' : ['e_pySolution','E','_e_py'],
'e_pyPrimary' : ['e_pySolution','E','_e_pyPrimary'],
'e_pySecondary' : ['e_pySolution','E','_e_pySecondary'],
'b_px' : ['e_pxSolution','F','_b_px'],
'b_pxPrimary' : ['e_pxSolution','F','_b_pxPrimary'],
'b_pxSecondary' : ['e_pxSolution','F','_b_pxSecondary'],
'b_py' : ['e_pySolution','F','_b_py'],
'b_pyPrimary' : ['e_pySolution','F','_b_pyPrimary'],
'b_pySecondary' : ['e_pySolution','F','_b_pySecondary']
}
def __init__(self,mesh,survey,**kwargs):
BaseMTFields.__init__(self,mesh,survey,**kwargs)
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)
#NOTE: For e_p?Deriv_u,
# v has to be u(2*nE) long for the not adjoint and nE long for adjoint.
# Returns nE long for not adjoint and 2*nE long for adjoint
def _e_pxDeriv_u(self, src, v, adjoint = False):
'''
Takes the derivative of e_px wrt u
'''
if adjoint:
# adjoint: returns a 2*nE long vector with zero's for py
return np.vstack((v,np.zeros_like(v)))
# Not adjoint: return only the px part of the vector
return v[:len(v)/2]
def _e_pyDeriv_u(self, src, v, adjoint = False):
'''
Takes the derivative of e_py wrt u
'''
if adjoint:
# adjoint: returns a 2*nE long vector with zero's for px
return np.vstack((np.zeros_like(v),v))
# Not adjoint: return only the px part of the vector
return v[len(v)/2::]
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.nF,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[:,0]
return b_pxPrimary
def _b_pyPrimary(self, e_pySolution, srcList):
b_pyPrimary = np.zeros([self.survey.mesh.nF,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._b_pxPrimary(eSolution, srcList) + self._b_pxSecondary(eSolution, srcList)
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 = 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.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)
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_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_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
SimPEG/MT/Problem1D/Probs.py
+242
View File
@@ -0,0 +1,242 @@
from SimPEG.EM.Utils import omega
from SimPEG import mkvc
from scipy.constants import mu_0
from SimPEG.MT.BaseMT import BaseMTProblem
from SimPEG.MT.SurveyMT import Survey, Data
from SimPEG.MT.FieldsMT import Fields1D_e
from SimPEG.MT.Utils.MT1Danalytic import getEHfields
import numpy as np
import multiprocessing, sys, time
class eForm_psField(BaseMTProblem):
"""
A MT problem soving a e formulation and primary/secondary fields decomposion.
Solves the equation
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
_fieldType = 'e_1d'
_eqLocs = 'EF'
_sigmaPrimary = None
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
self.fieldsPair = Fields1D_e
# self._sigmaPrimary = sigmaPrimary
@property
def sigmaPrimary(self):
return self._sigmaPrimary
@sigmaPrimary.setter
def sigmaPrimary(self, val):
# Note: TODO add logic for val, make sure it is the correct size.
self._sigmaPrimary = val
def getA(self, freq):
"""
Function to get the A matrix.
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
Msig = self.mesh.getFaceInnerProduct(self.curModel.sigma)
# Note: need to use the code above since in the 1D problem I want
# e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
# Possible that _fieldType and _eqLocs can fix this
# Mmui = self.MfMui
# Msig = self.MeSigma
C = self.mesh.nodalGrad
# Make A
A = C.T*Mmui*C + 1j*omega(freq)*Msig
# Either return full or only the inner part of A
return A
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
The derivative of A wrt sigma
"""
dsig_dm = self.curModel.sigmaDeriv
MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
#
u_src = u['e_1dSolution']
dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
if adjoint:
return 1j * omega(freq) * ( dMf_dsig.T * v )
# Note: output has to be nN/nF, not nC/nE.
# v should be nC
return 1j * omega(freq) * ( dMf_dsig * v )
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
:return: RHS for 1 polarizations, primary fields
"""
# Get sources for the frequncy(polarizations)
Src = self.survey.getSrcByFreq(freq)[0]
S_e = Src.S_e(self)
return -1j * omega(freq) * S_e
def getRHSDeriv_m(self, freq, v, adjoint=False):
"""
The derivative of the RHS wrt sigma
"""
Src = self.survey.getSrcByFreq(freq)[0]
S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
return -1j * omega(freq) * S_eDeriv
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
'''
# Set the current model
self.curModel = m
F = Fields1D_e(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
# Store the fields
Src = self.survey.getSrcByFreq(freq)[0]
# NOTE: only store the e_solution(secondary), all other components calculated in the fields object
F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
# Note curl e = -iwb so b = -curl e /iw
# b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
# F[Src, 'b_1d'] = b[:,1]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
class eForm_TotalField(BaseMTProblem):
"""
A MT problem solving a e formulation and a Total bondary domain decompostion.
Solves the equation:
Math:
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
_fieldType = 'e'
_eqLocs = 'EF'
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
def getA(self, freq, full=False):
"""
Function to get the A matrix.
:param float freq: Frequency
:param logic full: Return full A or the inner part
:rtype: scipy.sparse.csr_matrix
:return: A
"""
Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
Msig = self.mesh.getFaceInnerProduct(self.curModel.sigma)
# Note: need to use the code above since in the 1D problem I want
# e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
# Possible that _fieldType and _eqLocs can fix this
# Mmui = self.MfMui
# Msig = self.MeSigma
C = self.mesh.nodalGrad
# Make A
A = C.T*Mmui*C + 1j*omega(freq)*Msig
# Either return full or only the inner part of A
if full:
return A
else:
return A[1:-1,1:-1]
def getADeriv(self, freq, u, v, adjoint=False):
sig = self.curTModel
dsig_dm = self.curTModelDeriv
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u)
if adjoint:
return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
:return: RHS for both polarizations, primary fields
"""
# Get sources for the frequency
# NOTE: Need to use the source information, doesn't really apply in 1D
src = self.survey.getSrcByFreq(freq)
# Get the full A
A = self.getA(freq,full=True)
# Define the outer part of the solution matrix
Aio = A[1:-1,[0,-1]]
Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
Etot = (Ed + Eu)
sourceAmp = 1.0
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
## Note: The analytic solution is derived with e^iwt
eBC = np.r_[Etot[0],Etot[-1]]
# The right hand side
return -Aio*eBC, eBC
def getRHSderiv(self, freq, backSigma, u, v, adjoint=False):
raise NotImplementedError('getRHSDeriv not implemented yet')
return None
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
:param np.ndarray (nC,) m_back: Background conductivity model
'''
self.curModel = m
# RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
F = Fields1D_e(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs, e_o = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
e_i = Ainv * rhs
e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
# Store the fields
Src = self.survey.getSrcByFreq(freq)
# NOTE: only store e fields
F[Src, 'e_1dSolution'] = e[:,0]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
SimPEG/MT/Problem1D/__init__.py
+1
View File
@@ -0,0 +1 @@
from Probs import eForm_TotalField, eForm_psField
SimPEG/MT/Problem2D/Probs.py
View File
SimPEG/MT/Problem2D/__init__.py
+1
View File
@@ -0,0 +1 @@
pass
SimPEG/MT/Problem3D/Probs.py
+264
View File
@@ -0,0 +1,264 @@
from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from SimPEG.MT.BaseMT import BaseMTProblem
from SimPEG.MT.SurveyMT import Survey, Data
from SimPEG.MT.FieldsMT import Fields3D_e
import multiprocessing, sys, time
class eForm_ps(BaseMTProblem):
"""
A MT problem solving a e formulation and a primary/secondary fields decompostion.
Solves the equation:
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
_fieldType = 'e'
_eqLocs = 'FE'
fieldsPair = Fields3D_e
_sigmaPrimary = None
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
@property
def sigmaPrimary(self):
return self._sigmaPrimary
@sigmaPrimary.setter
def sigmaPrimary(self, val):
# Note: TODO add logic for val, make sure it is the correct size.
self._sigmaPrimary = val
def getA(self, freq):
"""
Function to get the A matrix.
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
Mmui = self.MfMui
Msig = self.MeSigma
C = self.mesh.edgeCurl
return C.T*Mmui*C + 1j*omega(freq)*Msig
def getADeriv_m(self, freq, u, v, adjoint=False):
# 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:
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
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
:return: RHS for both polarizations, primary fields
"""
# Get sources for the frequncy(polarizations)
Src = self.survey.getSrcByFreq(freq)[0]
S_e = Src.S_e(self)
return -1j * omega(freq) * S_e
def getRHSDeriv_m(self, freq, v, adjoint=False):
"""
The derivative of the RHS with respect to sigma
"""
Src = self.survey.getSrcByFreq(freq)[0]
S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
return -1j * omega(freq) * S_eDeriv
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
'''
# Set the current model
self.curModel = m
F = Fields3D_e(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
# Store the fields
Src = self.survey.getSrcByFreq(freq)[0]
# Store the fieldss
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]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
class eForm_Tp(BaseMTProblem):
"""
A MT problem solving a e formulation and a total/primary fields decompostion.
Solves the equation
"""
_fieldType = 'e'
_eqLocs = 'FE'
fieldsPair = Fields3D_e
# Set new properties
# Background model
@property
def backModel(self):
"""
Sets the model, and removes dependent mass matrices.
"""
return getattr(self, '_backModel', None)
@backModel.setter
def backModel(self, value):
if value is self.backModel:
return # it is the same!
self._backModel = Models.Model(value, self.mapping)
for prop in self.deleteTheseOnModelUpdate:
if hasattr(self, prop):
delattr(self, prop)
@property
def MeSigmaBack(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MeSigmaBack', None) is None:
sigma = self.curModel
sigmaBG = self.backModel
self._MeSigmaBack = self.mesh.getEdgeInnerProduct(sigmaBG)
return self._MeSigmaBack
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
def getA(self, freq):
"""
Function to get the A matrix.
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
mui = self.MfMui
sig = self.MeSigma
C = self.mesh.edgeCurl
return C.T*mui*C + 1j*omega(freq)*sig
def getADeriv(self, freq, u, v, adjoint=False):
sig = self.curTModel
dsig_dm = self.curTModelDeriv
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u)
if adjoint:
return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
def getRHS(self, freq, backSigma):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:param numpy.ndarray (nC,) backSigma: Background conductivity model
:rtype: numpy.ndarray (nE, 2)
:return: one RHS for both polarizations
"""
# Get sources for the frequency
src = self.survey.getSources(freq)
# Make sure that there is 2 polarizations.
# assert len()
# Get the background electric fields
from SimPEG.MT.Sources import homo1DModelSource
eBG_bp = homo1DModelSource(self.mesh,freq,backSigma)
MeBack = self.MeSigmaBack
# Set up the A system
mui = self.MfMui
C = self.mesh.edgeCurl
Abg = C.T*mui*C + 1j*omega(freq)*MeBack
return Abg*eBG_bp, eBG_bp
def getRHSderiv(self, freq, backSigma, u, v, adjoint=False):
raise NotImplementedError('getRHSDeriv not implemented yet')
return None
def fields(self, m, m_back):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
:param np.ndarray (nC,) m_back: Background conductivity model
'''
self.curModel = m
self.backModel = m_back
# RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
F = Fields3D_e(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs, e_p = self.getRHS(freq,m_back)
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
e = e_s
# Store the fields
Src = self.survey.getSources(freq)
# Store the fieldss
F[Src, 'e_px'] = e[:,0]
F[Src, 'e_py'] = e[:,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]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
SimPEG/MT/Problem3D/__init__.py
+1
View File
@@ -0,0 +1 @@
from Probs import eForm_ps, eForm_Tp
SimPEG/MT/Sources/__init__.py
+1
View File
@@ -0,0 +1 @@
from backgroundModelSources import *
SimPEG/MT/Sources/backgroundModelSources.py
+178
View File
@@ -0,0 +1,178 @@
import SimPEG as simpeg, numpy as np
def homo1DModelSource(mesh,freq,sigma_1d):
'''
Function that calculates and return background fields
:param Simpeg mesh object mesh: Holds information on the discretization
:param float freq: The frequency to solve at
:param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model.
:rtype: numpy.ndarray (mesh.nE,2)
:return: eBG_bp, E fields for the background model at both polarizations.
'''
# import
from SimPEG.MT.Utils import get1DEfields
# Get a 1d solution for a halfspace background
if mesh.dim == 1:
mesh1d = mesh
elif mesh.dim == 2:
mesh1d = simpeg.Mesh.TensorMesh([mesh.hy],np.array([mesh.x0[1]]))
elif mesh.dim == 3:
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
# # Note: Everything is using e^iwt
e0_1d = get1DEfields(mesh1d,sigma_1d,freq)
if mesh.dim == 1:
eBG_px = simpeg.mkvc(e0_1d,2)
eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents.
elif mesh.dim == 2:
ex_px = np.zeros(mesh.vnEx,dtype=complex)
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
for i in np.arange(mesh.vnEx[0]):
ex_px[i,:] = -e0_1d
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px))
# Setup y (north) polarization (_py)
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
# Assign the source to ey_py
for i in np.arange(mesh.vnEy[0]):
ey_py[i,:] = e0_1d
# ey_py[1:-1,1:-1,1:-1] = 0
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
elif mesh.dim == 3:
# Setup x (east) polarization (_x)
ex_px = np.zeros(mesh.vnEx,dtype=complex)
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
ez_px = np.zeros((mesh.nEz,1),dtype=complex)
# Assign the source to ex_x
for i in np.arange(mesh.vnEx[0]):
for j in np.arange(mesh.vnEx[1]):
ex_px[i,j,:] = -e0_1d
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
# Setup y (north) polarization (_py)
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
# Assign the source to ey_py
for i in np.arange(mesh.vnEy[0]):
for j in np.arange(mesh.vnEy[1]):
ey_py[i,j,:] = e0_1d
# ey_py[1:-1,1:-1,1:-1] = 0
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
# Return the electric fields
eBG_bp = np.hstack((eBG_px,eBG_py))
return eBG_bp
def analytic1DModelSource(mesh,freq,sigma_1d):
'''
Function that calculates and return background fields
:param Simpeg mesh object mesh: Holds information on the discretization
:param float freq: The frequency to solve at
:param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model.
:rtype: numpy.ndarray (mesh.nE,2)
:return: eBG_bp, E fields for the background model at both polarizations.
'''
# import
from SimPEG.MT.Utils import getEHfields
# Get a 1d solution for a halfspace background
if mesh.dim == 1:
mesh1d = mesh
elif mesh.dim == 2:
mesh1d = simpeg.Mesh.TensorMesh([mesh.hy],np.array([mesh.x0[1]]))
elif mesh.dim == 3:
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
# # Note: Everything is using e^iwt
Eu, Ed, _, _ = getEHfields(mesh1d,sigma_1d,freq,mesh.vectorNz)
# Make the fields into a dictionary of location and the fields
e0_1d = Eu+Ed
E1dFieldDict = dict(zip(mesh.vectorNz,e0_1d))
if mesh.dim == 1:
eBG_px = simpeg.mkvc(e0_1d,2)
eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents.
elif mesh.dim == 2:
ex_px = np.zeros(mesh.vnEx,dtype=complex)
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
for i in np.arange(mesh.vnEx[0]):
ex_px[i,:] = -e0_1d
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px))
# Setup y (north) polarization (_py)
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
# Assign the source to ey_py
for i in np.arange(mesh.vnEy[0]):
ey_py[i,:] = e0_1d
# ey_py[1:-1,1:-1,1:-1] = 0
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
elif mesh.dim == 3:
# Setup x (east) polarization (_x)
ex_px = -np.array([E1dFieldDict[i] for i in mesh.gridEx[:,2]]).reshape(-1,1)
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
ez_px = np.zeros((mesh.nEz,1),dtype=complex)
# Construct the full fields
eBG_px = np.vstack((ex_px,ey_px,ez_px))
# Setup y (north) polarization (_py)
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
ey_py = np.array([E1dFieldDict[i] for i in mesh.gridEy[:,2]]).reshape(-1,1)
ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
# Construct the full fields
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
# Return the electric fields
eBG_bp = np.hstack((eBG_px,eBG_py))
return eBG_bp
# def homo3DModelSource(mesh,model,freq):
# '''
# Function that estimates 1D analytic background fields from a 3D model.
# :param Simpeg mesh object mesh: Holds information on the discretization
# :param float freq: The frequency to solve at
# :param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model.
# :rtype: numpy.ndarray (mesh.nE,2)
# :return: eBG_bp, E fields for the background model at both polarizations.
# '''
# if mesh.dim < 3:
# raise IOError('Input mesh has to have 3 dimensions.')
# # Get the locations
# a = mesh.gridCC[:,0:2].copy()
# unixy = np.unique(a.view(a.dtype.descr * a.shape[1])).view(float).reshape(-1,2)
# uniz = np.unique(mesh.gridCC[:,2])
# # # Note: Everything is using e^iwt
# # Need to loop thourgh the xy locations, assess the model and calculate the fields at the phusdo cell centers.
# # Then interpolate the cc fields to the edges.
# e0_1d = get1DEfields(mesh1d,sigma_1d,freq)
# elif mesh.dim == 3:
# # Setup x (east) polarization (_x)
# ex_px = np.zeros(mesh.vnEx,dtype=complex)
# ey_px = np.zeros((mesh.nEy,1),dtype=complex)
# ez_px = np.zeros((mesh.nEz,1),dtype=complex)
# # Assign the source to ex_x
# for i in np.arange(mesh.vnEx[0]):
# for j in np.arange(mesh.vnEx[1]):
# ex_px[i,j,:] = -e0_1d
# eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
# # Setup y (north) polarization (_py)
# ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
# ey_py = np.zeros(mesh.vnEy, dtype='complex128')
# ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
# # Assign the source to ey_py
# for i in np.arange(mesh.vnEy[0]):
# for j in np.arange(mesh.vnEy[1]):
# ey_py[i,j,:] = e0_1d
# # ey_py[1:-1,1:-1,1:-1] = 0
# eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
# # Return the electric fields
# eBG_bp = np.hstack((eBG_px,eBG_py))
# return eBG_bp
SimPEG/MT/SrcMT.py
+206
View File
@@ -0,0 +1,206 @@
from SimPEG import Survey, Utils, Problem, Maps, np, sp, mkvc
from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from Sources import homo1DModelSource
from Utils import rec2ndarr
from SurveyMT import Rx
import sys
#################
### Sources ###
#################
class BaseMTSrc(FDEMBaseSrc):
'''
Sources for the MT problem.
Use the SimPEG BaseSrc, since the source fields share properties with the transmitters.
:param float freq: The frequency of the source
:param list rxList: A list of receivers associated with the source
'''
freq = None #: Frequency (float)
rxPair = Rx
def __init__(self, rxList, freq):
self.freq = float(freq)
Survey.BaseSrc.__init__(self, rxList)
# 1D sources
class polxy_1DhomotD(BaseMTSrc):
"""
MT source for both polarizations (x and y) for the total Domain. It calculates fields calculated based on conditions on the boundary of the domain.
"""
def __init__(self, rxList, freq):
BaseMTSrc.__init__(self, rxList, freq)
# TODO: need to add the primary fields calc and source terms into the problem.
# Need to implement such that it works for all dims.
class polxy_1Dprimary(BaseMTSrc):
"""
MT source for both polarizations (x and y) given a 1D primary models. It assigns fields calculated from the 1D model
as fields in the full space of the problem.
"""
def __init__(self, rxList, freq):
# assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
self.sigma1d = None
BaseMTSrc.__init__(self, rxList, freq)
# Hidden property of the ePrimary
self._ePrimary = None
def ePrimary(self,problem):
# Get primary fields for both polarizations
if self.sigma1d is None:
# Set the sigma1d as the 1st column in the background model
if len(problem._sigmaPrimary) == problem.mesh.nC:
if problem.mesh.dim == 1:
self.sigma1d = problem.mesh.r(problem._sigmaPrimary,'CC','CC','M')[:]
elif problem.mesh.dim == 3:
self.sigma1d = problem.mesh.r(problem._sigmaPrimary,'CC','CC','M')[0,0,:]
# Or as the 1D model that matches the vertical cell number
elif len(problem._sigmaPrimary) == problem.mesh.nCz:
self.sigma1d = problem._sigmaPrimary
if self._ePrimary is None:
self._ePrimary = homo1DModelSource(problem.mesh,self.freq,self.sigma1d)
return self._ePrimary
def bPrimary(self,problem):
# Project ePrimary to bPrimary
# Satisfies the primary(background) field conditions
if problem.mesh.dim == 1:
C = problem.mesh.nodalGrad
elif problem.mesh.dim == 3:
C = problem.mesh.edgeCurl
bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
return bBG_bp
def S_e(self,problem):
"""
Get the electrical field source
"""
e_p = self.ePrimary(problem)
Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
sigma_p = Map_sigma_p._transform(self.sigma1d)
# Make mass matrix
# Note: M(sig) - M(sig_p) = M(sig - sig_p)
# Need to deal with the edge/face discrepencies between 1d/2d/3d
if problem.mesh.dim == 1:
Mesigma = problem.mesh.getFaceInnerProduct(problem.curModel.sigma)
Mesigma_p = problem.mesh.getFaceInnerProduct(sigma_p)
if problem.mesh.dim == 2:
pass
if problem.mesh.dim == 3:
Mesigma = problem.MeSigma
Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p)
return (Mesigma - Mesigma_p) * e_p
def S_eDeriv_m(self, problem, v, adjoint = False):
'''
Get the derivative of S_e wrt to sigma (m)
'''
# Need to deal with
if problem.mesh.dim == 1:
# Need to use the faceInnerProduct
MsigmaDeriv = problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(self.ePrimary(problem)[:,1]) * problem.curModel.sigmaDeriv
# MsigmaDeriv = ( MsigmaDeriv * MsigmaDeriv.T)**2
if problem.mesh.dim == 2:
pass
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(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
else:
# v should be nC size
return MsigmaDeriv * v
class polxy_3Dprimary(BaseMTSrc):
"""
MT source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
as fields in the full space of the problem.
"""
def __init__(self, rxList, freq):
# assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
self.sigmaPrimary = None
BaseMTSrc.__init__(self, rxList, freq)
# Hidden property of the ePrimary
self._ePrimary = None
def ePrimary(self,problem):
# Get primary fields for both polarizations
self.sigmaPrimary = problem._sigmaPrimary
if self._ePrimary is None:
self._ePrimary = homo3DModelSource(problem.mesh,self.sigmaPrimary,self.freq)
return self._ePrimary
def bPrimary(self,problem):
# Project ePrimary to bPrimary
# Satisfies the primary(background) field conditions
if problem.mesh.dim == 1:
C = problem.mesh.nodalGrad
elif problem.mesh.dim == 3:
C = problem.mesh.edgeCurl
bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
return bBG_bp
def S_e(self,problem):
"""
Get the electrical field source
"""
e_p = self.ePrimary(problem)
Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
sigma_p = Map_sigma_p._transform(self.sigma1d)
# Make mass matrix
# Note: M(sig) - M(sig_p) = M(sig - sig_p)
# Need to deal with the edge/face discrepencies between 1d/2d/3d
if problem.mesh.dim == 1:
Mesigma = problem.mesh.getFaceInnerProduct(problem.curModel.sigma)
Mesigma_p = problem.mesh.getFaceInnerProduct(sigma_p)
if problem.mesh.dim == 2:
pass
if problem.mesh.dim == 3:
Mesigma = problem.MeSigma
Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p)
return (Mesigma - Mesigma_p) * e_p
def S_eDeriv_m(self, problem, v, adjoint = False):
'''
Get the derivative of S_e wrt to sigma (m)
'''
# Need to deal with
if problem.mesh.dim == 1:
# Need to use the faceInnerProduct
MsigmaDeriv = problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(self.ePrimary(problem)[:,1]) * problem.curModel.sigmaDeriv
# MsigmaDeriv = ( MsigmaDeriv * MsigmaDeriv.T)**2
if problem.mesh.dim == 2:
pass
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(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
else:
# v should be nC size
return MsigmaDeriv * v
SimPEG/MT/SurveyMT.py
+489
View File
@@ -0,0 +1,489 @@
from SimPEG import Survey as SimPEGsurvey, Utils, Problem, Maps, np, sp, mkvc
from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from Sources import homo1DModelSource
from Utils import rec2ndarr
import sys
#################
### Receivers ###
#################
class Rx(SimPEGsurvey.BaseRx):
knownRxTypes = {
# 3D impedance
'zxxr':['Z3D', 'real'],
'zxyr':['Z3D', 'real'],
'zyxr':['Z3D', 'real'],
'zyyr':['Z3D', 'real'],
'zxxi':['Z3D', 'imag'],
'zxyi':['Z3D', 'imag'],
'zyxi':['Z3D', 'imag'],
'zyyi':['Z3D', 'imag'],
# 2D impedance
# TODO:
# 1D impedance
'z1dr':['Z1D', 'real'],
'z1di':['Z1D', 'imag'],
# Tipper
'tzxr':['T3D','real'],
'tzxi':['T3D','imag'],
'tzyr':['T3D','real'],
'tzyi':['T3D','imag']
}
# TODO: Have locs as single or double coordinates for both or numerator and denominator separately, respectively.
def __init__(self, locs, rxType):
SimPEGsurvey.BaseRx.__init__(self, locs, rxType)
@property
def projField(self):
"""
Field Type projection (e.g. e b ...)
:param str fracPos: Position of the field in the data ratio
"""
if 'numerator' in fracPos:
return self.knownRxTypes[self.rxType][0][0]
elif 'denominator' in fracPos:
return self.knownRxTypes[self.rxType][1][0]
else:
raise Exception('{s} is an unknown option. Use numerator or denominator.')
@property
def projGLoc(self):
"""
Grid Location projection (e.g. Ex Fy ...)
:param str fracPos: Position of the field in the data ratio
"""
if 'numerator' in fracPos:
return self.knownRxTypes[self.rxType][0][1]
elif 'denominator' in fracPos:
return self.knownRxTypes[self.rxType][0][1]
else:
raise Exception('{s} is an unknown option. Use numerator or denominator.')
@property
def projType(self):
"""
Receiver type for projection.
"""
return self.knownRxTypes[self.rxType][0]
@property
def projComp(self):
"""Component projection (real/imag)"""
return self.knownRxTypes[self.rxType][1]
def projectFields(self, src, mesh, f):
'''
Project the fields and return the correct data.
'''
if self.projType is 'Z1D':
Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx')
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
ex = Pex*mkvc(f[src,'e_1d'],2)
bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
# Note: Has a minus sign in front, to comply with quadrant calculations.
# Can be derived from zyx case for the 3D case.
f_part_complex = -ex/bx
# elif self.projType is 'Z2D':
elif self.projType is 'Z3D':
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
# Make the complex data
if 'zxx' in self.rxType:
f_part_complex = ( ex_px*hy_py - ex_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zxy' in self.rxType:
f_part_complex = (-ex_px*hx_py + ex_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zyx' in self.rxType:
f_part_complex = ( ey_px*hy_py - ey_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zyy' in self.rxType:
f_part_complex = (-ey_px*hx_py + ey_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
elif self.projType is 'T3D':
if self.locs.ndim == 3:
horLoc = self.locs[:,:,0]
vertLoc = self.locs[:,:,1]
else:
horLoc = self.locs
vertLoc = self.locs
Pbx = mesh.getInterpolationMat(horLoc,'Fx')
Pby = mesh.getInterpolationMat(horLoc,'Fy')
Pbz = mesh.getInterpolationMat(vertLoc,'Fz')
bx_px = Pbx*f[src,'b_px']
by_px = Pby*f[src,'b_px']
bz_px = Pbz*f[src,'b_px']
bx_py = Pbx*f[src,'b_py']
by_py = Pby*f[src,'b_py']
bz_py = Pbz*f[src,'b_py']
if 'tzx' in self.rxType:
f_part_complex = (- by_px*bz_py + by_py*bz_px)/(bx_px*by_py - bx_py*by_px)
if 'tzy' in self.rxType:
f_part_complex = ( bx_px*bz_py - bx_py*bz_px)/(bx_px*by_py - bx_py*by_px)
else:
NotImplementedError('Projection of {:s} receiver type is not implemented.'.format(self.rxType))
# Get the real or imag component
real_or_imag = self.projComp
f_part = getattr(f_part_complex, real_or_imag)
# print f_part
return f_part
def projectFieldsDeriv(self, src, mesh, f, v, adjoint=False):
"""
The derivative of the projection wrt u
:param MTsrc src: MT source
:param TensorMesh mesh: Mesh defining the topology of the problem
:param MTfields f: MT fields object of the source
:param numpy.ndarray v: Random vector of size
"""
real_or_imag = self.projComp
if not adjoint:
if self.projType is 'Z1D':
Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx')
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
# ex = Pex*mkvc(f[src,'e_1d'],2)
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
dP_de = -mkvc(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v),2)
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 been implement for 2D impedance tensor')
elif self.projType is 'Z3D':
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
# The size of the diratives should be nD,nU
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: 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
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 = 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 ))
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 = 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 = 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 = 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 )
# Extract the real number for the real/imag components.
Pv = np.array(getattr(PDeriv_complex, real_or_imag))
elif adjoint:
# Note: The v vector is real and the return should be complex
if self.projType is 'Z1D':
Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx')
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
# ex = Pex*mkvc(f[src,'e_1d'],2)
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
dP_deTv = -mkvc(Pex.T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*v,2)
db_duv = Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v
dP_dbTv = mkvc(f._bDeriv_u(src,db_duv,adjoint=True),2)
PDeriv_real = np.sum(np.hstack((dP_deTv,dP_dbTv)),1)
elif self.projType is 'Z2D':
raise NotImplementedError('Has not be implement for 2D impedance tensor')
elif self.projType is 'Z3D':
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.
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)
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
# Update the input vector
# Define shortcuts
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(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:
Zij = sDiag(aHd*( sDiag(ahy_py)*aex_px - sDiag(ahy_px)*aex_py))
ZijN_uV = lambda x: aex_px_u(sDiag(ahy_py)*x) + ahy_py_u(sDiag(aex_px)*x) - ahy_px_u(sDiag(aex_py)*x) - aex_py_u(sDiag(ahy_px)*x)
elif 'zxy' in self.rxType:
Zij = sDiag(aHd*(-sDiag(ahx_py)*aex_px + sDiag(ahx_px)*aex_py))
ZijN_uV = lambda x:-aex_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(aex_px)*x) + ahx_px_u(sDiag(aex_py)*x) + aex_py_u(sDiag(ahx_px)*x)
elif 'zyx' in self.rxType:
Zij = sDiag(aHd*( sDiag(ahy_py)*aey_px - sDiag(ahy_px)*aey_py))
ZijN_uV = lambda x: aey_px_u(sDiag(ahy_py)*x) + ahy_py_u(sDiag(aey_px)*x) - ahy_px_u(sDiag(aey_py)*x) - aey_py_u(sDiag(ahy_px)*x)
elif 'zyy' in self.rxType:
Zij = sDiag(aHd*(-sDiag(ahx_py)*aey_px + sDiag(ahx_px)*aey_py))
ZijN_uV = lambda x:-aey_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(aey_px)*x) + ahx_px_u(sDiag(aey_py)*x) + aey_py_u(sDiag(ahx_px)*x)
# Calculate the complex derivative
PDeriv_real = ZijN_uV(aHd*v) - aHd_uV(Zij.T*aHd*v)#
# NOTE: Need to reshape the output to go from 2*nU array to a (nU,2) matrix for each polarization
# PDeriv_real = np.hstack((mkvc(PDeriv_real[:len(PDeriv_real)/2],2),mkvc(PDeriv_real[len(PDeriv_real)/2::],2)))
PDeriv_real = PDeriv_real.reshape((2,mesh.nE)).T
# Extract the data
if real_or_imag == 'imag':
Pv = 1j*PDeriv_real
elif real_or_imag == 'real':
Pv = PDeriv_real.astype(complex)
return Pv
#################
### Survey ###
#################
class Survey(SimPEGsurvey.BaseSurvey):
"""
Survey class for MT. Contains all the sources associated with the survey.
:param list srcList: List of sources associated with the survey
"""
import SrcMT
srcPair = SrcMT.BaseMTSrc
def __init__(self, srcList, **kwargs):
# Sort these by frequency
self.srcList = srcList
SimPEGsurvey.BaseSurvey.__init__(self, **kwargs)
_freqDict = {}
for src in srcList:
if src.freq not in _freqDict:
_freqDict[src.freq] = []
_freqDict[src.freq] += [src]
self._freqDict = _freqDict
self._freqs = sorted([f for f in self._freqDict])
@property
def freqs(self):
"""Frequencies"""
return self._freqs
@property
def nFreq(self):
"""Number of frequencies"""
return len(self._freqDict)
# TODO: Rename to getSources
def getSrcByFreq(self, freq):
"""Returns the sources associated with a specific frequency."""
assert freq in self._freqDict, "The requested frequency is not in this survey."
return self._freqDict[freq]
def projectFields(self, u):
data = Data(self)
for src in self.srcList:
sys.stdout.flush()
for rx in src.rxList:
data[src, rx] = rx.projectFields(src, self.mesh, u)
return data
def projectFieldsDeriv(self, u):
raise Exception('Use Transmitters to project fields deriv.')
#################
### Data ###
#################
class Data(SimPEGsurvey.Data):
'''
Data class for MTdata
:param SimPEG survey object survey:
:param v vector with data
'''
def __init__(self, survey, v=None):
# Pass the variables to the "parent" method
SimPEGsurvey.Data.__init__(self, survey, v)
# # Import data
# @classmethod
# def fromEDIFiles():
# pass
def toRecArray(self,returnType='RealImag'):
'''
Function that returns a numpy.recarray for a SimpegMT impedance data object.
:param str returnType: Switches between returning a rec array where the impedance is split to real and imaginary ('RealImag') or is a complex ('Complex')
'''
# Define the record fields
dtRI = [('freq',float),('x',float),('y',float),('z',float),('zxxr',float),('zxxi',float),('zxyr',float),('zxyi',float),
('zyxr',float),('zyxi',float),('zyyr',float),('zyyi',float),('tzxr',float),('tzxi',float),('tzyr',float),('tzyi',float)]
dtCP = [('freq',float),('x',float),('y',float),('z',float),('zxx',complex),('zxy',complex),('zyx',complex),('zyy',complex),('tzx',complex),('tzy',complex)]
impList = ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']
for src in self.survey.srcList:
# Temp array for all the receivers of the source.
# Note: needs to be written more generally, using diffterent rxTypes and not all the data at the locaitons
# Assume the same locs for all RX
locs = src.rxList[0].locs
if locs.shape[1] == 1:
locs = np.hstack((np.array([[0.0,0.0]]),locs))
elif locs.shape[1] == 2:
locs = np.hstack((np.array([[0.0]]),locs))
tArrRec = np.concatenate((src.freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],12))),axis=1).view(dtRI)
# np.array([(src.freq,rx.locs[0,0],rx.locs[0,1],rx.locs[0,2],np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ) for rx in src.rxList],dtype=dtRI)
# Get the type and the value for the DataMT object as a list
typeList = [[rx.rxType.replace('z1d','zyx'),self[src,rx]] for rx in src.rxList]
# Insert the values to the temp array
for nr,(key,val) in enumerate(typeList):
tArrRec[key] = mkvc(val,2)
# Masked array
mArrRec = np.ma.MaskedArray(rec2ndarr(tArrRec),mask=np.isnan(rec2ndarr(tArrRec))).view(dtype=tArrRec.dtype)
# Unique freq and loc of the masked array
uniFLmarr = np.unique(mArrRec[['freq','x','y','z']]).copy()
try:
outTemp = recFunc.stack_arrays((outTemp,mArrRec))
#outTemp = np.concatenate((outTemp,dataBlock),axis=0)
except NameError as e:
outTemp = mArrRec
if 'RealImag' in returnType:
outArr = outTemp
elif 'Complex' in returnType:
# Add the real and imaginary to a complex number
outArr = np.empty(outTemp.shape,dtype=dtCP)
for comp in ['freq','x','y','z']:
outArr[comp] = outTemp[comp].copy()
for comp in ['zxx','zxy','zyx','zyy','tzx','tzy']:
outArr[comp] = outTemp[comp+'r'].copy() + 1j*outTemp[comp+'i'].copy()
else:
raise NotImplementedError('{:s} is not implemented, as to be RealImag or Complex.')
# Return
return outArr
@classmethod
def fromRecArray(cls, recArray, srcType='primary'):
"""
Class method that reads in a numpy record array to MTdata object.
Only imports the impedance data.
"""
if srcType=='primary':
src = SrcMT.src_polxy_1Dprimary
elif srcType=='total':
src = SrcMT.src_polxy_1DhomotD
else:
raise NotImplementedError('{:s} is not a valid source type for MTdata')
# Find all the frequencies in recArray
uniFreq = np.unique(recArray['freq'])
srcList = []
dataList = []
for freq in uniFreq:
# Initiate rxList
rxList = []
# Find that data for freq
dFreq = recArray[recArray['freq'] == freq].copy()
# Find the impedance rxTypes in the recArray.
rxTypes = [ comp for comp in recArray.dtype.names if (len(comp)==4 or len(comp)==3) and 'z' in comp]
for rxType in rxTypes:
# Find index of not nan values in rxType
notNaNind = ~np.isnan(dFreq[rxType])
if np.any(notNaNind): # Make sure that there is any data to add.
locs = rec2ndarr(dFreq[['x','y','z']][notNaNind].copy())
if dFreq[rxType].dtype.name in 'complex128':
rxList.append(Rx(locs,rxType+'r'))
dataList.append(dFreq[rxType][notNaNind].real.copy())
rxList.append(Rx(locs,rxType+'i'))
dataList.append(dFreq[rxType][notNaNind].imag.copy())
else:
rxList.append(Rx(locs,rxType))
dataList.append(dFreq[rxType][notNaNind].copy())
srcList.append(src(rxList,freq))
# Make a survey
survey = Survey(srcList)
dataVec = np.hstack(dataList)
return cls(survey,dataVec)
SimPEG/MT/Utils/MT1Danalytic.py
+111
View File
@@ -0,0 +1,111 @@
# Analytic solution of EM fields due to a plane wave
import numpy as np, SimPEG as simpeg
def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
'''Analytic solution for MT 1D layered earth. Returns E and H fields.
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
:param numpy.array, vector sigma: Physical property of conductivity corresponding with the mesh.
:param float, freq: Frequency to calculate data at.
:param numpy array, vector zd: location to calculate EH fields at
:param bollean, scaleUD: scales the output to be 1 at the top, increases numeracal stability.
Assumes a halfspace with the same conductive as the last cell below.
'''
# Note add an error check for the mesh and sigma are the same size.
# Constants: Assume constant
mu = 4*np.pi*1e-7*np.ones((m1d.nC+1))
eps = 8.85*1e-12*np.ones((m1d.nC+1))
# Angular freq
w = 2*np.pi*freq
# Add the halfspace value to the property
sig = np.concatenate((np.array([sigma[0]]),sigma))
# Calculate the wave number
k = np.sqrt(eps*mu*w**2-1j*mu*sig*w)
# Initiate the propagation matrix, in the order down up.
UDp = np.zeros((2,m1d.nC+1),dtype=complex)
UDp[1,0] = 1. # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
# Loop over all the layers, starting at the bottom layer
for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
# Calculate
yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer
zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer
# Build the propagation matrix
# Convert fields to down/up going components in layer below current layer
Pj1 = np.array([[1,1],[yp1,-yp1]])
# Convert fields to down/up going components in current layer
Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
# Propagate down and up components through the current layer
elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]])
# The down and up component in current layer.
UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
if scaleUD:
UDp[:,lnr+1::-1] = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
# Calculate the fields
Ed = np.empty((zd.size,),dtype=complex)
Eu = np.empty((zd.size,),dtype=complex)
Hd = np.empty((zd.size,),dtype=complex)
Hu = np.empty((zd.size,),dtype=complex)
# Loop over the layers and calculate the fields
# In the halfspace below the mesh
dup = m1d.vectorNx[0]
dind = dup >= zd
Ed[dind] = UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
Eu[dind] = UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
Hd[dind] = (k[0]/(w*mu[0]))*UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
Hu[dind] = -(k[0]/(w*mu[0]))*UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]):
dind = np.logical_and(dup >= zd, zd > dlow)
Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind]))
Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind]))
Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind]))
# Return return the fields
return Ed, Eu, Hd, Hu
def getImpedance(m1d,sigma,freq):
"""Analytic solution for MT 1D layered earth. Returns the impedance at the surface.
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
:param numpy.array, vector sigma: Physical property corresponding with the mesh.
:param numpy.array, vector freq: Frequencies to calculate data at.
"""
# Define constants
mu0 = 4*np.pi*1e-7
eps0 = 8.85e-12
# Initiate the impedances
Z1d = np.empty(len(freq) , dtype='complex')
h = m1d.hx #vectorNx[:-1]
# Start the process
for nrFr, fr in enumerate(freq):
om = 2*np.pi*fr
Zall = np.empty(len(h)+1,dtype='complex')
# Calculate the impedance for the bottom layer
Zall[0] = (mu0*om)/np.sqrt(mu0*eps0*(om)**2 - 1j*mu0*sigma[0]*om)
for nr,hi in enumerate(h):
# Calculate the wave number
# print nr,sigma[nr]
k = np.sqrt(mu0*eps0*om**2 - 1j*mu0*sigma[nr]*om)
Z = (mu0*om)/k
Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))
#pdb.set_trace()
Z1d[nrFr] = Zall[-1]
return Z1d
SimPEG/MT/Utils/MT1Dsolutions.py
+45
View File
@@ -0,0 +1,45 @@
import numpy as np, SimPEG as simpeg
from MT1Danalytic import getEHfields
from scipy.constants import mu_0
def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
"""Function to get 1D electrical fields"""
# Get the gradient
G = m1d.nodalGrad
# Mass matrices
# Magnetic permeability
Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0))
# Conductivity
Msig = m1d.getFaceInnerProduct(sigma)
# Set up the solution matrix
A = G.T*Mmu*G + 1j*2.*np.pi*freq*Msig
# Define the inner part of the solution matrix
Aii = A[1:-1,1:-1]
# Define the outer part of the solution matrix
Aio = A[1:-1,[0,-1]]
# Set the boundary conditions
Ed, Eu, Hd, Hu = getEHfields(m1d,sigma,freq,m1d.vectorNx)
Etot = (Ed + Eu)
if sourceAmp is not None:
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
## Note: The analytic solution is derived with e^iwt
bc = np.r_[Etot[0],Etot[-1]]
# The right hand side
rhs = Aio*bc
# Solve the system
Aii_inv = simpeg.Solver(Aii)
eii = Aii_inv*rhs
# Assign the boundary conditions
e = np.r_[bc[0],eii,bc[1]]
# Return the electrical fields
return e
if __name__ == '__main__':
hz = [(100.,18)]
M = simpeg.Mesh.TensorMesh([hz],'C')
sig = np.zeros(M.nC) + 1e-8
sig[M.vectorCCx<=0] = sigHalf
SimPEG/MT/Utils/__init__.py
+4
View File
@@ -0,0 +1,4 @@
from MT1Dsolutions import * # Add the names of the functions
from MT1Danalytic import *
from dataUtils import *
from ediFilesUtils import *
SimPEG/MT/Utils/dataUtils.py
+245
View File
@@ -0,0 +1,245 @@
# Utils used for the data,
import numpy as np, matplotlib.pyplot as plt, sys
import SimPEG as simpeg
import numpy.lib.recfunctions as recFunc
from scipy.constants import mu_0
from scipy import interpolate as sciint
def getAppRes(MTdata):
# Make impedance
zList = []
for src in MTdata.survey.srcList:
zc = [src.freq]
for rx in src.rxList:
if 'i' in rx.rxType:
m=1j
else:
m = 1
zc.append(m*MTdata[src,rx])
zList.append(zc)
return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
def rotateData(MTdata,rotAngle):
'''
Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
'''
recData = MTdata.toRecArray('Complex')
impData = rec2ndarr(recData[['zxx','zxy','zyx','zyy']],complex)
# Make the rotation matrix
# c,s,zxx,zxy,zyx,zyy = sympy.symbols('c,s,zxx,zxy,zyx,zyy')
# rotM = sympy.Matrix([[c,-s],[s, c]])
# zM = sympy.Matrix([[zxx,zxy],[zyx,zyy]])
# rotM*zM*rotM.T
# [c*(c*zxx - s*zyx) - s*(c*zxy - s*zyy), c*(c*zxy - s*zyy) + s*(c*zxx - s*zyx)],
# [c*(c*zyx + s*zxx) - s*(c*zyy + s*zxy), c*(c*zyy + s*zxy) + s*(c*zyx + s*zxx)]])
s = np.sin(-np.deg2rad(rotAngle))
c = np.cos(-np.deg2rad(rotAngle))
rotMat = np.array([[c,-s],[s,c]])
rotData = (rotMat.dot(impData.reshape(-1,2,2).dot(rotMat.T))).transpose(1,0,2).reshape(-1,4)
outRec = recData.copy()
for nr,comp in enumerate(['zxx','zxy','zyx','zyy']):
outRec[comp] = rotData[:,nr]
from SimPEG import MT
return MT.Data.fromRecArray(outRec)
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
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)))
def makeAnalyticSolution(mesh,model,elev,freqs):
from SimPEG import MT
data1D = []
for freq in freqs:
anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
anaE = anaEd+anaEu
anaH = anaHd+anaHu
anaZ = anaE/anaH
# Add to the list
data1D.append((freq,0,0,elev,anaZ[0]))
dataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zyx',complex)])
return dataRec
def plotMT1DModelData(problem,models,symList=None):
from SimPEG import MT
# Setup the figure
fontSize = 15
fig = plt.figure(figsize=[9,7])
axM = fig.add_axes([0.075,.1,.25,.875])
axM.set_xlabel('Resistivity [Ohm*m]',fontsize=fontSize)
axM.set_xlim(1e-1,1e5)
axM.set_ylim(-10000,5000)
axM.set_ylabel('Depth [km]',fontsize=fontSize)
axR = fig.add_axes([0.42,.575,.5,.4])
axR.set_xscale('log')
axR.set_yscale('log')
axR.invert_xaxis()
# axR.set_xlabel('Frequency [Hz]')
axR.set_ylabel('Apparent resistivity [Ohm m]',fontsize=fontSize)
axP = fig.add_axes([0.42,.1,.5,.4])
axP.set_xscale('log')
axP.invert_xaxis()
axP.set_ylim(0,90)
axP.set_xlabel('Frequency [Hz]',fontsize=fontSize)
axP.set_ylabel('Apparent phase [deg]',fontsize=fontSize)
# if not symList:
# symList = ['x']*len(models)
import plotDataTypes as pDt
# Loop through the models.
modelList = [problem.survey.mtrue]
modelList.extend(models)
if False:
modelList = [problem.mapping.sigmaMap*mod for mod in modelList]
for nr, model in enumerate(modelList):
# Calculate the data
if nr==0:
data1D = problem.dataPair(problem.survey,problem.survey.dobs).toRecArray('Complex')
else:
data1D = problem.dataPair(problem.survey,problem.survey.dpred(model)).toRecArray('Complex')
# Plot the data and the model
colRat = nr/((len(modelList)-1.999)*1.)
if colRat > 1.:
col = 'k'
else:
col = plt.cm.seismic(1-colRat)
# The model - make the pts to plot
meshPts = np.concatenate((problem.mesh.gridN[0:1],np.kron(problem.mesh.gridN[1::],np.ones(2))[:-1]))
modelPts = np.kron(1./(problem.mapping.sigmaMap*model),np.ones(2,))
axM.semilogx(modelPts,meshPts,color=col)
## Data
# Appres
pDt.plotIsoStaImpedance(axR,np.array([0,0]),data1D,'zyx','res',pColor=col)
# Appphs
pDt.plotIsoStaImpedance(axP,np.array([0,0]),data1D,'zyx','phs',pColor=col)
try:
allData = np.concatenate((allData,simpeg.mkvc(data1D['zyx'],2)),1)
except:
allData = simpeg.mkvc(data1D['zyx'],2)
freq = simpeg.mkvc(data1D['freq'],2)
res, phs = appResPhs(freq,allData)
stdCol = 'gray'
axRtw = axR.twinx()
axRtw.set_ylabel('Std of log10',color=stdCol)
[(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
axPtw = axP.twinx()
axPtw.set_ylabel('Std ',color=stdCol)
[t.set_color(stdCol) for t in axPtw.get_yticklabels()]
axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
# Fix labels and ticks
yMtick = [l/1000 for l in axM.get_yticks().tolist()]
axM.set_yticklabels(yMtick)
[ l.set_rotation(90) for l in axM.get_yticklabels()]
[ l.set_rotation(90) for l in axR.get_yticklabels()]
[(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
[t.set_color(stdCol) for t in axPtw.get_yticklabels()]
for ax in [axM,axR,axP]:
ax.xaxis.set_tick_params(labelsize=fontSize)
ax.yaxis.set_tick_params(labelsize=fontSize)
return fig
def printTime():
import time
print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
from SimPEG import MT
# Find the unique locations
# Need to find the locations
recDataTemp = MTdata.toRecArray()
# Check if survey.std has been assigned.
## NEED TO: write this...
# Calculte and add the DET of the tensor to the recArray
if 'det' in rxType3D:
Zon = (recDataTemp['zxxr']+1j*recDataTemp['zxxi'])*(recDataTemp['zyyr']+1j*recDataTemp['zyyi'])
Zoff = (recDataTemp['zxyr']+1j*recDataTemp['zxyi'])*(recDataTemp['zyxr']+1j*recDataTemp['zyxi'])
det = np.sqrt(Zon.data - Zoff.data)
recData = recFunc.append_fields(recDataTemp,['zdetr','zdeti'],[det.real,det.imag] )
else:
recData = recDataTemp
uniLocs = rec2ndarr(np.unique(recData[['x','y','z']])).data
mtData1DList = []
if 'zxy' in rxType3D:
corr = -1 # Shift the data to comply with the quadtrature of the 1d problem
else:
corr = 1
for loc in uniLocs:
# Make the receiver list
rx1DList = []
for rxType in ['z1dr','z1di']:
rx1DList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
# Source list
locrecData = recData[np.sqrt(np.sum( (rec2ndarr(recData[['x','y','z']]).data - loc )**2,axis=1)) < 1e-5]
dat1DList = []
src1DList = []
for freq in locrecData['freq']:
src1DList.append(MT.SrcMT.src_polxy_1Dprimary(rx1DList,freq))
for comp in ['r','i']:
dat1DList.append( corr * locrecData[rxType3D+comp][locrecData['freq']== freq].data )
# Make the survey
sur1D = MT.Survey(src1DList)
# Make the data
dataVec = np.hstack(dat1DList)
dat1D = MT.Data(sur1D,dataVec)
sur1D.dobs = dataVec
# Need to take MTdata.survey.std and split it as well.
std=0.05
sur1D.std = np.abs(sur1D.dobs*std) #+ 0.01*np.linalg.norm(sur1D.dobs)
mtData1DList.append(dat1D)
# Return the the list of data.
return mtData1DList
def resampleMTdataAtFreq(MTdata,freqs):
"""
Function to resample MTdata at set of frequencies
"""
from SimPEG import MT
# Make a rec array
MTrec = MTdata.toRecArray().data
# Find unique locations
uniLoc = np.unique(MTrec[['x','y','z']])
uniFreq = MTdata.survey.freqs
# Get the comps
dNames = MTrec.dtype
# Loop over all the locations and interpolate
for loc in uniLoc:
# Find the index of the station
ind = np.sqrt(np.sum((rec2ndarr(MTrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
# Make a temporary recArray and interpolate all the components
tArrRec = np.concatenate((simpeg.mkvc(freqs,2),np.ones((len(freqs),1))*rec2ndarr(loc),np.nan*np.ones((len(freqs),12))),axis=1).view(dNames)
for comp in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
int1d = sciint.interp1d(MTrec[ind]['freq'],MTrec[ind][comp],bounds_error=False)
tArrRec[comp] = simpeg.mkvc(int1d(freqs),2)
# Join together
try:
outRecArr = recFunc.stack_arrays((outRecArr,tArrRec))
except NameError as e:
outRecArr = tArrRec
# Make the MTdata and return
return MT.Data.fromRecArray(outRecArr)
SimPEG/MT/Utils/ediFilesUtils.py
+175
View File
@@ -0,0 +1,175 @@
# Functions to import and export MT EDI files.
from SimPEG import mkvc
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from SimPEG.MT.Utils.dataUtils import rec2ndarr
# Import modules
import numpy as np
import os, sys, re
try:
import osr
except ImportError as e:
print 'Could not import osr, missing the gdal package'
pass
class EDIimporter:
"""
A class to import EDIfiles.
"""
_impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit
_impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)
# Properties
filesList = None
comps = None
# Hidden properties
_outEPSG = None
_2out = None
def __init__(self, EDIfilesList, compList=None, outEPSG=None):
# Set the fileList
self.filesList = EDIfilesList
# Set the components to import
if compList is None:
self.comps = ['ZXXR','ZXYR','ZYXR','ZYYR','ZXXI','ZXYI','ZYXI','ZYYI','ZXX.VAR','ZXY.VAR','ZYX.VAR','ZYY.VAR']
else:
self.comps = compList
if outEPSG is not None:
self._outEPSG = outEPSG
def __call__(self,comps=None):
if comps is None:
return self._data
return self._data[comps]
def importFiles(self):
"""
Function to import EDI files into a object.
"""
# Constants that are needed for convertion of units
# Temp lists
tmpStaList = []
tmpCompList = ['freq','x','y','z']
tmpCompList.extend(self.comps)
# Make the outarray
dtRI = [(compS.lower().replace('.',''),float) for compS in tmpCompList]
# Loop through all the files
for nrEDI, EDIfile in enumerate(self.filesList):
# Read the file into a list of the lines
with open(EDIfile,'r') as fid:
EDIlines = fid.readlines()
# Find the location
latD, longD, elevM = _findLatLong(EDIlines)
# Transfrom coordinates
transCoord = self._transfromPoints(longD,latD)
# Extract the name of the file (station)
EDIname = EDIfile.split(os.sep)[-1].split('.')[0]
# Arrange the data
staList = [EDIname, EDIfile, transCoord[0], transCoord[1], elevM[0]]
# Add to the station list
tmpStaList.extend(staList)
# Read the frequency data
freq = _findEDIcomp('>FREQ',EDIlines)
# Make the temporary rec array.
tArrRec = ( np.nan*np.ones( (len(freq),len(dtRI)) ) ).view(dtRI) #np.concatenate((freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],8))),axis=1).view(dtRI)
# Add data to the array
tArrRec['freq'] = mkvc(freq,2)
tArrRec['x'] = mkvc(np.ones((len(freq),1))*transCoord[0],2)
tArrRec['y'] = mkvc(np.ones((len(freq),1))*transCoord[1],2)
tArrRec['z'] = mkvc(np.ones((len(freq),1))*elevM[0],2)
for comp in self.comps:
# Deal with converting units of the impedance tensor
if 'Z' in comp:
unitConvert = self._impUnitEDI2SI
else:
unitConvert = 1
# Rotate the data since EDI x is *north, y *east but Simpeg uses x *east, y *north (* means internal reference frame)
key = [comp.lower().replace('.','').replace(s,t) for s,t in [['xx','yy'],['xy','yx'],['yx','xy'],['yy','xx']] if s in comp.lower()][0]
tArrRec[key] = mkvc(unitConvert*_findEDIcomp('>'+comp,EDIlines),2)
# Make a masked array
mArrRec = np.ma.MaskedArray(rec2ndarr(tArrRec),mask=np.isnan(rec2ndarr(tArrRec))).view(dtype=tArrRec.dtype)
try:
outTemp = recFunc.stack_arrays((outTemp,mArrRec))
except NameError as e:
outTemp = mArrRec
# Assign the data
self._data = outTemp
# % Assign the data to the obj
# nOutData=length(obj.data);
# obj.data(nOutData+1:nOutData+length(TEMP.data),:) = TEMP.data;
def _transfromPoints(self,longD,latD):
# Coordinates convertor
if self._2out is None:
src = osr.SpatialReference()
src.ImportFromEPSG(4326)
out = osr.SpatialReference()
if self._outEPSG is None:
# Find the UTM EPSG number
Nnr = 700 if latD < 0.0 else 600
utmZ = int(1+(longD+180.0)/6.0)
self._outEPSG = 32000 + Nnr + utmZ
out.ImportFromEPSG(self._outEPSG)
self._2out = osr.CoordinateTransformation(src,out)
# Return the transfrom
return self._2out.TransformPoint(longD,latD)
# Hidden functions
def _findLatLong(fileLines):
latDMS = np.array(fileLines[_findLine('LAT=',fileLines)[0]].split('=')[1].split()[0].split(':'),float)
longDMS = np.array(fileLines[_findLine('LONG=',fileLines)[0]].split('=')[1].split()[0].split(':'),float)
elevM = np.array([fileLines[_findLine('ELEV=',fileLines)[0]].split('=')[1].split()[0]],float)
# Convert to D.ddddd values
latS = np.sign(latDMS[0])
longS = np.sign(longDMS[0])
latD = latDMS[0] + latS*latDMS[1]/60 + latS*latDMS[2]/3600
longD = longDMS[0] + longS*longDMS[1]/60 + longS*longDMS[2]/3600
return latD, longD, elevM
def _findLine(comp,fileLines):
""" Find a line number in the file"""
# Line counter
c = 0
# List of indices for found lines
found = []
# Loop through all the lines
for line in fileLines:
if comp in line:
# Append if found
found.append(c)
# Increse the counter
c += 1
# Return the found indices
return found
def _findEDIcomp(comp,fileLines,dt=float):
"""
Extract the data vector.
Returns a list of the data.
"""
# Find the data
headLine, indHead = [(st,nr) for nr,st in enumerate(fileLines) if re.search(comp,st)][0]
# Extract the data
nrVec = int(headLine.split()[-1])
c = 0
dataList = []
while c < nrVec:
indHead += 1
dataList.extend(fileLines[indHead].split())
c = len(dataList)
return np.array(dataList,dt)
SimPEG/MT/Utils/plotDataTypes.py
+416
View File
@@ -0,0 +1,416 @@
from matplotlib import pyplot as plt, colors, numpy as np
def rec2nd(structArray):
""" Converts a structured/record array to ndarray to do operations on."""
return structArray.view((np.float,len(structArray.dtype.names)))
def plotIsoFreqNSimpedance(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=True,contour=True):
indUniFreq = np.where(freq==array['freq'])
x, y = array['x'][indUniFreq],array['y'][indUniFreq]
if par == 'abs':
zPlot = np.abs(array[flag][indUniFreq])
cmap = plt.get_cmap('OrRd_r')#seismic')
level = np.logspace(0,-5,31)
clevel = np.logspace(0,-4,5)
plotNorm = colors.LogNorm()
elif par == 'real':
zPlot = np.real(array[flag][indUniFreq])
cmap = plt.get_cmap('RdYlBu')
if cLevel:
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
else:
plotNorm = colors.Normalize()
elif par == 'imag':
zPlot = np.imag(array[flag][indUniFreq])
cmap = plt.get_cmap('RdYlBu')
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
if cLevel:
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
elif colorNorm=='Lin':
plotNorm = colors.Normalize()
if contour:
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
else:
uniX,uniY = np.unique(x),np.unique(y)
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
if colorbar:
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
ax.set_title(flag+' '+par,fontsize=8)
return cs
def plotIsoFreqNSDiff(ax,freq,arrayList,flag,par='abs',colorbar=True,cLevel=True,mask=None,contourLine=True,useLog=False):
indUniFreq0 = np.where(freq==arrayList[0]['freq'])
indUniFreq1 = np.where(freq==arrayList[1]['freq'])
seicmap = plt.get_cmap('RdYlBu')#seismic')
x, y = arrayList[0]['x'][indUniFreq0],arrayList[0]['y'][indUniFreq0]
if par == 'abs':
if useLog:
zPlot = (np.log10(np.abs(arrayList[0][flag][indUniFreq0])) - np.log10(np.abs(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs(arrayList[1][flag][indUniFreq1]))
else:
zPlot = (np.abs(arrayList[0][flag][indUniFreq0]) - np.abs(arrayList[1][flag][indUniFreq1]))/np.abs(arrayList[1][flag][indUniFreq1])
if mask:
maskInd = np.logical_or(np.abs(arrayList[0][flag][indUniFreq0])< 1e-3,np.abs(arrayList[1][flag][indUniFreq1]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
elif par == 'real':
if useLog:
zPlot = (np.log10(np.real(arrayList[0][flag][indUniFreq0])) -np.log10(np.real(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.real(arrayList[1][flag][indUniFreq1]))))
else:
zPlot = (np.real(arrayList[0][flag][indUniFreq0]) -np.real(arrayList[1][flag][indUniFreq1]))/np.abs((np.real(arrayList[1][flag][indUniFreq1])))
if mask:
maskInd = np.logical_or(np.abs(np.real(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.real(arrayList[1][flag][indUniFreq1])) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
elif par == 'imag':
if useLog:
zPlot = (np.log10(np.imag(arrayList[0][flag][indUniFreq0])) -np.log10(np.imag(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.imag(arrayList[1][flag][indUniFreq1]))))
else:
zPlot = (np.imag(arrayList[0][flag][indUniFreq0]) -np.imag(arrayList[1][flag][indUniFreq1]))/np.abs((np.imag(arrayList[1][flag][indUniFreq1])))
if mask:
maskInd = np.logical_or(np.abs(np.imag(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.imag(arrayList[1][flag][indUniFreq1])) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
cs = ax.tricontourf(x,y,zPlot*100,levels=level*100,cmap=seicmap,extend='both') #,norm=colors.SymLogNorm(1e-2,linscale=2))
if contourLine:
csl = ax.tricontour(x,y,zPlot*100,levels=clevel*100,colors='k')
plt.clabel(csl, fontsize=7, inline=1,fmt='%1.1e',inline_spacing=10)
if colorbar:
cb = plt.colorbar(cs,cax=ax.cax,ticks=clevel*100,format='%1.1e')
for t in cb.ax.get_yticklabels():
t.set_rotation(60)
t.set_fontsize(8)
ax.set_title(flag+' '+par,fontsize=8)
def plotIsoFreqNStipper(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=True,contour=True):
indUniFreq = np.where(freq==array['freq'])
x, y = array['x'][indUniFreq],array['y'][indUniFreq]
if par == 'abs':
cmap = plt.get_cmap('OrRd_r')#seismic')
zPlot = np.abs(array[flag][indUniFreq])
if cLevel:
level = np.logspace(-4,0,33)
clevel = np.logspace(-4,0,5)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
if colorNorm=='SymLog':
plotNorm = colors.LogNorm()
else:
plotNorm = colors.Normalize()
elif par == 'real':
cmap = plt.get_cmap('RdYlBu')
zPlot = np.real(array[flag][indUniFreq])
if cLevel:
level = np.concatenate((-np.logspace(0,-4,33),np.logspace(-4,0,33)))
clevel = np.concatenate((-np.logspace(0,-4,5),np.logspace(-4,0,5)))
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(1e-4,linscale=2)
else:
plotNorm = colors.Normalize()
elif par == 'imag':
cmap = plt.get_cmap('RdYlBu')
zPlot = np.imag(array[flag][indUniFreq])
if cLevel:
level = np.concatenate((-np.logspace(0,-4,33),np.logspace(-4,0,33)))
clevel = np.concatenate((-np.logspace(0,-4,5),np.logspace(-4,0,5)))
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(1e-4,linscale=2)
else:
plotNorm = colors.Normalize()
if contour:
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
else:
uniX,uniY = np.unique(x),np.unique(y)
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),levels=level,cmap=cmap,norm=plotNorm,edgecolors='k', linewidths=0.5)
if colorbar:
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
ax.set_title(flag+' '+par,fontsize=8)
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
appResFact = 1/(8*np.pi**2*10**(-7))
treshold = 1.0 # 1 meter
indUniSta = np.sqrt(np.sum((rec2nd(array[['x','y']])-loc)**2,axis=1)) < treshold
freq = array['freq'][indUniSta]
if par == 'abs':
zPlot = np.abs(array[flag][indUniSta])
elif par == 'real':
zPlot = np.real(array[flag][indUniSta])
elif par == 'imag':
zPlot = np.imag(array[flag][indUniSta])
elif par == 'res':
zPlot = (appResFact/freq)*np.abs(array[flag][indUniSta])**2
elif par == 'phs':
zPlot = np.arctan2(array[flag][indUniSta].imag,array[flag][indUniSta].real)*(180/np.pi)
if not pColor:
if 'xx' in flag:
lab = 'XX'
pColor = 'g'
elif 'xy' in flag:
lab = 'XY'
pColor = 'r'
elif 'yx' in flag:
lab = 'YX'
pColor = 'b'
elif 'yy' in flag:
lab = 'YY'
pColor = 'y'
ax.plot(freq,zPlot,color=pColor,marker=pSym,label=flag)
def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colorNorm='None',cLevel=None,contour=True):
indSect = np.where(sectDict.values()[0]==array[sectDict.keys()[0]])
# Define the plot axes
if 'x' in sectDict.keys()[0]:
x = array['y'][indSect]
else:
x = array['x'][indSect]
y = array['freq'][indSect]
if par == 'abs':
zPlot = np.abs(array[flag][indSect])
cmap = plt.get_cmap('OrRd_r')#seismic')
if cLevel:
level = np.logspace(0,-5,31,endpoint=True)
clevel = np.logspace(0,-4,5,endpoint=True)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100,endpoint=True)
clevel = np.linspace(zPlot.min(),zPlot.max(),10,endpoint=True)
elif par == 'ares':
zPlot = np.abs(array[flag][indSect])**2/(8*np.pi**2*10**(-7)*array['freq'][indSect])
cmap = plt.get_cmap('RdYlBu')#seismic)
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
else:
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
level = np.logspace(zMin,zMax,(zMax-zMin)*8+1,endpoint=True)
clevel = np.logspace(zMin,zMax,(zMax-zMin)*2+1,endpoint=True)
plotNorm = colors.LogNorm()
elif par == 'aphs':
zPlot = np.arctan2(array[flag][indSect].imag,array[flag][indSect].real)*(180/np.pi)
cmap = plt.get_cmap('RdYlBu')#seismic)
if cLevel:
zMax = cLevel[1]
zMin = cLevel[0]
else:
zMax = (np.ceil(zPlot).max())
zMin = (np.floor(zPlot).min())
level = np.arange(zMin,zMax+.1,1)
clevel = np.arange(zMin,zMax+.1,10)
plotNorm = colors.Normalize()
elif par == 'real':
zPlot = np.real(array[flag][indSect])
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
else:
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif par == 'imag':
zPlot = np.imag(array[flag][indSect])
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
else:
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif colorNorm=='Lin':
plotNorm = colors.Normalize()
elif colorNorm=='Log':
plotNorm = colors.LogNorm()
if contour:
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
else:
uniX,uniY = np.unique(x),np.unique(y)
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
if colorbar:
csB = plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
# csB.on_mappable_changed(cs)
ax.set_title(flag+' '+par,fontsize=8)
return cs, csB
return cs,None
def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=None,contour=True,mask=None,useLog=False):
def sortInArr(arr):
return np.sort(arr,order=['freq','x','y','z'])
# Find the index for the slice
indSect0 = np.where(sectDict.values()[0]==arrayList[0][sectDict.keys()[0]])
indSect1 = np.where(sectDict.values()[0]==arrayList[1][sectDict.keys()[0]])
# Extract and sort the mats
arr0 = sortInArr(arrayList[0][indSect0])
arr1 = sortInArr(arrayList[1][indSect1])
# Define the plot axes
if 'x' in sectDict.keys()[0]:
x0 = arr0['y']
x1 = arr1['y']
else:
x0 = arr0['x']
x1 = arr1['x']
y0 = arr0['freq']
y1 = arr1['freq']
if par == 'abs':
if useLog:
zPlot = (np.log10(np.abs(arr0[flag])) - np.log10(np.abs(arr1[flag])))/np.log10(np.abs(arr1[flag]))
else:
zPlot = (np.abs(arr0[flag]) - np.abs(arr1[flag]))/np.abs(arr1[flag])
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('RdYlBu')#seismic)
elif par == 'ares':
arF = 1/(8*np.pi**2*10**(-7))
if useLog:
zPlot = (np.log10((arF/arr0['freq'])*np.abs(arr0[flag])**2) - np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2))/np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2)
else:
zPlot = ((arF/arr0['freq'])*np.abs(arr0[flag])**2 - (arF/arr1['freq'])*np.abs(arr1[flag])**2)/((arF/arr1['freq'])*np.abs(arr1[flag])**2)
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral')#seismic)
elif par == 'aphs':
if useLog:
zPlot = (np.log10(np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi)) - np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi)) )/np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
else:
zPlot = ( np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi) - np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi) )/(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral')#seismic)
elif par == 'real':
if useLog:
zPlot = (np.log10(arr0[flag].real) - np.log10(arr1[flag].real))/np.log10(arr1[flag].real)
else:
zPlot = (arr0[flag].real - arr1[flag].real)/arr1[flag].real
if mask:
maskInd = np.logical_or(arr0[flag].real< 1e-3,arr1[flag].real < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral') #('Spectral')
elif par == 'imag':
if useLog:
zPlot = (np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag))/np.log10(arr1[flag].imag)
else:
zPlot = (arr0[flag].imag - arr1[flag].imag)/arr1[flag].imag
if mask:
maskInd = np.logical_or(arr0[flag].imag< 1e-3,arr1[flag].imag < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
else:
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
if colorNorm=='SymLog':
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif colorNorm=='Lin':
if cLevel:
level = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/50.)
clevel = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/10.)
else:
level = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/50.)
clevel = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/10.)
plotNorm = colors.Normalize()
elif colorNorm=='Log':
level = np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)
clevel = np.logspace(zMin,zMax,(zMax-zMin)*2+1,endpoint=True)
plotNorm = colors.LogNorm()
if contour:
cs = ax.tricontourf(x0,y0,zPlot*100,levels=level*100,cmap=cmap,norm=plotNorm,extend='both')#,extend='both')
else:
uniX,uniY = np.unique(x0),np.unique(y0)
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
if colorbar:
csB = plt.colorbar(cs,cax=ax.cax,ticks=clevel*100,format='%1.2e')
# csB.on_mappable_changed(cs)
ax.set_title(flag+' '+par + ' diff',fontsize=8)
return cs, csB
return cs,None
SimPEG/MT/Utils/srcUtils.py
+46
View File
@@ -0,0 +1,46 @@
import SimPEG as simpeg, numpy as np
def homo1DModelSource(mesh,freq,m_back):
'''
Function that calculates and return background fields for a 3D mesh and model.
The calculuations use 1D field solution for a vertical slice throught model (south-western most column),
which is assigned at the fields everywhere for the respective polarizations.2
:param Simpeg mesh object mesh: Holds information on the discretization
:param float freq: The frequency to solve at
:param np.array m_back: Background model of conductivity to base the calculations on.
:rtype: numpy.ndarray (mesh.nE,2)
:return: eBG_bp, E fields for the background model at both polarizations.
'''
# import
from SimPEG.MT.Utils import get1DEfields
# Get a 1d solution for a halfspace background
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
# Note: Everything is using e^iwt
e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq)
# Setup x (east) polarization (_x)
ex_px = np.zeros(mesh.vnEx,dtype=complex)
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
ez_px = np.zeros((mesh.nEz,1),dtype=complex)
# Assign the source to ex_x
for i in np.arange(mesh.vnEx[0]):
for j in np.arange(mesh.vnEx[1]):
ex_px[i,j,:] = -e0_1d
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
# Setup y (north) polarization (_py)
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
# Assign the source to ey_py
for i in np.arange(mesh.vnEy[0]):
for j in np.arange(mesh.vnEy[1]):
ey_py[i,j,:] = e0_1d
# ey_py[1:-1,1:-1,1:-1] = 0
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
# Return the electric fields
eBG_bp = np.hstack((eBG_px,eBG_py))
return eBG_bp
SimPEG/MT/__init__.py
+6
View File
@@ -0,0 +1,6 @@
import Utils
import Sources
from SurveyMT import Rx, Survey, Data
from FieldsMT import Fields1D_e, Fields3D_e
import Problem1D, Problem2D, Problem3D
import SrcMT
docs/mt/index.rst
+11
View File
@@ -0,0 +1,11 @@
SimPEG for Magnetotellurics
===========================
SimPEG (Simulation and Parameter Estimation in Geophysics) is a python
package for simulation and gradient based parameter estimation in the
context of geoscience applications.
simpegMT uses SimPEG as the framework for the forward and inverse
magnetotellurics geophysical problems.
tests/mt/__init__.py
+12
View File
@@ -0,0 +1,12 @@
import os
import glob
import unittest
if __name__ == '__main__':
test_file_strings = glob.glob('test_*.py')
module_strings = [str[0:len(str)-3] for str in test_file_strings]
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
in module_strings]
testSuite = unittest.TestSuite(suites)
unittest.TextTestRunner(verbosity=2).run(testSuite)
tests/mt/test_ApparentResistivityAnalytic.py
+48
View File
@@ -0,0 +1,48 @@
import unittest
from SimPEG import *
from SimPEG import MT
TOL = 1e-6
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi)
return app_res, app_phs
def appResNorm(sigmaHalf):
nFreq = 26
m1d = Mesh.TensorMesh([[(100,5,1.5),(100.,10),(100,5,1.5)]], x0=['C'])
sigma = np.zeros(m1d.nC) + sigmaHalf
sigma[m1d.gridCC[:]>200] = 1e-8
# Calculate the analytic fields
freqs = np.logspace(4,-4,nFreq)
Z = []
for freq in freqs:
Ed, Eu, Hd, Hu = MT.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
Z.append((Ed + Eu)/(Hd + Hu))
Zarr = np.concatenate(Z)
app_r, app_p = appResPhs(freqs,Zarr)
return np.linalg.norm(np.abs(app_r - np.ones(nFreq)/sigmaHalf)) / np.log10(sigmaHalf)
class TestAnalytics(unittest.TestCase):
def setUp(self):
pass
def test_appRes2en1(self):self.assertLess(appResNorm(2e-1), TOL)
def test_appRes2en2(self):self.assertLess(appResNorm(2e-2), TOL)
def test_appRes2en3(self):self.assertLess(appResNorm(2e-3), TOL)
def test_appRes2en4(self):self.assertLess(appResNorm(2e-4), TOL)
def test_appRes2en5(self):self.assertLess(appResNorm(2e-5), TOL)
def test_appRes2en6(self):self.assertLess(appResNorm(2e-6), TOL)
if __name__ == '__main__':
unittest.main()
tests/mt/test_Problem1D_againstAnalyticHalfspace.py
+162
View File
@@ -0,0 +1,162 @@
import unittest
import SimPEG as simpeg
from SimPEG import MT
from SimPEG.Utils import meshTensor
import numpy as np
# Define the tolerances
TOLr = 5e-2
TOLp = 5e-2
def setupSurvey(sigmaHalf,tD=True):
# Frequency
nFreq = 33
freqs = np.logspace(3,-3,nFreq)
# Make the mesh
ct = 5
air = meshTensor([(ct,25,1.3)])
# coreT0 = meshTensor([(ct,15,1.2)])
# coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
core = np.concatenate( ( np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
bot = meshTensor([(core[0],10,-1.3)])
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
# Make the model
sigma = np.zeros(m1d.nC) + sigmaHalf
sigma[m1d.gridCC > 0 ] = 1e-8
rxList = []
for rxType in ['z1dr','z1di']:
rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
# Source list
srcList =[]
if tD:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq))
else:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
survey = MT.Survey(srcList)
return survey, sigma, m1d
def getAppResPhs(MTdata):
# Make impedance
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
return app_res, app_phs
zList = []
for src in MTdata.survey.srcList:
zc = [src.freq]
for rx in src.rxList:
if 'i' in rx.rxType:
m=1j
else:
m = 1
zc.append(m*MTdata[src,rx])
zList.append(zc)
return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
def appRes_TotalFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf)
problem = MT.Problem1D.eForm_TotalField(mesh)
problem.pair(survey)
# Get the fields
fields = problem.fields(sigma)
# Project the data
data = survey.projectFields(fields)
# Calculate the app res and phs
app_r = np.array(getAppResPhs(data))[:,0]
return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf)
def appPhs_TotalFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf)
problem = MT.Problem1D.eForm_TotalField(mesh)
problem.pair(survey)
# Get the fields
fields = problem.fields(sigma)
# Project the data
data = survey.projectFields(fields)
# Calculate the app phs
app_p = np.array(getAppResPhs(data))[:,1]
return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
def appRes_psFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf,False)
problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
problem.pair(survey)
# Get the fields
fields = problem.fields(sigma)
# Project the data
data = survey.projectFields(fields)
# Calculate the app res and phs
app_r = np.array(getAppResPhs(data))[:,0]
return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf)
def appPhs_psFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf,False)
problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
problem.pair(survey)
# Get the fields
fields = problem.fields(sigma)
# Project the data
data = survey.projectFields(fields)
# Calculate the app phs
app_p = np.array(getAppResPhs(data))[:,1]
return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
class TestAnalytics(unittest.TestCase):
def setUp(self):
pass
# Total Fields
# def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr)
# def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp)
# def test_appRes2en2(self):self.assertLess(appRes_TotalFieldNorm(2e-2), TOLr)
# def test_appPhs2en2(self):self.assertLess(appPhs_TotalFieldNorm(2e-2), TOLp)
# def test_appRes2en3(self):self.assertLess(appRes_TotalFieldNorm(2e-3), TOLr)
# def test_appPhs2en3(self):self.assertLess(appPhs_TotalFieldNorm(2e-3), TOLp)
# def test_appRes2en4(self):self.assertLess(appRes_TotalFieldNorm(2e-4), TOLr)
# def test_appPhs2en4(self):self.assertLess(appPhs_TotalFieldNorm(2e-4), TOLp)
# def test_appRes2en5(self):self.assertLess(appRes_TotalFieldNorm(2e-5), TOLr)
# def test_appPhs2en5(self):self.assertLess(appPhs_TotalFieldNorm(2e-5), TOLp)
# def test_appRes2en6(self):self.assertLess(appRes_TotalFieldNorm(2e-6), TOLr)
# def test_appPhs2en6(self):self.assertLess(appPhs_TotalFieldNorm(2e-6), TOLp)
# Primary/secondary
def test_appRes2en2_ps(self):self.assertLess(appRes_psFieldNorm(2e-2), TOLr)
def test_appPhs2en2_ps(self):self.assertLess(appPhs_psFieldNorm(2e-2), TOLp)
if __name__ == '__main__':
unittest.main()
tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py
+135
View File
@@ -0,0 +1,135 @@
import unittest
import SimPEG as simpeg
from SimPEG import MT
from SimPEG.Utils import meshTensor
import numpy as np
# Define the tolerances
TOLr = 5e-2
TOLp = 5e-2
def setupSurvey(sigmaHalf,tD=True):
# Frequency
nFreq = 33
freqs = np.logspace(3,-3,nFreq)
# Make the mesh
ct = 5
air = meshTensor([(ct,25,1.3)])
# coreT0 = meshTensor([(ct,15,1.2)])
# coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
core = np.concatenate( ( np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
bot = meshTensor([(core[0],15,-1.3)])
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
# Make the model
sigma = np.zeros(m1d.nC) + sigmaHalf
sigma[m1d.gridCC > 0 ] = 1e-8
sigmaBack = sigma.copy()
# Add structure
shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600)
deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000)
sigma[shallow] = 1
sigma[deep] = 0.1
rxList = []
for rxType in ['z1dr','z1di']:
rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
# Source list
srcList =[]
if tD:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq))
else:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
survey = MT.Survey(srcList)
return survey, sigma, m1d
def getAppResPhs(MTdata):
# Make impedance
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
return app_res, app_phs
zList = []
for src in MTdata.survey.srcList:
zc = [src.freq]
for rx in src.rxList:
if 'i' in rx.rxType:
m=1j
else:
m = 1
zc.append(m*MTdata[src,rx])
zList.append(zc)
return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
def calculateAnalyticSolution(srcList,mesh,model):
surveyAna = MT.Survey(srcList)
data1D = MT.Data(surveyAna)
for src in surveyAna.srcList:
elev = src.rxList[0].locs[0]
anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev)
anaE = anaEd+anaEu
anaH = anaHd+anaHu
# Scale the solution
# anaE = (anaEtemp/anaEtemp[-1])#.conj()
# anaH = (anaHtemp/anaEtemp[-1])#.conj()
anaZ = anaE/anaH
for rx in src.rxList:
data1D[src,rx] = getattr(anaZ, rx.projComp)
return data1D
def dataMis_AnalyticTotalDomain(sigmaHalf):
# Make the survey
# Total domain solution
surveyTD, sigma, mesh = setupSurvey(sigmaHalf)
problemTD = MT.Problem1D.eForm_TotalField(mesh)
problemTD.pair(surveyTD)
# Analytic data
dataAnaObj = calculateAnalyticSolution(surveyTD.srcList,mesh,sigma)
# dataTDObj = MT.DataMT.DataMT(surveyTD, surveyTD.dpred(sigma))
dataTD = surveyTD.dpred(sigma)
dataAna = simpeg.mkvc(dataAnaObj)
return np.all((dataTD - dataAna)/dataAna < 2.)
# surveyTD.dtrue = -simpeg.mkvc(dataAna,2)
# surveyTD.dobs = -simpeg.mkvc(dataAna,2)
# surveyTD.Wd = np.ones(surveyTD.dtrue.shape) #/(np.abs(surveyTD.dtrue)*0.01)
# # Setup the data misfit
# dmis = simpeg.DataMisfit.l2_DataMisfit(surveyTD)
# dmis.Wd = surveyTD.Wd
# return dmis.eval(sigma)
def dataMis_AnalyticPrimarySecondary(sigmaHalf):
# Make the survey
# Primary secondary
surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,tD=False)
problemPS = MT.Problem1D.eForm_psField(mesh)
problemPS.sigmaPrimary = sigmaPS
problemPS.pair(surveyPS)
# Analytic data
dataAnaObj = calculateAnalyticSolution(surveyPS.srcList,mesh,sigmaPS)
dataPS = surveyPS.dpred(sigmaPS)
dataAna = simpeg.mkvc(dataAnaObj)
return np.all((dataPS - dataAna)/dataAna < 2.)
class TestNumericVsAnalytics(unittest.TestCase):
def setUp(self):
pass
# Total Fields
# def test_appRes2en2(self):self.assertTrue(dataMis_AnalyticTotalDomain(2e-2))
# Primary/secondary
def test_appRes2en2_ps(self):self.assertTrue(dataMis_AnalyticPrimarySecondary(2e-2))
if __name__ == '__main__':
unittest.main()
tests/mt/test_Problem3D_againstAnalytic.py
+290
View File
@@ -0,0 +1,290 @@
# Test functions
from glob import glob
import numpy as np, sys, os, time, scipy, subprocess
import SimPEG as simpeg
import unittest
from SimPEG import MT
from SimPEG.Utils import meshTensor
from scipy.constants import mu_0
TOLr = 5e-2
TOL = 1e-4
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
CONDUCTIVITY = 1e1
MU = mu_0
freq = [1e-1, 2e-1]
addrandoms = True
def getInputs():
"""
Function that returns Mesh, freqs, rx_loc, elev.
"""
# 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,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
## 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]]) #,[-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
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()
# Model
ccM = M.gridCC
# conds = [1e-2]
groundInd = ccM[:,2] < elev
sig = np.zeros(M.nC) + 1e-8
sig[groundInd] = conds
# Set the background, not the same as the model
sigBG = np.zeros(M.nC) + 1e-8
sigBG[groundInd] = conds
return (M, freqs, sig, sigBG, rx_loc)
def twoLayer(conds):
''' Returns a 2 layer model based on the conductivity values given'''
M, freqs, rx_loc, elev = getInputs()
# Model
ccM = M.gridCC
groundInd = ccM[:,2] < elev
botInd = ccM[:,2] < -3000
sig = np.zeros(M.nC) + 1e-8
sig[groundInd] = conds[1]
sig[botInd] = conds[0]
# Set the background, not the same as the model
sigBG = np.zeros(M.nC) + 1e-8
sigBG[groundInd] = conds[1]
return (M, freqs, sig, sigBG, rx_loc)
def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False,expMap=True):
M,freqs,sig,sigBG,rx_loc = inputSetup
# Make a receiver list
rxList = []
if comp == 'All':
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi',]:
rxList.append(MT.Rx(rx_loc,rxType))
else:
rxList.append(MT.Rx(rx_loc,comp))
# Source list
srcList =[]
if singleFreq:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,singleFreq))
else:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
# Survey MT
survey = MT.Survey(srcList)
## Setup the problem object
sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
if expMap:
problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) )
problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
problem.curModel = np.log(sig)
else:
problem = MT.Problem3D.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 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(MT.Rx(rx_loc,rxType))
else:
rxList.append(MT.Rx(rx_loc,comp))
# Source list
srcList =[]
if singleFreq:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,singleFreq))
else:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
# Survey MT
survey = MT.SurveyMT.SurveyMT(srcList)
## Setup the problem object
sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
if expMap:
problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) )
problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
problem.curModel = np.log(sig)
else:
problem = MT.Problem3D.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):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
return app_res, app_phs
recData = MTdata.toRecArray('Complex')
return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx'])
def JvecAdjointTest(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,)
# print problem.PropMap.PropModel.nP
w = np.random.rand(problem.mesh.nC,)
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 DerivJvecTest(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 DerivProjfieldsTest(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)
# Initate things for the derivs Test
src = survey.srcList[0]
rx = src.rxList[0]
u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))
f0 = problem.fieldsPair(survey.mesh,survey)
# u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))
f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x
f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y
def fun(u):
f = problem.fieldsPair(survey.mesh,survey)
f[src,'e_pxSolution'] = u[:len(u)/2]
f[src,'e_pySolution'] = u[len(u)/2::]
return rx.projectFields(src,survey.mesh,f), lambda t: rx.projectFieldsDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2))
return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)
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),expMap=expMap)
print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf)
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 < .4)
else:
return np.all(np.abs(app_rpxy[1,:] + np.ones(survey.nFreq)*135) / 135 < .4)
class TestAnalytics(unittest.TestCase):
def setUp(self):
pass
# # 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))
# Do a derivative test
def test_derivProj1(self):self.assertTrue(DerivProjfieldsTest(halfSpace(1e-2)))
# Do a derivative test of Jvec
# def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
# def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))
# def test_derivJvec_zxyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyr',.1))
# def test_derivJvec_zxyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyi',.1))
# def test_derivJvec_zyxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxr',.1))
# def test_derivJvec_zyxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxi',.1))
# def test_derivJvec_zyyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyr',.1))
# def test_derivJvec_zyyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyi',.1))
def test_derivJvec_All(self):self.assertTrue(DerivJvecTest(random(1e-2),'All',.1))
# Test the adjoint of Jvec and Jtvec
# def test_JvecAdjoint_zxxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxr',.1))
# def test_JvecAdjoint_zxxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxi',.1))
# def test_JvecAdjoint_zxyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyr',.1))
# def test_JvecAdjoint_zxyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyi',.1))
# def test_JvecAdjoint_zyxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxr',.1))
# def test_JvecAdjoint_zyxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxi',.1))
# def test_JvecAdjoint_zyyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyr',.1))
# def test_JvecAdjoint_zyyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyi',.1))
def test_JvecAdjoint_All(self):self.assertTrue(JvecAdjointTest(random(1e-2),'All',.1))
if __name__ == '__main__':
unittest.main()