mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 09:58:23 +08:00
Changed the name space to correspond with FDEM.
This commit is contained in:
GudniRos
+5
-34
@@ -1,8 +1,7 @@
|
||||
from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
|
||||
from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem
|
||||
from SurveyMT import SurveyMT
|
||||
from DataMT import DataMT
|
||||
from FieldsMT import FieldsMT
|
||||
from SurveyMT import Survey, Data
|
||||
from FieldsMT import BaseMTFields
|
||||
|
||||
|
||||
class BaseMTProblem(BaseFDEMProblem):
|
||||
@@ -11,37 +10,9 @@ class BaseMTProblem(BaseFDEMProblem):
|
||||
BaseFDEMProblem.__init__(self, mesh, **kwargs)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
# Set the default pairs of the problem
|
||||
surveyPair = SurveyMT
|
||||
dataPair = DataMT
|
||||
fieldsPair = FieldsMT
|
||||
|
||||
# # Pickleing support methods
|
||||
# def __getstate__(self):
|
||||
# '''
|
||||
# Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
|
||||
|
||||
# Used when doing:
|
||||
# pickle.dump(pickleFile,object)
|
||||
# '''
|
||||
# odict = self.__dict__.copy()
|
||||
# # Remove fields that are not needed
|
||||
# del odict['hook']
|
||||
# del odict['setKwargs']
|
||||
# # Return the dict
|
||||
# return odict
|
||||
|
||||
# def __setstate__(self,odict):
|
||||
# '''
|
||||
# Function that sets a pickle dictionary in to an object.
|
||||
|
||||
# Used when doing:
|
||||
# object = pickle.load(pickleFile)
|
||||
# '''
|
||||
# # Update the dict
|
||||
# self.__dict__.update(odict)
|
||||
# # Re-hook the methods to the object
|
||||
# Utils.codeutils.hook(self,Utils.codeutils.hook)
|
||||
# Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
|
||||
surveyPair = Survey
|
||||
dataPair = Data
|
||||
fieldsPair = BaseMTFields
|
||||
|
||||
# Set the solver
|
||||
Solver = SimpegSolver
|
||||
|
||||
@@ -1,128 +0,0 @@
|
||||
from SimPEG import Survey, Utils, Problem, np, sp, mkvc
|
||||
from simpegMT.Utils import rec2ndarr
|
||||
import simpegMT
|
||||
from scipy.constants import mu_0
|
||||
import sys
|
||||
from numpy.lib import recfunctions as recFunc
|
||||
|
||||
|
||||
############
|
||||
### Data ###
|
||||
############
|
||||
class DataMT(Survey.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
|
||||
Survey.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 = simpegMT.SurveyMT.srcMT_polxy_1Dprimary
|
||||
elif srcType=='total':
|
||||
src = sdsimpegMT.SurveyMT.srcMT_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(simpegMT.SurveyMT.RxMT(locs,rxType+'r'))
|
||||
dataList.append(dFreq[rxType][notNaNind].real.copy())
|
||||
rxList.append(simpegMT.SurveyMT.RxMT(locs,rxType+'i'))
|
||||
dataList.append(dFreq[rxType][notNaNind].imag.copy())
|
||||
else:
|
||||
rxList.append(simpegMT.SurveyMT.RxMT(locs,rxType))
|
||||
dataList.append(dFreq[rxType][notNaNind].copy())
|
||||
srcList.append(src(rxList,freq))
|
||||
|
||||
# Make a survey
|
||||
survey = simpegMT.SurveyMT.SurveyMT(srcList)
|
||||
dataVec = np.hstack(dataList)
|
||||
return cls(survey,dataVec)
|
||||
@@ -1,7 +1,7 @@
|
||||
# Test script to use simpegMT platform to forward model synthetic data.
|
||||
|
||||
# Import
|
||||
import simpegMT as simpegmt, SimPEG as simpeg
|
||||
# Import
|
||||
import simpegMT as simpegmt, SimPEG as simpeg
|
||||
import numpy as np
|
||||
|
||||
# Make a mesh
|
||||
@@ -23,13 +23,13 @@ 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(simpegmt.SurveyMT.RxMT(simpeg.mkvc(loc,2).T,rxType))
|
||||
rxList.append(simpegmt.RxMT(simpeg.mkvc(loc,2).T,rxType))
|
||||
# Source list
|
||||
srcList =[]
|
||||
for freq in np.logspace(3,-3,7):
|
||||
srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList))
|
||||
# Survey MT
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
# Survey MT
|
||||
survey = simpegmt.Survey(srcList)
|
||||
|
||||
## Setup the problem object
|
||||
problem = simpegmt.ProblemMT.MTProblem(M)
|
||||
|
||||
@@ -7,13 +7,13 @@ from SimPEG.EM.Utils import omega
|
||||
##############
|
||||
### Fields ###
|
||||
##############
|
||||
class FieldsMT(Problem.Fields):
|
||||
class BaseMTFields(Problem.Fields):
|
||||
"""Field Storage for a MT survey."""
|
||||
knownFields = {}
|
||||
dtype = complex
|
||||
|
||||
|
||||
class FieldsMT_1D(FieldsMT):
|
||||
class Fields1D_e(BaseMTFields):
|
||||
"""
|
||||
Fields storage for the 1D MT solution.
|
||||
"""
|
||||
@@ -28,7 +28,7 @@ class FieldsMT_1D(FieldsMT):
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsMT.__init__(self,mesh,survey,**kwargs)
|
||||
BaseMTFields.__init__(self,mesh,survey,**kwargs)
|
||||
|
||||
def _ePrimary(self, eSolution, srcList):
|
||||
ePrimary = np.zeros_like(eSolution)
|
||||
@@ -120,7 +120,7 @@ class FieldsMT_1D(FieldsMT):
|
||||
"""
|
||||
return None
|
||||
|
||||
class FieldsMT_3D(FieldsMT):
|
||||
class Fields3D_e(BaseMTFields):
|
||||
"""
|
||||
Fields storage for the 3D MT solution.
|
||||
"""
|
||||
@@ -144,7 +144,7 @@ class FieldsMT_3D(FieldsMT):
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsMT.__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)
|
||||
|
||||
@@ -2,9 +2,8 @@ from SimPEG.EM.Utils import omega
|
||||
from SimPEG import mkvc
|
||||
from scipy.constants import mu_0
|
||||
from simpegMT.BaseMT import BaseMTProblem
|
||||
from simpegMT.SurveyMT import SurveyMT
|
||||
from simpegMT.FieldsMT import FieldsMT_1D
|
||||
from simpegMT.DataMT import DataMT
|
||||
from simpegMT.SurveyMT import Survey, Data
|
||||
from simpegMT.FieldsMT import Fields1D_e
|
||||
from simpegMT.Utils.MT1Danalytic import getEHfields
|
||||
import numpy as np
|
||||
import multiprocessing, sys, time
|
||||
@@ -25,7 +24,7 @@ class eForm_psField(BaseMTProblem):
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseMTProblem.__init__(self, mesh, **kwargs)
|
||||
self.fieldsPair = FieldsMT_1D
|
||||
self.fieldsPair = Fields1D_e
|
||||
# self._sigmaPrimary = sigmaPrimary
|
||||
|
||||
@property
|
||||
@@ -105,7 +104,7 @@ class eForm_psField(BaseMTProblem):
|
||||
# Set the current model
|
||||
self.curModel = m
|
||||
|
||||
F = FieldsMT_1D(self.mesh, self.survey)
|
||||
F = Fields1D_e(self.mesh, self.survey)
|
||||
for freq in self.survey.freqs:
|
||||
if self.verbose:
|
||||
startTime = time.time()
|
||||
@@ -222,7 +221,7 @@ class eForm_TotalField(BaseMTProblem):
|
||||
self.curModel = m
|
||||
# RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
|
||||
|
||||
F = FieldsMT_1D(self.mesh, self.survey)
|
||||
F = Fields1D_e(self.mesh, self.survey)
|
||||
for freq in self.survey.freqs:
|
||||
if self.verbose:
|
||||
startTime = time.time()
|
||||
|
||||
@@ -2,9 +2,8 @@ from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as Sim
|
||||
from SimPEG.EM.Utils import omega
|
||||
from scipy.constants import mu_0
|
||||
from simpegMT.BaseMT import BaseMTProblem
|
||||
from simpegMT.SurveyMT import SurveyMT
|
||||
from simpegMT.FieldsMT import FieldsMT_3D
|
||||
from simpegMT.DataMT import DataMT
|
||||
from simpegMT.SurveyMT import Survey, Data
|
||||
from simpegMT.FieldsMT import Fields3D_e
|
||||
import multiprocessing, sys, time
|
||||
|
||||
|
||||
@@ -22,7 +21,7 @@ class eForm_ps(BaseMTProblem):
|
||||
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
|
||||
_fieldType = 'e'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = FieldsMT_3D
|
||||
fieldsPair = Fields3D_e
|
||||
_sigmaPrimary = None
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
@@ -110,7 +109,7 @@ class eForm_ps(BaseMTProblem):
|
||||
# Set the current model
|
||||
self.curModel = m
|
||||
|
||||
F = FieldsMT_3D(self.mesh, self.survey)
|
||||
F = Fields3D_e(self.mesh, self.survey)
|
||||
for freq in self.survey.freqs:
|
||||
if self.verbose:
|
||||
startTime = time.time()
|
||||
@@ -144,7 +143,7 @@ class eForm_Tp(BaseMTProblem):
|
||||
|
||||
_fieldType = 'e'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = FieldsMT_3D
|
||||
fieldsPair = Fields3D_e
|
||||
|
||||
# Set new properties
|
||||
# Background model
|
||||
@@ -238,7 +237,7 @@ class eForm_Tp(BaseMTProblem):
|
||||
self.backModel = m_back
|
||||
# RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
|
||||
|
||||
F = FieldsMT_3D(self.mesh, self.survey)
|
||||
F = Fields3D_e(self.mesh, self.survey)
|
||||
for freq in self.survey.freqs:
|
||||
if self.verbose:
|
||||
startTime = time.time()
|
||||
|
||||
@@ -1 +0,0 @@
|
||||
import 1D, 2D, 3D
|
||||
@@ -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
|
||||
+138
-214
@@ -1,16 +1,17 @@
|
||||
from SimPEG import Survey, Utils, Problem, Maps, np, sp, mkvc
|
||||
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
|
||||
import sys
|
||||
from numpy.lib import recfunctions as recFunc
|
||||
from DataMT import DataMT
|
||||
from simpegMT.Sources import homo1DModelSource
|
||||
from Sources import homo1DModelSource
|
||||
from Utils import rec2ndarr
|
||||
|
||||
import sys
|
||||
|
||||
#################
|
||||
### Receivers ###
|
||||
#################
|
||||
|
||||
class RxMT(Survey.BaseRx):
|
||||
class Rx(SimPEGsurvey.BaseRx):
|
||||
|
||||
knownRxTypes = {
|
||||
# 3D impedance
|
||||
@@ -35,7 +36,7 @@ class RxMT(Survey.BaseRx):
|
||||
}
|
||||
# TODO: Have locs as single or double coordinates for both or numerator and denominator separately, respectively.
|
||||
def __init__(self, locs, rxType):
|
||||
Survey.BaseRx.__init__(self, locs, rxType)
|
||||
SimPEGsurvey.BaseRx.__init__(self, locs, rxType)
|
||||
|
||||
@property
|
||||
def projField(self):
|
||||
@@ -64,6 +65,7 @@ class RxMT(Survey.BaseRx):
|
||||
return self.knownRxTypes[self.rxType][0][1]
|
||||
else:
|
||||
raise Exception('{s} is an unknown option. Use numerator or denominator.')
|
||||
|
||||
@property
|
||||
def projType(self):
|
||||
"""
|
||||
@@ -310,222 +312,23 @@ class RxMT(Survey.BaseRx):
|
||||
|
||||
return Pv
|
||||
|
||||
|
||||
###############
|
||||
### Sources ###
|
||||
###############
|
||||
|
||||
class srcMT(FDEMBaseSrc): # Survey.BaseSrc):
|
||||
'''
|
||||
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 = RxMT
|
||||
|
||||
|
||||
def __init__(self, rxList, freq):
|
||||
|
||||
self.freq = float(freq)
|
||||
Survey.BaseSrc.__init__(self, rxList)
|
||||
|
||||
# 1D sources
|
||||
class srcMT_polxy_1DhomotD(srcMT):
|
||||
"""
|
||||
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):
|
||||
srcMT.__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 srcMT_polxy_1Dprimary(srcMT):
|
||||
"""
|
||||
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
|
||||
srcMT.__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 srcMT_polxy_3Dprimary(srcMT):
|
||||
"""
|
||||
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
|
||||
srcMT.__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
|
||||
|
||||
##############
|
||||
### Survey ###
|
||||
##############
|
||||
class SurveyMT(Survey.BaseSurvey):
|
||||
#################
|
||||
### 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
|
||||
|
||||
"""
|
||||
|
||||
srcPair = srcMT
|
||||
import SrcMT
|
||||
srcPair = SrcMT.BaseMTSrc
|
||||
|
||||
def __init__(self, srcList, **kwargs):
|
||||
# Sort these by frequency
|
||||
self.srcList = srcList
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
SimPEGsurvey.BaseSurvey.__init__(self, **kwargs)
|
||||
|
||||
_freqDict = {}
|
||||
for src in srcList:
|
||||
@@ -553,7 +356,7 @@ class SurveyMT(Survey.BaseSurvey):
|
||||
return self._freqDict[freq]
|
||||
|
||||
def projectFields(self, u):
|
||||
data = DataMT(self)
|
||||
data = Data(self)
|
||||
for src in self.srcList:
|
||||
sys.stdout.flush()
|
||||
for rx in src.rxList:
|
||||
@@ -563,3 +366,124 @@ class SurveyMT(Survey.BaseSurvey):
|
||||
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)
|
||||
|
||||
|
||||
@@ -28,17 +28,17 @@ def setupSurvey(sigmaHalf,tD=True):
|
||||
|
||||
rxList = []
|
||||
for rxType in ['z1dr','z1di']:
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))
|
||||
rxList.append(simpegmt.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
|
||||
# Source list
|
||||
srcList =[]
|
||||
if tD:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))
|
||||
srcList.append(simpegmt.SrcMT.src_polxy_1DhomotD(rxList,freq))
|
||||
else:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
|
||||
srcList.append(simpegmt.SrcMT.src_polxy_1Dprimary(rxList,freq))
|
||||
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
survey = simpegmt.Survey(srcList)
|
||||
return survey, sigma, m1d
|
||||
|
||||
def getAppResPhs(MTdata):
|
||||
|
||||
@@ -34,17 +34,17 @@ def setupSurvey(sigmaHalf,tD=True):
|
||||
|
||||
rxList = []
|
||||
for rxType in ['z1dr','z1di']:
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))
|
||||
rxList.append(simpegmt.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
|
||||
# Source list
|
||||
srcList =[]
|
||||
if tD:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))
|
||||
srcList.append(simpegmt.SrcMT.src_polxy_1DhomotD(rxList,freq))
|
||||
else:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
|
||||
srcList.append(simpegmt.SrcMT.src_polxy_1Dprimary(rxList,freq))
|
||||
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
survey = simpegmt.Survey(srcList)
|
||||
return survey, sigma, m1d
|
||||
|
||||
def getAppResPhs(MTdata):
|
||||
@@ -66,8 +66,8 @@ def getAppResPhs(MTdata):
|
||||
return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
|
||||
|
||||
def calculateAnalyticSolution(srcList,mesh,model):
|
||||
surveyAna = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
data1D = simpegmt.DataMT.DataMT(surveyAna)
|
||||
surveyAna = simpegmt.Survey(srcList)
|
||||
data1D = simpegmt.Data(surveyAna)
|
||||
for src in surveyAna.srcList:
|
||||
elev = src.rxList[0].locs[0]
|
||||
anaEd, anaEu, anaHd, anaHu = simpegmt.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev)
|
||||
|
||||
@@ -88,19 +88,19 @@ def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False,expMap=True
|
||||
rxList = []
|
||||
if comp == 'All':
|
||||
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi',]:
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,rxType))
|
||||
rxList.append(simpegmt.Rx(rx_loc,rxType))
|
||||
else:
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,comp))
|
||||
rxList.append(simpegmt.Rx(rx_loc,comp))
|
||||
# Source list
|
||||
srcList =[]
|
||||
|
||||
if singleFreq:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,singleFreq))
|
||||
srcList.append(simpegmt.SrcMT.polxy_1Dprimary(rxList,singleFreq))
|
||||
else:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
|
||||
srcList.append(simpegmt.SrcMT.polxy_1Dprimary(rxList,freq))
|
||||
# Survey MT
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
survey = simpegmt.Survey(srcList)
|
||||
|
||||
## Setup the problem object
|
||||
sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
|
||||
@@ -129,17 +129,17 @@ def setupSimpegMTfwd_eForm_ps_multiRx(inputSetup,comp='All',singleFreq=False,exp
|
||||
rxList = []
|
||||
if comp == 'All':
|
||||
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi',]:
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,rxType))
|
||||
rxList.append(simpegmt.Rx(rx_loc,rxType))
|
||||
else:
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,comp))
|
||||
rxList.append(simpegmt.Rx(rx_loc,comp))
|
||||
# Source list
|
||||
srcList =[]
|
||||
|
||||
if singleFreq:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,singleFreq))
|
||||
srcList.append(simpegmt.SrcMT.polxy_1Dprimary(rxList,singleFreq))
|
||||
else:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
|
||||
srcList.append(simpegmt.SrcMT.polxy_1Dprimary(rxList,freq))
|
||||
# Survey MT
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
|
||||
|
||||
@@ -40,7 +40,7 @@ def rotateData(MTdata,rotAngle):
|
||||
outRec = recData.copy()
|
||||
for nr,comp in enumerate(['zxx','zxy','zyx','zyy']):
|
||||
outRec[comp] = rotData[:,nr]
|
||||
return simpegmt.DataMT.DataMT.fromRecArray(outRec)
|
||||
return simpegmt.Data.fromRecArray(outRec)
|
||||
|
||||
|
||||
def appResPhs(freq,z):
|
||||
@@ -182,22 +182,22 @@ def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
|
||||
# Make the receiver list
|
||||
rx1DList = []
|
||||
for rxType in ['z1dr','z1di']:
|
||||
rx1DList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(loc,2).T,rxType))
|
||||
rx1DList.append(simpegmt.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(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rx1DList,freq))
|
||||
src1DList.append(simpegmt.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 = simpegmt.SurveyMT.SurveyMT(src1DList)
|
||||
sur1D = simpegmt.Survey(src1DList)
|
||||
|
||||
# Make the data
|
||||
dataVec = np.hstack(dat1DList)
|
||||
dat1D = simpegmt.DataMT.DataMT(sur1D,dataVec)
|
||||
dat1D = simpegmt.Data(sur1D,dataVec)
|
||||
sur1D.dobs = dataVec
|
||||
# Need to take MTdata.survey.std and split it as well.
|
||||
std=0.05
|
||||
@@ -238,4 +238,4 @@ def resampleMTdataAtFreq(MTdata,freqs):
|
||||
outRecArr = tArrRec
|
||||
|
||||
# Make the MTdata and return
|
||||
return simpegmt.DataMT.DataMT.fromRecArray(outRecArr)
|
||||
return simpegmt.Data.fromRecArray(outRecArr)
|
||||
@@ -1,8 +1,6 @@
|
||||
# from EM import *
|
||||
import Utils
|
||||
# import Tests
|
||||
import Sources
|
||||
# from BaseMT import SurveyMT, RxMT, srcMT, DataMT, FieldsMT
|
||||
# import BaseMT
|
||||
import SurveyMT, DataMT, FieldsMT
|
||||
from SurveyMT import Rx, Survey, Data
|
||||
from FieldsMT import Fields1D_e, Fields3D_e
|
||||
import Problem1D, Problem2D, Problem3D
|
||||
import SrcMT
|
||||
Reference in New Issue
Block a user