diff --git a/.travis.yml b/.travis.yml index 258c4b7b..efa1f473 100644 --- a/.travis.yml +++ b/.travis.yml @@ -19,6 +19,7 @@ env: - TEST_DIR=tests/em/fdem/inverse/derivs - TEST_DIR=tests/em/tdem - TEST_DIR=tests/flow + - TEST_DIR=tests/mt - TEST_DIR=tests/examples - TEST_DIR=tests/em/fdem/inverse/adjoint - TEST_DIR=tests/em/fdem/forward diff --git a/SimPEG/MT/BaseMT.py b/SimPEG/MT/BaseMT.py new file mode 100644 index 00000000..7ebf66d3 --- /dev/null +++ b/SimPEG/MT/BaseMT.py @@ -0,0 +1,116 @@ +from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np +from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem +from SurveyMT import Survey, Data +from FieldsMT import BaseMTFields + + +class BaseMTProblem(BaseFDEMProblem): + + def __init__(self, mesh, **kwargs): + BaseFDEMProblem.__init__(self, mesh, **kwargs) + Utils.setKwargs(self, **kwargs) + # Set the default pairs of the problem + surveyPair = Survey + dataPair = Data + fieldsPair = BaseMTFields + + # Set the solver + Solver = SimpegSolver + solverOpts = {} + + verbose = False + # Notes: + # Use the forward and devs from BaseFDEMProblem + # Might need to add more stuff here. + + def Jvec(self, m, v, u=None): + """ + Function to calculate the data sensitivities dD/dm times a vector. + + :param numpy.ndarray (nC, 1) - conductive model + :param numpy.ndarray (nC, 1) - random vector + :param MTfields object (optional) - MT fields object, if not given it is calculated + :rtype: MTdata object + :return: Data sensitivities wrt m + """ + + # Calculate the fields + if u is None: + u = self.fields(m) + # Set current model + self.curModel = m + # Initiate the Jv object + Jv = self.dataPair(self.survey) + + # Loop all the frequenies + for freq in self.survey.freqs: + dA_du = self.getA(freq) # + + dA_duI = self.Solver(dA_du, **self.solverOpts) + + for src in self.survey.getSrcByFreq(freq): + # We need fDeriv_m = df/du*du/dm + df/dm + # Construct du/dm, it requires a solve + # NOTE: need to account for the 2 polarizations in the derivatives. + u_src = u[src,:] + # dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations. + dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns. + dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns. + if dRHS_dm is None: + du_dm = dA_duI * ( -dA_dm ) + else: + du_dm = dA_duI * ( -dA_dm + dRHS_dm ) + # Calculate the projection derivatives + for rx in src.rxList: + # Get the projection derivative + # v should be of size 2*nE (for 2 polarizations) + PDeriv_u = lambda t: rx.projectFieldsDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m + Jv[src, rx] = PDeriv_u(mkvc(du_dm)) + # Return the vectorized sensitivities + return mkvc(Jv) + + def Jtvec(self, m, v, u=None): + if u is None: + u = self.fields(m) + + self.curModel = m + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + Jtv = np.zeros(m.size) + + for freq in self.survey.freqs: + AT = self.getA(freq).T + + ATinv = self.Solver(AT, **self.solverOpts) + + for src in self.survey.getSrcByFreq(freq): + ftype = self._fieldType + 'Solution' + u_src = u[src, :] + + for rx in src.rxList: + # Get the adjoint projectFieldsDeriv + # PTv needs to be nE, + PTv = rx.projectFieldsDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m + # Get the + dA_duIT = ATinv * PTv + dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True) + dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True) + # Make du_dmT + if dRHS_dmT is None: + du_dmT = -dA_dmT + else: + du_dmT = -dA_dmT + dRHS_dmT + # Select the correct component + # du_dmT needs to be of size nC, + real_or_imag = rx.projComp + if real_or_imag == 'real': + Jtv += du_dmT.real + elif real_or_imag == 'imag': + Jtv += -du_dmT.real + else: + raise Exception('Must be real or imag') + + return Jtv \ No newline at end of file diff --git a/SimPEG/MT/Examples/simple3DfowardProblem.py b/SimPEG/MT/Examples/simple3DfowardProblem.py new file mode 100644 index 00000000..feabae28 --- /dev/null +++ b/SimPEG/MT/Examples/simple3DfowardProblem.py @@ -0,0 +1,41 @@ +# Test script to use SimPEG.MT platform to forward model synthetic data. + +# Import +import SimPEG as simpeg +from SimPEG import MT +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 +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(MT.RxMT(simpeg.mkvc(loc,2).T,rxType)) +# Source list +srcList =[] +for freq in np.logspace(3,-3,7): + srcList.append(MT.SurveyMT.srcMT(freq,rxList)) +# Survey MT +survey = MT.Survey(srcList) + +## Setup the problem object +problem = MT.ProblemMT.MTProblem(M) +problem.pair(survey) + +fields = problem.fields(sig,sigBG) +mtData = survey.projectFields(fields) + diff --git a/SimPEG/MT/FieldsMT.py b/SimPEG/MT/FieldsMT.py new file mode 100644 index 00000000..87ea109d --- /dev/null +++ b/SimPEG/MT/FieldsMT.py @@ -0,0 +1,348 @@ +from SimPEG import Survey, Utils, Problem, np, sp, mkvc +from scipy.constants import mu_0 +import sys +from numpy.lib import recfunctions as recFunc +from SimPEG.EM.Utils import omega + +############## +### Fields ### +############## +class BaseMTFields(Problem.Fields): + """Field Storage for a MT survey.""" + knownFields = {} + dtype = complex + + +class Fields1D_e(BaseMTFields): + """ + Fields storage for the 1D MT solution. + """ + knownFields = {'e_1dSolution':'F'} + aliasFields = { + 'e_1d' : ['e_1dSolution','F','_e'], + 'e_1dPrimary' : ['e_1dSolution','F','_ePrimary'], + 'e_1dSecondary' : ['e_1dSolution','F','_eSecondary'], + 'b_1d' : ['e_1dSolution','E','_b'], + 'b_1dPrimary' : ['e_1dSolution','E','_bPrimary'], + 'b_1dSecondary' : ['e_1dSolution','E','_bSecondary'] + } + + def __init__(self,mesh,survey,**kwargs): + BaseMTFields.__init__(self,mesh,survey,**kwargs) + + def _ePrimary(self, eSolution, srcList): + ePrimary = np.zeros_like(eSolution) + for i, src in enumerate(srcList): + ep = src.ePrimary(self.survey.prob) + if ep is not None: + ePrimary[:,i] = ep[:,-1] + return ePrimary + + def _eSecondary(self, eSolution, srcList): + return eSolution + + def _e(self, eSolution, srcList): + return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) + + def _eDeriv_u(self, src, v, adjoint = False): + return v + + def _eDeriv_m(self, src, v, adjoint = False): + # assuming primary does not depend on the model + return None + + def _bPrimary(self, eSolution, srcList): + bPrimary = np.zeros([self.survey.mesh.nE,eSolution.shape[1]], dtype = complex) + for i, src in enumerate(srcList): + bp = src.bPrimary(self.survey.prob) + if bp is not None: + bPrimary[:,i] += bp[:,-1] + return bPrimary + + def _bSecondary(self, eSolution, srcList): + C = self.mesh.nodalGrad + b = (C * eSolution) + for i, src in enumerate(srcList): + b[:,i] *= - 1./(1j*omega(src.freq)) + # There is no magnetic source in the MT problem + # S_m, _ = src.eval(self.survey.prob) + # if S_m is not None: + # b[:,i] += 1./(1j*omega(src.freq)) * S_m + return b + + def _b(self, eSolution, srcList): + return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList) + + def _bSecondaryDeriv_u(self, src, v, adjoint = False): + C = self.mesh.nodalGrad + if adjoint: + return - 1./(1j*omega(src.freq)) * (C.T * v) + return - 1./(1j*omega(src.freq)) * (C * v) + + def _bSecondaryDeriv_m(self, src, v, adjoint = False): + # Doesn't depend on m + # _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint) + # S_eDeriv = S_eDeriv(v) + # if S_eDeriv is not None: + # return 1./(1j * omega(src.freq)) * S_eDeriv + return None + + def _bDeriv_u(self, src, v, adjoint=False): + # Primary does not depend on u + return self._bSecondaryDeriv_u(src, v, adjoint) + + def _bDeriv_m(self, src, v, adjoint=False): + # Assuming the primary does not depend on the model + return self._bSecondaryDeriv_m(src, v, adjoint) + + def _fDeriv_u(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt u. + + :param MTsrc src: MT source + :param numpy.ndarray v: random vector of f_sol.size + This function stacks the fields derivatives appropriately + + return a vector of size (nreEle+nrbEle) + """ + + de_du = v #Utils.spdiag(np.ones((self.nF,))) + db_du = self._bDeriv_u(src, v, adjoint) + # Return the stack + # This doesn't work... + return np.vstack((de_du,db_du)) + + def _fDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt m. + + This function stacks the fields derivatives appropriately + """ + return None + +class Fields3D_e(BaseMTFields): + """ + Fields storage for the 3D MT solution. + """ + # Define the known the alias fields + # Assume that the solution of e on the E. + ## NOTE: Need to make this more general, to allow for other solutions formats. + knownFields = {'e_pxSolution':'E','e_pySolution':'E'} + aliasFields = { + 'e_px' : ['e_pxSolution','E','_e_px'], + 'e_pxPrimary' : ['e_pxSolution','E','_e_pxPrimary'], + 'e_pxSecondary' : ['e_pxSolution','E','_e_pxSecondary'], + 'e_py' : ['e_pySolution','E','_e_py'], + 'e_pyPrimary' : ['e_pySolution','E','_e_pyPrimary'], + 'e_pySecondary' : ['e_pySolution','E','_e_pySecondary'], + 'b_px' : ['e_pxSolution','F','_b_px'], + 'b_pxPrimary' : ['e_pxSolution','F','_b_pxPrimary'], + 'b_pxSecondary' : ['e_pxSolution','F','_b_pxSecondary'], + 'b_py' : ['e_pySolution','F','_b_py'], + 'b_pyPrimary' : ['e_pySolution','F','_b_pyPrimary'], + 'b_pySecondary' : ['e_pySolution','F','_b_pySecondary'] + } + + def __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) + for i, src in enumerate(srcList): + ep = src.ePrimary(self.survey.prob) + if ep is not None: + e_pxPrimary[:,i] = ep[:,0] + return e_pxPrimary + + def _e_pyPrimary(self, e_pySolution, srcList): + e_pyPrimary = np.zeros_like(e_pySolution) + for i, src in enumerate(srcList): + ep = src.ePrimary(self.survey.prob) + if ep is not None: + e_pyPrimary[:,i] = ep[:,1] + return e_pyPrimary + + def _e_pxSecondary(self, e_pxSolution, srcList): + return e_pxSolution + + def _e_pySecondary(self, e_pySolution, srcList): + return e_pySolution + + def _e_px(self, e_pxSolution, srcList): + return self._e_pxPrimary(e_pxSolution,srcList) + self._e_pxSecondary(e_pxSolution,srcList) + + def _e_py(self, e_pySolution, srcList): + return self._e_pyPrimary(e_pySolution,srcList) + self._e_pySecondary(e_pySolution,srcList) + + #NOTE: For e_p?Deriv_u, + # v has to be u(2*nE) long for the not adjoint and nE long for adjoint. + # Returns nE long for not adjoint and 2*nE long for adjoint + def _e_pxDeriv_u(self, src, v, adjoint = False): + ''' + Takes the derivative of e_px wrt u + ''' + if adjoint: + # adjoint: returns a 2*nE long vector with zero's for py + return np.vstack((v,np.zeros_like(v))) + # Not adjoint: return only the px part of the vector + return v[:len(v)/2] + + def _e_pyDeriv_u(self, src, v, adjoint = False): + ''' + Takes the derivative of e_py wrt u + ''' + if adjoint: + # adjoint: returns a 2*nE long vector with zero's for px + return np.vstack((np.zeros_like(v),v)) + # Not adjoint: return only the px part of the vector + return v[len(v)/2::] + + def _e_pxDeriv_m(self, src, v, adjoint = False): + # assuming primary does not depend on the model + return None + def _e_pyDeriv_m(self, src, v, adjoint = False): + # assuming primary does not depend on the model + return None + + def _b_pxPrimary(self, e_pxSolution, srcList): + b_pxPrimary = np.zeros([self.survey.mesh.nF,e_pxSolution.shape[1]], dtype = complex) + for i, src in enumerate(srcList): + bp = src.bPrimary(self.survey.prob) + if bp is not None: + b_pxPrimary[:,i] += bp[:,0] + return b_pxPrimary + + def _b_pyPrimary(self, e_pySolution, srcList): + b_pyPrimary = np.zeros([self.survey.mesh.nF,e_pySolution.shape[1]], dtype = complex) + for i, src in enumerate(srcList): + bp = src.bPrimary(self.survey.prob) + if bp is not None: + b_pyPrimary[:,i] += bp[:,1] + return b_pyPrimary + + def _b_pxSecondary(self, e_pxSolution, srcList): + C = self.mesh.edgeCurl + b = (C * e_pxSolution) + for i, src in enumerate(srcList): + b[:,i] *= - 1./(1j*omega(src.freq)) + # There is no magnetic source in the MT problem + # S_m, _ = src.eval(self.survey.prob) + # if S_m is not None: + # b[:,i] += 1./(1j*omega(src.freq)) * S_m + return b + + def _b_pySecondary(self, e_pySolution, srcList): + C = self.mesh.edgeCurl + b = (C * e_pySolution) + for i, src in enumerate(srcList): + b[:,i] *= - 1./(1j*omega(src.freq)) + # There is no magnetic source in the MT problem + # S_m, _ = src.eval(self.survey.prob) + # if S_m is not None: + # b[:,i] += 1./(1j*omega(src.freq)) * S_m + return b + + def _b_px(self, eSolution, srcList): + return self._b_pxPrimary(eSolution, srcList) + self._b_pxSecondary(eSolution, srcList) + + def _b_py(self, eSolution, srcList): + return self._b_pyPrimary(eSolution, srcList) + self._b_pySecondary(eSolution, srcList) + + # NOTE: v needs to be length 2*nE to account for both polarizations + def _b_pxSecondaryDeriv_u(self, src, v, adjoint = False): + # C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,0]]) + C = sp.hstack((self.mesh.edgeCurl,Utils.spzeros(self.mesh.nF,self.mesh.nE))) # This works for adjoint = None + if adjoint: + return - 1./(1j*omega(src.freq)) * (C.T * v) + return - 1./(1j*omega(src.freq)) * (C * v) + + def _b_pySecondaryDeriv_u(self, src, v, adjoint = False): + # C = sp.kron(self.mesh.edgeCurl,[[0,0],[0,1]]) + C = sp.hstack((Utils.spzeros(self.mesh.nF,self.mesh.nE),self.mesh.edgeCurl)) # This works for adjoint = None + if adjoint: + return - 1./(1j*omega(src.freq)) * (C.T * v) + return - 1./(1j*omega(src.freq)) * (C * v) + + def _b_pxSecondaryDeriv_m(self, src, v, adjoint = False): + # Doesn't depend on m + # _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint) + # S_eDeriv = S_eDeriv(v) + # if S_eDeriv is not None: + # return 1./(1j * omega(src.freq)) * S_eDeriv + return None + + def _b_pySecondaryDeriv_m(self, src, v, adjoint = False): + # Doesn't depend on m + # _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint) + # S_eDeriv = S_eDeriv(v) + # if S_eDeriv is not None: + # return 1./(1j * omega(src.freq)) * S_eDeriv + return None + + def _b_pxDeriv_u(self, src, v, adjoint=False): + # Primary does not depend on u + return self._b_pxSecondaryDeriv_u(src, v, adjoint) + + def _b_pyDeriv_u(self, src, v, adjoint=False): + # Primary does not depend on u + return self._b_pySecondaryDeriv_u(src, v, adjoint) + + def _b_pxDeriv_m(self, src, v, adjoint=False): + # Assuming the primary does not depend on the model + return self._b_pxSecondaryDeriv_m(src, v, adjoint) + + def _b_pyDeriv_m(self, src, v, adjoint=False): + # Assuming the primary does not depend on the model + return self._b_pySecondaryDeriv_m(src, v, adjoint) + + def _f_pxDeriv_u(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt u. + + :param MTsrc src: MT source + :param numpy.ndarray v: random vector of f_sol.size + This function stacks the fields derivatives appropriately + + return a vector of size (nreEle+nrbEle) + """ + + de_du = v #Utils.spdiag(np.ones((self.nF,))) + db_du = self._b_pxDeriv_u(src, v, adjoint) + # Return the stack + # This doesn't work... + return np.vstack((de_du,db_du)) + + def _f_pyDeriv_u(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt u. + + :param MTsrc src: MT source + :param numpy.ndarray v: random vector of f_sol.size + This function stacks the fields derivatives appropriately + + return a vector of size (nreEle+nrbEle) + """ + + de_du = v #Utils.spdiag(np.ones((self.nF,))) + db_du = self._b_pyDeriv_u(src, v, adjoint) + # Return the stack + # This doesn't work... + return np.vstack((de_du,db_du)) + + def _f_pxDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt m. + + This function stacks the fields derivatives appropriately + """ + # The fields have no dependance to the model. + return None + + def _f_pyDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt m. + + This function stacks the fields derivatives appropriately + """ + # The fields have no dependance to the model. + return None \ No newline at end of file diff --git a/SimPEG/MT/Problem1D/Probs.py b/SimPEG/MT/Problem1D/Probs.py new file mode 100644 index 00000000..54d7c91b --- /dev/null +++ b/SimPEG/MT/Problem1D/Probs.py @@ -0,0 +1,242 @@ +from SimPEG.EM.Utils import omega +from SimPEG import mkvc +from scipy.constants import mu_0 +from SimPEG.MT.BaseMT import BaseMTProblem +from SimPEG.MT.SurveyMT import Survey, Data +from SimPEG.MT.FieldsMT import Fields1D_e +from SimPEG.MT.Utils.MT1Danalytic import getEHfields +import numpy as np +import multiprocessing, sys, time + + +class eForm_psField(BaseMTProblem): + """ + A MT problem soving a e formulation and primary/secondary fields decomposion. + + Solves the equation + + """ + # From FDEMproblem: Used to project the fields. Currently not used for MTproblem. + _fieldType = 'e_1d' + _eqLocs = 'EF' + _sigmaPrimary = None + + + def __init__(self, mesh, **kwargs): + BaseMTProblem.__init__(self, mesh, **kwargs) + self.fieldsPair = Fields1D_e + # self._sigmaPrimary = sigmaPrimary + + @property + def sigmaPrimary(self): + return self._sigmaPrimary + @sigmaPrimary.setter + def sigmaPrimary(self, val): + # Note: TODO add logic for val, make sure it is the correct size. + self._sigmaPrimary = val + + def getA(self, freq): + """ + Function to get the A matrix. + + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A + """ + + Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0) + 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 + # Either return full or only the inner part of A + return A + + def getADeriv_m(self, freq, u, v, adjoint=False): + """ + The derivative of A wrt sigma + """ + + dsig_dm = self.curModel.sigmaDeriv + MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0) + # + u_src = u['e_1dSolution'] + dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv + if adjoint: + return 1j * omega(freq) * ( dMf_dsig.T * v ) + # Note: output has to be nN/nF, not nC/nE. + # v should be nC + return 1j * omega(freq) * ( dMf_dsig * v ) + + def getRHS(self, freq): + """ + Function to return the right hand side for the system. + :param float freq: Frequency + :rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1) + :return: RHS for 1 polarizations, primary fields + """ + + # 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_m(self, freq, v, adjoint=False): + """ + The derivative of the RHS wrt sigma + """ + + Src = self.survey.getSrcByFreq(freq)[0] + S_eDeriv = Src.S_eDeriv_m(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 + ''' + # Set the current model + self.curModel = m + + F = Fields1D_e(self.mesh, self.survey) + for freq in self.survey.freqs: + if self.verbose: + startTime = time.time() + print 'Starting work for {:.3e}'.format(freq) + sys.stdout.flush() + A = self.getA(freq) + rhs = self.getRHS(freq) + Ainv = self.Solver(A, **self.solverOpts) + e_s = Ainv * rhs + + # Store the fields + Src = self.survey.getSrcByFreq(freq)[0] + # NOTE: only store the e_solution(secondary), all other components calculated in the fields object + F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d + + # Note curl e = -iwb so b = -curl e /iw + # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) ) + # F[Src, 'b_1d'] = b[:,1] + if self.verbose: + print 'Ran for {:f} seconds'.format(time.time()-startTime) + sys.stdout.flush() + return F + +class eForm_TotalField(BaseMTProblem): + """ + A MT problem solving a e formulation and a Total bondary domain decompostion. + + Solves the equation: + + Math: + + + """ + + # From FDEMproblem: Used to project the fields. Currently not used for MTproblem. + _fieldType = 'e' + _eqLocs = 'EF' + + + def __init__(self, mesh, **kwargs): + BaseMTProblem.__init__(self, mesh, **kwargs) + + def getA(self, freq, full=False): + """ + Function to get the A matrix. + + :param float freq: Frequency + :param logic full: Return full A or the inner part + :rtype: scipy.sparse.csr_matrix + :return: A + """ + + Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0) + 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 + # Either return full or only the inner part of A + if full: + return A + else: + return A[1:-1,1:-1] + + def getADeriv(self, freq, u, v, adjoint=False): + sig = self.curTModel + dsig_dm = self.curTModelDeriv + dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, 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): + """ + Function to return the right hand side for the system. + :param float freq: Frequency + :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2) + :return: RHS for both polarizations, primary fields + """ + # Get sources for the frequency + # NOTE: Need to use the source information, doesn't really apply in 1D + 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.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 + ## Note: The analytic solution is derived with e^iwt + eBC = np.r_[Etot[0],Etot[-1]] + # The right hand side + + return -Aio*eBC, eBC + + def getRHSderiv(self, freq, backSigma, u, v, adjoint=False): + raise NotImplementedError('getRHSDeriv not implemented yet') + return None + + 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 + ''' + self.curModel = m + # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields + + F = Fields1D_e(self.mesh, self.survey) + for freq in self.survey.freqs: + if self.verbose: + startTime = time.time() + print 'Starting work for {:.3e}'.format(freq) + sys.stdout.flush() + A = self.getA(freq) + rhs, e_o = self.getRHS(freq) + Ainv = self.Solver(A, **self.solverOpts) + e_i = Ainv * rhs + e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2) + # Store the fields + Src = self.survey.getSrcByFreq(freq) + # NOTE: only store e fields + F[Src, 'e_1dSolution'] = e[:,0] + if self.verbose: + print 'Ran for {:f} seconds'.format(time.time()-startTime) + sys.stdout.flush() + return F diff --git a/SimPEG/MT/Problem1D/__init__.py b/SimPEG/MT/Problem1D/__init__.py new file mode 100644 index 00000000..b0a77250 --- /dev/null +++ b/SimPEG/MT/Problem1D/__init__.py @@ -0,0 +1 @@ +from Probs import eForm_TotalField, eForm_psField \ No newline at end of file diff --git a/SimPEG/MT/Problem2D/Probs.py b/SimPEG/MT/Problem2D/Probs.py new file mode 100644 index 00000000..e69de29b diff --git a/SimPEG/MT/Problem2D/__init__.py b/SimPEG/MT/Problem2D/__init__.py new file mode 100644 index 00000000..fc80254b --- /dev/null +++ b/SimPEG/MT/Problem2D/__init__.py @@ -0,0 +1 @@ +pass \ No newline at end of file diff --git a/SimPEG/MT/Problem3D/Probs.py b/SimPEG/MT/Problem3D/Probs.py new file mode 100644 index 00000000..64866cf1 --- /dev/null +++ b/SimPEG/MT/Problem3D/Probs.py @@ -0,0 +1,264 @@ +from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver +from SimPEG.EM.Utils import omega +from scipy.constants import mu_0 +from SimPEG.MT.BaseMT import BaseMTProblem +from SimPEG.MT.SurveyMT import Survey, Data +from SimPEG.MT.FieldsMT import Fields3D_e +import multiprocessing, sys, time + + + +class eForm_ps(BaseMTProblem): + """ + A MT problem solving a e formulation and a primary/secondary fields decompostion. + + Solves the equation: + + + + """ + + # From FDEMproblem: Used to project the fields. Currently not used for MTproblem. + _fieldType = 'e' + _eqLocs = 'FE' + fieldsPair = Fields3D_e + _sigmaPrimary = None + + def __init__(self, mesh, **kwargs): + BaseMTProblem.__init__(self, mesh, **kwargs) + + @property + def sigmaPrimary(self): + return self._sigmaPrimary + @sigmaPrimary.setter + def sigmaPrimary(self, val): + # Note: TODO add logic for val, make sure it is the correct size. + self._sigmaPrimary = val + + def getA(self, freq): + """ + Function to get the A matrix. + + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A + """ + Mmui = self.MfMui + Msig = self.MeSigma + C = self.mesh.edgeCurl + + return C.T*Mmui*C + 1j*omega(freq)*Msig + + def getADeriv_m(self, freq, u, v, adjoint=False): + + # Nee to account for both the polarizations + # dMe_dsig = (self.MeSigmaDeriv( u['e_pxSolution'] ) + self.MeSigmaDeriv( u['e_pySolution'] )) + # dMe_dsig = (self.MeSigmaDeriv( u['e_pxSolution'] + u['e_pySolution'] )) + + # # dMe_dsig = self.MeSigmaDeriv( u ) + # if adjoint: + # return 1j * omega(freq) * ( dMe_dsig.T * v ) # As in simpegEM + + # return 1j * omega(freq) * ( dMe_dsig * v ) # As in simpegEM + + # This considers both polarizations and returns a nE,2 matrix for each polarization + if adjoint: + dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v + else: + # Need a nE,2 matrix to be returned + dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) )) + return 1j * omega(freq) * dMe_dsigV + # Stacking them + + # if adjoint: + # dMe_dsig = sp.vstack((self.MeSigmaDeriv( u['e_pxSolution'] ), self.MeSigmaDeriv( u['e_pxSolution'] ) )).T + # # self.MeSigmaDeriv(u['e_pySolution'] ).T*v,2) )) + # else: + # dMe_dsig = sp.vstack((self.MeSigmaDeriv( u['e_pxSolution'] ), self.MeSigmaDeriv( u['e_pxSolution'] ) )) + # return 1j * omega(freq) * dMe_dsig*v + + + def getRHS(self, freq): + """ + Function to return the right hand side for the system. + :param float freq: Frequency + :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2) + :return: RHS for both polarizations, primary fields + """ + + # 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_m(self, freq, v, adjoint=False): + """ + The derivative of the RHS with respect to sigma + """ + + Src = self.survey.getSrcByFreq(freq)[0] + S_eDeriv = Src.S_eDeriv_m(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 + ''' + # Set the current model + self.curModel = m + + F = Fields3D_e(self.mesh, self.survey) + for freq in self.survey.freqs: + if self.verbose: + startTime = time.time() + print 'Starting work for {:.3e}'.format(freq) + sys.stdout.flush() + A = self.getA(freq) + rhs = self.getRHS(freq) + Ainv = self.Solver(A, **self.solverOpts) + e_s = Ainv * rhs + + # Store the fields + Src = self.survey.getSrcByFreq(freq)[0] + # Store the fieldss + F[Src, 'e_pxSolution'] = e_s[:,0] + F[Src, 'e_pySolution'] = e_s[:,1] + # Note curl e = -iwb so b = -curl/iw + # b = -( self.mesh.edgeCurl * e )/( 1j*omega(freq) ) + # F[Src, 'b_px'] = b[:,0] + # F[Src, 'b_py'] = b[:,1] + if self.verbose: + 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 + """ + + _fieldType = 'e' + _eqLocs = 'FE' + fieldsPair = Fields3D_e + + # Set new properties + # 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 MeSigmaBack(self): + #TODO: hardcoded to sigma as the model + if getattr(self, '_MeSigmaBack', None) is None: + sigma = self.curModel + sigmaBG = self.backModel + self._MeSigmaBack = self.mesh.getEdgeInnerProduct(sigmaBG) + return self._MeSigmaBack + + def __init__(self, mesh, **kwargs): + BaseMTProblem.__init__(self, mesh, **kwargs) + + def getA(self, freq): + """ + Function to get the A matrix. + + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A + """ + mui = self.MfMui + sig = self.MeSigma + C = self.mesh.edgeCurl + + return C.T*mui*C + 1j*omega(freq)*sig + + def getADeriv(self, freq, u, v, adjoint=False): + sig = self.curTModel + dsig_dm = self.curTModelDeriv + dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, 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): + """ + 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) + :return: one RHS for both polarizations + """ + # 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 SimPEG.MT.Sources import homo1DModelSource + eBG_bp = homo1DModelSource(self.mesh,freq,backSigma) + MeBack = self.MeSigmaBack + # Set up the A system + mui = self.MfMui + C = self.mesh.edgeCurl + Abg = C.T*mui*C + 1j*omega(freq)*MeBack + + return Abg*eBG_bp, eBG_bp + + def getRHSderiv(self, freq, backSigma, u, v, adjoint=False): + raise NotImplementedError('getRHSDeriv not implemented yet') + return None + + def fields(self, m, m_back): + ''' + 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 + ''' + self.curModel = m + self.backModel = m_back + # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields + + F = Fields3D_e(self.mesh, self.survey) + for freq in self.survey.freqs: + if self.verbose: + startTime = time.time() + print 'Starting work for {:.3e}'.format(freq) + sys.stdout.flush() + A = self.getA(freq) + rhs, e_p = self.getRHS(freq,m_back) + Ainv = self.Solver(A, **self.solverOpts) + e_s = Ainv * rhs + e = e_s + # Store the fields + Src = self.survey.getSources(freq) + # Store the fieldss + F[Src, 'e_px'] = e[:,0] + F[Src, 'e_py'] = e[:,1] + # Note curl e = -iwb so b = -curl/iw + b = -( self.mesh.edgeCurl * e )/( 1j*omega(freq) ) + F[Src, 'b_px'] = b[:,0] + F[Src, 'b_py'] = b[:,1] + if self.verbose: + print 'Ran for {:f} seconds'.format(time.time()-startTime) + sys.stdout.flush() + return F + diff --git a/SimPEG/MT/Problem3D/__init__.py b/SimPEG/MT/Problem3D/__init__.py new file mode 100644 index 00000000..58c374ca --- /dev/null +++ b/SimPEG/MT/Problem3D/__init__.py @@ -0,0 +1 @@ +from Probs import eForm_ps, eForm_Tp \ No newline at end of file diff --git a/SimPEG/MT/Sources/__init__.py b/SimPEG/MT/Sources/__init__.py new file mode 100644 index 00000000..7b4a53bc --- /dev/null +++ b/SimPEG/MT/Sources/__init__.py @@ -0,0 +1 @@ +from backgroundModelSources import * \ No newline at end of file diff --git a/SimPEG/MT/Sources/backgroundModelSources.py b/SimPEG/MT/Sources/backgroundModelSources.py new file mode 100644 index 00000000..24b36bfe --- /dev/null +++ b/SimPEG/MT/Sources/backgroundModelSources.py @@ -0,0 +1,178 @@ +import SimPEG as simpeg, numpy as np + +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 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 SimPEG.MT.Utils import get1DEfields + # Get a 1d solution for a halfspace background + 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) + if mesh.dim == 1: + eBG_px = simpeg.mkvc(e0_1d,2) + eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents. + elif mesh.dim == 2: + ex_px = np.zeros(mesh.vnEx,dtype=complex) + ey_px = np.zeros((mesh.nEy,1),dtype=complex) + for i in np.arange(mesh.vnEx[0]): + ex_px[i,:] = -e0_1d + eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px)) + # Setup y (north) polarization (_py) + ex_py = np.zeros((mesh.nEx,1), dtype='complex128') + ey_py = np.zeros(mesh.vnEy, dtype='complex128') + # Assign the source to ey_py + for i in np.arange(mesh.vnEy[0]): + ey_py[i,:] = e0_1d + # ey_py[1:-1,1:-1,1:-1] = 0 + eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + elif mesh.dim == 3: + # 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 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(ex_px,2),ey_px,ez_px)) + # Setup y (north) polarization (_py) + 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 np.arange(mesh.vnEy[0]): + for j in np.arange(mesh.vnEy[1]): + ey_py[i,j,:] = e0_1d + # ey_py[1:-1,1:-1,1:-1] = 0 + eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + + # Return the electric fields + eBG_bp = np.hstack((eBG_px,eBG_py)) + return eBG_bp + +def analytic1DModelSource(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 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 SimPEG.MT.Utils import getEHfields + # Get a 1d solution for a halfspace background + 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 + Eu, Ed, _, _ = getEHfields(mesh1d,sigma_1d,freq,mesh.vectorNz) + # Make the fields into a dictionary of location and the fields + e0_1d = Eu+Ed + E1dFieldDict = dict(zip(mesh.vectorNz,e0_1d)) + if mesh.dim == 1: + eBG_px = simpeg.mkvc(e0_1d,2) + eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents. + elif mesh.dim == 2: + ex_px = np.zeros(mesh.vnEx,dtype=complex) + ey_px = np.zeros((mesh.nEy,1),dtype=complex) + for i in np.arange(mesh.vnEx[0]): + ex_px[i,:] = -e0_1d + eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px)) + # Setup y (north) polarization (_py) + ex_py = np.zeros((mesh.nEx,1), dtype='complex128') + ey_py = np.zeros(mesh.vnEy, dtype='complex128') + # Assign the source to ey_py + for i in np.arange(mesh.vnEy[0]): + ey_py[i,:] = e0_1d + # ey_py[1:-1,1:-1,1:-1] = 0 + eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + elif mesh.dim == 3: + # Setup x (east) polarization (_x) + ex_px = -np.array([E1dFieldDict[i] for i in mesh.gridEx[:,2]]).reshape(-1,1) + ey_px = np.zeros((mesh.nEy,1),dtype=complex) + ez_px = np.zeros((mesh.nEz,1),dtype=complex) + # Construct the full fields + eBG_px = np.vstack((ex_px,ey_px,ez_px)) + # Setup y (north) polarization (_py) + ex_py = np.zeros((mesh.nEx,1), dtype='complex128') + ey_py = np.array([E1dFieldDict[i] for i in mesh.gridEy[:,2]]).reshape(-1,1) + ez_py = np.zeros((mesh.nEz,1), dtype='complex128') + # Construct the full fields + eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + + # Return the electric fields + eBG_bp = np.hstack((eBG_px,eBG_py)) + return eBG_bp + +# def homo3DModelSource(mesh,model,freq): +# ''' +# Function that estimates 1D analytic background fields from a 3D model. + +# :param Simpeg mesh object mesh: Holds information on the discretization +# :param float freq: The frequency to solve at +# :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. + +# ''' + +# if mesh.dim < 3: +# raise IOError('Input mesh has to have 3 dimensions.') + + +# # Get the locations +# a = mesh.gridCC[:,0:2].copy() +# unixy = np.unique(a.view(a.dtype.descr * a.shape[1])).view(float).reshape(-1,2) +# uniz = np.unique(mesh.gridCC[:,2]) +# # # Note: Everything is using e^iwt +# # Need to loop thourgh the xy locations, assess the model and calculate the fields at the phusdo cell centers. +# # Then interpolate the cc fields to the edges. + +# e0_1d = get1DEfields(mesh1d,sigma_1d,freq) + +# elif mesh.dim == 3: +# # 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 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(ex_px,2),ey_px,ez_px)) +# # Setup y (north) polarization (_py) +# 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 np.arange(mesh.vnEy[0]): +# for j in np.arange(mesh.vnEy[1]): +# ey_py[i,j,:] = e0_1d +# # ey_py[1:-1,1:-1,1:-1] = 0 +# eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + +# # Return the electric fields +# eBG_bp = np.hstack((eBG_px,eBG_py)) +# return eBG_bp diff --git a/SimPEG/MT/SrcMT.py b/SimPEG/MT/SrcMT.py new file mode 100644 index 00000000..e00dbd52 --- /dev/null +++ b/SimPEG/MT/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/SimPEG/MT/SurveyMT.py b/SimPEG/MT/SurveyMT.py new file mode 100644 index 00000000..b5ede631 --- /dev/null +++ b/SimPEG/MT/SurveyMT.py @@ -0,0 +1,489 @@ +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 +from numpy.lib import recfunctions as recFunc +from Sources import homo1DModelSource +from Utils import rec2ndarr + +import sys + +################# +### Receivers ### +################# +class Rx(SimPEGsurvey.BaseRx): + + knownRxTypes = { + # 3D impedance + 'zxxr':['Z3D', 'real'], + 'zxyr':['Z3D', 'real'], + 'zyxr':['Z3D', 'real'], + 'zyyr':['Z3D', 'real'], + 'zxxi':['Z3D', 'imag'], + 'zxyi':['Z3D', 'imag'], + 'zyxi':['Z3D', 'imag'], + 'zyyi':['Z3D', 'imag'], + # 2D impedance + # TODO: + # 1D impedance + 'z1dr':['Z1D', 'real'], + 'z1di':['Z1D', 'imag'], + # Tipper + 'tzxr':['T3D','real'], + 'tzxi':['T3D','imag'], + 'tzyr':['T3D','real'], + 'tzyi':['T3D','imag'] + } + # TODO: Have locs as single or double coordinates for both or numerator and denominator separately, respectively. + def __init__(self, locs, rxType): + SimPEGsurvey.BaseRx.__init__(self, locs, rxType) + + @property + def projField(self): + """ + 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 ...) + :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 projType(self): + """ + Receiver type for projection. + + """ + return self.knownRxTypes[self.rxType][0] + + @property + def projComp(self): + """Component projection (real/imag)""" + return self.knownRxTypes[self.rxType][1] + + def projectFields(self, src, mesh, f): + ''' + Project the fields and return the correct data. + ''' + + if self.projType is 'Z1D': + Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx') + Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex') + ex = Pex*mkvc(f[src,'e_1d'],2) + bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0 + # Note: Has a minus sign in front, to comply with quadrant calculations. + # Can be derived from zyx case for the 3D case. + 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(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'] + ey_px = Pey*f[src,'e_px'] + ex_py = Pex*f[src,'e_py'] + ey_py = Pey*f[src,'e_py'] + hx_px = Pbx*f[src,'b_px']/mu_0 + hy_px = Pby*f[src,'b_px']/mu_0 + hx_py = Pbx*f[src,'b_py']/mu_0 + 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) + 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) + elif self.projType is 'T3D': + 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'] + bx_py = Pbx*f[src,'b_py'] + by_py = Pby*f[src,'b_py'] + bz_py = Pbz*f[src,'b_py'] + 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) + + else: + NotImplementedError('Projection of {:s} receiver type is not implemented.'.format(self.rxType)) + # Get the real or imag component + real_or_imag = self.projComp + f_part = getattr(f_part_complex, real_or_imag) + # print f_part + return f_part + + def projectFieldsDeriv(self, src, mesh, f, v, adjoint=False): + """ + The derivative of the projection wrt u + + :param MTsrc src: MT source + :param TensorMesh mesh: Mesh defining the topology of the problem + :param MTfields f: MT fields object of the source + :param numpy.ndarray v: Random vector of size + """ + + real_or_imag = self.projComp + + if not adjoint: + if self.projType is 'Z1D': + Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx') + Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex') + # ex = Pex*mkvc(f[src,'e_1d'],2) + # bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0 + dP_de = -mkvc(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v),2) + dP_db = mkvc( Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pbx*f._bDeriv_u(src,v)/mu_0),2) + PDeriv_complex = np.sum(np.hstack((dP_de,dP_db)),1) + elif self.projType is 'Z2D': + raise NotImplementedError('Has not been implement for 2D impedance tensor') + 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(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'] + ey_px = Pey*f[src,'e_px'] + ex_py = Pex*f[src,'e_py'] + ey_py = Pey*f[src,'e_py'] + hx_px = Pbx*f[src,'b_px']/mu_0 + hy_px = Pby*f[src,'b_px']/mu_0 + hx_py = Pbx*f[src,'b_py']/mu_0 + hy_py = Pby*f[src,'b_py']/mu_0 + # Derivatives as lambda functions + # The size of the diratives should be nD,nU + ex_px_u = lambda vec: Pex*f._e_pxDeriv_u(src,vec) + ey_px_u = lambda vec: Pey*f._e_pxDeriv_u(src,vec) + ex_py_u = lambda vec: Pex*f._e_pyDeriv_u(src,vec) + ey_py_u = lambda vec: Pey*f._e_pyDeriv_u(src,vec) + # NOTE: Think b_p?Deriv_u should return a 2*nF size matrix + hx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec)/mu_0 + hy_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec)/mu_0 + hx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec)/mu_0 + hy_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec)/mu_0 + # Update the input vector + sDiag = lambda t: Utils.sdiag(mkvc(t,2)) + # Define the components of the derivative + Hd = sDiag(1./(sDiag(hx_px)*hy_py - sDiag(hx_py)*hy_px)) + Hd_uV = sDiag(hy_py)*hx_px_u(v) + sDiag(hx_px)*hy_py_u(v) - sDiag(hx_py)*hy_px_u(v) - sDiag(hy_px)*hx_py_u(v) + # Calculate components + if 'zxx' in self.rxType: + Zij = sDiag(Hd*( sDiag(ex_px)*hy_py - sDiag(ex_py)*hy_px )) + ZijN_uV = sDiag(hy_py)*ex_px_u(v) + sDiag(ex_px)*hy_py_u(v) - sDiag(ex_py)*hy_px_u(v) - sDiag(hy_px)*ex_py_u(v) + elif 'zxy' in self.rxType: + Zij = sDiag(Hd*(-sDiag(ex_px)*hx_py + sDiag(ex_py)*hx_px )) + ZijN_uV = -sDiag(hx_py)*ex_px_u(v) - sDiag(ex_px)*hx_py_u(v) + sDiag(ex_py)*hx_px_u(v) + sDiag(hx_px)*ex_py_u(v) + elif 'zyx' in self.rxType: + Zij = sDiag(Hd*( sDiag(ey_px)*hy_py - sDiag(ey_py)*hy_px )) + ZijN_uV = sDiag(hy_py)*ey_px_u(v) + sDiag(ey_px)*hy_py_u(v) - sDiag(ey_py)*hy_px_u(v) - sDiag(hy_px)*ey_py_u(v) + elif 'zyy' in self.rxType: + Zij = sDiag(Hd*(-sDiag(ey_px)*hx_py + sDiag(ey_py)*hx_px )) + ZijN_uV = -sDiag(hx_py)*ey_px_u(v) - sDiag(ey_px)*hx_py_u(v) + sDiag(ey_py)*hx_px_u(v) + sDiag(hx_px)*ey_py_u(v) + + # Calculate the complex derivative + PDeriv_complex = Hd * (ZijN_uV - Zij * Hd_uV ) + # Extract the real number for the real/imag components. + Pv = np.array(getattr(PDeriv_complex, real_or_imag)) + elif adjoint: + # Note: The v vector is real and the return should be complex + if self.projType is 'Z1D': + Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx') + Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex') + # ex = Pex*mkvc(f[src,'e_1d'],2) + # bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0 + dP_deTv = -mkvc(Pex.T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*v,2) + db_duv = Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v + dP_dbTv = mkvc(f._bDeriv_u(src,db_duv,adjoint=True),2) + PDeriv_real = np.sum(np.hstack((dP_deTv,dP_dbTv)),1) + elif self.projType is 'Z2D': + raise NotImplementedError('Has not be implement for 2D impedance tensor') + 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(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. + aex_px = mkvc(mkvc(f[src,'e_px'],2).T*Pex.T) + aey_px = mkvc(mkvc(f[src,'e_px'],2).T*Pey.T) + aex_py = mkvc(mkvc(f[src,'e_py'],2).T*Pex.T) + aey_py = mkvc(mkvc(f[src,'e_py'],2).T*Pey.T) + ahx_px = mkvc(mkvc(f[src,'b_px'],2).T/mu_0*Pbx.T) + ahy_px = mkvc(mkvc(f[src,'b_px'],2).T/mu_0*Pby.T) + ahx_py = mkvc(mkvc(f[src,'b_py'],2).T/mu_0*Pbx.T) + ahy_py = mkvc(mkvc(f[src,'b_py'],2).T/mu_0*Pby.T) + # Derivatives as lambda functions + aex_px_u = lambda vec: f._e_pxDeriv_u(src,Pex.T*vec,adjoint=True) + aey_px_u = lambda vec: f._e_pxDeriv_u(src,Pey.T*vec,adjoint=True) + aex_py_u = lambda vec: f._e_pyDeriv_u(src,Pex.T*vec,adjoint=True) + aey_py_u = lambda vec: f._e_pyDeriv_u(src,Pey.T*vec,adjoint=True) + ahx_px_u = lambda vec: f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0 + ahy_px_u = lambda vec: f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0 + ahx_py_u = lambda vec: f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0 + ahy_py_u = lambda vec: f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0 + + # Update the input vector + # Define shortcuts + sDiag = lambda t: Utils.sdiag(mkvc(t,2)) + sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2)) + # Define the components of the derivative + aHd = sDiag(1./(sDiag(ahx_px)*ahy_py - sDiag(ahx_py)*ahy_px)) + aHd_uV = lambda x: ahx_px_u(sDiag(ahy_py)*x) + ahx_px_u(sDiag(ahy_py)*x) - ahy_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(ahy_px)*x) + # Need to fix this to reflect the adjoint + if 'zxx' in self.rxType: + Zij = sDiag(aHd*( sDiag(ahy_py)*aex_px - sDiag(ahy_px)*aex_py)) + ZijN_uV = lambda x: aex_px_u(sDiag(ahy_py)*x) + ahy_py_u(sDiag(aex_px)*x) - ahy_px_u(sDiag(aex_py)*x) - aex_py_u(sDiag(ahy_px)*x) + elif 'zxy' in self.rxType: + Zij = sDiag(aHd*(-sDiag(ahx_py)*aex_px + sDiag(ahx_px)*aex_py)) + ZijN_uV = lambda x:-aex_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(aex_px)*x) + ahx_px_u(sDiag(aex_py)*x) + aex_py_u(sDiag(ahx_px)*x) + elif 'zyx' in self.rxType: + Zij = sDiag(aHd*( sDiag(ahy_py)*aey_px - sDiag(ahy_px)*aey_py)) + ZijN_uV = lambda x: aey_px_u(sDiag(ahy_py)*x) + ahy_py_u(sDiag(aey_px)*x) - ahy_px_u(sDiag(aey_py)*x) - aey_py_u(sDiag(ahy_px)*x) + elif 'zyy' in self.rxType: + Zij = sDiag(aHd*(-sDiag(ahx_py)*aey_px + sDiag(ahx_px)*aey_py)) + ZijN_uV = lambda x:-aey_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(aey_px)*x) + ahx_px_u(sDiag(aey_py)*x) + aey_py_u(sDiag(ahx_px)*x) + + # Calculate the complex derivative + PDeriv_real = ZijN_uV(aHd*v) - aHd_uV(Zij.T*aHd*v)# + # NOTE: Need to reshape the output to go from 2*nU array to a (nU,2) matrix for each polarization + # PDeriv_real = np.hstack((mkvc(PDeriv_real[:len(PDeriv_real)/2],2),mkvc(PDeriv_real[len(PDeriv_real)/2::],2))) + PDeriv_real = PDeriv_real.reshape((2,mesh.nE)).T + # Extract the data + if real_or_imag == 'imag': + Pv = 1j*PDeriv_real + elif real_or_imag == 'real': + Pv = PDeriv_real.astype(complex) + + + return Pv + +################# +### 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 + + """ + import SrcMT + srcPair = SrcMT.BaseMTSrc + + def __init__(self, srcList, **kwargs): + # Sort these by frequency + self.srcList = srcList + SimPEGsurvey.BaseSurvey.__init__(self, **kwargs) + + _freqDict = {} + 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]) + + @property + def freqs(self): + """Frequencies""" + return self._freqs + + @property + def nFreq(self): + """Number of frequencies""" + return len(self._freqDict) + + # TODO: Rename to getSources + 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] + + def projectFields(self, u): + data = Data(self) + for src in self.srcList: + sys.stdout.flush() + for rx in src.rxList: + data[src, rx] = rx.projectFields(src, self.mesh, u) + return data + + 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/SimPEG/MT/Utils/MT1Danalytic.py b/SimPEG/MT/Utils/MT1Danalytic.py new file mode 100644 index 00000000..cf0847d9 --- /dev/null +++ b/SimPEG/MT/Utils/MT1Danalytic.py @@ -0,0 +1,111 @@ +# Analytic solution of EM fields due to a plane wave + +import numpy as np, SimPEG as simpeg + +def getEHfields(m1d,sigma,freq,zd,scaleUD=True): + '''Analytic solution for MT 1D layered earth. Returns E and H fields. + + :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information. + :param numpy.array, vector sigma: Physical property of conductivity corresponding with the mesh. + :param float, freq: Frequency to calculate data at. + :param numpy array, vector zd: location to calculate EH fields at + :param bollean, scaleUD: scales the output to be 1 at the top, increases numeracal stability. + + Assumes a halfspace with the same conductive as the last cell below. + + ''' + # Note add an error check for the mesh and sigma are the same size. + + # Constants: Assume constant + mu = 4*np.pi*1e-7*np.ones((m1d.nC+1)) + eps = 8.85*1e-12*np.ones((m1d.nC+1)) + # Angular freq + w = 2*np.pi*freq + # Add the halfspace value to the property + sig = np.concatenate((np.array([sigma[0]]),sigma)) + # Calculate the wave number + k = np.sqrt(eps*mu*w**2-1j*mu*sig*w) + + # Initiate the propagation matrix, in the order down up. + UDp = np.zeros((2,m1d.nC+1),dtype=complex) + UDp[1,0] = 1. # Set the wave amplitude as 1 into the half-space at the bottom of the mesh + # Loop over all the layers, starting at the bottom layer + for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer + # Calculate + yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer + zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer + # Build the propagation matrix + + # Convert fields to down/up going components in layer below current layer + Pj1 = np.array([[1,1],[yp1,-yp1]]) + # Convert fields to down/up going components in current layer + Pjinv = 1./2*np.array([[1,zp],[1,-zp]]) + # Propagate down and up components through the current layer + elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]]) + + # The down and up component in current layer. + UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr]) + + if scaleUD: + UDp[:,lnr+1::-1] = UDp[:,lnr+1::-1]/UDp[1,lnr+1] + + # Calculate the fields + Ed = np.empty((zd.size,),dtype=complex) + Eu = np.empty((zd.size,),dtype=complex) + Hd = np.empty((zd.size,),dtype=complex) + Hu = np.empty((zd.size,),dtype=complex) + + # Loop over the layers and calculate the fields + # In the halfspace below the mesh + dup = m1d.vectorNx[0] + dind = dup >= zd + Ed[dind] = UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind])) + Eu[dind] = UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind])) + Hd[dind] = (k[0]/(w*mu[0]))*UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind])) + Hu[dind] = -(k[0]/(w*mu[0]))*UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind])) + for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]): + dind = np.logical_and(dup >= zd, zd > dlow) + Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind])) + Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind])) + Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind])) + Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind])) + + # Return return the fields + return Ed, Eu, Hd, Hu + +def getImpedance(m1d,sigma,freq): + """Analytic solution for MT 1D layered earth. Returns the impedance at the surface. + + :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information. + :param numpy.array, vector sigma: Physical property corresponding with the mesh. + :param numpy.array, vector freq: Frequencies to calculate data at. + + + """ + + # Define constants + mu0 = 4*np.pi*1e-7 + eps0 = 8.85e-12 + + # Initiate the impedances + Z1d = np.empty(len(freq) , dtype='complex') + h = m1d.hx #vectorNx[:-1] + # Start the process + for nrFr, fr in enumerate(freq): + om = 2*np.pi*fr + Zall = np.empty(len(h)+1,dtype='complex') + # Calculate the impedance for the bottom layer + Zall[0] = (mu0*om)/np.sqrt(mu0*eps0*(om)**2 - 1j*mu0*sigma[0]*om) + + for nr,hi in enumerate(h): + # Calculate the wave number + # print nr,sigma[nr] + k = np.sqrt(mu0*eps0*om**2 - 1j*mu0*sigma[nr]*om) + Z = (mu0*om)/k + + Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi))) + + #pdb.set_trace() + Z1d[nrFr] = Zall[-1] + + return Z1d diff --git a/SimPEG/MT/Utils/MT1Dsolutions.py b/SimPEG/MT/Utils/MT1Dsolutions.py new file mode 100644 index 00000000..aae83162 --- /dev/null +++ b/SimPEG/MT/Utils/MT1Dsolutions.py @@ -0,0 +1,45 @@ +import numpy as np, SimPEG as simpeg +from MT1Danalytic import getEHfields +from scipy.constants import mu_0 + +def get1DEfields(m1d,sigma,freq,sourceAmp=1.0): + """Function to get 1D electrical fields""" + + # Get the gradient + G = m1d.nodalGrad + # Mass matrices + # Magnetic permeability + Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0)) + # Conductivity + Msig = m1d.getFaceInnerProduct(sigma) + # Set up the solution matrix + 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 + Aio = A[1:-1,[0,-1]] + + # Set the boundary conditions + Ed, Eu, Hd, Hu = getEHfields(m1d,sigma,freq,m1d.vectorNx) + Etot = (Ed + Eu) + if sourceAmp is not None: + Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top + ## Note: The analytic solution is derived with e^iwt + bc = np.r_[Etot[0],Etot[-1]] + # The right hand side + rhs = Aio*bc + # Solve the system + Aii_inv = simpeg.Solver(Aii) + eii = Aii_inv*rhs + # Assign the boundary conditions + e = np.r_[bc[0],eii,bc[1]] + # Return the electrical fields + return e + + +if __name__ == '__main__': + + hz = [(100.,18)] + M = simpeg.Mesh.TensorMesh([hz],'C') + sig = np.zeros(M.nC) + 1e-8 + sig[M.vectorCCx<=0] = sigHalf diff --git a/SimPEG/MT/Utils/__init__.py b/SimPEG/MT/Utils/__init__.py new file mode 100644 index 00000000..b683f8b4 --- /dev/null +++ b/SimPEG/MT/Utils/__init__.py @@ -0,0 +1,4 @@ +from MT1Dsolutions import * # Add the names of the functions +from MT1Danalytic import * +from dataUtils import * +from ediFilesUtils import * diff --git a/SimPEG/MT/Utils/dataUtils.py b/SimPEG/MT/Utils/dataUtils.py new file mode 100644 index 00000000..b1947cb8 --- /dev/null +++ b/SimPEG/MT/Utils/dataUtils.py @@ -0,0 +1,245 @@ +# Utils used for the data, +import numpy as np, matplotlib.pyplot as plt, sys +import SimPEG as simpeg +import numpy.lib.recfunctions as recFunc +from scipy.constants import mu_0 +from scipy import interpolate as sciint + +def getAppRes(MTdata): + # Make impedance + zList = [] + for src in MTdata.survey.srcList: + zc = [src.freq] + for rx in src.rxList: + if 'i' in rx.rxType: + m=1j + else: + m = 1 + zc.append(m*MTdata[src,rx]) + zList.append(zc) + return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))] + +def rotateData(MTdata,rotAngle): + ''' + Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation) + ''' + recData = MTdata.toRecArray('Complex') + impData = rec2ndarr(recData[['zxx','zxy','zyx','zyy']],complex) + # Make the rotation matrix + # c,s,zxx,zxy,zyx,zyy = sympy.symbols('c,s,zxx,zxy,zyx,zyy') + # rotM = sympy.Matrix([[c,-s],[s, c]]) + # zM = sympy.Matrix([[zxx,zxy],[zyx,zyy]]) + # rotM*zM*rotM.T + # [c*(c*zxx - s*zyx) - s*(c*zxy - s*zyy), c*(c*zxy - s*zyy) + s*(c*zxx - s*zyx)], + # [c*(c*zyx + s*zxx) - s*(c*zyy + s*zxy), c*(c*zyy + s*zxy) + s*(c*zyx + s*zxx)]]) + s = np.sin(-np.deg2rad(rotAngle)) + c = np.cos(-np.deg2rad(rotAngle)) + rotMat = np.array([[c,-s],[s,c]]) + rotData = (rotMat.dot(impData.reshape(-1,2,2).dot(rotMat.T))).transpose(1,0,2).reshape(-1,4) + outRec = recData.copy() + for nr,comp in enumerate(['zxx','zxy','zyx','zyy']): + outRec[comp] = rotData[:,nr] + + from SimPEG import MT + return MT.Data.fromRecArray(outRec) + + +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) + return app_res, app_phs + +def skindepth(rho,freq): + ''' Function to calculate the skindepth of EM waves''' + return np.sqrt( (rho*((1/(freq * mu_0 * np.pi ))))) + +def rec2ndarr(x,dt=float): + return x.view((dt, len(x.dtype.names))) + +def makeAnalyticSolution(mesh,model,elev,freqs): + from SimPEG import MT + data1D = [] + for freq in freqs: + anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev) + anaE = anaEd+anaEu + anaH = anaHd+anaHu + + anaZ = anaE/anaH + # Add to the list + data1D.append((freq,0,0,elev,anaZ[0])) + dataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zyx',complex)]) + return dataRec + +def plotMT1DModelData(problem,models,symList=None): + from SimPEG import MT + # Setup the figure + fontSize = 15 + + fig = plt.figure(figsize=[9,7]) + axM = fig.add_axes([0.075,.1,.25,.875]) + axM.set_xlabel('Resistivity [Ohm*m]',fontsize=fontSize) + axM.set_xlim(1e-1,1e5) + axM.set_ylim(-10000,5000) + axM.set_ylabel('Depth [km]',fontsize=fontSize) + axR = fig.add_axes([0.42,.575,.5,.4]) + axR.set_xscale('log') + axR.set_yscale('log') + axR.invert_xaxis() + # axR.set_xlabel('Frequency [Hz]') + axR.set_ylabel('Apparent resistivity [Ohm m]',fontsize=fontSize) + + axP = fig.add_axes([0.42,.1,.5,.4]) + axP.set_xscale('log') + axP.invert_xaxis() + axP.set_ylim(0,90) + axP.set_xlabel('Frequency [Hz]',fontsize=fontSize) + axP.set_ylabel('Apparent phase [deg]',fontsize=fontSize) + + # if not symList: + # symList = ['x']*len(models) + import plotDataTypes as pDt + # Loop through the models. + modelList = [problem.survey.mtrue] + modelList.extend(models) + if False: + modelList = [problem.mapping.sigmaMap*mod for mod in modelList] + for nr, model in enumerate(modelList): + # Calculate the data + if nr==0: + data1D = problem.dataPair(problem.survey,problem.survey.dobs).toRecArray('Complex') + else: + data1D = problem.dataPair(problem.survey,problem.survey.dpred(model)).toRecArray('Complex') + # Plot the data and the model + colRat = nr/((len(modelList)-1.999)*1.) + if colRat > 1.: + col = 'k' + else: + col = plt.cm.seismic(1-colRat) + # The model - make the pts to plot + meshPts = np.concatenate((problem.mesh.gridN[0:1],np.kron(problem.mesh.gridN[1::],np.ones(2))[:-1])) + modelPts = np.kron(1./(problem.mapping.sigmaMap*model),np.ones(2,)) + axM.semilogx(modelPts,meshPts,color=col) + + ## Data + # Appres + pDt.plotIsoStaImpedance(axR,np.array([0,0]),data1D,'zyx','res',pColor=col) + # Appphs + pDt.plotIsoStaImpedance(axP,np.array([0,0]),data1D,'zyx','phs',pColor=col) + try: + allData = np.concatenate((allData,simpeg.mkvc(data1D['zyx'],2)),1) + except: + allData = simpeg.mkvc(data1D['zyx'],2) + freq = simpeg.mkvc(data1D['freq'],2) + res, phs = appResPhs(freq,allData) + + stdCol = 'gray' + axRtw = axR.twinx() + axRtw.set_ylabel('Std of log10',color=stdCol) + [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()] + axPtw = axP.twinx() + axPtw.set_ylabel('Std ',color=stdCol) + [t.set_color(stdCol) for t in axPtw.get_yticklabels()] + axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol) + axPtw.plot(freq, np.std(phs,1),'--',color=stdCol) + + # Fix labels and ticks + + yMtick = [l/1000 for l in axM.get_yticks().tolist()] + axM.set_yticklabels(yMtick) + [ l.set_rotation(90) for l in axM.get_yticklabels()] + [ l.set_rotation(90) for l in axR.get_yticklabels()] + [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()] + [t.set_color(stdCol) for t in axPtw.get_yticklabels()] + for ax in [axM,axR,axP]: + ax.xaxis.set_tick_params(labelsize=fontSize) + ax.yaxis.set_tick_params(labelsize=fontSize) + return fig + +def printTime(): + import time + print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime()) + +def convert3Dto1Dobject(MTdata,rxType3D='zyx'): + from SimPEG import MT + # Find the unique locations + # Need to find the locations + recDataTemp = MTdata.toRecArray() + # Check if survey.std has been assigned. + ## NEED TO: write this... + # Calculte and add the DET of the tensor to the recArray + if 'det' in rxType3D: + Zon = (recDataTemp['zxxr']+1j*recDataTemp['zxxi'])*(recDataTemp['zyyr']+1j*recDataTemp['zyyi']) + Zoff = (recDataTemp['zxyr']+1j*recDataTemp['zxyi'])*(recDataTemp['zyxr']+1j*recDataTemp['zyxi']) + det = np.sqrt(Zon.data - Zoff.data) + recData = recFunc.append_fields(recDataTemp,['zdetr','zdeti'],[det.real,det.imag] ) + else: + recData = recDataTemp + + uniLocs = rec2ndarr(np.unique(recData[['x','y','z']])).data + mtData1DList = [] + if 'zxy' in rxType3D: + corr = -1 # Shift the data to comply with the quadtrature of the 1d problem + else: + corr = 1 + for loc in uniLocs: + # Make the receiver list + rx1DList = [] + for rxType in ['z1dr','z1di']: + rx1DList.append(MT.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(MT.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 = MT.Survey(src1DList) + + # Make the data + dataVec = np.hstack(dat1DList) + dat1D = MT.Data(sur1D,dataVec) + sur1D.dobs = dataVec + # Need to take MTdata.survey.std and split it as well. + std=0.05 + sur1D.std = np.abs(sur1D.dobs*std) #+ 0.01*np.linalg.norm(sur1D.dobs) + mtData1DList.append(dat1D) + + # Return the the list of data. + return mtData1DList + +def resampleMTdataAtFreq(MTdata,freqs): + """ + Function to resample MTdata at set of frequencies + + """ + from SimPEG import MT + # Make a rec array + MTrec = MTdata.toRecArray().data + + # Find unique locations + uniLoc = np.unique(MTrec[['x','y','z']]) + uniFreq = MTdata.survey.freqs + # Get the comps + dNames = MTrec.dtype + + # Loop over all the locations and interpolate + for loc in uniLoc: + # Find the index of the station + ind = np.sqrt(np.sum((rec2ndarr(MTrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy + # Make a temporary recArray and interpolate all the components + tArrRec = np.concatenate((simpeg.mkvc(freqs,2),np.ones((len(freqs),1))*rec2ndarr(loc),np.nan*np.ones((len(freqs),12))),axis=1).view(dNames) + for comp in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']: + int1d = sciint.interp1d(MTrec[ind]['freq'],MTrec[ind][comp],bounds_error=False) + tArrRec[comp] = simpeg.mkvc(int1d(freqs),2) + + # Join together + try: + outRecArr = recFunc.stack_arrays((outRecArr,tArrRec)) + except NameError as e: + outRecArr = tArrRec + + # Make the MTdata and return + return MT.Data.fromRecArray(outRecArr) diff --git a/SimPEG/MT/Utils/ediFilesUtils.py b/SimPEG/MT/Utils/ediFilesUtils.py new file mode 100644 index 00000000..55818a87 --- /dev/null +++ b/SimPEG/MT/Utils/ediFilesUtils.py @@ -0,0 +1,175 @@ +# Functions to import and export MT EDI files. +from SimPEG import mkvc +from scipy.constants import mu_0 +from numpy.lib import recfunctions as recFunc +from SimPEG.MT.Utils.dataUtils import rec2ndarr + +# Import modules +import numpy as np +import os, sys, re +try: + import osr +except ImportError as e: + print 'Could not import osr, missing the gdal package' + pass + +class EDIimporter: + """ + A class to import EDIfiles. + + """ + _impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit + _impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI) + + # Properties + filesList = None + comps = None + + # Hidden properties + _outEPSG = None + _2out = None + + + def __init__(self, EDIfilesList, compList=None, outEPSG=None): + + # Set the fileList + self.filesList = EDIfilesList + # Set the components to import + if compList is None: + self.comps = ['ZXXR','ZXYR','ZYXR','ZYYR','ZXXI','ZXYI','ZYXI','ZYYI','ZXX.VAR','ZXY.VAR','ZYX.VAR','ZYY.VAR'] + else: + self.comps = compList + if outEPSG is not None: + self._outEPSG = outEPSG + + def __call__(self,comps=None): + + if comps is None: + return self._data + + return self._data[comps] + + def importFiles(self): + """ + Function to import EDI files into a object. + + + """ + + # Constants that are needed for convertion of units + + # Temp lists + tmpStaList = [] + + tmpCompList = ['freq','x','y','z'] + tmpCompList.extend(self.comps) + # Make the outarray + dtRI = [(compS.lower().replace('.',''),float) for compS in tmpCompList] + # Loop through all the files + for nrEDI, EDIfile in enumerate(self.filesList): + # Read the file into a list of the lines + with open(EDIfile,'r') as fid: + EDIlines = fid.readlines() + # Find the location + latD, longD, elevM = _findLatLong(EDIlines) + # Transfrom coordinates + transCoord = self._transfromPoints(longD,latD) + # Extract the name of the file (station) + EDIname = EDIfile.split(os.sep)[-1].split('.')[0] + # Arrange the data + staList = [EDIname, EDIfile, transCoord[0], transCoord[1], elevM[0]] + # Add to the station list + tmpStaList.extend(staList) + + # Read the frequency data + freq = _findEDIcomp('>FREQ',EDIlines) + # Make the temporary rec array. + tArrRec = ( np.nan*np.ones( (len(freq),len(dtRI)) ) ).view(dtRI) #np.concatenate((freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],8))),axis=1).view(dtRI) + # Add data to the array + tArrRec['freq'] = mkvc(freq,2) + tArrRec['x'] = mkvc(np.ones((len(freq),1))*transCoord[0],2) + tArrRec['y'] = mkvc(np.ones((len(freq),1))*transCoord[1],2) + tArrRec['z'] = mkvc(np.ones((len(freq),1))*elevM[0],2) + for comp in self.comps: + # Deal with converting units of the impedance tensor + if 'Z' in comp: + unitConvert = self._impUnitEDI2SI + else: + unitConvert = 1 + # Rotate the data since EDI x is *north, y *east but Simpeg uses x *east, y *north (* means internal reference frame) + key = [comp.lower().replace('.','').replace(s,t) for s,t in [['xx','yy'],['xy','yx'],['yx','xy'],['yy','xx']] if s in comp.lower()][0] + tArrRec[key] = mkvc(unitConvert*_findEDIcomp('>'+comp,EDIlines),2) + # Make a masked array + mArrRec = np.ma.MaskedArray(rec2ndarr(tArrRec),mask=np.isnan(rec2ndarr(tArrRec))).view(dtype=tArrRec.dtype) + try: + outTemp = recFunc.stack_arrays((outTemp,mArrRec)) + except NameError as e: + outTemp = mArrRec + + # Assign the data + self._data = outTemp + + # % Assign the data to the obj + # nOutData=length(obj.data); + # obj.data(nOutData+1:nOutData+length(TEMP.data),:) = TEMP.data; + def _transfromPoints(self,longD,latD): + # Coordinates convertor + if self._2out is None: + src = osr.SpatialReference() + src.ImportFromEPSG(4326) + out = osr.SpatialReference() + if self._outEPSG is None: + # Find the UTM EPSG number + Nnr = 700 if latD < 0.0 else 600 + utmZ = int(1+(longD+180.0)/6.0) + self._outEPSG = 32000 + Nnr + utmZ + out.ImportFromEPSG(self._outEPSG) + self._2out = osr.CoordinateTransformation(src,out) + # Return the transfrom + return self._2out.TransformPoint(longD,latD) + +# Hidden functions +def _findLatLong(fileLines): + latDMS = np.array(fileLines[_findLine('LAT=',fileLines)[0]].split('=')[1].split()[0].split(':'),float) + longDMS = np.array(fileLines[_findLine('LONG=',fileLines)[0]].split('=')[1].split()[0].split(':'),float) + elevM = np.array([fileLines[_findLine('ELEV=',fileLines)[0]].split('=')[1].split()[0]],float) + # Convert to D.ddddd values + latS = np.sign(latDMS[0]) + longS = np.sign(longDMS[0]) + latD = latDMS[0] + latS*latDMS[1]/60 + latS*latDMS[2]/3600 + longD = longDMS[0] + longS*longDMS[1]/60 + longS*longDMS[2]/3600 + return latD, longD, elevM + +def _findLine(comp,fileLines): + """ Find a line number in the file""" + # Line counter + c = 0 + # List of indices for found lines + found = [] + # Loop through all the lines + for line in fileLines: + if comp in line: + # Append if found + found.append(c) + # Increse the counter + c += 1 + # Return the found indices + return found + +def _findEDIcomp(comp,fileLines,dt=float): + """ + Extract the data vector. + + Returns a list of the data. + """ + # Find the data + headLine, indHead = [(st,nr) for nr,st in enumerate(fileLines) if re.search(comp,st)][0] + # Extract the data + nrVec = int(headLine.split()[-1]) + c = 0 + dataList = [] + while c < nrVec: + indHead += 1 + dataList.extend(fileLines[indHead].split()) + c = len(dataList) + return np.array(dataList,dt) diff --git a/SimPEG/MT/Utils/plotDataTypes.py b/SimPEG/MT/Utils/plotDataTypes.py new file mode 100644 index 00000000..c0180f00 --- /dev/null +++ b/SimPEG/MT/Utils/plotDataTypes.py @@ -0,0 +1,416 @@ +from matplotlib import pyplot as plt, colors, numpy as np + + +def rec2nd(structArray): + """ Converts a structured/record array to ndarray to do operations on.""" + return structArray.view((np.float,len(structArray.dtype.names))) + +def plotIsoFreqNSimpedance(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=True,contour=True): + + indUniFreq = np.where(freq==array['freq']) + + + x, y = array['x'][indUniFreq],array['y'][indUniFreq] + if par == 'abs': + zPlot = np.abs(array[flag][indUniFreq]) + cmap = plt.get_cmap('OrRd_r')#seismic') + level = np.logspace(0,-5,31) + clevel = np.logspace(0,-4,5) + plotNorm = colors.LogNorm() + elif par == 'real': + zPlot = np.real(array[flag][indUniFreq]) + cmap = plt.get_cmap('RdYlBu') + if cLevel: + level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31))) + clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5))) + else: + level = np.linspace(zPlot.min(),zPlot.max(),100) + clevel = np.linspace(zPlot.min(),zPlot.max(),10) + if colorNorm=='SymLog': + plotNorm = colors.SymLogNorm(1e-10,linscale=2) + else: + plotNorm = colors.Normalize() + elif par == 'imag': + zPlot = np.imag(array[flag][indUniFreq]) + cmap = plt.get_cmap('RdYlBu') + level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31))) + clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5))) + plotNorm = colors.SymLogNorm(1e-10,linscale=2) + if cLevel: + level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31))) + clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5))) + else: + level = np.linspace(zPlot.min(),zPlot.max(),100) + clevel = np.linspace(zPlot.min(),zPlot.max(),10) + if colorNorm=='SymLog': + plotNorm = colors.SymLogNorm(1e-10,linscale=2) + elif colorNorm=='Lin': + plotNorm = colors.Normalize() + if contour: + cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both') + else: + uniX,uniY = np.unique(x),np.unique(y) + X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25)) + cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm) + if colorbar: + plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e') + ax.set_title(flag+' '+par,fontsize=8) + return cs + +def plotIsoFreqNSDiff(ax,freq,arrayList,flag,par='abs',colorbar=True,cLevel=True,mask=None,contourLine=True,useLog=False): + + indUniFreq0 = np.where(freq==arrayList[0]['freq']) + indUniFreq1 = np.where(freq==arrayList[1]['freq']) + seicmap = plt.get_cmap('RdYlBu')#seismic') + x, y = arrayList[0]['x'][indUniFreq0],arrayList[0]['y'][indUniFreq0] + if par == 'abs': + if useLog: + zPlot = (np.log10(np.abs(arrayList[0][flag][indUniFreq0])) - np.log10(np.abs(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs(arrayList[1][flag][indUniFreq1])) + else: + zPlot = (np.abs(arrayList[0][flag][indUniFreq0]) - np.abs(arrayList[1][flag][indUniFreq1]))/np.abs(arrayList[1][flag][indUniFreq1]) + if mask: + maskInd = np.logical_or(np.abs(arrayList[0][flag][indUniFreq0])< 1e-3,np.abs(arrayList[1][flag][indUniFreq1]) < 1e-3) + zPlot = np.ma.array(zPlot) + zPlot[maskInd] = mask + if cLevel: + level = np.arange(-200,201,10) + clevel = np.arange(-200,201,25) + else: + level = np.linspace(zPlot.min(),zPlot.max(),100) + clevel = np.linspace(zPlot.min(),zPlot.max(),10) + elif par == 'real': + if useLog: + zPlot = (np.log10(np.real(arrayList[0][flag][indUniFreq0])) -np.log10(np.real(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.real(arrayList[1][flag][indUniFreq1])))) + else: + zPlot = (np.real(arrayList[0][flag][indUniFreq0]) -np.real(arrayList[1][flag][indUniFreq1]))/np.abs((np.real(arrayList[1][flag][indUniFreq1]))) + if mask: + maskInd = np.logical_or(np.abs(np.real(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.real(arrayList[1][flag][indUniFreq1])) < 1e-3) + zPlot = np.ma.array(zPlot) + zPlot[maskInd] = mask + if cLevel: + level = np.arange(-200,201,10) + clevel = np.arange(-200,201,25) + else: + level = np.linspace(zPlot.min(),zPlot.max(),100) + clevel = np.linspace(zPlot.min(),zPlot.max(),10) + elif par == 'imag': + if useLog: + zPlot = (np.log10(np.imag(arrayList[0][flag][indUniFreq0])) -np.log10(np.imag(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.imag(arrayList[1][flag][indUniFreq1])))) + else: + zPlot = (np.imag(arrayList[0][flag][indUniFreq0]) -np.imag(arrayList[1][flag][indUniFreq1]))/np.abs((np.imag(arrayList[1][flag][indUniFreq1]))) + if mask: + maskInd = np.logical_or(np.abs(np.imag(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.imag(arrayList[1][flag][indUniFreq1])) < 1e-3) + zPlot = np.ma.array(zPlot) + zPlot[maskInd] = mask + if cLevel: + level = np.arange(-200,201,10) + clevel = np.arange(-200,201,25) + else: + level = np.linspace(zPlot.min(),zPlot.max(),100) + clevel = np.linspace(zPlot.min(),zPlot.max(),10) + cs = ax.tricontourf(x,y,zPlot*100,levels=level*100,cmap=seicmap,extend='both') #,norm=colors.SymLogNorm(1e-2,linscale=2)) + if contourLine: + csl = ax.tricontour(x,y,zPlot*100,levels=clevel*100,colors='k') + plt.clabel(csl, fontsize=7, inline=1,fmt='%1.1e',inline_spacing=10) + if colorbar: + cb = plt.colorbar(cs,cax=ax.cax,ticks=clevel*100,format='%1.1e') + for t in cb.ax.get_yticklabels(): + t.set_rotation(60) + t.set_fontsize(8) + + ax.set_title(flag+' '+par,fontsize=8) + +def plotIsoFreqNStipper(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=True,contour=True): + + indUniFreq = np.where(freq==array['freq']) + + x, y = array['x'][indUniFreq],array['y'][indUniFreq] + if par == 'abs': + cmap = plt.get_cmap('OrRd_r')#seismic') + zPlot = np.abs(array[flag][indUniFreq]) + if cLevel: + level = np.logspace(-4,0,33) + clevel = np.logspace(-4,0,5) + else: + level = np.linspace(zPlot.min(),zPlot.max(),100) + clevel = np.linspace(zPlot.min(),zPlot.max(),10) + if colorNorm=='SymLog': + plotNorm = colors.LogNorm() + else: + plotNorm = colors.Normalize() + elif par == 'real': + cmap = plt.get_cmap('RdYlBu') + zPlot = np.real(array[flag][indUniFreq]) + if cLevel: + level = np.concatenate((-np.logspace(0,-4,33),np.logspace(-4,0,33))) + clevel = np.concatenate((-np.logspace(0,-4,5),np.logspace(-4,0,5))) + else: + level = np.linspace(zPlot.min(),zPlot.max(),100) + clevel = np.linspace(zPlot.min(),zPlot.max(),10) + if colorNorm=='SymLog': + plotNorm = colors.SymLogNorm(1e-4,linscale=2) + else: + plotNorm = colors.Normalize() + elif par == 'imag': + cmap = plt.get_cmap('RdYlBu') + zPlot = np.imag(array[flag][indUniFreq]) + if cLevel: + level = np.concatenate((-np.logspace(0,-4,33),np.logspace(-4,0,33))) + clevel = np.concatenate((-np.logspace(0,-4,5),np.logspace(-4,0,5))) + else: + level = np.linspace(zPlot.min(),zPlot.max(),100) + clevel = np.linspace(zPlot.min(),zPlot.max(),10) + if colorNorm=='SymLog': + plotNorm = colors.SymLogNorm(1e-4,linscale=2) + else: + plotNorm = colors.Normalize() + if contour: + cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both') + else: + uniX,uniY = np.unique(x),np.unique(y) + X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25)) + cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),levels=level,cmap=cmap,norm=plotNorm,edgecolors='k', linewidths=0.5) + if colorbar: + plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e') + ax.set_title(flag+' '+par,fontsize=8) + +def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None): + + appResFact = 1/(8*np.pi**2*10**(-7)) + treshold = 1.0 # 1 meter + indUniSta = np.sqrt(np.sum((rec2nd(array[['x','y']])-loc)**2,axis=1)) < treshold + freq = array['freq'][indUniSta] + + if par == 'abs': + zPlot = np.abs(array[flag][indUniSta]) + elif par == 'real': + zPlot = np.real(array[flag][indUniSta]) + elif par == 'imag': + zPlot = np.imag(array[flag][indUniSta]) + elif par == 'res': + zPlot = (appResFact/freq)*np.abs(array[flag][indUniSta])**2 + elif par == 'phs': + zPlot = np.arctan2(array[flag][indUniSta].imag,array[flag][indUniSta].real)*(180/np.pi) + + if not pColor: + if 'xx' in flag: + lab = 'XX' + pColor = 'g' + elif 'xy' in flag: + lab = 'XY' + pColor = 'r' + elif 'yx' in flag: + lab = 'YX' + pColor = 'b' + elif 'yy' in flag: + lab = 'YY' + pColor = 'y' + + ax.plot(freq,zPlot,color=pColor,marker=pSym,label=flag) + + +def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colorNorm='None',cLevel=None,contour=True): + + indSect = np.where(sectDict.values()[0]==array[sectDict.keys()[0]]) + + # Define the plot axes + if 'x' in sectDict.keys()[0]: + x = array['y'][indSect] + else: + x = array['x'][indSect] + y = array['freq'][indSect] + + if par == 'abs': + zPlot = np.abs(array[flag][indSect]) + cmap = plt.get_cmap('OrRd_r')#seismic') + if cLevel: + level = np.logspace(0,-5,31,endpoint=True) + clevel = np.logspace(0,-4,5,endpoint=True) + else: + level = np.linspace(zPlot.min(),zPlot.max(),100,endpoint=True) + clevel = np.linspace(zPlot.min(),zPlot.max(),10,endpoint=True) + + elif par == 'ares': + zPlot = np.abs(array[flag][indSect])**2/(8*np.pi**2*10**(-7)*array['freq'][indSect]) + cmap = plt.get_cmap('RdYlBu')#seismic) + if cLevel: + zMax = np.log10(cLevel[1]) + zMin = np.log10(cLevel[0]) + else: + zMax = (np.ceil(np.log10(np.abs(zPlot).max()))) + zMin = (np.floor(np.log10(np.abs(zPlot).min()))) + level = np.logspace(zMin,zMax,(zMax-zMin)*8+1,endpoint=True) + clevel = np.logspace(zMin,zMax,(zMax-zMin)*2+1,endpoint=True) + plotNorm = colors.LogNorm() + + elif par == 'aphs': + zPlot = np.arctan2(array[flag][indSect].imag,array[flag][indSect].real)*(180/np.pi) + cmap = plt.get_cmap('RdYlBu')#seismic) + if cLevel: + zMax = cLevel[1] + zMin = cLevel[0] + else: + zMax = (np.ceil(zPlot).max()) + zMin = (np.floor(zPlot).min()) + level = np.arange(zMin,zMax+.1,1) + clevel = np.arange(zMin,zMax+.1,10) + plotNorm = colors.Normalize() + + elif par == 'real': + zPlot = np.real(array[flag][indSect]) + cmap = plt.get_cmap('Spectral') #('RdYlBu') + if cLevel: + zMax = np.log10(cLevel[1]) + zMin = np.log10(cLevel[0]) + else: + zMax = (np.ceil(np.log10(np.abs(zPlot).max()))) + zMin = (np.floor(np.log10(np.abs(zPlot).min()))) + level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True))) + clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True))) + plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1) + elif par == 'imag': + zPlot = np.imag(array[flag][indSect]) + cmap = plt.get_cmap('Spectral') #('RdYlBu') + + if cLevel: + zMax = np.log10(cLevel[1]) + zMin = np.log10(cLevel[0]) + else: + zMax = (np.ceil(np.log10(np.abs(zPlot).max()))) + zMin = (np.floor(np.log10(np.abs(zPlot).min()))) + level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True))) + clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True))) + plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1) + + if colorNorm=='SymLog': + plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1) + elif colorNorm=='Lin': + plotNorm = colors.Normalize() + elif colorNorm=='Log': + plotNorm = colors.LogNorm() + if contour: + cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both') + else: + uniX,uniY = np.unique(x),np.unique(y) + X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25)) + cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm) + if colorbar: + csB = plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e') + # csB.on_mappable_changed(cs) + ax.set_title(flag+' '+par,fontsize=8) + return cs, csB + return cs,None + + +def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=None,contour=True,mask=None,useLog=False): + + def sortInArr(arr): + return np.sort(arr,order=['freq','x','y','z']) + # Find the index for the slice + indSect0 = np.where(sectDict.values()[0]==arrayList[0][sectDict.keys()[0]]) + indSect1 = np.where(sectDict.values()[0]==arrayList[1][sectDict.keys()[0]]) + # Extract and sort the mats + arr0 = sortInArr(arrayList[0][indSect0]) + arr1 = sortInArr(arrayList[1][indSect1]) + + # Define the plot axes + if 'x' in sectDict.keys()[0]: + x0 = arr0['y'] + x1 = arr1['y'] + else: + x0 = arr0['x'] + x1 = arr1['x'] + y0 = arr0['freq'] + y1 = arr1['freq'] + + + if par == 'abs': + if useLog: + zPlot = (np.log10(np.abs(arr0[flag])) - np.log10(np.abs(arr1[flag])))/np.log10(np.abs(arr1[flag])) + else: + zPlot = (np.abs(arr0[flag]) - np.abs(arr1[flag]))/np.abs(arr1[flag]) + if mask: + maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3) + zPlot = np.ma.array(zPlot) + zPlot[maskInd] = mask + cmap = plt.get_cmap('RdYlBu')#seismic) + elif par == 'ares': + arF = 1/(8*np.pi**2*10**(-7)) + if useLog: + zPlot = (np.log10((arF/arr0['freq'])*np.abs(arr0[flag])**2) - np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2))/np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2) + else: + zPlot = ((arF/arr0['freq'])*np.abs(arr0[flag])**2 - (arF/arr1['freq'])*np.abs(arr1[flag])**2)/((arF/arr1['freq'])*np.abs(arr1[flag])**2) + if mask: + maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3) + zPlot = np.ma.array(zPlot) + zPlot[maskInd] = mask + cmap = plt.get_cmap('Spectral')#seismic) + + elif par == 'aphs': + if useLog: + zPlot = (np.log10(np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi)) - np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi)) )/np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi)) + else: + zPlot = ( np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi) - np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi) )/(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi)) + if mask: + maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3) + zPlot = np.ma.array(zPlot) + zPlot[maskInd] = mask + cmap = plt.get_cmap('Spectral')#seismic) + elif par == 'real': + if useLog: + zPlot = (np.log10(arr0[flag].real) - np.log10(arr1[flag].real))/np.log10(arr1[flag].real) + else: + zPlot = (arr0[flag].real - arr1[flag].real)/arr1[flag].real + if mask: + maskInd = np.logical_or(arr0[flag].real< 1e-3,arr1[flag].real < 1e-3) + zPlot = np.ma.array(zPlot) + zPlot[maskInd] = mask + cmap = plt.get_cmap('Spectral') #('Spectral') + + elif par == 'imag': + if useLog: + zPlot = (np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag))/np.log10(arr1[flag].imag) + else: + zPlot = (arr0[flag].imag - arr1[flag].imag)/arr1[flag].imag + if mask: + maskInd = np.logical_or(arr0[flag].imag< 1e-3,arr1[flag].imag < 1e-3) + zPlot = np.ma.array(zPlot) + zPlot[maskInd] = mask + cmap = plt.get_cmap('Spectral') #('RdYlBu') + + if cLevel: + zMax = np.log10(cLevel[1]) + zMin = np.log10(cLevel[0]) + else: + zMax = (np.ceil(np.log10(np.abs(zPlot).max()))) + zMin = (np.floor(np.log10(np.abs(zPlot).min()))) + + + if colorNorm=='SymLog': + level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True))) + clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True))) + plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1) + elif colorNorm=='Lin': + if cLevel: + level = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/50.) + clevel = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/10.) + else: + level = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/50.) + clevel = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/10.) + plotNorm = colors.Normalize() + elif colorNorm=='Log': + level = np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True) + clevel = np.logspace(zMin,zMax,(zMax-zMin)*2+1,endpoint=True) + plotNorm = colors.LogNorm() + if contour: + cs = ax.tricontourf(x0,y0,zPlot*100,levels=level*100,cmap=cmap,norm=plotNorm,extend='both')#,extend='both') + else: + uniX,uniY = np.unique(x0),np.unique(y0) + X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25)) + cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm) + if colorbar: + csB = plt.colorbar(cs,cax=ax.cax,ticks=clevel*100,format='%1.2e') + # csB.on_mappable_changed(cs) + ax.set_title(flag+' '+par + ' diff',fontsize=8) + return cs, csB + return cs,None diff --git a/SimPEG/MT/Utils/srcUtils.py b/SimPEG/MT/Utils/srcUtils.py new file mode 100644 index 00000000..a98cc5b2 --- /dev/null +++ b/SimPEG/MT/Utils/srcUtils.py @@ -0,0 +1,46 @@ +import SimPEG as simpeg, numpy as np + +def homo1DModelSource(mesh,freq,m_back): + ''' + Function that calculates and return background fields for a 3D mesh and model. + The calculuations use 1D field solution for a vertical slice throught model (south-western most column), + which is assigned at the fields everywhere for the respective polarizations.2 + + :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. + :rtype: numpy.ndarray (mesh.nE,2) + :return: eBG_bp, E fields for the background model at both polarizations. + + ''' + + # import + from SimPEG.MT.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) + # 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 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(ex_px,2),ey_px,ez_px)) + # Setup y (north) polarization (_py) + 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 np.arange(mesh.vnEy[0]): + for j in np.arange(mesh.vnEy[1]): + ey_py[i,j,:] = e0_1d + # ey_py[1:-1,1:-1,1:-1] = 0 + eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + + # Return the electric fields + eBG_bp = np.hstack((eBG_px,eBG_py)) + return eBG_bp diff --git a/SimPEG/MT/__init__.py b/SimPEG/MT/__init__.py new file mode 100644 index 00000000..dfbd9a15 --- /dev/null +++ b/SimPEG/MT/__init__.py @@ -0,0 +1,6 @@ +import Utils +import Sources +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 diff --git a/docs/mt/index.rst b/docs/mt/index.rst new file mode 100644 index 00000000..8e3a2248 --- /dev/null +++ b/docs/mt/index.rst @@ -0,0 +1,11 @@ + + +SimPEG for Magnetotellurics +=========================== + +SimPEG (Simulation and Parameter Estimation in Geophysics) is a python +package for simulation and gradient based parameter estimation in the +context of geoscience applications. + +simpegMT uses SimPEG as the framework for the forward and inverse +magnetotellurics geophysical problems. diff --git a/tests/mt/__init__.py b/tests/mt/__init__.py new file mode 100644 index 00000000..420388ef --- /dev/null +++ b/tests/mt/__init__.py @@ -0,0 +1,12 @@ +import os +import glob +import unittest + +if __name__ == '__main__': + test_file_strings = glob.glob('test_*.py') + module_strings = [str[0:len(str)-3] for str in test_file_strings] + suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str + in module_strings] + testSuite = unittest.TestSuite(suites) + + unittest.TextTestRunner(verbosity=2).run(testSuite) diff --git a/tests/mt/test_ApparentResistivityAnalytic.py b/tests/mt/test_ApparentResistivityAnalytic.py new file mode 100644 index 00000000..2a3b1ba9 --- /dev/null +++ b/tests/mt/test_ApparentResistivityAnalytic.py @@ -0,0 +1,48 @@ +import unittest +from SimPEG import * +from SimPEG import MT + +TOL = 1e-6 + +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) + return app_res, app_phs + +def appResNorm(sigmaHalf): + nFreq = 26 + + m1d = Mesh.TensorMesh([[(100,5,1.5),(100.,10),(100,5,1.5)]], x0=['C']) + sigma = np.zeros(m1d.nC) + sigmaHalf + sigma[m1d.gridCC[:]>200] = 1e-8 + + # Calculate the analytic fields + freqs = np.logspace(4,-4,nFreq) + Z = [] + for freq in freqs: + Ed, Eu, Hd, Hu = MT.Utils.getEHfields(m1d,sigma,freq,np.array([200])) + Z.append((Ed + Eu)/(Hd + Hu)) + + Zarr = np.concatenate(Z) + + app_r, app_p = appResPhs(freqs,Zarr) + + return np.linalg.norm(np.abs(app_r - np.ones(nFreq)/sigmaHalf)) / np.log10(sigmaHalf) + + +class TestAnalytics(unittest.TestCase): + + def setUp(self): + pass + def test_appRes2en1(self):self.assertLess(appResNorm(2e-1), TOL) + def test_appRes2en2(self):self.assertLess(appResNorm(2e-2), TOL) + def test_appRes2en3(self):self.assertLess(appResNorm(2e-3), TOL) + def test_appRes2en4(self):self.assertLess(appResNorm(2e-4), TOL) + def test_appRes2en5(self):self.assertLess(appResNorm(2e-5), TOL) + def test_appRes2en6(self):self.assertLess(appResNorm(2e-6), TOL) + + + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mt/test_Problem1D_againstAnalyticHalfspace.py b/tests/mt/test_Problem1D_againstAnalyticHalfspace.py new file mode 100644 index 00000000..66f5f41d --- /dev/null +++ b/tests/mt/test_Problem1D_againstAnalyticHalfspace.py @@ -0,0 +1,162 @@ +import unittest +import SimPEG as simpeg +from SimPEG import MT +from SimPEG.Utils import meshTensor +import numpy as np +# Define the tolerances +TOLr = 5e-2 +TOLp = 5e-2 + + +def setupSurvey(sigmaHalf,tD=True): + + # Frequency + nFreq = 33 + freqs = np.logspace(3,-3,nFreq) + # Make the mesh + ct = 5 + air = meshTensor([(ct,25,1.3)]) + # coreT0 = meshTensor([(ct,15,1.2)]) + # coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,))) + core = np.concatenate( ( np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) ) + bot = meshTensor([(core[0],10,-1.3)]) + x0 = -np.array([np.sum(np.concatenate((core,bot)))]) + m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0) + # Make the model + sigma = np.zeros(m1d.nC) + sigmaHalf + sigma[m1d.gridCC > 0 ] = 1e-8 + + rxList = [] + for rxType in ['z1dr','z1di']: + rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType)) + # Source list + srcList =[] + if tD: + for freq in freqs: + srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq)) + else: + for freq in freqs: + srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq)) + + survey = MT.Survey(srcList) + return survey, sigma, m1d + +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) + return app_res, app_phs + zList = [] + for src in MTdata.survey.srcList: + zc = [src.freq] + for rx in src.rxList: + if 'i' in rx.rxType: + m=1j + else: + m = 1 + zc.append(m*MTdata[src,rx]) + zList.append(zc) + return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))] + +def appRes_TotalFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = setupSurvey(sigmaHalf) + problem = MT.Problem1D.eForm_TotalField(mesh) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.projectFields(fields) + + # Calculate the app res and phs + app_r = np.array(getAppResPhs(data))[:,0] + + return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf) + +def appPhs_TotalFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = setupSurvey(sigmaHalf) + problem = MT.Problem1D.eForm_TotalField(mesh) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.projectFields(fields) + + # Calculate the app phs + app_p = np.array(getAppResPhs(data))[:,1] + + return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45) + +def appRes_psFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = setupSurvey(sigmaHalf,False) + problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.projectFields(fields) + + # Calculate the app res and phs + app_r = np.array(getAppResPhs(data))[:,0] + + return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf) + +def appPhs_psFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = setupSurvey(sigmaHalf,False) + problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.projectFields(fields) + + # Calculate the app phs + app_p = np.array(getAppResPhs(data))[:,1] + + return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45) + +class TestAnalytics(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.assertLess(appRes_TotalFieldNorm(2e-2), TOLr) + # def test_appPhs2en2(self):self.assertLess(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) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py b/tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py new file mode 100644 index 00000000..2a2e82fd --- /dev/null +++ b/tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py @@ -0,0 +1,135 @@ +import unittest +import SimPEG as simpeg +from SimPEG import MT +from SimPEG.Utils import meshTensor +import numpy as np +# Define the tolerances +TOLr = 5e-2 +TOLp = 5e-2 + + +def setupSurvey(sigmaHalf,tD=True): + + # Frequency + nFreq = 33 + freqs = np.logspace(3,-3,nFreq) + # Make the mesh + ct = 5 + air = meshTensor([(ct,25,1.3)]) + # coreT0 = meshTensor([(ct,15,1.2)]) + # coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,))) + core = np.concatenate( ( np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) ) + bot = meshTensor([(core[0],15,-1.3)]) + x0 = -np.array([np.sum(np.concatenate((core,bot)))]) + m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0) + # Make the model + sigma = np.zeros(m1d.nC) + sigmaHalf + sigma[m1d.gridCC > 0 ] = 1e-8 + sigmaBack = sigma.copy() + # Add structure + shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600) + deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000) + sigma[shallow] = 1 + sigma[deep] = 0.1 + + rxList = [] + for rxType in ['z1dr','z1di']: + rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType)) + # Source list + srcList =[] + if tD: + for freq in freqs: + srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq)) + else: + for freq in freqs: + srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq)) + + survey = MT.Survey(srcList) + return survey, sigma, m1d + +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) + return app_res, app_phs + zList = [] + for src in MTdata.survey.srcList: + zc = [src.freq] + for rx in src.rxList: + if 'i' in rx.rxType: + m=1j + else: + m = 1 + zc.append(m*MTdata[src,rx]) + zList.append(zc) + return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))] + +def calculateAnalyticSolution(srcList,mesh,model): + surveyAna = MT.Survey(srcList) + data1D = MT.Data(surveyAna) + for src in surveyAna.srcList: + elev = src.rxList[0].locs[0] + anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev) + anaE = anaEd+anaEu + anaH = anaHd+anaHu + # Scale the solution + # anaE = (anaEtemp/anaEtemp[-1])#.conj() + # anaH = (anaHtemp/anaEtemp[-1])#.conj() + anaZ = anaE/anaH + for rx in src.rxList: + data1D[src,rx] = getattr(anaZ, rx.projComp) + return data1D + +def dataMis_AnalyticTotalDomain(sigmaHalf): + + # Make the survey + + # Total domain solution + surveyTD, sigma, mesh = setupSurvey(sigmaHalf) + problemTD = MT.Problem1D.eForm_TotalField(mesh) + problemTD.pair(surveyTD) + # Analytic data + dataAnaObj = calculateAnalyticSolution(surveyTD.srcList,mesh,sigma) + # dataTDObj = MT.DataMT.DataMT(surveyTD, surveyTD.dpred(sigma)) + dataTD = surveyTD.dpred(sigma) + dataAna = simpeg.mkvc(dataAnaObj) + return np.all((dataTD - dataAna)/dataAna < 2.) + # surveyTD.dtrue = -simpeg.mkvc(dataAna,2) + # surveyTD.dobs = -simpeg.mkvc(dataAna,2) + # surveyTD.Wd = np.ones(surveyTD.dtrue.shape) #/(np.abs(surveyTD.dtrue)*0.01) + # # Setup the data misfit + # dmis = simpeg.DataMisfit.l2_DataMisfit(surveyTD) + # dmis.Wd = surveyTD.Wd + # return dmis.eval(sigma) + + +def dataMis_AnalyticPrimarySecondary(sigmaHalf): + + # Make the survey + # Primary secondary + surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,tD=False) + problemPS = MT.Problem1D.eForm_psField(mesh) + problemPS.sigmaPrimary = sigmaPS + problemPS.pair(surveyPS) + # Analytic data + dataAnaObj = calculateAnalyticSolution(surveyPS.srcList,mesh,sigmaPS) + + dataPS = surveyPS.dpred(sigmaPS) + dataAna = simpeg.mkvc(dataAnaObj) + return np.all((dataPS - dataAna)/dataAna < 2.) + + + +class TestNumericVsAnalytics(unittest.TestCase): + + def setUp(self): + pass + # Total Fields + # def test_appRes2en2(self):self.assertTrue(dataMis_AnalyticTotalDomain(2e-2)) + + # Primary/secondary + def test_appRes2en2_ps(self):self.assertTrue(dataMis_AnalyticPrimarySecondary(2e-2)) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mt/test_Problem3D_againstAnalytic.py b/tests/mt/test_Problem3D_againstAnalytic.py new file mode 100644 index 00000000..5a07d39a --- /dev/null +++ b/tests/mt/test_Problem3D_againstAnalytic.py @@ -0,0 +1,290 @@ +# Test functions +from glob import glob +import numpy as np, sys, os, time, scipy, subprocess +import SimPEG as simpeg +import unittest +from SimPEG import MT +from SimPEG.Utils import meshTensor +from scipy.constants import mu_0 + +TOLr = 5e-2 +TOL = 1e-4 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = [1e-1, 2e-1] +addrandoms = True + + +def getInputs(): + """ + Function that returns Mesh, freqs, rx_loc, elev. + """ + # 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]) + # M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model + M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(500,8,-1.3),(500.,8),(500,8,1.3)]], x0=['C','C','C'])# Setup the model + # Set the frequencies + freqs = np.logspace(1,-3,5) + elev = 0 + + ## Setup the the survey object + # Receiver locations + rx_x, rx_y = np.meshgrid(np.arange(-3000,3001,500),np.arange(-1000,1001,500)) + rx_loc = np.array([[0, 0, 0]]) #,[-100,-100,0],[100,100,0]]) #np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1)))) + + return M, freqs, rx_loc, elev + + +def random(conds): + ''' Returns a halfspace model based on the inputs''' + M, freqs, rx_loc, elev = getInputs() + + # Backround + sigBG = np.ones(M.nC)*conds + # Add randomness to the model (10% of the value). + sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 ) + + return (M, freqs, sig, sigBG, rx_loc) + +def halfSpace(conds): + ''' Returns a halfspace model based on the inputs''' + M, freqs, rx_loc, elev = getInputs() + + # Model + ccM = M.gridCC + # conds = [1e-2] + groundInd = ccM[:,2] < elev + sig = np.zeros(M.nC) + 1e-8 + sig[groundInd] = conds + # Set the background, not the same as the model + sigBG = np.zeros(M.nC) + 1e-8 + sigBG[groundInd] = conds + + return (M, freqs, sig, sigBG, rx_loc) + +def twoLayer(conds): + ''' Returns a 2 layer model based on the conductivity values given''' + M, freqs, rx_loc, elev = getInputs() + + # Model + ccM = M.gridCC + groundInd = ccM[:,2] < elev + botInd = ccM[:,2] < -3000 + sig = np.zeros(M.nC) + 1e-8 + sig[groundInd] = conds[1] + sig[botInd] = conds[0] + # Set the background, not the same as the model + sigBG = np.zeros(M.nC) + 1e-8 + sigBG[groundInd] = conds[1] + + + return (M, freqs, sig, sigBG, rx_loc) + +def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False,expMap=True): + M,freqs,sig,sigBG,rx_loc = inputSetup + # Make a receiver list + rxList = [] + if comp == 'All': + for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi',]: + rxList.append(MT.Rx(rx_loc,rxType)) + else: + rxList.append(MT.Rx(rx_loc,comp)) + # Source list + srcList =[] + + if singleFreq: + srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,singleFreq)) + else: + for freq in freqs: + srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq)) + # Survey MT + survey = MT.Survey(srcList) + + ## Setup the problem object + sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:] + if expMap: + problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) ) + problem.mapping = simpeg.Maps.ExpMap(problem.mesh) + problem.curModel = np.log(sig) + else: + problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= sigma1d) + problem.curModel = sig + problem.pair(survey) + problem.verbose = False + try: + from pymatsolver import MumpsSolver + problem.Solver = MumpsSolver + except: + pass + + return (survey, problem) + +def setupSimpegMTfwd_eForm_ps_multiRx(inputSetup,comp='All',singleFreq=False,expMap=True): + M,freqs,sig,sigBG,rx_loc = inputSetup + # Add to the receiver list + rx_loc = np.vstack((rx_loc,np.array([[-100,100,0]])))# ,[-100,-100,0],[100,-100,0],[100,100,0]]))) + # Make a receiver list + rxList = [] + if comp == 'All': + for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi',]: + rxList.append(MT.Rx(rx_loc,rxType)) + else: + rxList.append(MT.Rx(rx_loc,comp)) + # Source list + srcList =[] + + if singleFreq: + srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,singleFreq)) + else: + for freq in freqs: + srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq)) + # Survey MT + survey = MT.SurveyMT.SurveyMT(srcList) + + ## Setup the problem object + sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:] + if expMap: + problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) ) + problem.mapping = simpeg.Maps.ExpMap(problem.mesh) + problem.curModel = np.log(sig) + else: + problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= sigma1d) + problem.curModel = sig + problem.pair(survey) + problem.verbose = False + try: + from pymatsolver import MumpsSolver + problem.Solver = MumpsSolver + except: + pass + + return (survey, problem) +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) + return app_res, app_phs + recData = MTdata.toRecArray('Complex') + return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx']) + +def JvecAdjointTest(inputSetup,comp='All',freq=False): + (M, freqs, sig, sigBG, rx_loc) = inputSetup + survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=freq) + print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs)) + + m = sig + u = problem.fields(m) + + v = np.random.rand(survey.nD,) + # print problem.PropMap.PropModel.nP + w = np.random.rand(problem.mesh.nC,) + + vJw = v.ravel().dot(problem.Jvec(m, w, u)) + wJtv = w.ravel().dot(problem.Jtvec(m, v, u)) + tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) + print ' vJw wJtv vJw - wJtv tol abs(vJw - wJtv) < tol' + print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol + return np.abs(vJw - wJtv) < tol + +# Test the Jvec derivative +def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True): + (M, freqs, sig, sigBG, rx_loc) = inputSetup + survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp=comp,singleFreq=freq,expMap=expMap) + print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs) + # problem.mapping = simpeg.Maps.ExpMap(problem.mesh) + # problem.sigmaPrimary = np.log(sigBG) + x0 = np.log(sigBG) + # cond = sig[0] + # x0 = np.log(np.ones(problem.mesh.nC)*cond) + # problem.sigmaPrimary = x0 + # if True: + # x0 = x0 + np.random.randn(problem.mesh.nC)*cond*1e-1 + survey = problem.survey + def fun(x): + return survey.dpred(x), lambda x: problem.Jvec(x0, x) + return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR) + + +def DerivProjfieldsTest(inputSetup,comp='All',freq=False): + + survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp,freq) + print 'Derivative test of data projection for eFormulation primary/secondary\n\n' + # problem.mapping = simpeg.Maps.ExpMap(problem.mesh) + # Initate things for the derivs Test + src = survey.srcList[0] + rx = src.rxList[0] + + u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j + u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j + u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2))) + f0 = problem.fieldsPair(survey.mesh,survey) + # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2))) + f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x + f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y + + def fun(u): + f = problem.fieldsPair(survey.mesh,survey) + f[src,'e_pxSolution'] = u[:len(u)/2] + f[src,'e_pySolution'] = u[len(u)/2::] + return rx.projectFields(src,survey.mesh,f), lambda t: rx.projectFieldsDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2)) + + return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR) + + +def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False): + if appR: + label = 'resistivity' + else: + label = 'phase' + # Make the survey and the problem + survey, problem = setupSimpegMTfwd_eForm_ps(halfSpace(sigmaHalf),expMap=expMap) + print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf) + + data = problem.dataPair(survey,survey.dpred(problem.curModel)) + # Calculate the app phs + app_rpxy, app_rpyx = np.array(getAppResPhs(data)) + if appR: + return np.all(np.abs(app_rpxy[0,:] - np.ones(survey.nFreq)/sigmaHalf) * sigmaHalf < .4) + else: + return np.all(np.abs(app_rpxy[1,:] + np.ones(survey.nFreq)*135) / 135 < .4) + +class TestAnalytics(unittest.TestCase): + + def setUp(self): + pass + # # Test apparent resistivity and phase + def test_appRes1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2)) + def test_appPhs1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False)) + + def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3)) + def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False)) + + # Do a derivative test + def test_derivProj1(self):self.assertTrue(DerivProjfieldsTest(halfSpace(1e-2))) + + # Do a derivative test of Jvec + # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1)) + # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1)) + # def test_derivJvec_zxyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyr',.1)) + # def test_derivJvec_zxyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyi',.1)) + # def test_derivJvec_zyxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxr',.1)) + # def test_derivJvec_zyxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxi',.1)) + # def test_derivJvec_zyyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyr',.1)) + # def test_derivJvec_zyyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyi',.1)) + def test_derivJvec_All(self):self.assertTrue(DerivJvecTest(random(1e-2),'All',.1)) + + # Test the adjoint of Jvec and Jtvec + # def test_JvecAdjoint_zxxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxr',.1)) + # def test_JvecAdjoint_zxxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxi',.1)) + # def test_JvecAdjoint_zxyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyr',.1)) + # def test_JvecAdjoint_zxyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyi',.1)) + # def test_JvecAdjoint_zyxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxr',.1)) + # def test_JvecAdjoint_zyxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxi',.1)) + # def test_JvecAdjoint_zyyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyr',.1)) + # def test_JvecAdjoint_zyyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyi',.1)) + def test_JvecAdjoint_All(self):self.assertTrue(JvecAdjointTest(random(1e-2),'All',.1)) + +if __name__ == '__main__': + unittest.main()