diff --git a/simpegMT/BaseMT.py b/simpegMT/BaseMT.py index f04fc5f1..7ebf66d3 100644 --- a/simpegMT/BaseMT.py +++ b/simpegMT/BaseMT.py @@ -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 diff --git a/simpegMT/DataMT.py b/simpegMT/DataMT.py deleted file mode 100644 index 6ad911f6..00000000 --- a/simpegMT/DataMT.py +++ /dev/null @@ -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) \ No newline at end of file diff --git a/simpegMT/Examples/simple3DfowardProblem.py b/simpegMT/Examples/simple3DfowardProblem.py index d4f73281..e03ee09d 100644 --- a/simpegMT/Examples/simple3DfowardProblem.py +++ b/simpegMT/Examples/simple3DfowardProblem.py @@ -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) diff --git a/simpegMT/FieldsMT.py b/simpegMT/FieldsMT.py index 206d982f..87ea109d 100644 --- a/simpegMT/FieldsMT.py +++ b/simpegMT/FieldsMT.py @@ -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) diff --git a/simpegMT/Problem1D/Probs.py b/simpegMT/Problem1D/Probs.py index a2c36bed..9b0f74a7 100644 --- a/simpegMT/Problem1D/Probs.py +++ b/simpegMT/Problem1D/Probs.py @@ -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() diff --git a/simpegMT/Problem3D/Probs.py b/simpegMT/Problem3D/Probs.py index 5bbafc62..fbe92f00 100644 --- a/simpegMT/Problem3D/Probs.py +++ b/simpegMT/Problem3D/Probs.py @@ -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() diff --git a/simpegMT/Problems/__init__.py b/simpegMT/Problems/__init__.py deleted file mode 100644 index ac5c4623..00000000 --- a/simpegMT/Problems/__init__.py +++ /dev/null @@ -1 +0,0 @@ -import 1D, 2D, 3D \ No newline at end of file diff --git a/simpegMT/SrcMT.py b/simpegMT/SrcMT.py new file mode 100644 index 00000000..e00dbd52 --- /dev/null +++ b/simpegMT/SrcMT.py @@ -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 diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 703030a7..b5ede631 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -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) + diff --git a/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py b/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py index 86d67538..e6f56800 100644 --- a/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py +++ b/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py @@ -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): diff --git a/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py b/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py index b91f8b30..a72a5998 100644 --- a/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py +++ b/simpegMT/Tests/test_Problem1D_totalDvsPSvsAnalytic.py @@ -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) diff --git a/simpegMT/Tests/test_Problem3D_againstAnalytic.py b/simpegMT/Tests/test_Problem3D_againstAnalytic.py index 37ff9dc1..2aabe2f3 100644 --- a/simpegMT/Tests/test_Problem3D_againstAnalytic.py +++ b/simpegMT/Tests/test_Problem3D_againstAnalytic.py @@ -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) diff --git a/simpegMT/Utils/dataUtils.py b/simpegMT/Utils/dataUtils.py index a2cd89cd..f0c96017 100644 --- a/simpegMT/Utils/dataUtils.py +++ b/simpegMT/Utils/dataUtils.py @@ -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) \ No newline at end of file + return simpegmt.Data.fromRecArray(outRecArr) \ No newline at end of file diff --git a/simpegMT/__init__.py b/simpegMT/__init__.py index 129d68d6..dfbd9a15 100644 --- a/simpegMT/__init__.py +++ b/simpegMT/__init__.py @@ -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 \ No newline at end of file