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

Updated MT-->NSEM for all the classes.

Browse Source 
Changed Vertical1DMap --> SurjectVertical1D
This commit is contained in:
GudniRos
2016-05-10 14:41:20 -07:00
parent a96e9e08d7
commit 1ab67b5790
24 changed files with 151 additions and 151 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/MT/__init__.py
-5
View File
@@ -1,5 +0,0 @@
import Utils
from SurveyMT import Rx, Survey, Data
from FieldsMT import Fields1D_e, Fields3D_e
import Problem1D, Problem2D, Problem3D
import SrcMT
SimPEG/MT/FieldsMT.py → SimPEG/NSEM/FieldsNSEM.py
+13 -13
View File
@@ -7,15 +7,15 @@ from SimPEG.EM.Utils import omega
##############
### Fields ###
##############
class BaseMTFields(Problem.Fields):
"""Field Storage for a MT survey."""
class BaseNSEMFields(Problem.Fields):
"""Field Storage for a NSEM survey."""
knownFields = {}
dtype = complex
class Fields1D_e(BaseMTFields):
class Fields1D_e(BaseNSEMFields):
"""
Fields storage for the 1D MT solution.
Fields storage for the 1D NSEM solution.
"""
knownFields = {'e_1dSolution':'F'}
aliasFields = {
@@ -28,7 +28,7 @@ class Fields1D_e(BaseMTFields):
}
def __init__(self,mesh,survey,**kwargs):
BaseMTFields.__init__(self,mesh,survey,**kwargs)
BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
def _ePrimary(self, eSolution, srcList):
ePrimary = np.zeros_like(eSolution)
@@ -99,7 +99,7 @@ class Fields1D_e(BaseMTFields):
"""
Derivative of the fields object wrt u.
:param MTsrc src: MT source
:param NSEMsrc src: NSEM source
:param numpy.ndarray v: random vector of f_sol.size
This function stacks the fields derivatives appropriately
@@ -120,9 +120,9 @@ class Fields1D_e(BaseMTFields):
"""
return None
class Fields3D_e(BaseMTFields):
class Fields3D_e(BaseNSEMFields):
"""
Fields storage for the 3D MT solution. Labels polarizations by px and py.
Fields storage for the 3D NSEM solution. Labels polarizations by px and py.
:param SimPEG object mesh: The solution mesh
:param SimPEG object survey: A survey object
@@ -147,7 +147,7 @@ class Fields3D_e(BaseMTFields):
}
def __init__(self,mesh,survey,**kwargs):
BaseMTFields.__init__(self,mesh,survey,**kwargs)
BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
def _e_pxPrimary(self, e_pxSolution, srcList):
e_pxPrimary = np.zeros_like(e_pxSolution)
@@ -228,7 +228,7 @@ class Fields3D_e(BaseMTFields):
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
# There is no magnetic source in the NSEM problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
# b[:,i] += 1./(1j*omega(src.freq)) * S_m
@@ -239,7 +239,7 @@ class Fields3D_e(BaseMTFields):
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
# There is no magnetic source in the NSEM problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
# b[:,i] += 1./(1j*omega(src.freq)) * S_m
@@ -302,7 +302,7 @@ class Fields3D_e(BaseMTFields):
"""
Derivative of the fields object wrt u.
:param MTsrc src: MT source
:param NSEMsrc src: NSEM source
:param numpy.ndarray v: random vector of f_sol.size
This function stacks the fields derivatives appropriately
@@ -319,7 +319,7 @@ class Fields3D_e(BaseMTFields):
"""
Derivative of the fields object wrt u.
:param MTsrc src: MT source
:param NSEMsrc src: NSEM source
:param numpy.ndarray v: random vector of f_sol.size
This function stacks the fields derivatives appropriately
SimPEG/MT/BaseMT.py → SimPEG/NSEM/NSEM.py
+8 -8
View File
@@ -1,10 +1,10 @@
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
from SurveyNSEM import Survey, Data
from FieldsNSEM import BaseNSEMFields
class BaseMTProblem(BaseFDEMProblem):
class BaseNSEMProblem(BaseFDEMProblem):
"""
Base class for all Natural source problems.
"""
@@ -15,7 +15,7 @@ class BaseMTProblem(BaseFDEMProblem):
# Set the default pairs of the problem
surveyPair = Survey
dataPair = Data
fieldsPair = BaseMTFields
fieldsPair = BaseNSEMFields
# Set the solver
Solver = SimpegSolver
@@ -33,8 +33,8 @@ class BaseMTProblem(BaseFDEMProblem):
:param numpy.ndarray m (nC, 1) - conductive model
:param numpy.ndarray v (nC, 1) - random vector
:param MTfields object (optional) - MT fields object, if not given it is calculated
:rtype: MTdata object
:param NSEMfields object (optional) - NSEM fields object, if not given it is calculated
:rtype: NSEMdata object
:return: Data sensitivities wrt m
"""
@@ -80,8 +80,8 @@ class BaseMTProblem(BaseFDEMProblem):
:param numpy.ndarray m (nC, 1) - conductive model
:param numpy.ndarray v (nD, 1) - vector
:param MTfields object u (optional) - MT fields object, if not given it is calculated
:rtype: MTdata object
:param NSEMfields object u (optional) - NSEM fields object, if not given it is calculated
:rtype: NSEMdata object
:return: Data sensitivities wrt m
"""
SimPEG/MT/Problem1D/Probs.py → SimPEG/NSEM/Problem1D/Probs.py
+12 -12
View File
@@ -1,17 +1,17 @@
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
from SimPEG.NSEM.NSEM import BaseNSEMProblem
from SimPEG.NSEM.SurveyNSEM import Survey, Data
from SimPEG.NSEM.FieldsNSEM import Fields1D_e
from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields
import numpy as np
import multiprocessing, sys, time
class eForm_psField(BaseMTProblem):
class eForm_psField(BaseNSEMProblem):
"""
A MT problem soving a e formulation and primary/secondary fields decomposion.
A NSEM problem soving a e formulation and primary/secondary fields decomposion.
By eliminating the magnetic flux density using
@@ -30,14 +30,14 @@ class eForm_psField(BaseMTProblem):
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
# From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
_fieldType = 'e_1d'
_eqLocs = 'EF'
_sigmaPrimary = None
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
BaseNSEMProblem.__init__(self, mesh, **kwargs)
self.fieldsPair = Fields1D_e
# self._sigmaPrimary = sigmaPrimary
@property
@@ -167,9 +167,9 @@ class eForm_psField(BaseMTProblem):
# Fields class corresponding to the fields
# Update Jvec and Jtvec to include all the derivatives components
# Other things ...
class eForm_TotalField(BaseMTProblem):
class eForm_TotalField(BaseNSEMProblem):
"""
A MT problem solving a e formulation and a Total bondary domain decompostion.
A NSEM problem solving a e formulation and a Total bondary domain decompostion.
Solves the equation:
@@ -178,13 +178,13 @@ class eForm_TotalField(BaseMTProblem):
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
# From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
_fieldType = 'e'
_eqLocs = 'EF'
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
BaseNSEMProblem.__init__(self, mesh, **kwargs)
@property
def MeMui(self):
"""
SimPEG/MT/Problem1D/__init__.py → SimPEG/NSEM/Problem1D/__init__.py
View File
SimPEG/MT/Problem2D/Probs.py → SimPEG/NSEM/Problem2D/Probs.py
View File
SimPEG/MT/Problem2D/__init__.py → SimPEG/NSEM/Problem2D/__init__.py
View File
SimPEG/MT/Problem3D/Probs.py → SimPEG/NSEM/Problem3D/Probs.py
+7 -7
View File
@@ -1,16 +1,16 @@
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
from SimPEG.NSEM.NSEM import BaseNSEMProblem
from SimPEG.NSEM.SurveyNSEM import Survey, Data
from SimPEG.NSEM.FieldsNSEM import Fields3D_e
import multiprocessing, sys, time
class eForm_ps(BaseMTProblem):
class eForm_ps(BaseNSEMProblem):
"""
A MT problem solving a e formulation and a primary/secondary fields decompostion.
A NSEM problem solving a e formulation and a primary/secondary fields decompostion.
By eliminating the magnetic flux density using
@@ -29,14 +29,14 @@ class eForm_ps(BaseMTProblem):
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
# From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
_fieldType = 'e'
_eqLocs = 'FE'
fieldsPair = Fields3D_e
_sigmaPrimary = None
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
BaseNSEMProblem.__init__(self, mesh, **kwargs)
@property
def sigmaPrimary(self):
SimPEG/MT/Problem3D/__init__.py → SimPEG/NSEM/Problem3D/__init__.py
View File
SimPEG/MT/SrcMT.py → SimPEG/NSEM/SrcNSEM.py
+13 -13
View File
@@ -11,9 +11,9 @@ import sys
### Sources ###
#################
class BaseMTSrc(FDEMBaseSrc):
class BaseNSEMSrc(FDEMBaseSrc):
'''
Sources for the MT problem.
Sources for the NSEM problem.
Use the SimPEG BaseSrc, since the source fields share properties with the transmitters.
:param float freq: The frequency of the source
@@ -29,28 +29,28 @@ class BaseMTSrc(FDEMBaseSrc):
FDEMBaseSrc.__init__(self, rxList)
# 1D sources
class polxy_1DhomotD(BaseMTSrc):
class polxy_1DhomotD(BaseNSEMSrc):
"""
MT source for both polarizations (x and y) for the total Domain.
NSEM 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)
BaseNSEMSrc.__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):
class polxy_1Dprimary(BaseNSEMSrc):
"""
MT source for both polarizations (x and y) given a 1D primary models.
NSEM 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)
BaseNSEMSrc.__init__(self, rxList, freq)
# Hidden property of the ePrimary
self._ePrimary = None
@@ -86,7 +86,7 @@ class polxy_1Dprimary(BaseMTSrc):
Get the electrical field source
"""
e_p = self.ePrimary(problem)
Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
sigma_p = Map_sigma_p._transform(self.sigma1d)
# Make mass matrix
# Note: M(sig) - M(sig_p) = M(sig - sig_p)
@@ -128,15 +128,15 @@ class polxy_1Dprimary(BaseMTSrc):
# v should be nC size
return MsigmaDeriv * v
class polxy_3Dprimary(BaseMTSrc):
class polxy_3Dprimary(BaseNSEMSrc):
"""
MT source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
NSEM 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)
BaseNSEMSrc.__init__(self, rxList, freq)
# Hidden property of the ePrimary
self._ePrimary = None
@@ -163,7 +163,7 @@ class polxy_3Dprimary(BaseMTSrc):
Get the electrical field source
"""
e_p = self.ePrimary(problem)
Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
sigma_p = Map_sigma_p._transform(self.sigma1d)
# Make mass matrix
# Note: M(sig) - M(sig_p) = M(sig - sig_p)
SimPEG/MT/SurveyMT.py → SimPEG/NSEM/SurveyNSEM.py
+14 -14
View File
@@ -4,7 +4,7 @@ from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from Utils import rec2ndarr
import SrcMT
import SrcNSEM
import sys
#################
@@ -63,9 +63,9 @@ class Rx(SimPEGsurvey.BaseRx):
'''
Project the fields to natural source data.
:param SrcMT src: The source of the fields to project
:param SrcNSEM src: The source of the fields to project
:param SimPEG.Mesh mesh:
:param FieldsMT f: Natural source fields object to project
:param FieldsNSEM f: Natural source fields object to project
'''
## NOTE: Assumes that e is on t
@@ -143,9 +143,9 @@ class Rx(SimPEGsurvey.BaseRx):
"""
The derivative of the projection wrt u
:param MTsrc src: MT source
:param NSEMsrc src: NSEM source
:param TensorMesh mesh: Mesh defining the topology of the problem
:param MTfields f: MT fields object of the source
:param NSEMfields f: NSEM fields object of the source
:param numpy.ndarray v: Random vector of size
"""
@@ -390,12 +390,12 @@ class Rx(SimPEGsurvey.BaseRx):
#################
class Survey(SimPEGsurvey.BaseSurvey):
"""
Survey class for MT. Contains all the sources associated with the survey.
Survey class for NSEM. Contains all the sources associated with the survey.
:param list srcList: List of sources associated with the survey
"""
srcPair = SrcMT.BaseMTSrc
srcPair = SrcNSEM.BaseNSEMSrc
def __init__(self, srcList, **kwargs):
# Sort these by frequency
@@ -443,7 +443,7 @@ class Survey(SimPEGsurvey.BaseSurvey):
#################
class Data(SimPEGsurvey.Data):
'''
Data class for MTdata. Stores the data vector indexed by the survey.
Data class for NSEMdata. Stores the data vector indexed by the survey.
:param SimPEG survey object survey:
:param v vector of the data in order matching of the survey
@@ -461,7 +461,7 @@ class Data(SimPEGsurvey.Data):
def toRecArray(self,returnType='RealImag'):
'''
Function that returns a numpy.recarray for a SimpegMT impedance data object.
Function that returns a numpy.recarray for a SimpegNSEM 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')
@@ -483,7 +483,7 @@ class Data(SimPEGsurvey.Data):
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
# Get the type and the value for the DataNSEM 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):
@@ -517,17 +517,17 @@ class Data(SimPEGsurvey.Data):
@classmethod
def fromRecArray(cls, recArray, srcType='primary'):
"""
Class method that reads in a numpy record array to MTdata object.
Class method that reads in a numpy record array to NSEMdata object.
Only imports the impedance data.
"""
if srcType=='primary':
src = SrcMT.polxy_1Dprimary
src = SrcNSEM.polxy_1Dprimary
elif srcType=='total':
src = SrcMT.polxy_1DhomotD
src = SrcNSEM.polxy_1DhomotD
else:
raise NotImplementedError('{:s} is not a valid source type for MTdata')
raise NotImplementedError('{:s} is not a valid source type for NSEMdata')
# Find all the frequencies in recArray
uniFreq = np.unique(recArray['freq'])
SimPEG/MT/Utils/MT1Danalytic.py → SimPEG/NSEM/Utils/MT1Danalytic.py
View File
SimPEG/MT/Utils/MT1Dsolutions.py → SimPEG/NSEM/Utils/MT1Dsolutions.py
View File
SimPEG/MT/Utils/__init__.py → SimPEG/NSEM/Utils/__init__.py
View File
SimPEG/MT/Utils/dataUtils.py → SimPEG/NSEM/Utils/dataUtils.py
+29 -29
View File
@@ -5,25 +5,25 @@ import numpy.lib.recfunctions as recFunc
from scipy.constants import mu_0
from scipy import interpolate as sciint
def getAppRes(MTdata):
def getAppRes(NSEMdata):
# Make impedance
zList = []
for src in MTdata.survey.srcList:
for src in NSEMdata.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])
zc.append(m*NSEMdata[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):
def rotateData(NSEMdata,rotAngle):
'''
Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
'''
recData = MTdata.toRecArray('Complex')
recData = NSEMdata.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')
@@ -40,8 +40,8 @@ def rotateData(MTdata,rotAngle):
for nr,comp in enumerate(['zxx','zxy','zyx','zyy']):
outRec[comp] = rotData[:,nr]
from SimPEG import MT
return MT.Data.fromRecArray(outRec)
from SimPEG import NSEM
return NSEM.Data.fromRecArray(outRec)
def appResPhs(freq,z):
@@ -57,10 +57,10 @@ def rec2ndarr(x,dt=float):
return x.view((dt, len(x.dtype.names)))
def makeAnalyticSolution(mesh,model,elev,freqs):
from SimPEG import MT
from SimPEG import NSEM
data1D = []
for freq in freqs:
anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
anaEd, anaEu, anaHd, anaHu = NSEM.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
anaE = anaEd+anaEu
anaH = anaHd+anaHu
@@ -71,7 +71,7 @@ def makeAnalyticSolution(mesh,model,elev,freqs):
return dataRec
def plotMT1DModelData(problem,models,symList=None):
from SimPEG import MT
from SimPEG import NSEM
# Setup the figure
fontSize = 15
@@ -214,11 +214,11 @@ 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
def convert3Dto1Dobject(NSEMdata,rxType3D='zyx'):
from SimPEG import NSEM
# Find the unique locations
# Need to find the locations
recDataTemp = MTdata.toRecArray()
recDataTemp = NSEMdata.toRecArray()
# Check if survey.std has been assigned.
## NEED TO: write this...
# Calculte and add the DET of the tensor to the recArray
@@ -240,24 +240,24 @@ def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
# Make the receiver list
rx1DList = []
for rxType in ['z1dr','z1di']:
rx1DList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
rx1DList.append(NSEM.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))
src1DList.append(NSEM.SrcNSEM.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)
sur1D = NSEM.Survey(src1DList)
# Make the data
dataVec = np.hstack(dat1DList)
dat1D = MT.Data(sur1D,dataVec)
dat1D = NSEM.Data(sur1D,dataVec)
sur1D.dobs = dataVec
# Need to take MTdata.survey.std and split it as well.
# Need to take NSEMdata.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)
@@ -265,29 +265,29 @@ def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
# Return the the list of data.
return mtData1DList
def resampleMTdataAtFreq(MTdata,freqs):
def resampleNSEMdataAtFreq(NSEMdata,freqs):
"""
Function to resample MTdata at set of frequencies
Function to resample NSEMdata at set of frequencies
"""
from SimPEG import MT
from SimPEG import NSEM
# Make a rec array
MTrec = MTdata.toRecArray().data
NSEMrec = NSEMdata.toRecArray().data
# Find unique locations
uniLoc = np.unique(MTrec[['x','y','z']])
uniFreq = MTdata.survey.freqs
uniLoc = np.unique(NSEMrec[['x','y','z']])
uniFreq = NSEMdata.survey.freqs
# Get the comps
dNames = MTrec.dtype
dNames = NSEMrec.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
ind = np.sqrt(np.sum((rec2ndarr(NSEMrec[['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)
int1d = sciint.interp1d(NSEMrec[ind]['freq'],NSEMrec[ind][comp],bounds_error=False)
tArrRec[comp] = simpeg.mkvc(int1d(freqs),2)
# Join together
@@ -296,5 +296,5 @@ def resampleMTdataAtFreq(MTdata,freqs):
except NameError as e:
outRecArr = tArrRec
# Make the MTdata and return
return MT.Data.fromRecArray(outRecArr)
# Make the NSEMdata and return
return NSEM.Data.fromRecArray(outRecArr)
SimPEG/MT/Utils/ediFilesUtils.py → SimPEG/NSEM/Utils/ediFilesUtils.py
+1 -1
View File
@@ -2,7 +2,7 @@
from SimPEG import mkvc
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from SimPEG.MT.Utils.dataUtils import rec2ndarr
from SimPEG.NSEM.Utils.dataUtils import rec2ndarr
# Import modules
import numpy as np
SimPEG/MT/Utils/plotDataTypes.py → SimPEG/NSEM/Utils/plotDataTypes.py
View File
SimPEG/MT/Utils/sourceUtils.py → SimPEG/NSEM/Utils/sourceUtils.py
+2 -2
View File
@@ -12,7 +12,7 @@ def homo1DModelSource(mesh,freq,sigma_1d):
'''
# import
from SimPEG.MT.Utils import get1DEfields
from SimPEG.NSEM.Utils import get1DEfields
# Get a 1d solution for a halfspace background
if mesh.dim == 1:
mesh1d = mesh
@@ -77,7 +77,7 @@ def analytic1DModelSource(mesh,freq,sigma_1d):
'''
# import
from SimPEG.MT.Utils import getEHfields
from SimPEG.NSEM.Utils import getEHfields
# Get a 1d solution for a halfspace background
if mesh.dim == 1:
mesh1d = mesh
SimPEG/MT/Utils/srcUtils.py → SimPEG/NSEM/Utils/srcUtils.py
+1 -1
View File
@@ -15,7 +15,7 @@ def homo1DModelSource(mesh,freq,m_back):
'''
# import
from SimPEG.MT.Utils import get1DEfields
from SimPEG.NSEM.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
SimPEG/NSEM/__init__.py
+5
View File
@@ -0,0 +1,5 @@
import Utils
from SurveyNSEM import Rx, Survey, Data
from FieldsNSEM import Fields1D_e, Fields3D_e
import Problem1D, Problem2D, Problem3D
import SrcNSEM
tests/mt/test_ApparentResistivityAnalytic.py
+2 -2
View File
@@ -1,6 +1,6 @@
import unittest
from SimPEG import *
from SimPEG import MT
from SimPEG import NSEM
TOL = 1e-6
@@ -20,7 +20,7 @@ def appResNorm(sigmaHalf):
freqs = np.logspace(4,-4,nFreq)
Z = []
for freq in freqs:
Ed, Eu, Hd, Hu = MT.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
Ed, Eu, Hd, Hu = NSEM.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
Z.append((Ed + Eu)/(Hd + Hu))
Zarr = np.concatenate(Z)
tests/mt/test_Problem1D_againstAnalyticHalfspace.py
+12 -12
View File
@@ -1,6 +1,6 @@
import unittest
import SimPEG as simpeg
from SimPEG import MT
from SimPEG import NSEM
from SimPEG.Utils import meshTensor
import numpy as np
# Define the tolerances
@@ -28,34 +28,34 @@ def setupSurvey(sigmaHalf,tD=True):
rxList = []
for rxType in ['z1dr','z1di']:
rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
rxList.append(NSEM.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))
srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq))
else:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
survey = MT.Survey(srcList)
survey = NSEM.Survey(srcList)
return survey, sigma, m1d
def getAppResPhs(MTdata):
def getAppResPhs(NSEMdata):
# 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:
for src in NSEMdata.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])
zc.append(m*NSEMdata[src,rx])
zList.append(zc)
return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
@@ -63,7 +63,7 @@ def appRes_TotalFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf)
problem = MT.Problem1D.eForm_TotalField(mesh)
problem = NSEM.Problem1D.eForm_TotalField(mesh)
problem.pair(survey)
# Get the fields
@@ -81,7 +81,7 @@ def appPhs_TotalFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf)
problem = MT.Problem1D.eForm_TotalField(mesh)
problem = NSEM.Problem1D.eForm_TotalField(mesh)
problem.pair(survey)
# Get the fields
@@ -99,7 +99,7 @@ def appRes_psFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf,False)
problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
problem = NSEM.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
problem.pair(survey)
# Get the fields
@@ -117,7 +117,7 @@ def appPhs_psFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf,False)
problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
problem = NSEM.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
problem.pair(survey)
# Get the fields
tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py
+14 -14
View File
@@ -1,6 +1,6 @@
import unittest
import SimPEG as simpeg
from SimPEG import MT
from SimPEG import NSEM
from SimPEG.Utils import meshTensor
import numpy as np
# Define the tolerances
@@ -34,43 +34,43 @@ def setupSurvey(sigmaHalf,tD=True):
rxList = []
for rxType in ['z1dr','z1di']:
rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
rxList.append(NSEM.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))
srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq))
else:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
survey = MT.Survey(srcList)
survey = NSEM.Survey(srcList)
return survey, sigma, m1d
def getAppResPhs(MTdata):
def getAppResPhs(NSEMdata):
# 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:
for src in NSEMdata.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])
zc.append(m*NSEMdata[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)
surveyAna = NSEM.Survey(srcList)
data1D = NSEM.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)
anaEd, anaEu, anaHd, anaHu = NSEM.Utils.NSEM1Danalytic.getEHfields(mesh,model,src.freq,elev)
anaE = anaEd+anaEu
anaH = anaHd+anaHu
# Scale the solution
@@ -87,11 +87,11 @@ def dataMis_AnalyticTotalDomain(sigmaHalf):
# Total domain solution
surveyTD, sigma, mesh = setupSurvey(sigmaHalf)
problemTD = MT.Problem1D.eForm_TotalField(mesh)
problemTD = NSEM.Problem1D.eForm_TotalField(mesh)
problemTD.pair(surveyTD)
# Analytic data
dataAnaObj = calculateAnalyticSolution(surveyTD.srcList,mesh,sigma)
# dataTDObj = MT.DataMT.DataMT(surveyTD, surveyTD.dpred(sigma))
# dataTDObj = NSEM.DataNSEM.DataNSEM(surveyTD, surveyTD.dpred(sigma))
dataTD = surveyTD.dpred(sigma)
dataAna = simpeg.mkvc(dataAnaObj)
return np.all((dataTD - dataAna)/dataAna < 2.)
@@ -109,7 +109,7 @@ def dataMis_AnalyticPrimarySecondary(sigmaHalf):
# Make the survey
# Primary secondary
surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,tD=False)
problemPS = MT.Problem1D.eForm_psField(mesh)
problemPS = NSEM.Problem1D.eForm_psField(mesh)
problemPS.sigmaPrimary = sigmaPS
problemPS.pair(surveyPS)
# Analytic data
tests/mt/test_Problem3D_againstAnalytic.py
+18 -18
View File
@@ -3,7 +3,7 @@ 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 import NSEM
from SimPEG.Utils import meshTensor
from scipy.constants import mu_0
@@ -98,40 +98,40 @@ def twoLayer(conds):
def setupSimpegMTfwd_eForm_ps(inputSetup,comp='Imp',singleFreq=False,expMap=True):
def setupSimpegNSEMfwd_eForm_ps(inputSetup,comp='Imp',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','tzxr','tzxi','tzyr','tzyi']:
rxList.append(MT.Rx(rx_loc,rxType))
rxList.append(NSEM.Rx(rx_loc,rxType))
elif comp == 'Imp':
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
rxList.append(MT.Rx(rx_loc,rxType))
rxList.append(NSEM.Rx(rx_loc,rxType))
elif comp == 'Tip':
for rxType in ['tzxr','tzxi','tzyr','tzyi']:
rxList.append(MT.Rx(rx_loc,rxType))
rxList.append(NSEM.Rx(rx_loc,rxType))
else:
rxList.append(MT.Rx(rx_loc,comp))
rxList.append(NSEM.Rx(rx_loc,comp))
# Source list
srcList =[]
if singleFreq:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,singleFreq))
srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,singleFreq))
else:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
# Survey MT
survey = MT.Survey(srcList)
srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
# Survey NSEM
survey = NSEM.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 = NSEM.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 = NSEM.Problem3D.eForm_ps(M,sigmaPrimary= sigma1d)
problem.curModel = sig
problem.pair(survey)
problem.verbose = False
@@ -143,18 +143,18 @@ def setupSimpegMTfwd_eForm_ps(inputSetup,comp='Imp',singleFreq=False,expMap=True
return (survey, problem)
def getAppResPhs(MTdata):
def getAppResPhs(NSEMdata):
# 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')
recData = NSEMdata.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)
survey, problem = setupSimpegNSEMfwd_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
@@ -174,7 +174,7 @@ def JvecAdjointTest(inputSetup,comp='All',freq=False):
# 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)
survey, problem = setupSimpegNSEMfwd_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)
@@ -191,7 +191,7 @@ def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True):
def DerivProjfieldsTest(inputSetup,comp='All',freq=False):
survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp,freq)
survey, problem = setupSimpegNSEMfwd_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
@@ -220,7 +220,7 @@ def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False):
else:
label = 'phase'
# Make the survey and the problem
survey, problem = setupSimpegMTfwd_eForm_ps(halfSpace(sigmaHalf),expMap=expMap)
survey, problem = setupSimpegNSEMfwd_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))