wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-03 14:02:31 +08:00
Code Issues Packages Projects Releases Wiki Activity

Implementing 3D derivatives

Browse Source
This commit is contained in:
GudniRos
2015-10-14 11:07:49 -07:00
parent 5eccffffb5
commit 1654c1c8b5
4 changed files with 46 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegMT/BaseMT.py
-1
View File
@@ -43,7 +43,6 @@ class BaseMTProblem(BaseFDEMProblem):
# Re-hook the methods to the object
Utils.codeutils.hook(self,Utils.codeutils.hook)
Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
self.
# Set the solver
Solver = SimpegSolver
simpegMT/DataMT.py
+5 -4
View File
@@ -35,8 +35,9 @@ class DataMT(Survey.Data):
'''
# 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)]
dtCP = [('freq',float),('x',float),('y',float),('z',float),('zxx',complex),('zxy',complex),('zyx',complex),('zyy',complex)]
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.
@@ -47,7 +48,7 @@ class DataMT(Survey.Data):
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],8))),axis=1).view(dtRI)
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]
@@ -72,7 +73,7 @@ class DataMT(Survey.Data):
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']:
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.')
simpegMT/SurveyMT.py
+32 -15
View File
@@ -26,8 +26,7 @@ class RxMT(Survey.BaseRx):
# TODO:
# 1D impedance
'z1dr':['Z1D', 'real'],
'z1di':['Z1D', 'imag']
#TODO: Add tipper fractions as well. Bz/B(x|y)
'z1di':['Z1D', 'imag'],
# Tipper
'tzxr':['T3D','real'],
'tzxi':['T3D','imag'],
@@ -93,11 +92,17 @@ class RxMT(Survey.BaseRx):
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(self.locs,'Ex')
Pey = mesh.getInterpolationMat(self.locs,'Ey')
Pbx = mesh.getInterpolationMat(self.locs,'Fx')
Pby = mesh.getInterpolationMat(self.locs,'Fy')
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']
@@ -110,17 +115,23 @@ class RxMT(Survey.BaseRx):
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)
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)
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)
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)
f_part_complex = (-ey_px*hx_py + ey_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
elif self.projType is 'T3D':
Pbx = mesh.getInterpolationMat(self.locs,'Fx')
Pby = mesh.getInterpolationMat(self.locs,'Fy')
Pbz = mesh.getInterpolationMat(self.locs,'Fz')
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']
@@ -130,7 +141,7 @@ class RxMT(Survey.BaseRx):
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)
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))
@@ -238,7 +249,13 @@ class srcMT_polxy_1Dprimary(srcMT):
def ePrimary(self,problem):
# Get primary fields for both polarizations
self.sigma1d = problem._sigmaPrimary
if self.sigma1d is None:
# Set the sigma1d as the 1st column in the background model
if len(problem._sigmaPrimary) == problem.mesh.nC:
self.sigma1d = problem.mesh.r(problem._sigmaPrimary,'CC','CC','M')[0,0,:]
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
simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py
+9 -29
View File
@@ -42,7 +42,7 @@ def setupSurvey(sigmaHalf,tD=True):
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))
else:
for freq in freqs:
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigmaBack))
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
survey = simpegmt.SurveyMT.SurveyMT(srcList)
return survey, sigma, m1d
@@ -108,19 +108,16 @@ def dataMis_AnalyticPrimarySecondary(sigmaHalf):
# Make the survey
# Primary secondary
surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,False)
problemPS = simpegmt.ProblemMT1D.eForm_psField(mesh,sigmaPS)
surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,tD=False)
problemPS = simpegmt.ProblemMT1D.eForm_psField(mesh)
problemPS.sigmaPrimary = sigmaPS
problemPS.pair(surveyPS)
# Analytic data
dataAna = calculateAnalyticSolution(surveyPS.srcList,mesh,sigma)
dataAnaObj = calculateAnalyticSolution(surveyPS.srcList,mesh,sigmaPS)
surveyPS.dtrue = dataAna
# Project the data
data = surveyPS.dpred(sigmaPS)
# Setup the data misfit
dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
return dmis.eval(sigma)
dataPS = surveyPS.dpred(sigmaPS)
dataAna = simpeg.mkvc(dataAnaObj)
return np.all((dataPS - dataAna)/dataAna < 2.)
@@ -129,27 +126,10 @@ class TestNumericVsAnalytics(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.assertTrue(dataMis_AnalyticTotalDomain(2e-2))
# def test_appPhs2en2(self):self.assert(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)
def test_appRes2en2_ps(self):self.assertTrue(dataMis_AnalyticPrimarySecondary(2e-2))
if __name__ == '__main__':
unittest.main()