From 9b650041d1fa34b68b2d0c003997679d5635dce8 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Tue, 17 Feb 2015 00:09:24 -0800 Subject: [PATCH] Working on MT problem and survey. Fixed bugs and completed example --- simpegMT/Examples/simple3DfowardProblem.py | 33 ++++++++++- simpegMT/ProblemMT.py | 65 ++++++++++++++-------- simpegMT/Sources/backgroundModelSources.py | 24 ++++---- simpegMT/SurveyMT.py | 51 +++++++++++++---- 4 files changed, 126 insertions(+), 47 deletions(-) diff --git a/simpegMT/Examples/simple3DfowardProblem.py b/simpegMT/Examples/simple3DfowardProblem.py index ee0f0bad..1fe8842e 100644 --- a/simpegMT/Examples/simple3DfowardProblem.py +++ b/simpegMT/Examples/simple3DfowardProblem.py @@ -1,9 +1,36 @@ # Test script to use simpegMT platform to forward model synthetic data. # Import -import simpegMT as simpegMT, SimPEG as simpeg - +import simpegMT as simpegmt, SimPEG as simpeg +import numpy as np +# Make a mesh +M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360]) +# Setup the model +conds = [1e-2,1] +sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,[-1000,-1000,-400],[1000,1000,-200],conds) +sig[M.gridCC[:,2]>0] = 1e-8 +sig[M.gridCC[:,2]<-600] = 1e-1 +sigBG = np.zeros(M.nC) + conds[0] +sigBG[M.gridCC[:,2]>0] = 1e-8 +## Setup the the survey object +# Receiver locations +rx_x, rx_y = np.meshgrid(np.arange(-500,501,50),np.arange(-500,501,50)) +rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),np.zeros((np.prod(rx_x.shape),1)))) # Make a receiver list -for \ No newline at end of file +rxList = [] +for loc in rx_loc: + for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']: + rxList.append(simpegmt.SurveyMT.RxMT(loc,rxType)) +# Source list +srcList =[] +for freq in np.logspace(3,-1,21): + srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList)) +# Survey MT +survey = simpegmt.SurveyMT.SurveyMT(srcList) + +## Setup the problem object +problem = simpegmt.ProblemMT.MTProblem(M) +problem.pair(survey) + diff --git a/simpegMT/ProblemMT.py b/simpegMT/ProblemMT.py index 04fabe30..1d6f51ce 100644 --- a/simpegMT/ProblemMT.py +++ b/simpegMT/ProblemMT.py @@ -1,4 +1,4 @@ -from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver +from SimPEG import Survey, Problem, Utils, Models, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 from SurveyMT import SurveyMT, FieldsMT @@ -47,15 +47,6 @@ class MTProblem(Problem.BaseProblem): self._MeSigma = self.mesh.getEdgeInnerProduct(sigma) return self._MeSigma - # 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): @@ -70,38 +61,68 @@ class MTProblem(Problem.BaseProblem): @property def curTModel(self): if getattr(self, '_curTModel', None) is None: - self._curTModel = self.mapping.transform(self.curModel) + self._curTModel = self.mapping*self.curModel return self._curTModel @property def curTModelDeriv(self): if getattr(self, '_curTModelDeriv', None) is None: - self._curTModelDeriv = self.mapping.transformDeriv(self.curModel) + self._curTModelDeriv = self.mapping*self.curModel return self._curTModelDeriv - def fields(self, m): + + # Background model + @property + def backModel(self): + """ + Sets the model, and removes dependent mass matrices. + """ + return getattr(self, '_backModel', None) + + @backModel.setter + def backModel(self, value): + if value is self.backModel: + return # it is the same! + self._backModel = Models.Model(value, self.mapping) + for prop in self.deleteTheseOnModelUpdate: + if hasattr(self, prop): + delattr(self, prop) + + @property + def MeSigmaBG(self): + #TODO: hardcoded to sigma as the model + if getattr(self, '_MeSigmaBG', None) is None: + sigmaBG = self.backModel + self._MeSigmaBG = self.mesh.getEdgeInnerProduct(sigmaBG) + return self._MeSigmaBG + + def fields(self, m, m_back): ''' Function to calculate all the fields for the model m. - :param np.ndarray (nC,) m: f + :param np.ndarray (nC,) m: Conductivity model + :param np.ndarray (nC,) m_back: Background conductivity model ''' self.curModel = m - RHS, CalcFields = self.getRHS(freq), self.calcFields + self.backModel = m_back + # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields F = FieldsMT(self.mesh, self.survey) for freq in self.survey.freqs: A = self.getA(freq) - rhs = RHS(freq) + rhs = self.getRHS(freq,m_back) Ainv = self.Solver(A, **self.solverOpts) e = Ainv * rhs # Store the fields Src = self.survey.getSources(freq) - # Stroe the fields - F[Src, 'e'] = e - F[Src, 'b'] = self.mesh.edgeCurl * e - + # Store the fields + F[Src, 'e_px'] = e[:,0] + F[Src, 'e_py'] = e[:,1] + b = self.mesh.edgeCurl * e + F[Src, 'b_px'] = b[:,0] + F[Src, 'b_py'] = b[:,1] return F @@ -157,8 +178,8 @@ class MTProblem(Problem.BaseProblem): # assert len() # Get the background electric fields from simpegMT.Sources import homo1DModelSource - eBG_bp = home1DModelSource(self.mesh,freq,backSigma) - Abg = getAbg(freq) + eBG_bp = homo1DModelSource(self.mesh,freq,backSigma) + Abg = self.getAbg(freq) return Abg*eBG_bp diff --git a/simpegMT/Sources/backgroundModelSources.py b/simpegMT/Sources/backgroundModelSources.py index f74b4ad1..072a35d5 100644 --- a/simpegMT/Sources/backgroundModelSources.py +++ b/simpegMT/Sources/backgroundModelSources.py @@ -1,36 +1,40 @@ -def homo1DModelSource(mesh,freq,bgMod): +import SimPEG as simpeg, numpy as np +from simpegMT.Utils import get1DEfields + +def homo1DModelSource(mesh,freq,m_back): ''' Function that calculates and return background fields :param Simpeg mesh object mesh: Holds information on the discretization :param float freq: The frequency to solve at - :param np.array bgMod: Background model to base the calculations on. + :param np.array m_back: Background model of conductivity to base the calculations on. :rtype: numpy.ndarray (mesh.nE,2) :return: eBG_bp, E fields for the background model at both polarizations. ''' # import + import SimPEG as simpeg from simpegMT.Utils import get1DEfields # Get a 1d solution for a halfspace background mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]])) - e0_1d = get1DEfields(mesh1d,mesh.r(sigBG,'CC','CC','M')[0,0,:],freq) + e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq) # Setup x (east) polarization (_x) 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(mesh.vnEx[0]): - for j in arange(mesh.vnEx[1]): + for i in np.arange(mesh.vnEx[0]): + for j in np.arange(mesh.vnEx[1]): ex_px[i,j,:] = e0_1d 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(mesh.nEx, dtype='complex128') - ey_py = np.zeros((mesh.vnEy), dtype='complex128') - ez_py = np.zeros(mesh.nEz, dtype='complex128') + ex_py = np.zeros((mesh.nEx,1), dtype='complex128') + ey_py = np.zeros(mesh.vnEy, dtype='complex128') + ez_py = np.zeros((mesh.nEz,1), dtype='complex128') # Assign the source to ey_py - for i in arange(mesh.vnEy[0]): - for j in arange(mesh.vnEy[1]): + for i in np.arange(mesh.vnEy[0]): + for j in np.arange(mesh.vnEy[1]): ey_py[i,j,:] = e0_1d eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 829d99b7..8604ed62 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -1,4 +1,5 @@ from SimPEG import Survey, Utils, Problem, np, sp +from scipy.constants import mu_0 class RxMT(Survey.BaseRx): @@ -65,14 +66,37 @@ class RxMT(Survey.BaseRx): ''' 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 projection + Pex = self.getP(mesh,'Ex') + Pey = self.getP(mesh,'Ey') + Pbx = self.getP(mesh,'Fx') + Pby = self.getP(mesh,'Fy') + # Get the fields at location + ex_px = Pex*u[src,'e_px'] + ey_px = Pey*u[src,'e_px'] + ex_py = Pex*u[src,'e_py'] + ey_py = Pey*u[src,'e_py'] + hx_px = Pbx*u[src,'b_px']/mu_0 + hy_px = Pby*u[src,'b_px']/mu_0 + hx_py = Pbx*u[src,'b_py']/mu_0 + hy_py = Pby*u[src,'b_py']/mu_0 + if 'zxx' in self.rxType: + 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) + elif 'zyx' in self.rxType: + 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) + + # 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 f_part = getattr(f_part_complex, real_or_imag) @@ -100,6 +124,7 @@ class RxMT(Survey.BaseRx): # Call this Source or polarization or something...? +# Note: Might need to add tests to make sure that both polarization have the same rxList. class srcMT(Survey.BaseTx): ''' Sources for the MT problem. @@ -107,24 +132,25 @@ class srcMT(Survey.BaseTx): :param float freq: The frequency of the source :param list rxList: A list of receivers associated with the source + :param str srcPol: The polarization of the source ''' freq = None #: Frequency (float) rxPair = RxMT - knownTxTypes = ['ORTPOL'] # ORThogonal POLarization + knownTxTypes = ['pol_xy','pol_x','pol_y'] # ORThogonal POLarization - def __init__(self, freq, rxList): # remove txType? hardcode to one thing. always polarizations + def __init__(self, freq, rxList, srcPol = 'pol_xy'): # remove txType? hardcode to one thing. always polarizations self.freq = float(freq) - Survey.BaseTx.__init__(self, None, 'ORTPOL', rxList) + Survey.BaseTx.__init__(self, None, srcPol, rxList) # Survey.BaseTx.__init__(self, loc, 'polarization', rxList) class FieldsMT(Problem.Fields): """Fancy Field Storage for a MT survey.""" - knownFields = {'b': 'F', 'e': 'E'} + knownFields = {'b_px': 'F','b_py': 'F', 'e_px': 'E','e_py': 'E'} dtype = complex @@ -141,6 +167,7 @@ class SurveyMT(Survey.BaseSurvey): def __init__(self, srcList, **kwargs): # Sort these by frequency self.srcList = srcList + self.txList = self.srcList # Hack - make txList index srcList, since it is used in the backend. Survey.BaseSurvey.__init__(self, **kwargs) _freqDict = {}