From 096a72fdb41392d1a6e2e7f1ba09f48758f08d69 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Tue, 26 Jan 2016 12:07:42 -0800 Subject: [PATCH] Tipper derivatives and adjoint working. --- SimPEG/MT/SurveyMT.py | 78 +++++++++++----------- tests/mt/test_Problem3D_againstAnalytic.py | 77 ++++++++------------- 2 files changed, 68 insertions(+), 87 deletions(-) diff --git a/SimPEG/MT/SurveyMT.py b/SimPEG/MT/SurveyMT.py index c8cae777..b86ac6c2 100644 --- a/SimPEG/MT/SurveyMT.py +++ b/SimPEG/MT/SurveyMT.py @@ -230,7 +230,7 @@ class Rx(SimPEGsurvey.BaseRx): # Calculate the complex derivative PDeriv_complex = Hd * (ZijN_uV - Zij * Hd_uV ) - elif self.projTyep is 'T3D': + elif self.projType is 'T3D': if self.locs.ndim == 3: eFLocs = self.locs[:,:,0] bFLocs = self.locs[:,:,1] @@ -244,31 +244,31 @@ class Rx(SimPEGsurvey.BaseRx): # Get the fields at location # px: x-polaration and py: y-polaration. - hx_px = Pbx*f[src,'b_px']/mu_0 - hy_px = Pby*f[src,'b_px']/mu_0 - hz_px = Pbz*f[src,'b_px']/mu_0 - hx_py = Pbx*f[src,'b_py']/mu_0 - hy_py = Pby*f[src,'b_py']/mu_0 - hz_py = Pbz*f[src,'b_py']/mu_0 + 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'] # Derivatives as lambda functions # 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 - hz_px_u = lambda vec: Pbz*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 - hz_py_u = lambda vec: Pbz*f._b_pyDeriv_u(src,vec)/mu_0 + bx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec) + by_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec) + bz_px_u = lambda vec: Pbz*f._b_pxDeriv_u(src,vec) + bx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec) + by_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec) + bz_py_u = lambda vec: Pbz*f._b_pyDeriv_u(src,vec) # 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) - if 'zxx' in self.rxType: - Tij = sDiag(Hd*( sDiag(hy_px)*hz_py + sDiag(hy_py)*hz_px )) - TijN_uV = sDiag(hy_px)*hz_py_u(v) + sDiag(hz_py)*hy_px_u(v) + sDiag(hy_py)*hz_px_u(v) + sDiag(hz_px)*hy_py_u(v) - elif 'zxy' in self.rxType: - Tij = sDiag(Hd*( sDiag(hx_px)*hz_py + sDiag(hx_py)*hz_px )) - TijN_uV = sDiag(hz_py)*hx_px_u(v) - sDiag(hx_px)*hz_py_u(v) + sDiag(hx_py)*hz_px_u(v) + sDiag(hz_px)*hx_py_u(v) + Hd = sDiag(1./(sDiag(bx_px)*by_py - sDiag(bx_py)*by_px)) + Hd_uV = sDiag(by_py)*bx_px_u(v) + sDiag(bx_px)*by_py_u(v) - sDiag(bx_py)*by_px_u(v) - sDiag(by_px)*bx_py_u(v) + if 'tzx' in self.rxType: + Tij = sDiag(Hd*( - sDiag(by_px)*bz_py + sDiag(by_py)*bz_px )) + TijN_uV = -sDiag(by_px)*bz_py_u(v) - sDiag(bz_py)*by_px_u(v) + sDiag(by_py)*bz_px_u(v) + sDiag(bz_px)*by_py_u(v) + elif 'tzy' in self.rxType: + Tij = sDiag(Hd*( sDiag(bx_px)*bz_py - sDiag(bx_py)*bz_px )) + TijN_uV = sDiag(bz_py)*bx_px_u(v) + sDiag(bx_px)*bz_py_u(v) - sDiag(bx_py)*bz_px_u(v) - sDiag(bz_px)*bx_py_u(v) # Calculate the complex derivative PDeriv_complex = Hd * (TijN_uV - Tij * Hd_uV ) @@ -357,34 +357,34 @@ class Rx(SimPEGsurvey.BaseRx): Pbz = mesh.getInterpolationMat(bFLocs,'Fz') # Get the fields at location # px: x-polaration and py: y-polaration. - 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) - ahz_px = mkvc(mkvc(f[src,'b_px'],2).T/mu_0*Pbz.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) - ahz_py = mkvc(mkvc(f[src,'b_py'],2).T/mu_0*Pbz.T) + abx_px = mkvc(mkvc(f[src,'b_px'],2).T*Pbx.T) + aby_px = mkvc(mkvc(f[src,'b_px'],2).T*Pby.T) + abz_px = mkvc(mkvc(f[src,'b_px'],2).T*Pbz.T) + abx_py = mkvc(mkvc(f[src,'b_py'],2).T*Pbx.T) + aby_py = mkvc(mkvc(f[src,'b_py'],2).T*Pby.T) + abz_py = mkvc(mkvc(f[src,'b_py'],2).T*Pbz.T) # Derivatives as lambda functions - 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 - ahz_px_u = lambda vec: f._b_pxDeriv_u(src,Pbz.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 - ahz_py_u = lambda vec: f._b_pyDeriv_u(src,Pbz.T*vec,adjoint=True)/mu_0 + abx_px_u = lambda vec: f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True) + aby_px_u = lambda vec: f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True) + abz_px_u = lambda vec: f._b_pxDeriv_u(src,Pbz.T*vec,adjoint=True) + abx_py_u = lambda vec: f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True) + aby_py_u = lambda vec: f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True) + abz_py_u = lambda vec: f._b_pyDeriv_u(src,Pbz.T*vec,adjoint=True) # 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) + aHd = sDiag(1./(sDiag(abx_px)*aby_py - sDiag(abx_py)*aby_px)) + aHd_uV = lambda x: abx_px_u(sDiag(aby_py)*x) + abx_px_u(sDiag(aby_py)*x) - aby_px_u(sDiag(abx_py)*x) - abx_py_u(sDiag(aby_px)*x) # Need to fix this to reflect the adjoint if 'tzx' in self.rxType: - Tij = sDiag(aHd*( sDiag(ahz_py)*ahy_px - sDiag(ahz_px)*ahy_py)) - TijN_uV = lambda x: ahz_px_u(sDiag(ahy_px)*x) + ahy_px_u(sDiag(ahz_px)*x) + ahy_py_u(sDiag(ahz_px)*x) + ahz_px_u(sDiag(ahy_py)*x) + Tij = sDiag(aHd*( -sDiag(abz_py)*aby_px + sDiag(abz_px)*aby_py)) + TijN_uV = lambda x: -abz_py_u(sDiag(aby_px)*x) - aby_px_u(sDiag(abz_py)*x) + aby_py_u(sDiag(abz_px)*x) + abz_px_u(sDiag(aby_py)*x) elif 'tzy' in self.rxType: - Tij = sDiag(aHd*( sDiag(ahz_py)*ahx_px + sDiag(ahz_px)*ahx_py)) - TijN_uV = lambda x: ahx_px_u(sDiag(ahz_py)*x) + ahz_py_u(sDiag(ahx_px)*x) + ahx_py_u(sDiag(ahz_px)*x) + ahz_px_u(sDiag(ahx_py)*x) + Tij = sDiag(aHd*( sDiag(abz_py)*abx_px - sDiag(abz_px)*abx_py)) + TijN_uV = lambda x: abx_px_u(sDiag(abz_py)*x) + abz_py_u(sDiag(abx_px)*x) - abx_py_u(sDiag(abz_px)*x) - abz_px_u(sDiag(abx_py)*x) # Calculate the complex derivative PDeriv_real = TijN_uV(aHd*v) - aHd_uV(Tij.T*aHd*v)# # NOTE: Need to reshape the output to go from 2*nU array to a (nU,2) matrix for each polarization diff --git a/tests/mt/test_Problem3D_againstAnalytic.py b/tests/mt/test_Problem3D_againstAnalytic.py index 5a07d39a..1cef276d 100644 --- a/tests/mt/test_Problem3D_againstAnalytic.py +++ b/tests/mt/test_Problem3D_againstAnalytic.py @@ -30,12 +30,11 @@ def getInputs(): ## 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)))) + rx_x, rx_y = np.meshgrid(np.arange(-1000,1001,500),np.arange(-1000,1001,500)) + rx_loc = 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() @@ -63,6 +62,22 @@ def halfSpace(conds): return (M, freqs, sig, sigBG, rx_loc) +def blockInhalfSpace(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 = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds) + sig[~groundInd] = 1e-8 + # 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 twoLayer(conds): ''' Returns a 2 layer model based on the conductivity values given''' M, freqs, rx_loc, elev = getInputs() @@ -81,13 +96,21 @@ def twoLayer(conds): return (M, freqs, sig, sigBG, rx_loc) -def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False,expMap=True): + + +def setupSimpegMTfwd_eForm_ps(inputSetup,comp='Imp',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)) + for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']: + rxList.append(MT.Rx(rx_loc,rxType)) + elif comp == 'Imp': + for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']: + rxList.append(MT.Rx(rx_loc,rxType)) + elif comp == 'Tip': + for rxType in ['tzxr','tzxi','tzyr','tzyi']: + rxList.append(MT.Rx(rx_loc,rxType)) else: rxList.append(MT.Rx(rx_loc,comp)) # Source list @@ -120,46 +143,6 @@ def setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=False,expMap=True 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): @@ -206,7 +189,6 @@ def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True): 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) @@ -232,7 +214,6 @@ def DerivProjfieldsTest(inputSetup,comp='All',freq=False): 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'