diff --git a/simpegMT/ProblemMT.py b/simpegMT/ProblemMT.py index 7762086d..159c26eb 100644 --- a/simpegMT/ProblemMT.py +++ b/simpegMT/ProblemMT.py @@ -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 diff --git a/simpegMT/Sources/__init__.py b/simpegMT/Sources/__init__.py new file mode 100644 index 00000000..7b4a53bc --- /dev/null +++ b/simpegMT/Sources/__init__.py @@ -0,0 +1 @@ +from backgroundModelSources import * \ No newline at end of file diff --git a/simpegMT/Sources/backgroundModelSources.py b/simpegMT/Sources/backgroundModelSources.py index 4703dd96..d5e97750 100644 --- a/simpegMT/Sources/backgroundModelSources.py +++ b/simpegMT/Sources/backgroundModelSources.py @@ -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] - \ No newline at end of file + # Return the electric fields + return np.hstack((eBG_px,eBG_py)) \ No newline at end of file diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 4072a4ad..80a77d53 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -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): diff --git a/simpegMT/__init__.py b/simpegMT/__init__.py index e75af9a5..fe5cf927 100644 --- a/simpegMT/__init__.py +++ b/simpegMT/__init__.py @@ -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