mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-29 22:51:44 +08:00
Updating Survey and Problem MT, work in progress. Code not tested.
This commit is contained in:
Gudni Karl Rosenkjaer
+27
-35
@@ -49,6 +49,13 @@ class MTProblem(Problem.BaseProblem):
|
||||
|
||||
# TODO:
|
||||
# MeSigmaBG
|
||||
@property
|
||||
def MeSigmaBG(self):
|
||||
#TODO: hardcoded to sigma as the model
|
||||
if getattr(self, '_MeSigma', None) is None:
|
||||
sigmaBG = self.curTModelBG
|
||||
self._MeSigmaBG = self.mesh.getEdgeInnerProduct(sigmaBG)
|
||||
return self._MeSigmaBG
|
||||
|
||||
@property
|
||||
def MeSigmaI(self):
|
||||
@@ -73,8 +80,13 @@ class MTProblem(Problem.BaseProblem):
|
||||
return self._curTModelDeriv
|
||||
|
||||
def fields(self, m):
|
||||
'''
|
||||
Function to calculate all the fields for the model m.
|
||||
|
||||
:param np.ndarray (nC,) m: f
|
||||
'''
|
||||
self.curModel = m
|
||||
RHS, CalcFields = self.getRHS, self.calcFields
|
||||
RHS, CalcFields = self.getRHS(freq), self.calcFields
|
||||
|
||||
F = FieldsMT(self.mesh, self.survey)
|
||||
|
||||
@@ -82,8 +94,9 @@ class MTProblem(Problem.BaseProblem):
|
||||
A = self.getA(freq)
|
||||
rhs = RHS(freq)
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
e = Ainv * rhs # is this e?
|
||||
e = Ainv * rhs
|
||||
|
||||
# Store the fields
|
||||
Src = self.survey.getSources(freq)
|
||||
# Stroe the fields
|
||||
F[Src, 'e'] = e
|
||||
@@ -105,6 +118,7 @@ class MTProblem(Problem.BaseProblem):
|
||||
C = self.mesh.edgeCurl
|
||||
|
||||
return C.T*mui*C + 1j*omega(freq)*sig
|
||||
|
||||
def getAbg(self, freq):
|
||||
"""
|
||||
Function to get the A matrix for the background model.
|
||||
@@ -129,45 +143,23 @@ class MTProblem(Problem.BaseProblem):
|
||||
|
||||
return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
|
||||
|
||||
def getRHS(self, freq):
|
||||
def getRHS(self, freq, backSigma):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
:param numpy.ndarray (nC,) backSigma: Background conductivity model
|
||||
:rtype: numpy.ndarray (nE, 2)
|
||||
:return: one RHS for both polarizations
|
||||
"""
|
||||
raise NotImplementedError()
|
||||
# Get sources for the frequency
|
||||
src = self.survey.getSources(freq)
|
||||
# Make sure that there is 2 polarizations.
|
||||
assert
|
||||
# Get the background electric fields
|
||||
from simpegMT.Sources import homo1DModelSource
|
||||
eBG_bp = home1DModelSource(self.mesh,freq,backSigma)
|
||||
Abg = getAbg(freq)
|
||||
|
||||
"""
|
||||
Put this in MT.Sources.EldadsSource
|
||||
|
||||
"""
|
||||
|
||||
Txs = self.survey.getTransmitters(freq)
|
||||
|
||||
# assert that only one Tx/src?
|
||||
# Create the two polarizations at this freq and return np array (nE,2).
|
||||
|
||||
# solve analytic.... get p1 p2
|
||||
|
||||
# Abg * [p1,p2] = rhs
|
||||
|
||||
rhs = range(len(Txs))
|
||||
for i, tx in enumerate(Txs):
|
||||
if tx.txType == 'VMD': # EH source.
|
||||
src = Sources.MagneticDipoleVectorPotential # this is where you would put multiple types of boundary conditions.
|
||||
else:
|
||||
raise NotImplemented('%s txType is not implemented' % tx.txType)
|
||||
SRCx = src(tx.loc, self.mesh.gridEx, 'x')
|
||||
SRCy = src(tx.loc, self.mesh.gridEy, 'y')
|
||||
SRCz = src(tx.loc, self.mesh.gridEz, 'z')
|
||||
rhs[i] = np.concatenate((SRCx, SRCy, SRCz))
|
||||
|
||||
a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F')
|
||||
mui = self.MfMui
|
||||
C = self.mesh.edgeCurl
|
||||
|
||||
j_s = C.T*mui*C*a
|
||||
return -1j*omega(freq)*j_s
|
||||
return Abg*eBG_bp
|
||||
|
||||
##################################################################
|
||||
# Inversion stuff
|
||||
|
||||
@@ -0,0 +1 @@
|
||||
from backgroundModelSources import *
|
||||
@@ -1,31 +1,32 @@
|
||||
def homo1DModelSource(mesh,bgMod):
|
||||
'''
|
||||
Function that calculates and return backround fields
|
||||
def homo1DModelSource(mesh,freq,bgMod):
|
||||
'''
|
||||
Function that calculates and return backround fields
|
||||
|
||||
'''
|
||||
'''
|
||||
|
||||
# import
|
||||
# import
|
||||
from simpegMT.Utils import get1DEfields
|
||||
# Get a 1d solution for a halfspace background
|
||||
mesh1d = simpeg.Mesh.TensorMesh([M.hz],np.array([M.x0[2]]))
|
||||
e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)
|
||||
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
|
||||
e0_1d = get1DEfields(mesh1d,mesh.r(sigBG,'CC','CC','M')[0,0,:],freq)
|
||||
# Setup x (east) polarization (_x)
|
||||
ex_px = np.zeros(M.vnEx,dtype=complex)
|
||||
ey_px = np.zeros((M.nEy,1),dtype=complex)
|
||||
ez_px = np.zeros((M.nEz,1),dtype=complex)
|
||||
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 arange(M.vnEx[0]):
|
||||
for j in arange(M.vnEx[2]):
|
||||
for i in arange(mesh.vnEx[0]):
|
||||
for j in arange(mesh.vnEx[1]):
|
||||
ex_px[i,j,:] = e0_1d
|
||||
eBG_px = np.vstack((simpeg.Utils.mkvc(M.r(ex_px,'Ex','Ex','V'),2),ey_px,ez_px))
|
||||
eBG_px = np.vstack((simpeg.Utils.mkvc(mesh.r(ex_px,'Ex','Ex','V'),2),ey_px,ez_px))
|
||||
# Setup y (north) polarization (_py)
|
||||
ex_py = np.zeros(M.nEx, dtype='complex128')
|
||||
ey_py = np.zeros((M.vnEy), dtype='complex128')
|
||||
ez_py = np.zeros(M.nEz, dtype='complex128')
|
||||
# Assign the source to ey_py
|
||||
for i in arange(M.vnEy[0]):
|
||||
for j in arange(M.vnEy[1]):
|
||||
ey_py[i,j,:] = e0_1d
|
||||
ex_py = np.zeros(mesh.nEx, dtype='complex128')
|
||||
ey_py = np.zeros((mesh.vnEy), dtype='complex128')
|
||||
ez_py = np.zeros(mesh.nEz, dtype='complex128')
|
||||
# Assign the source to ey_py
|
||||
for i in arange(mesh.vnEy[0]):
|
||||
for j in arange(mesh.vnEy[1]):
|
||||
ey_py[i,j,:] = e0_1d
|
||||
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
|
||||
|
||||
eBG_py = np.r_[ex_py,simpeg.Utils.mkvc(ey_py),ez_py]
|
||||
|
||||
# Return the electric fields
|
||||
return np.hstack((eBG_px,eBG_py))
|
||||
+78
-42
@@ -1,11 +1,17 @@
|
||||
from SimPEG import Survey, Utils, np, sp
|
||||
|
||||
class RxFDEM(Survey.BaseRx):
|
||||
class RxMT(Survey.BaseRx):
|
||||
|
||||
knownRxTypes = {
|
||||
'exr':['e', 'Ex', 'real'],
|
||||
'eyr':['e', 'Ey', 'real'],
|
||||
'ezr':['e', 'Ez', 'real'],
|
||||
'zxxr':[['e', 'Ex'],['b','Fx'], 'real'],
|
||||
'zxyr':[['e', 'Ex'],['b','Fy'], 'real'],
|
||||
'zyxr':[['e', 'Ey'],['b','Fx'], 'real'],
|
||||
'zyyr':[['e', 'Ey'],['b','Fy'], 'real'],
|
||||
'zxxi':[['e', 'Ex'],['b','Fx'], 'imag'],
|
||||
'zxyi':[['e', 'Ex'],['b','Fy'], 'imag'],
|
||||
'zyxi':[['e', 'Ey'],['b','Fx'], 'imag'],
|
||||
'zyyi':[['e', 'Ey'],['b','Fy'], 'imag'],
|
||||
|
||||
'exi':['e', 'Ex', 'imag'],
|
||||
'eyi':['e', 'Ey', 'imag'],
|
||||
'ezi':['e', 'Ez', 'imag'],
|
||||
@@ -22,29 +28,56 @@ class RxFDEM(Survey.BaseRx):
|
||||
Survey.BaseRx.__init__(self, locs, rxType)
|
||||
|
||||
@property
|
||||
def projField(self):
|
||||
"""Field Type projection (e.g. e b ...)"""
|
||||
return self.knownRxTypes[self.rxType][0]
|
||||
def projField(self,fracPos):
|
||||
"""
|
||||
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 ...)"""
|
||||
return self.knownRxTypes[self.rxType][1]
|
||||
def projGLoc(self,fracPos):
|
||||
"""
|
||||
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 projComp(self):
|
||||
"""Component projection (real/imag)"""
|
||||
return self.knownRxTypes[self.rxType][2]
|
||||
|
||||
def projectFields(self, tx, mesh, u):
|
||||
P = self.getP(mesh)
|
||||
u_part_complex = u[tx, self.projField]
|
||||
def projectFields(self, src, mesh, u):
|
||||
'''
|
||||
Project the fields and return the
|
||||
'''
|
||||
# Get the numerator information
|
||||
P_num = self.getP(mesh,self.projGLoc('numerator'))
|
||||
u_num_complex = u[src, self.projField('numerator')]
|
||||
# Get the denominator information
|
||||
P_den = self.getP(mesh,self.projGLoc('denominator'))
|
||||
u_den_complex = u[src, self.projField('denominator')]
|
||||
# Calculate the fraction
|
||||
f_part_complex = (P_num*u_num_complex)/(P_den*u_den_complex)
|
||||
# get the real or imag component
|
||||
real_or_imag = self.projComp
|
||||
u_part = getattr(u_part_complex, real_or_imag)
|
||||
return P*u_part
|
||||
f_part = getattr(f_part_complex, real_or_imag)
|
||||
return f_part
|
||||
|
||||
def projectFieldsDeriv(self, tx, mesh, u, v, adjoint=False):
|
||||
def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False):
|
||||
P = self.getP(mesh)
|
||||
|
||||
if not adjoint:
|
||||
@@ -66,13 +99,16 @@ class RxFDEM(Survey.BaseRx):
|
||||
|
||||
|
||||
# Call this Source or polarization or something...?
|
||||
class TxFDEM(Survey.BaseTx):
|
||||
class srcMT(Survey.BaseTx):
|
||||
'''
|
||||
Sources for the MT problem.
|
||||
'''
|
||||
|
||||
freq = None #: Frequency (float)
|
||||
|
||||
rxPair = RxFDEM
|
||||
rxPair = RxMT
|
||||
|
||||
knownTxTypes = ['VMD'] # Polarization...
|
||||
knownSrcTypes = ['ORTPOL'] # ORThogonal POLarization
|
||||
|
||||
def __init__(self, loc, txType, freq, rxList): # remove txType? hardcode to one thing. always polarizations
|
||||
self.freq = float(freq)
|
||||
@@ -81,29 +117,29 @@ class TxFDEM(Survey.BaseTx):
|
||||
|
||||
|
||||
|
||||
class FieldsFDEM(Survey.Fields):
|
||||
"""Fancy Field Storage for a FDEM survey."""
|
||||
class FieldsMT(Survey.Fields):
|
||||
"""Fancy Field Storage for a MT survey."""
|
||||
knownFields = {'b': 'F', 'e': 'E'}
|
||||
dtype = complex
|
||||
|
||||
|
||||
class SurveyFDEM(Survey.BaseSurvey):
|
||||
class SurveyMT(Survey.BaseSurvey):
|
||||
"""
|
||||
docstring for SurveyFDEM
|
||||
docstring for SurveyMT
|
||||
"""
|
||||
|
||||
txPair = TxFDEM
|
||||
srcPair = srcMT
|
||||
|
||||
def __init__(self, txList, **kwargs):
|
||||
def __init__(self, srcList, **kwargs):
|
||||
# Sort these by frequency
|
||||
self.txList = txList
|
||||
self.srcList = srcList
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
|
||||
_freqDict = {}
|
||||
for tx in txList:
|
||||
if tx.freq not in _freqDict:
|
||||
_freqDict[tx.freq] = []
|
||||
_freqDict[tx.freq] += [tx]
|
||||
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])
|
||||
@@ -117,26 +153,26 @@ class SurveyFDEM(Survey.BaseSurvey):
|
||||
def nFreq(self):
|
||||
"""Number of frequencies"""
|
||||
return len(self._freqDict)
|
||||
|
||||
@property
|
||||
def nTxByFreq(self):
|
||||
if getattr(self, '_nTxByFreq', None) is None:
|
||||
self._nTxByFreq = {}
|
||||
for freq in self.freqs:
|
||||
self._nTxByFreq[freq] = len(self.getTransmitters(freq))
|
||||
return self._nTxByFreq
|
||||
# Don't need this
|
||||
# @property
|
||||
# def nTxByFreq(self):
|
||||
# if getattr(self, '_nTxByFreq', None) is None:
|
||||
# self._nTxByFreq = {}
|
||||
# for freq in self.freqs:
|
||||
# self._nTxByFreq[freq] = len(self.getTransmitters(freq))
|
||||
# return self._nTxByFreq
|
||||
|
||||
# TODO: Rename to getSources
|
||||
def getTransmitters(self, freq):
|
||||
def getSources(self, freq):
|
||||
"""Returns the transmitters 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 = Survey.Data(self)
|
||||
for tx in self.txList:
|
||||
for rx in tx.rxList:
|
||||
data[tx, rx] = rx.projectFields(tx, self.mesh, u)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.projectFields(src, self.mesh, u)
|
||||
return data
|
||||
|
||||
def projectFieldsDeriv(self, u):
|
||||
|
||||
@@ -1,6 +1,6 @@
|
||||
# from EM import *
|
||||
import Utils
|
||||
# import FDEM
|
||||
# import Sources
|
||||
# import PromblemMT
|
||||
# import SurveyMT
|
||||
# import Tests
|
||||
import Sources
|
||||
import ProblemMT
|
||||
import SurveyMT
|
||||
|
||||
Reference in New Issue
Block a user