diff --git a/simpegMT/BaseMT.py b/simpegMT/BaseMT.py index be5ab8c5..12261a01 100644 --- a/simpegMT/BaseMT.py +++ b/simpegMT/BaseMT.py @@ -2,8 +2,7 @@ from simpegEM.FDEM import BaseFDEMProblem from SurveyMT import SurveyMT from DataMT import DataMT from FieldsMT import FieldsMT -from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc -import numpy as np +from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np class BaseMTProblem(BaseFDEMProblem): @@ -75,24 +74,27 @@ class BaseMTProblem(BaseFDEMProblem): # 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 - ftype = self._fieldType + 'Solution' - u_src = u[src, ftype] - dA_dm = self.getADeriv_m(freq, u_src, v) - dRHS_dm = self.getRHSDeriv_m(freq, v) + # 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 ) + du_dm = dA_duI * ( -dA_dm ) else: - du_dm = dA_duI * ( - dA_dm + dRHS_dm ) + du_dm = dA_duI * ( -dA_dm + dRHS_dm ) # Calculate the projection derivatives for rx in src.rxList: # Get the projection derivative - PDeriv_u = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m - Jv[src, rx] = PDeriv_u(du_dm) + # v should be of size nE,2 (each column for 2 polarizations) + PDeriv_u = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # 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) @@ -110,25 +112,28 @@ class BaseMTProblem(BaseFDEMProblem): 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, ftype] + u_src = u[src, :] for rx in src.rxList: # Get the adjoint projectFieldsDeriv - PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m + # 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, dA_duIT, adjoint=True) - dRHS_dmT = self.getRHSDeriv_m(freq, dA_duIT, adjoint=True) + 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 diff --git a/simpegMT/FieldsMT.py b/simpegMT/FieldsMT.py index 505c0510..3d2b23c5 100644 --- a/simpegMT/FieldsMT.py +++ b/simpegMT/FieldsMT.py @@ -228,23 +228,18 @@ class FieldsMT_3D(FieldsMT): def _b_px(self, eSolution, srcList): return self._b_pxPrimary(eSolution, srcList) + self._b_pxSecondary(eSolution, srcList) - def _b_pxSecondaryDeriv_u(self, src, v, adjoint = False): - C = self.mesh.edgeCurl - if adjoint: - return - 1./(1j*omega(src.freq)) * (C.T * v) - return - 1./(1j*omega(src.freq)) * (C * v) - 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 = self.mesh.edgeCurl + C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,1]]) 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 = self.mesh.edgeCurl + C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,1]]) if adjoint: return - 1./(1j*omega(src.freq)) * (C.T * v) return - 1./(1j*omega(src.freq)) * (C * v) diff --git a/simpegMT/ProblemMT1D/Problems.py b/simpegMT/ProblemMT1D/Problems.py index 8b1a8085..d7c2cdb0 100644 --- a/simpegMT/ProblemMT1D/Problems.py +++ b/simpegMT/ProblemMT1D/Problems.py @@ -92,7 +92,7 @@ class eForm_psField(BaseMTProblem): """ Src = self.survey.getSrcByFreq(freq)[0] - S_eDeriv = Src.S_eDeriv(self, v, adjoint) + S_eDeriv = Src.S_eDeriv_m(self, v, adjoint) return -1j * omega(freq) * S_eDeriv def fields(self, m): diff --git a/simpegMT/ProblemMT3D/Problems.py b/simpegMT/ProblemMT3D/Problems.py index 691bbd33..2ae646b1 100644 --- a/simpegMT/ProblemMT3D/Problems.py +++ b/simpegMT/ProblemMT3D/Problems.py @@ -1,4 +1,4 @@ -from SimPEG import Survey, Problem, Utils, Models, np, sp, SolverLU as SimpegSolver +from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver from simpegEM.Utils.EMUtils import omega from scipy.constants import mu_0 from simpegMT.BaseMT import BaseMTProblem @@ -20,7 +20,7 @@ class eForm_ps(BaseMTProblem): """ # From FDEMproblem: Used to project the fields. Currently not used for MTproblem. - _fieldType = 'e_px' + _fieldType = 'e' _eqLocs = 'FE' fieldsPair = FieldsMT_3D _sigmaPrimary = None @@ -52,14 +52,32 @@ class eForm_ps(BaseMTProblem): def getADeriv_m(self, freq, u, v, adjoint=False): - dMe_dsig = self.MeSigmaDeriv( u ) + # 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: - return 1j * omega(freq) * ( dMe_dsig.T * v ) # As in simpegEM - # return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) ) + 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 - return 1j * omega(freq) * ( dMe_dsig * v ) # As in simpegEM - # return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) ) def getRHS(self, freq): """ @@ -80,7 +98,7 @@ class eForm_ps(BaseMTProblem): """ Src = self.survey.getSrcByFreq(freq)[0] - S_eDeriv = Src.S_eDeriv(self, v, adjoint) + S_eDeriv = Src.S_eDeriv_m(self, v, adjoint) return -1j * omega(freq) * S_eDeriv def fields(self, m): diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index d91a41c8..a99e9212 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -173,7 +173,7 @@ class RxMT(Survey.BaseRx): 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 be implement for 2D impedance tensor') + raise NotImplementedError('Has not been implement for 2D impedance tensor') elif self.projType is 'Z3D': if self.locs.ndim == 3: eFLocs = self.locs[:,:,0] @@ -197,14 +197,15 @@ class RxMT(Survey.BaseRx): hx_py = Pbx*mkvc(f[src,'b_py']/mu_0,2) hy_py = Pby*mkvc(f[src,'b_py']/mu_0,2) # Derivatives as lambda functions - 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) - 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 + ex_px_u = lambda vec: sp.hstack((Pex,Pex))*f._e_pxDeriv_u(src,vec) + ey_px_u = lambda vec: sp.hstack((Pey,Pey))*f._e_pxDeriv_u(src,vec) + ex_py_u = lambda vec: sp.hstack((Pex,Pex))*f._e_pyDeriv_u(src,vec) + ey_py_u = lambda vec: sp.hstack((Pey,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: sp.hstack((Pbx,Pbx))*f._b_pxDeriv_u(src,vec)/mu_0 + hy_px_u = lambda vec: sp.hstack((Pby,Pby))*f._b_pxDeriv_u(src,vec)/mu_0 + hx_py_u = lambda vec: sp.hstack((Pbx,Pbx))*f._b_pyDeriv_u(src,vec)/mu_0 + hy_py_u = lambda vec: sp.hstack((Pby,Pby))*f._b_pyDeriv_u(src,vec)/mu_0 # Update the input vector v = mkvc(v,2) # Make v into a column vector @@ -226,7 +227,7 @@ class RxMT(Survey.BaseRx): ZijN_uV = -ey_px_u(hx_py*v) + ey_py*hx_px_u(v) - ey_px*hx_py_u(v) +ey_py_u(hx_px*v) # Calculate the complex derivative - PDeriv_complex = ZijN_uV*Hd - Zij * (Hd_uV*Hd) + PDeriv_complex = ZijN_uV*Hd.toarray() - Zij * (Hd_uV*Hd.toarray()) # Extract the real number for the real/imag components. Pv = np.array(getattr(PDeriv_complex, real_or_imag)) @@ -266,14 +267,14 @@ class RxMT(Survey.BaseRx): ahx_py = mkvc(f[src,'b_py'],2).T/mu_0*Pbx.T ahy_py = 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 + aex_px_u = lambda vec: f._e_pxDeriv_u(src,sp.hstack((Pex,Pex)).T*vec,adjoint=True) + aey_px_u = lambda vec: f._e_pxDeriv_u(src,sp.hstack((Pey,Pey)).T*vec,adjoint=True) + aex_py_u = lambda vec: f._e_pyDeriv_u(src,sp.hstack((Pex,Pex)).T*vec,adjoint=True) + aey_py_u = lambda vec: f._e_pyDeriv_u(src,sp.hstack((Pey,Pey)).T*vec,adjoint=True) + ahx_px_u = lambda vec: f._b_pxDeriv_u(src,sp.hstack((Pbx,Pbx)).T*vec,adjoint=True)/mu_0 + ahy_px_u = lambda vec: f._b_pxDeriv_u(src,sp.hstack((Pby,Pby)).T*vec,adjoint=True)/mu_0 + ahx_py_u = lambda vec: f._b_pyDeriv_u(src,sp.hstack((Pbx,Pbx)).T*vec,adjoint=True)/mu_0 + ahy_py_u = lambda vec: f._b_pyDeriv_u(src,sp.hstack((Pby,Pby)).T*vec,adjoint=True)/mu_0 # Update the input vector v = mkvc(v,2) # Make v into a column vector @@ -376,7 +377,7 @@ class srcMT_polxy_1Dprimary(srcMT): C = problem.mesh.nodalGrad elif problem.mesh.dim == 3: C = problem.mesh.edgeCurl - bBG_bp = (- C * self.ePrimary(problem) )/( 1j*omega(self.freq) ) + bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) )) return bBG_bp def S_e(self,problem): @@ -399,7 +400,7 @@ class srcMT_polxy_1Dprimary(srcMT): Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p) return (Mesigma - Mesigma_p) * e_p - def S_eDeriv(self, problem, v, adjoint = False): + def S_eDeriv_m(self, problem, v, adjoint = False): ''' Get the derivative of S_e wrt to sigma (m) ''' @@ -413,7 +414,12 @@ class srcMT_polxy_1Dprimary(srcMT): 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(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 diff --git a/simpegMT/Tests/test_Problem3D_againstAnalytic.py b/simpegMT/Tests/test_Problem3D_againstAnalytic.py index 347ec03c..cca91ef8 100644 --- a/simpegMT/Tests/test_Problem3D_againstAnalytic.py +++ b/simpegMT/Tests/test_Problem3D_againstAnalytic.py @@ -23,7 +23,8 @@ def getInputs(): """ # 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,10,-1.3),(1000.,2),(1000,10,1.3)]], x0=['C','C','C'])# Setup the model + # 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 @@ -36,6 +37,17 @@ def getInputs(): 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() @@ -70,7 +82,7 @@ def twoLayer(conds): return (M, freqs, sig, sigBG, rx_loc) -def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False): +def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False,expMap=True): M,freqs,sig,sigBG,rx_loc = inputSetup # Make a receiver list rxList = [] @@ -81,9 +93,9 @@ def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False): rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,comp)) # Source list srcList =[] - sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:] + if singleFreq: - srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freqs[2])) + srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,singleFreq)) else: for freq in freqs: srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq)) @@ -91,16 +103,21 @@ def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False): survey = simpegmt.SurveyMT.SurveyMT(srcList) ## Setup the problem object - problem = simpegmt.ProblemMT3D.eForm_ps(M,sigmaPrimary=sigma1d) + sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:] + if expMap: + problem = simpegmt.ProblemMT3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) ) + problem.mapping = simpeg.Maps.ExpMap(problem.mesh) + problem.curModel = np.log(sig) + else: + problem = simpegmt.ProblemMT3D.eForm_ps(M,sigmaPrimary= sigma1d) + problem.curModel = sig + problem.pair(survey) problem.verbose = False try: from pymatsolver import MumpsSolver problem.Solver = MumpsSolver except: pass - problem.pair(survey) - problem.curMod = sig - problem.mapping = simpeg.Maps.ExpMap(problem.mesh) return (survey, problem) @@ -113,34 +130,49 @@ def getAppResPhs(MTdata): recData = MTdata.toRecArray('Complex') return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx']) -def adjointTest(inputSetup): - - survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup) - print 'Adjoint test of eForm primary/secondary\n' - - m = problem.curMod - - # if addrandoms is True: - # m = m + np.random.randn(problem.mesh.nC)*CONDUCTIVITY*1e-1 +def test_JvecAdjoint(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) + v = np.random.rand(survey.nD,) # print problem.PropMap.PropModel.nP - w = np.random.rand(problem.mesh.nC) + w = np.random.rand(problem.mesh.nC,) - vJw = v.dot(problem.Jvec(m, w, u)) - wJtv = w.dot(problem.Jtvec(m, v, u)) + 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 test_DerivJvec(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 derivProjfields(inputSetup): - survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup) +def test_DerivProjfields(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) # Define a src and rx src = survey.srcList[-1] @@ -158,39 +190,58 @@ def derivProjfields(inputSetup): return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR) -def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True): +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)) + 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.curMod)) + 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 < .35) + 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 < .35) + return np.all(np.abs(app_rpxy[1,:] + np.ones(survey.nFreq)*135) / 135 < .4) class TestAnalytics(unittest.TestCase): def setUp(self): pass - # def test_appRes2en1(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(2e-1)) - 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)) + # # 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)) + # 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_deriv1(self):self.assertTrue(derivProjfields(halfSpace(1e-3))) + def test_derivProj1(self):self.assertTrue(test_DerivProjfields(halfSpace(1e-2))) + + # Do a derivative test of Jvec + # def test_derivJvec_zxxr(self):self.assertTrue(test_DerivJvec(random(1e-2),'zxxr',1)) + # def test_derivJvec_zxxi(self):self.assertTrue(test_DerivJvec(random(1e-2),'zxxi',1)) + def test_derivJvec_zxyr(self):self.assertTrue(test_DerivJvec(random(1e-2),'zxyr',.1)) + # def test_derivJvec_zxyi(self):self.assertTrue(test_DerivJvec(random(1e-2),'zxyi',1)) + def test_derivJvec_zyxr(self):self.assertTrue(test_DerivJvec(random(1e-2),'zyxr',1)) + # def test_derivJvec_zyxi(self):self.assertTrue(test_DerivJvec(random(1e-2),'zyxi',1)) + # def test_derivJvec_zyyr(self):self.assertTrue(test_DerivJvec(random(1e-2),'zyyr',1)) + # def test_derivJvec_zyyi(self):self.assertTrue(test_DerivJvec(random(1e-2),'zyyi',1)) + # def test_derivJvec_All(self):self.assertTrue(test_DerivJvec(random(1e-2),'All',1)) # Test the adjoint of Jvec and Jtvec - def test_adjointDeriv1(self):self.assertTrue(adjointTest(halfSpace(1e-3))) + # def test_JvecAdjoint_zxxr(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zxxr',1)) + # def test_JvecAdjoint_zxxi(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zxxi',1)) + def test_JvecAdjoint_zxyr(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zxyr',.1)) + # def test_JvecAdjoint_zxyi(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zxyi',1)) + def test_JvecAdjoint_zyxr(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zyxr',1)) + # def test_JvecAdjoint_zyxi(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zyxi',1)) + # def test_JvecAdjoint_zyyr(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zyyr',1)) + # def test_JvecAdjoint_zyyi(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'zyyi',1)) + # def test_JvecAdjoint_All(self):self.assertTrue(test_JvecAdjoint(random(1e-2),'All',1)) if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/simpegMT/Utils/dataUtils.py b/simpegMT/Utils/dataUtils.py index 1c0dd742..f401bc59 100644 --- a/simpegMT/Utils/dataUtils.py +++ b/simpegMT/Utils/dataUtils.py @@ -3,6 +3,8 @@ import numpy as np, matplotlib.pyplot as plt, sys import SimPEG as simpeg import simpegMT as simpegmt import numpy.lib.recfunctions as recFunc +from scipy.constants import mu_0 + def getAppRes(MTdata): # Make impedance zList = [] @@ -23,6 +25,10 @@ def appResPhs(freq,z): 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)))