diff --git a/simpegMT/BaseMT.py b/simpegMT/BaseMT.py index 9ec143ec..973147fc 100644 --- a/simpegMT/BaseMT.py +++ b/simpegMT/BaseMT.py @@ -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 diff --git a/simpegMT/DataMT.py b/simpegMT/DataMT.py index 144d79a3..03b88045 100644 --- a/simpegMT/DataMT.py +++ b/simpegMT/DataMT.py @@ -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.') diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 5a0f7abb..766e0472 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -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 diff --git a/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py b/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py index 850ecf40..5ab70fd4 100644 --- a/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py +++ b/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py @@ -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() \ No newline at end of file