From 4bcc9e08502c99a732646b215adff5f1f54d61ec Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 28 Jan 2016 14:00:47 -0800 Subject: [PATCH] - name cleanup in Jvec and Jtvec (use u instead of f to be consistent with the rest of SimPEG) - example 1D inversion for FDEM --- SimPEG/EM/FDEM/FDEM.py | 30 +++--- SimPEG/Examples/EM_FDEM_1D_Inversion.py | 116 ++++++++++++++++++++++++ 2 files changed, 131 insertions(+), 15 deletions(-) create mode 100644 SimPEG/Examples/EM_FDEM_1D_Inversion.py diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index f2167fd8..cb6adc12 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -53,13 +53,13 @@ class BaseFDEMProblem(BaseEMProblem): return F - def Jvec(self, m, v, f=None): + def Jvec(self, m, v, u=None): """ Sensitivity times a vector """ - if f is None: - f = self.fields(m) + if u is None: + u = self.fields(m) self.curModel = m @@ -71,33 +71,33 @@ class BaseFDEMProblem(BaseEMProblem): for src in self.survey.getSrcByFreq(freq): ftype = self._fieldType + 'Solution' - u_src = f[src, ftype] + u_src = u[src, ftype] dA_dm = self.getADeriv_m(freq, u_src, v) dRHS_dm = self.getRHSDeriv_m(freq, src, v) du_dm = Ainv * ( - dA_dm + dRHS_dm ) for rx in src.rxList: - df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None) + df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None) df_dudu_dm = df_duFun(src, du_dm, adjoint=False) - df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) + df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None) df_dm = df_dmFun(src, v, adjoint=False) Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex) - P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m + P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m Jv[src, rx] = P(Df_Dm) return Utils.mkvc(Jv) - def Jtvec(self, m, v, f=None): + def Jtvec(self, m, v, u=None): """ Sensitivity transpose times a vector """ - if f is None: - f = self.fields(m) + if u is None: + u = self.fields(m) self.curModel = m @@ -113,12 +113,12 @@ class BaseFDEMProblem(BaseEMProblem): for src in self.survey.getSrcByFreq(freq): ftype = self._fieldType + 'Solution' - u_src = f[src, ftype] + u_src = u[src, ftype] for rx in src.rxList: - PTv = rx.projectFieldsDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt u, need possibility wrt m + PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m - df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None) + df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None) df_duT = df_duTFun(src, PTv, adjoint=True) ATinvdf_duT = ATinv * df_duT @@ -127,7 +127,7 @@ class BaseFDEMProblem(BaseEMProblem): dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT - df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) + df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None) dfT_dm = df_dmFun(src, PTv, adjoint=True) du_dmT += dfT_dm @@ -140,7 +140,7 @@ class BaseFDEMProblem(BaseEMProblem): else: raise Exception('Must be real or imag') - return Jtv + return Utils.mkvc(Jtv) def getSourceTerm(self, freq): """ diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py new file mode 100644 index 00000000..a6478016 --- /dev/null +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -0,0 +1,116 @@ +from SimPEG import * +import SimPEG.EM as EM +from scipy.constants import mu_0 + + +def run(plotIt=True): + """ + EM: FDEM: 1D: Inversion + ======================= + + Here we will create and run a FDEM 1D inversion. + + """ + + cs, ncx, ncz, npad = 5., 25, 15, 15 + hx = [(cs,ncx), (cs,npad,1.3)] + hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] + mesh = Mesh.CylMesh([hx,1,hz], '00C') + + layerz = -100. + + active = mesh.vectorCCz<0. + layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=layerz) + actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) + mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap + sig_half = 2e-2 + sig_air = 1e-8 + sig_layer = 1e-2 + sigma = np.ones(mesh.nCz)*sig_air + sigma[active] = sig_half + sigma[layer] = sig_layer + mtrue = np.log(sigma[active]) + + if plotIt: + import matplotlib.pyplot as plt + fig, ax = plt.subplots(1,1, figsize = (3, 6)) + plt.semilogx(sigma[active], mesh.vectorCCz[active]) + ax.set_ylim(-500, 0) + ax.set_xlim(1e-3, 1e-1) + ax.set_xlabel('Conductivity (S/m)', fontsize = 14) + ax.set_ylabel('Depth (m)', fontsize = 14) + ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5) + + + rxOffset=10. + bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi') + + freqs = np.logspace(1,3,10) + srcLoc = np.array([0., 0., 10.]) + + srcList = [] + [srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs] + + survey = EM.FDEM.Survey(srcList) + prb = EM.FDEM.Problem_b(mesh, mapping=mapping) + + try: + from pymatsolver import MumpsSolver + prb.Solver = MumpsSolver + except ImportError, e: + prb.Solver = SolverLU + + prb.pair(survey) + + std = 0.05 + survey.makeSyntheticData(mtrue, std) + + survey.std = std + survey.eps = np.linalg.norm(survey.dtrue)*1e-5 + + if plotIt: + import matplotlib.pyplot as plt + fig, ax = plt.subplots(1,1, figsize = (6, 6)) + ax.semilogx(freqs,survey.dtrue[:freqs.size], 'b.-') + ax.semilogx(freqs,survey.dobs[:freqs.size], 'r.-') + ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16) + ax.set_xlabel('Time (s)', fontsize = 14) + ax.set_ylabel('$B_z$ (T)', fontsize = 16) + ax.set_xlabel('Time (s)', fontsize = 14) + ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5) + + dmisfit = DataMisfit.l2_DataMisfit(survey) + regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]]) + reg = Regularization.Tikhonov(regMesh) + opt = Optimization.InexactGaussNewton(maxIter = 6) + invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt) + + # Create an inversion object + beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2) + betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0) + inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaest]) + m0 = np.log(np.ones(mtrue.size)*sig_half) + reg.alpha_s = 1e-3 + reg.alpha_x = 1. + prb.counter = opt.counter = Utils.Counter() + opt.LSshorten = 0.5 + opt.remember('xc') + + mopt = inv.run(m0) + + if plotIt: + import matplotlib.pyplot as plt + fig, ax = plt.subplots(1,1, figsize = (3, 6)) + plt.semilogx(sigma[active], mesh.vectorCCz[active]) + plt.semilogx(np.exp(mopt), mesh.vectorCCz[active]) + ax.set_ylim(-500, 0) + ax.set_xlim(1e-3, 1e-1) + ax.set_xlabel('Conductivity (S/m)', fontsize = 14) + ax.set_ylabel('Depth (m)', fontsize = 14) + ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5) + plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'],loc='best') + plt.show() + + +if __name__ == '__main__': + run()