From 422911a95f871b455628923b2400325ea84aa5d7 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Thu, 11 Jun 2015 16:26:11 -0700 Subject: [PATCH] Fixed 1D test and current code to work, where the src in the 1D problem is partly implemented --- simpegMT/ProblemMT1D/Problems.py | 28 ++-- simpegMT/ProblemMT3D/Problems.py | 129 +++++++++--------- simpegMT/Sources/backgroundModelSources.py | 18 ++- simpegMT/SurveyMT.py | 85 ++++++++++-- ...test_Problem1D_againstAnalyticHalfspace.py | 2 +- .../Tests/test_Problem3D_againstAnalytic.py | 9 +- simpegMT/Utils/MT1Dsolutions.py | 4 +- 7 files changed, 174 insertions(+), 101 deletions(-) diff --git a/simpegMT/ProblemMT1D/Problems.py b/simpegMT/ProblemMT1D/Problems.py index 35fee655..74c62292 100644 --- a/simpegMT/ProblemMT1D/Problems.py +++ b/simpegMT/ProblemMT1D/Problems.py @@ -10,8 +10,11 @@ import numpy as np import multiprocessing, sys, time +# class eForm_ps(BaseMTProblem): + + class eForm_TotalField(BaseMTProblem): - """ + """ A MT problem solving a e formulation and a primary/secondary fields decompostion. Solves the equation: @@ -21,8 +24,8 @@ class eForm_TotalField(BaseMTProblem): # From FDEMproblem: Used to project the fields. Currently not used for MTproblem. _fieldType = 'e' - _eqLocs = 'FE' - + _eqLocs = 'EF' + def __init__(self, mesh, **kwargs): BaseMTProblem.__init__(self, mesh, **kwargs) @@ -36,8 +39,14 @@ class eForm_TotalField(BaseMTProblem): :rtype: scipy.sparse.csr_matrix :return: A """ + Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0) - Msig = self.mesh.getFaceInnerProduct(self.curModel) + Msig = self.mesh.getFaceInnerProduct(self.curModel.sigma) + # Note: need to use the code above since in the 1D problem I want + # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this + # Possible that _fieldType and _eqLocs can fix this + # Mmui = self.MfMui + # Msig = self.MeSigma C = self.mesh.nodalGrad # Make A A = C.T*Mmui*C + 1j*omega(freq)*Msig @@ -66,12 +75,12 @@ class eForm_TotalField(BaseMTProblem): """ # Get sources for the frequency # NOTE: Need to use the source information, doesn't really apply in 1D - src = self.survey.getSources(freq) + src = self.survey.getSrcByFreq(freq) # Get the full A A = self.getA(freq,full=True) # Define the outer part of the solution matrix Aio = A[1:-1,[0,-1]] - Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel,freq,self.mesh.vectorNx) + Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx) Etot = (Ed + Eu) sourceAmp = 1.0 Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top @@ -104,12 +113,12 @@ class eForm_TotalField(BaseMTProblem): A = self.getA(freq) rhs, e_o = self.getRHS(freq) Ainv = self.Solver(A, **self.solverOpts) - e_i = Ainv * rhs + e_i = Ainv * rhs e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2) # Store the fields - Src = self.survey.getSources(freq) + Src = self.survey.getSrcByFreq(freq) # Store the fields - # NOTE: only store + # NOTE: only store F[Src, 'e_1d'] = e # F[Src, 'e_py'] = 0*e[:,0] # Note curl e = -iwb so b = -curl e /iw @@ -120,4 +129,3 @@ class eForm_TotalField(BaseMTProblem): print 'Ran for {:f} seconds'.format(time.time()-startTime) sys.stdout.flush() return F - \ No newline at end of file diff --git a/simpegMT/ProblemMT3D/Problems.py b/simpegMT/ProblemMT3D/Problems.py index efec3d64..b067ebba 100644 --- a/simpegMT/ProblemMT3D/Problems.py +++ b/simpegMT/ProblemMT3D/Problems.py @@ -4,49 +4,55 @@ from scipy.constants import mu_0 from simpegMT.BaseMT import BaseMTProblem from simpegMT.SurveyMT import SurveyMT from simpegMT.FieldsMT import FieldsMT -from simpegMT.DataMT import DataMT +from simpegMT.DataMT import DataMT import multiprocessing, sys, time class eForm_ps(BaseMTProblem): - """ + """ A MT problem solving a e formulation and a primary/secondary fields decompostion. - Solves the equation + Solves the equation: + + + """ # From FDEMproblem: Used to project the fields. Currently not used for MTproblem. _fieldType = 'e' _eqLocs = 'FE' - - # Set new properties + # Need to add the src .... + + + # Set new properties # Background model - @property - def backModel(self): - """ - Sets the model, and removes dependent mass matrices. - """ - return getattr(self, '_backModel', None) + # Shouldn't need the commented block. + # @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) + # @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 MeDeltaSigma(self): - #TODO: hardcoded to sigma as the model - if getattr(self, '_MeDeltaSigma', None) is None: - sigma = self.curModel - sigmaBG = self.backModel - self._MeDeltaSigma = self.mesh.getEdgeInnerProduct(sigma - sigmaBG) - return self._MeDeltaSigma + # @property + # def MeDeltaSigma(self): + # #TODO: hardcoded to sigma as the model + # if getattr(self, '_MeDeltaSigma', None) is None: + # sigma = self.curModel + # sigmaBG = self.backModel + # self._MeDeltaSigma = self.mesh.getEdgeInnerProduct(sigma - sigmaBG) + # return self._MeDeltaSigma def __init__(self, mesh, **kwargs): BaseMTProblem.__init__(self, mesh, **kwargs) @@ -59,56 +65,52 @@ class eForm_ps(BaseMTProblem): :rtype: scipy.sparse.csr_matrix :return: A """ - mui = self.MfMui - sig = self.MeSigma + Mmui = self.MfMui + Msig = self.MeSigma C = self.mesh.edgeCurl - return C.T*mui*C + 1j*omega(freq)*sig + return C.T*Mmui*C + 1j*omega(freq)*Msig def getADeriv(self, freq, u, v, adjoint=False): - sig = self.curTModel - dsig_dm = self.curTModelDeriv - dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u) + + dsig_dm = self.curModel.sigmaDeriv + dMe_dsig = self.MeSimgaDeriv( v=u) if adjoint: return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) ) return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) ) - def getRHS(self, freq, backSigma): + def getRHS(self, freq): """ Function to return the right hand side for the system. :param float freq: Frequency - :param numpy.ndarray (nC,) backSigma: Background conductivity model :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2) :return: RHS for both polarizations, primary fields """ - # Get sources for the frequency - src = self.survey.getSources(freq) - # Make sure that there is 2 polarizations. - # assert len() - # Get the background electric fields - from simpegMT.Sources import homo1DModelSource - eBG_bp = homo1DModelSource(self.mesh,freq,backSigma) - deltM = self.MeDeltaSigma - Abg = -1j*omega(freq)*deltM - return Abg*eBG_bp, eBG_bp + # Get sources for the frequncy(polarizations) + Src = self.survey.getSrcByFreq(freq)[0] + S_e = Src.S_e(self) + return -1j * omega(freq) * S_e - def getRHSderiv(self, freq, backSigma, u, v, adjoint=False): - raise NotImplementedError('getRHSDeriv not implemented yet') - return None + def getRHSderiv(self, freq, u, v, adjoint=False): + """ + The derivative of the RHS with respect to sigma + """ - def fields(self, m, m_back): + Src = self.survey.getSrcByFreq(freq)[0] + S_eDeriv = Src.S_eDeriv(self, v, adjoint) + return -1j * omega(freq) * S_eDeriv + + def fields(self, m): ''' Function to calculate all the fields for the model m. :param np.ndarray (nC,) m: Conductivity model - :param np.ndarray (nC,) m_back: Background conductivity model ''' + # Set the current model self.curModel = m - 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: @@ -117,12 +119,15 @@ class eForm_ps(BaseMTProblem): print 'Starting work for {:.3e}'.format(freq) sys.stdout.flush() A = self.getA(freq) - rhs, e_p = self.getRHS(freq,m_back) + rhs = self.getRHS(freq) Ainv = self.Solver(A, **self.solverOpts) - e_s = Ainv * rhs - e = e_p + e_s + e_s = Ainv * rhs + # Store the fields - Src = self.survey.getSources(freq) + Src = self.survey.getSrcByFreq(freq)[0] + # Calculate total e + + e = Src.ePrimary(self) + e_s # Store the fieldss F[Src, 'e_px'] = e[:,0] F[Src, 'e_py'] = e[:,1] @@ -134,9 +139,9 @@ class eForm_ps(BaseMTProblem): print 'Ran for {:f} seconds'.format(time.time()-startTime) sys.stdout.flush() return F - + class eForm_Tp(BaseMTProblem): - """ + """ A MT problem solving a e formulation and a total/primary fields decompostion. Solves the equation @@ -146,7 +151,7 @@ class eForm_Tp(BaseMTProblem): _eqLocs = 'FE' fieldsPair = FieldsMT - # Set new properties + # Set new properties # Background model @property def backModel(self): @@ -210,7 +215,7 @@ class eForm_Tp(BaseMTProblem): """ # Get sources for the frequency src = self.survey.getSources(freq) - # Make sure that there is 2 polarizations. + # Make sure that there is 2 polarizations. # assert len() # Get the background electric fields from simpegMT.Sources import homo1DModelSource @@ -246,7 +251,7 @@ class eForm_Tp(BaseMTProblem): A = self.getA(freq) rhs, e_p = self.getRHS(freq,m_back) Ainv = self.Solver(A, **self.solverOpts) - e_s = Ainv * rhs + e_s = Ainv * rhs e = e_s # Store the fields Src = self.survey.getSources(freq) @@ -261,4 +266,4 @@ class eForm_Tp(BaseMTProblem): print 'Ran for {:f} seconds'.format(time.time()-startTime) sys.stdout.flush() return F - + diff --git a/simpegMT/Sources/backgroundModelSources.py b/simpegMT/Sources/backgroundModelSources.py index 74d21964..261cdbe5 100644 --- a/simpegMT/Sources/backgroundModelSources.py +++ b/simpegMT/Sources/backgroundModelSources.py @@ -1,23 +1,27 @@ import SimPEG as simpeg, numpy as np -def homo1DModelSource(mesh,freq,m_back): +def homo1DModelSource(mesh,freq,sigma_1d): ''' 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 m_back: Background model of conductivity to base the calculations on. + :param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model. :rtype: numpy.ndarray (mesh.nE,2) :return: eBG_bp, E fields for the background model at both polarizations. ''' - # import from simpegMT.Utils import get1DEfields # Get a 1d solution for a halfspace background - mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]])) - # Note: Everything is using e^iwt - e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq) + if mesh.dim == 1: + mesh1d = mesh + elif mesh.dim == 2: + mesh1d = simpeg.Mesh.TensorMesh([mesh.hy],np.array([mesh.x0[1]])) + elif mesh.dim == 3: + mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]])) + # # Note: Everything is using e^iwt + e0_1d = get1DEfields(mesh1d,sigma_1d,freq) # Setup x (east) polarization (_x) ex_px = np.zeros(mesh.vnEx,dtype=complex) ey_px = np.zeros((mesh.nEy,1),dtype=complex) @@ -32,7 +36,7 @@ def homo1DModelSource(mesh,freq,m_back): 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 np.arange(mesh.vnEy[0]): for j in np.arange(mesh.vnEy[1]): ey_py[i,j,:] = e0_1d diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 96bc7723..dfe174a4 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -1,8 +1,11 @@ -from SimPEG import Survey, Utils, Problem, np, sp, mkvc +from SimPEG import Survey, Utils, Problem, Maps, np, sp, mkvc +from simpegEM.FDEM.SurveyFDEM import SrcFDEM +from simpegEM.Utils.EMUtils 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 ################# ### Receivers ### ################# @@ -45,7 +48,7 @@ class RxMT(Survey.BaseRx): """ 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] @@ -59,7 +62,7 @@ class RxMT(Survey.BaseRx): """ 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] @@ -74,7 +77,7 @@ class RxMT(Survey.BaseRx): """ return self.knownRxTypes[self.rxType][0] - + @property def projComp(self): """Component projection (real/imag)""" @@ -82,12 +85,12 @@ class RxMT(Survey.BaseRx): def projectFields(self, src, mesh, u): ''' - Project the fields and return the + Project the fields and return the ''' if self.projType is 'Z1D': Pex = mesh.getInterpolationMat(self.locs,'Fx') - Pbx = mesh.getInterpolationMat(self.locs,'Ex') + Pbx = mesh.getInterpolationMat(self.locs,'Ex') ex = Pex*mkvc(u[src,'e_1d'],2) bx = Pbx*mkvc(u[src,'b_1d'],2)/mu_0 f_part_complex = ex/bx @@ -144,30 +147,82 @@ class RxMT(Survey.BaseRx): return Pv -# Note: Might need to add tests to make sure that both polarization have the same rxList. +# Note: Might need to add tests to make sure that both polarization have the same rxList. ############### ### Sources ### ############### class srcMT(Survey.BaseSrc): ''' - Sources for the MT problem. + 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 - :param str srcPol: The polarization of the source ''' freq = None #: Frequency (float) - rxPair = RxMT - knownSrcTypes = ['pol_xy','pol_x','pol_y'] # ORThogonal POLarization - - def __init__(self, freq, rxList, srcPol = 'pol_xy'): # remove rxType? hardcode to one thing. always polarizations + def __init__(self, rxList, freq): self.freq = float(freq) - Survey.BaseSrc.__init__(self, None, srcPol, rxList) + 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, sigma1d): + 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 = sigma1d + srcMT.__init__(self, rxList, freq) + + + + def ePrimary(self,problem): + # Get primary fields for both polarizations + eBG_bp = homo1DModelSource(problem.mesh,self.freq,self.sigma1d) + return eBG_bp + + def bPrimary(self,problem): + # Project ePrimary to bPrimary + # Satisfies the primary(background) field conditions + bBG_bp = (- self.mesh.edgeCurl * self.ePrimary )/( 1j*omega(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) + Mesigma = problem.MeSigma + Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p) + return (Mesigma - Mesigma_p) * e_p + + def S_eDeriv(self, problem, v, adjoint = False): + MesigmaDeriv = problem.MeSigmaDeriv(self.ePrimary(problem)) + if adjoint: + return MesigmaDeriv.T * v + else: + return MesigmaDeriv * v ############## @@ -208,7 +263,7 @@ class SurveyMT(Survey.BaseSurvey): return len(self._freqDict) # TODO: Rename to getSources - def getSources(self, freq): + def getSrcByFreq(self, freq): """Returns the sources associated with a specific frequency.""" assert freq in self._freqDict, "The requested frequency is not in this survey." return self._freqDict[freq] diff --git a/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py b/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py index 04d491ee..8aa54880 100644 --- a/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py +++ b/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py @@ -32,7 +32,7 @@ def setupSurvey(sigmaHalf): # Source list srcList =[] for freq in freqs: - srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList)) + srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq)) survey = simpegmt.SurveyMT.SurveyMT(srcList) return survey, sigma, m1d diff --git a/simpegMT/Tests/test_Problem3D_againstAnalytic.py b/simpegMT/Tests/test_Problem3D_againstAnalytic.py index 8b459ad7..57af1427 100644 --- a/simpegMT/Tests/test_Problem3D_againstAnalytic.py +++ b/simpegMT/Tests/test_Problem3D_againstAnalytic.py @@ -71,8 +71,9 @@ def runSimpegMTfwd_eForm_ps(inputsProblem): rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,rxType)) # Source list srcList =[] + sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:] for freq in freqs: - srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList)) + srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigma1d)) # Survey MT survey = simpegmt.SurveyMT.SurveyMT(srcList) @@ -83,7 +84,7 @@ def runSimpegMTfwd_eForm_ps(inputsProblem): problem.Solver = MumpsSolver problem.pair(survey) - fields = problem.fields(sig,sigBG) + fields = problem.fields(sig) mtData = survey.projectFields(fields) return (survey, problem, fields, mtData) @@ -93,7 +94,7 @@ def getAppResPhs(MTdata): # Make impedance def appResPhs(freq,z): app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2 - app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi) + app_phs = np.arctan2(z.imag,z.real)*(180/np.pi) return app_res, app_phs recData = MTdata.toRecArray('Complex') return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx']) @@ -107,7 +108,7 @@ def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True): if appR: return np.linalg.norm(np.abs(app_rpxy[0,:] - np.ones(survey.nFreq)/sigmaHalf) * sigmaHalf) else: - return np.linalg.norm(np.abs(app_rpxy[1,:] - np.ones(survey.nFreq)/135) * 135) + return np.linalg.norm(np.abs(app_rpxy[1,:] + np.ones(survey.nFreq)*135) / 135) class TestAnalytics(unittest.TestCase): diff --git a/simpegMT/Utils/MT1Dsolutions.py b/simpegMT/Utils/MT1Dsolutions.py index 1de508d9..aae83162 100644 --- a/simpegMT/Utils/MT1Dsolutions.py +++ b/simpegMT/Utils/MT1Dsolutions.py @@ -13,7 +13,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0): # Conductivity Msig = m1d.getFaceInnerProduct(sigma) # Set up the solution matrix - A = G.T*Mmu*G - 1j*2.*np.pi*freq*Msig + A = G.T*Mmu*G + 1j*2.*np.pi*freq*Msig # Define the inner part of the solution matrix Aii = A[1:-1,1:-1] # Define the outer part of the solution matrix @@ -27,7 +27,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0): ## Note: The analytic solution is derived with e^iwt bc = np.r_[Etot[0],Etot[-1]] # The right hand side - rhs = -Aio*bc + rhs = Aio*bc # Solve the system Aii_inv = simpeg.Solver(Aii) eii = Aii_inv*rhs