From 6e21d3223086f883fc6410d2ef93a14fed13dbc4 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 29 Apr 2014 10:52:15 -0700 Subject: [PATCH 1/7] test Combinations of rxs and txs for J and Jt --- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 5 +- simpegEM/Tests/test_TDEM_combos.py | 77 ++++++++++++++++++++++ 2 files changed, 79 insertions(+), 3 deletions(-) create mode 100644 simpegEM/Tests/test_TDEM_combos.py diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 071e0de5..36176944 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -22,9 +22,8 @@ class TDEM_bDerivTests(unittest.TestCase): [Maps.ExpMap, Maps.Vertical1DMap, activeMap]) rxOffset = 40. - rxTypes = 'bx,bz' - rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), rxType) for rxType in rxTypes.split(',')] - tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', rxs) + rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz') + tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx]) survey = EM.TDEM.SurveyTDEM([tx]) diff --git a/simpegEM/Tests/test_TDEM_combos.py b/simpegEM/Tests/test_TDEM_combos.py new file mode 100644 index 00000000..f089d182 --- /dev/null +++ b/simpegEM/Tests/test_TDEM_combos.py @@ -0,0 +1,77 @@ +import unittest +from SimPEG import * +import simpegEM as EM + +plotIt = False + +def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1): + cs = 5. + ncx = 20 + ncy = 6 + npad = 20 + hx = [(cs,ncx), (cs,npad,1.3)] + hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] + mesh = Mesh.CylMesh([hx,1,hy], '00C') + + active = mesh.vectorCCz<0. + activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) + mapping = Maps.ComboMap(mesh, [Maps.ExpMap, Maps.Vertical1DMap, activeMap]) + + rxOffset = 40. + + txs = [] + for ii in range(nTx): + rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')] + txs += [EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', rxs)] + + survey = EM.TDEM.SurveyTDEM(txs) + + prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) + # prb.timeSteps = [1e-5] + prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] + # prb.timeSteps = [(1e-05, 100)] + + sigma = np.ones(mesh.nCz)*1e-8 + sigma[mesh.vectorCCz<0] = 1e-1 + sigma = np.log(sigma[active]) + + prb.pair(survey) + return prb, mesh, sigma + +def testJvec(prb, mesh, sigma): + prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) + d_sig = 10*np.random.rand(prb.mapping.nP) + derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] + return Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20) + +def testAdjoint(prb, mesh, sigma): + m = np.random.rand(prb.mapping.nP) + d = np.random.rand(prb.survey.nD) + + V1 = d.dot(prb.Jvec(sigma, m)) + V2 = m.dot(prb.Jtvec(sigma, d)) + print 'AdjointTest', V1, V2 + return np.abs(V1-V2)/np.abs(V1), 1e-6 + +class TDEM_bDerivTests(unittest.TestCase): + + def test_Jvec_bx(self): self.assertTrue(testJvec(*getProb(rxTypes='bx'))) + def test_Adjoint_bx(self): self.assertLess(*testAdjoint(*getProb(rxTypes='bx'))) + + def test_Jvec_bxbz(self): self.assertTrue(testJvec(*getProb(rxTypes='bx,bz'))) + def test_Adjoint_bxbz(self): self.assertLess(*testAdjoint(*getProb(rxTypes='bx,bz'))) + + def test_Jvec_bxbz_2tx(self): self.assertTrue(testJvec(*getProb(rxTypes='bx,bz',nTx=2))) + def test_Adjoint_bxbz_2tx(self): self.assertLess(*testAdjoint(*getProb(rxTypes='bx,bz',nTx=2))) + + def test_Jvec_bxbzbz(self): self.assertTrue(testJvec(*getProb(rxTypes='bx,bz,bz'))) + def test_Adjoint_bxbzbz(self): self.assertLess(*testAdjoint(*getProb(rxTypes='bx,bz,bz'))) + + # def test_Jvec_ey(self): self.assertTrue(testJvec(*getProb(rxTypes='ey'))) + # def test_Adjoint_ey(self): self.assertLess(*testAdjoint(*getProb(rxTypes='ey'))) + + + +if __name__ == '__main__': + unittest.main() From c77a0279c8b4f90ed27a3f9e2bc8e020340383a2 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 29 Apr 2014 11:36:35 -0700 Subject: [PATCH 2/7] Documentation updates. --- docs/api_FDEM.rst | 1 - docs/api_TDEM.rst | 5 +- simpegEM/TDEM/SurveyTDEM.py | 126 ------------------ simpegEM/TDEM/TDEM_b.py | 110 +++++++++++++-- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 14 +- .../Tests/test_TDEM_b_MultiTx_DerivAdjoint.py | 2 +- 6 files changed, 109 insertions(+), 149 deletions(-) diff --git a/docs/api_FDEM.rst b/docs/api_FDEM.rst index faa93cd0..2c063216 100644 --- a/docs/api_FDEM.rst +++ b/docs/api_FDEM.rst @@ -107,4 +107,3 @@ FDEM Survey :show-inheritance: :members: :undoc-members: - diff --git a/docs/api_TDEM.rst b/docs/api_TDEM.rst index 6dbd5446..a77ed227 100644 --- a/docs/api_TDEM.rst +++ b/docs/api_TDEM.rst @@ -55,13 +55,12 @@ TDEM - B formulation :show-inheritance: :members: :undoc-members: - :inherited-members: Field Storage ============= -.. automodule:: simpegEM.TDEM.FieldsTDEM +.. autoclass:: simpegEM.TDEM.SurveyTDEM.FieldsTDEM :show-inheritance: :members: :undoc-members: @@ -71,7 +70,7 @@ Field Storage TDEM Survey Classes =================== -.. automodule:: simpegEM.TDEM.SurveyTDEM +.. autoclass:: simpegEM.TDEM.SurveyTDEM.SurveyTDEM :show-inheritance: :members: :undoc-members: diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 3df191cd..34bbfe42 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -135,129 +135,3 @@ class SurveyTDEM(Survey.BaseSurvey): return f - -# class SurveyTDEM1D(BaseSurvey): -# """ -# docstring for SurveyTDEM1D -# """ - -# txLoc = None #: txLoc -# txType = None #: txType -# rxLoc = None #: rxLoc -# rxType = None #: rxType -# timeCh = None #: timeCh -# nTx = 1 #: Number of transmitters - -# @property -# def nTimeCh(self): -# """Number of time channels""" -# return self.timeCh.size - -# def __init__(self, **kwargs): -# BaseSurvey.__init__(self, **kwargs) -# Utils.setKwargs(self, **kwargs) - -# def projectFields(self, u): -# #TODO: this is hardcoded to 1Tx -# return self.Qrx.dot(u.b[:,:,0].T).T - -# def projectFieldsAdjoint(self, d): -# # TODO: make the following self.nTimeCh -# d = d.reshape((self.prob.nT, self.nTx), order='F') -# #TODO: *Qtime.T need to multiply by a time projection. (outside for loop??) -# ii = 0 -# F = FieldsTDEM(self.prob.mesh, self.nTx, self.prob.nT, 'b') -# for ii in range(self.prob.nT): -# b = self.Qrx.T*d[ii,:] -# F.set_b(b, ii) -# F.set_e(np.zeros((self.prob.mesh.nE,self.nTx)), ii) -# return F - -# #################################################### -# # Interpolation Matrices -# #################################################### - -# @property -# def Qrx(self): -# if self._Qrx is None: -# if self.rxType == 'bz': -# locType = 'Fz' -# self._Qrx = self.prob.mesh.getInterpolationMat(self.rxLoc, locType=locType) -# return self._Qrx -# _Qrx = None - - -# class FieldsTDEM_OLD(object): -# """docstring for FieldsTDEM""" - -# phi0 = None #: Initial electric potential -# A0 = None #: Initial magnetic vector potential -# e0 = None #: Initial electric field -# b0 = None #: Initial magnetic flux density -# j0 = None #: Initial current density -# h0 = None #: Initial magnetic field - -# phi = None #: Electric potential -# A = None #: Magnetic vector potential -# e = None #: Electric field -# b = None #: Magnetic flux density -# j = None #: Current density -# h = None #: Magnetic field - -# def __init__(self, mesh, nTx, nT, store='b'): - -# self.nT = nT #: Number of times -# self.nTx = nTx #: Number of transmitters -# self.mesh = mesh - -# def update(self, newFields, tInd): -# self.set_b(newFields['b'], tInd) -# self.set_e(newFields['e'], tInd) - -# def fieldVec(self): -# u = np.ndarray((0, self.nTx)) -# for i in range(self.nT): -# u = np.r_[u, self.get_b(i), self.get_e(i)] -# if self.nTx == 1: -# u = u.flatten() -# return u - -# #################################################### -# # Get Methods -# #################################################### - -# def get_b(self, ind): -# if ind == -1: -# return self.b0 -# else: -# return self.b[ind,:,:] - -# def get_e(self, ind): -# if ind == -1: -# return self.e0 -# else: -# return self.e[ind,:,:] - -# #################################################### -# # Set Methods -# #################################################### - -# def set_b(self, b, ind): -# if self.b is None: -# self.b = np.zeros((self.nT, np.sum(self.mesh.nF), self.nTx)) -# self.b[:] = np.nan -# if len(b.shape) == 1: -# b = b[:, np.newaxis] -# self.b[ind,:,:] = b - -# def set_e(self, e, ind): -# if self.e is None: -# self.e = np.zeros((self.nT, np.sum(self.mesh.nE), self.nTx)) -# self.e[:] = np.nan -# if len(e.shape) == 1: -# e = e[:, np.newaxis] -# self.e[ind,:,:] = e - - -# def __contains__(self, key): -# return key in self.children diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 207baf8a..65695d00 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -18,7 +18,7 @@ class ProblemTDEM_b(BaseTDEMProblem): def __init__(self, mesh, mapping=None, **kwargs): BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs) - solType = 'b' + solType = 'b' #: Type of the solution, in this case the 'b' field surveyPair = SurveyTDEM @@ -47,16 +47,42 @@ class ProblemTDEM_b(BaseTDEMProblem): #################################################### def Jvec(self, m, v, u=None): - if u is None: - u = self.fields(m) + """ + :param numpy.array m: Conductivity model + :param numpy.ndarray v: vector (model object) + :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m + :rtype: numpy.ndarray + :return: w (data object) + + Multiplying \\\(\\\mathbf{J}\\\) onto a vector can be broken into three steps + + * Compute \\\(\\\\vec{p} = \\\mathbf{G}v\\\) + * Solve \\\(\\\hat{\\\mathbf{A}} \\\\vec{y} = \\\\vec{p}\\\) + * Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\) + + """ + u = u or self.fields(m) p = self.Gvec(m, v, u) y = self.solveAh(m, p) Jv = self.survey.projectFieldsDeriv(u, v=y) return - mkvc(Jv) def Jtvec(self, m, v, u=None): - if u is None: - u = self.fields(m) + """ + :param numpy.array m: Conductivity model + :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object) + :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m + :rtype: numpy.ndarray + :return: w (model object) + + Multiplying \\\(\\\mathbf{J}^\\\\top\\\) onto a vector can be broken into three steps + + * Compute \\\(\\\\vec{p} = \\\mathbf{Q}^\\\\top \\\\vec{v}\\\) + * Solve \\\(\\\hat{\\\mathbf{A}}^\\\\top \\\\vec{y} = \\\\vec{p}\\\) + * Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\) + + """ + u = u or self.fields(m) if not isinstance(v, self.dataPair): v = self.dataPair(self.survey, v) @@ -76,8 +102,7 @@ class ProblemTDEM_b(BaseTDEMProblem): Multiply G by a vector """ - if u is None: - u = self.fields(m) + u = u or self.fields(m) # Note: Fields has shape (nF/E, nTx, nT+1) # However, p will only really fill (:,:,1:nT+1) @@ -104,8 +129,7 @@ class ProblemTDEM_b(BaseTDEMProblem): Multiply G.T by a vector """ - if u is None: - u = self.fields(m) + u = u or self.fields(m) nTx, nE = self.survey.nTx, self.mesh.nE tmp = np.zeros(nE) # Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go. @@ -120,6 +144,38 @@ class ProblemTDEM_b(BaseTDEMProblem): return p def solveAh(self, m, p): + """ + :param numpy.array m: Conductivity model + :param simpegEM.TDEM.FieldsTDEM p: Fields object + :rtype: simpegEM.TDEM.FieldsTDEM + :return: y + + Solve the block-matrix system \\\(\\\hat{A} \\\hat{y} = \\\hat{p}\\\): + + .. math:: + \mathbf{\hat{A}} = \left[ + \\begin{array}{cccc} + A & 0 & & \\\\ + B & A & & \\\\ + & \ddots & \ddots & \\\\ + & & B & A + \end{array} + \\right] \\\\ + \mathbf{A} = + \left[ + \\begin{array}{cc} + \\frac{1}{\delta t} \MfMui & \MfMui\dcurl \\\\ + \dcurl^\\top \MfMui & -\MeSig + \end{array} + \\right] \\\\ + \mathbf{B} = + \left[ + \\begin{array}{cc} + -\\frac{1}{\delta t} \MfMui & 0 \\\\ + 0 & 0 + \end{array} + \\right] \\\\ + """ def AhRHS(tInd, y): rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1] + p[:,'b',tInd+1] @@ -139,6 +195,38 @@ class ProblemTDEM_b(BaseTDEMProblem): return self.forward(m, AhRHS, AhCalcFields) def solveAht(self, m, p): + """ + :param numpy.array m: Conductivity model + :param simpegEM.TDEM.FieldsTDEM p: Fields object + :rtype: simpegEM.TDEM.FieldsTDEM + :return: y + + Solve the block-matrix system \\\(\\\hat{A}^\\\\top \\\hat{y} = \\\hat{p}\\\): + + .. math:: + \mathbf{\hat{A}}^\\top = \left[ + \\begin{array}{cccc} + A & B & & \\\\ + & \ddots & \ddots & \\\\ + & & A & B \\\\ + & & 0 & A + \end{array} + \\right] \\\\ + \mathbf{A} = + \left[ + \\begin{array}{cc} + \\frac{1}{\delta t} \MfMui & \MfMui\dcurl \\\\ + \dcurl^\\top \MfMui & -\MeSig + \end{array} + \\right] \\\\ + \mathbf{B} = + \left[ + \\begin{array}{cc} + -\\frac{1}{\delta t} \MfMui & 0 \\\\ + 0 & 0 + \end{array} + \\right] \\\\ + """ # Mini Example: # @@ -184,7 +272,7 @@ class ProblemTDEM_b(BaseTDEMProblem): # Functions for tests #################################################### - def AhVec(self, m, vec): + def _AhVec(self, m, vec): """ :param numpy.array m: Conductivity model :param simpegEM.TDEM.FieldsTDEM vec: Fields object @@ -229,7 +317,7 @@ class ProblemTDEM_b(BaseTDEMProblem): f[:,'e',i] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i] - self.MeSigma*vec[:,'e',i] return f - def AhtVec(self, m, vec): + def _AhtVec(self, m, vec): """ :param numpy.array m: Conductivity model :param simpegEM.TDEM.FieldsTDEM vec: Fields object diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 36176944..db1a3dfb 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -48,7 +48,7 @@ class TDEM_bDerivTests(unittest.TestCase): sigma = self.sigma u = prb.fields(sigma) - Ahu = prb.AhVec(sigma, u) + Ahu = prb._AhVec(sigma, u) V1 = Ahu[:,'b',1] V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0] @@ -87,7 +87,7 @@ class TDEM_bDerivTests(unittest.TestCase): f = prb.fields(sigma) u1 = A*f.tovec() - u2 = prb.AhVec(sigma,f).tovec() + u2 = prb._AhVec(sigma,f).tovec() self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12) @@ -130,7 +130,7 @@ class TDEM_bDerivTests(unittest.TestCase): for i in range(prb.nT): f[:,'e', i] = np.random.rand(mesh.nE, 1) - Ahf = prb.AhVec(sigma, f) + Ahf = prb._AhVec(sigma, f) f_test = prb.solveAh(sigma, Ahf) u1 = f.tovec() @@ -149,7 +149,7 @@ class TDEM_bDerivTests(unittest.TestCase): dm = 1000*np.random.rand(self.prb.mapping.nP) h = 0.01 - derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] + derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] print '\ntest_DerivG' passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) self.assertTrue(passed) @@ -221,8 +221,8 @@ class TDEM_bDerivTests(unittest.TestCase): f2[:,'b',i] = np.random.rand(mesh.nF, 1) f2[:,'e',i] = np.random.rand(mesh.nE, 1) - V1 = f2.tovec().dot(prb.AhVec(sigma, f1).tovec()) - V2 = f1.tovec().dot(prb.AhtVec(sigma, f2).tovec()) + V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec()) + V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec()) self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) # def test_solveAhtVsAhtVec(self): @@ -236,7 +236,7 @@ class TDEM_bDerivTests(unittest.TestCase): # f1[:,'e',i] = np.random.rand(mesh.nE, 1) # f2 = prb.solveAht(sigma, f1) - # f3 = prb.AhtVec(sigma, f2) + # f3 = prb._AhtVec(sigma, f2) # if True: # import matplotlib.pyplot as plt diff --git a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py index 477952e1..1d7ea8f2 100644 --- a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py @@ -53,7 +53,7 @@ class TDEM_bDerivTests(unittest.TestCase): dm = 1000*np.random.rand(self.prb.mapping.nP) h = 0.01 - derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] + derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] print '\ntest_DerivG' passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) self.assertTrue(passed) From ce2ab576695555a1514b0c8b1e83879f46bbda3d Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 29 Apr 2014 11:42:53 -0700 Subject: [PATCH 3/7] Rearrange methods. --- simpegEM/TDEM/BaseTDEM.py | 57 +++++++++++++++++----- simpegEM/TDEM/TDEM_b.py | 99 ++++++--------------------------------- 2 files changed, 60 insertions(+), 96 deletions(-) diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index 8b3c681f..11193cf4 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -16,17 +16,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): surveyPair = SurveyTDEM - def calcFields(self, sol, tInd): - - if self.solType == 'b': - b = sol - e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - # Todo: implement non-zero js - else: - raise NotImplementedError('solType "%s" is not implemented in CalcFields.' % self.solType) - - return {'b':b, 'e':e} - def fields(self, m): self.curModel = m # Create a fields storage object @@ -73,3 +62,49 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): F[:,:,tInd+1] = CalcFields(sol, tInd) return F + def Jvec(self, m, v, u=None): + """ + :param numpy.array m: Conductivity model + :param numpy.ndarray v: vector (model object) + :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m + :rtype: numpy.ndarray + :return: w (data object) + + Multiplying \\\(\\\mathbf{J}\\\) onto a vector can be broken into three steps + + * Compute \\\(\\\\vec{p} = \\\mathbf{G}v\\\) + * Solve \\\(\\\hat{\\\mathbf{A}} \\\\vec{y} = \\\\vec{p}\\\) + * Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\) + + """ + u = u or self.fields(m) + p = self.Gvec(m, v, u) + y = self.solveAh(m, p) + Jv = self.survey.projectFieldsDeriv(u, v=y) + return - mkvc(Jv) + + def Jtvec(self, m, v, u=None): + """ + :param numpy.array m: Conductivity model + :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object) + :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m + :rtype: numpy.ndarray + :return: w (model object) + + Multiplying \\\(\\\mathbf{J}^\\\\top\\\) onto a vector can be broken into three steps + + * Compute \\\(\\\\vec{p} = \\\mathbf{Q}^\\\\top \\\\vec{v}\\\) + * Solve \\\(\\\hat{\\\mathbf{A}}^\\\\top \\\\vec{y} = \\\\vec{p}\\\) + * Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\) + + """ + u = u or self.fields(m) + + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True) + y = self.solveAht(m, p) + w = self.Gtvec(m, y, u) + return - mkvc(w) + diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 65695d00..fa69b090 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -41,57 +41,22 @@ class ProblemTDEM_b(BaseTDEMProblem): RHS = (1.0/dt)*self.MfMui*B_n return RHS + def calcFields(self, sol, tInd): + + if self.solType == 'b': + b = sol + e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b + # Todo: implement non-zero js + else: + raise NotImplementedError('solType "%s" is not implemented in CalcFields.' % self.solType) + + return {'b':b, 'e':e} + #################################################### # Derivatives #################################################### - def Jvec(self, m, v, u=None): - """ - :param numpy.array m: Conductivity model - :param numpy.ndarray v: vector (model object) - :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m - :rtype: numpy.ndarray - :return: w (data object) - - Multiplying \\\(\\\mathbf{J}\\\) onto a vector can be broken into three steps - - * Compute \\\(\\\\vec{p} = \\\mathbf{G}v\\\) - * Solve \\\(\\\hat{\\\mathbf{A}} \\\\vec{y} = \\\\vec{p}\\\) - * Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\) - - """ - u = u or self.fields(m) - p = self.Gvec(m, v, u) - y = self.solveAh(m, p) - Jv = self.survey.projectFieldsDeriv(u, v=y) - return - mkvc(Jv) - - def Jtvec(self, m, v, u=None): - """ - :param numpy.array m: Conductivity model - :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object) - :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m - :rtype: numpy.ndarray - :return: w (model object) - - Multiplying \\\(\\\mathbf{J}^\\\\top\\\) onto a vector can be broken into three steps - - * Compute \\\(\\\\vec{p} = \\\mathbf{Q}^\\\\top \\\\vec{v}\\\) - * Solve \\\(\\\hat{\\\mathbf{A}}^\\\\top \\\\vec{y} = \\\\vec{p}\\\) - * Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\) - - """ - u = u or self.fields(m) - - if not isinstance(v, self.dataPair): - v = self.dataPair(self.survey, v) - - p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True) - y = self.solveAht(m, p) - w = self.Gtvec(m, y, u) - return - mkvc(w) - def Gvec(self, m, vec, u=None): """ :param numpy.array m: Conductivity model @@ -236,10 +201,10 @@ class ProblemTDEM_b(BaseTDEMProblem): # ^ # fLoc 0 1 2 3 # |-----|-----|-----| - # tInd 0 1 2 / / - # / __/ + # tInd 0 1 2 + # / ___/ # 2 (tInd=2 uses fields 3 and would use 4 but it doesn't exist) - # / __/ + # / ___/ # 1 (tInd=1 uses fields 2 and 3) def AhtRHS(tInd, y): @@ -359,39 +324,3 @@ class ProblemTDEM_b(BaseTDEMProblem): f[:,'b', i] = b f[:,'e', i] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i] - self.MeSigma*vec[:,'e',i] return f - - - -if __name__ == '__main__': - from SimPEG import * - import simpegEM as EM - from simpegEM.Utils.Ana import hzAnalyticDipoleT - from scipy.constants import mu_0 - import matplotlib.pyplot as plt - - cs, ncx, ncz, npad = 5., 20, 6, 20 - 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') - mapping = Maps.Vertical1DMap(mesh) - - opts = {'txLoc':0., - 'txType':'VMD_MVP', - 'rxLoc':np.r_[150., 0., 0.], - 'rxType':'bz', - 'timeCh':np.logspace(-4,-2,20), - } - survey = EM.TDEM.SurveyTDEM1D(**opts) - - prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) - # prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150]) - # prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10]) - prb.timeSteps = [(1e-5, 10)] - prb.pair(survey) - m = np.random.rand(mesh.nCz) - - print survey.dpred(m) - - - - From a834b6974450e85583541b96779ded479d85df61 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 30 Apr 2014 09:11:17 -0700 Subject: [PATCH 4/7] Choose vertical down source for VMD_MVP --- simpegEM/TDEM/TDEM_b.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index fa69b090..360697cc 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -317,10 +317,10 @@ class ProblemTDEM_b(BaseTDEMProblem): """ self.curModel = m f = FieldsTDEM(self.mesh, self.survey) - for i in range(1,self.nT+1): - b = 1.0/self.timeSteps[i-1]*self.MfMui*vec[:,'b',i] + self.MfMui*self.mesh.edgeCurl*vec[:,'e',i] - if i < self.nT: - b = b - 1.0/self.timeSteps[i]*self.MfMui*vec[:,'b',i+1] - f[:,'b', i] = b - f[:,'e', i] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i] - self.MeSigma*vec[:,'e',i] + for i in range(self.nT): + b = 1.0/self.timeSteps[i]*self.MfMui*vec[:,'b',i+1] + self.MfMui*self.mesh.edgeCurl*vec[:,'e',i+1] + if i < self.nT-1: + b = b - 1.0/self.timeSteps[i+1]*self.MfMui*vec[:,'b',i+2] + f[:,'b', i+1] = b + f[:,'e', i+1] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i+1] - self.MeSigma*vec[:,'e',i+1] return f From 8b30a1fe6a2570436ca392cbf92dade1b531f217 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 30 Apr 2014 18:29:41 -0700 Subject: [PATCH 5/7] MVP choose down. --- simpegEM/Utils/Sources/magneticDipole.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpegEM/Utils/Sources/magneticDipole.py b/simpegEM/Utils/Sources/magneticDipole.py index 94381aaf..162af089 100644 --- a/simpegEM/Utils/Sources/magneticDipole.py +++ b/simpegEM/Utils/Sources/magneticDipole.py @@ -36,7 +36,7 @@ def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0. dR = obsLoc - txLoc[i, np.newaxis].repeat(nEdges, axis=0) mCr = np.cross(m, dR) r = np.sqrt((dR**2).sum(axis=1)) - A[:, i] = -(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3) + A[:, i] = +(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3) if nTx == 1: return A.flatten() return A From f0e0cade97ccf9fbe779d306efa0f1dc7886513a Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 30 Apr 2014 19:17:32 -0700 Subject: [PATCH 6/7] FDEM sign switch due to VMD, tdem renaming of 'test*' to dotest so travis doesn't run them.. --- simpegEM/Tests/test_FDEM_analytics.py | 4 ++-- simpegEM/Tests/test_TDEM_combos.py | 25 ++++++++++++------------- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/simpegEM/Tests/test_FDEM_analytics.py b/simpegEM/Tests/test_FDEM_analytics.py index d9d7408d..5772be29 100644 --- a/simpegEM/Tests/test_FDEM_analytics.py +++ b/simpegEM/Tests/test_FDEM_analytics.py @@ -3,7 +3,7 @@ from SimPEG import * import simpegEM as EM from scipy.constants import mu_0 -plotIt = False +plotIt = True class FDEM_analyticTests(unittest.TestCase): @@ -55,7 +55,7 @@ class FDEM_analyticTests(unittest.TestCase): an = EM.Utils.Ana.FEM.hzAnalyticDipoleF(x, self.Tx0.freq, self.sig) - diff = np.log10(np.abs(P*np.imag(u[self.Tx0, 'b']) - np.abs(mu_0*np.imag(an)))) + diff = np.log10(np.abs(P*np.imag(u[self.Tx0, 'b']) - mu_0*np.imag(an))) if plotIt: import matplotlib.pyplot as plt diff --git a/simpegEM/Tests/test_TDEM_combos.py b/simpegEM/Tests/test_TDEM_combos.py index f089d182..c63fbcdd 100644 --- a/simpegEM/Tests/test_TDEM_combos.py +++ b/simpegEM/Tests/test_TDEM_combos.py @@ -38,14 +38,14 @@ def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1): prb.pair(survey) return prb, mesh, sigma -def testJvec(prb, mesh, sigma): +def dotestJvec(prb, mesh, sigma): prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) d_sig = 10*np.random.rand(prb.mapping.nP) derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] return Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20) -def testAdjoint(prb, mesh, sigma): +def dotestAdjoint(prb, mesh, sigma): m = np.random.rand(prb.mapping.nP) d = np.random.rand(prb.survey.nD) @@ -56,21 +56,20 @@ def testAdjoint(prb, mesh, sigma): class TDEM_bDerivTests(unittest.TestCase): - def test_Jvec_bx(self): self.assertTrue(testJvec(*getProb(rxTypes='bx'))) - def test_Adjoint_bx(self): self.assertLess(*testAdjoint(*getProb(rxTypes='bx'))) + def test_Jvec_bx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx'))) + def test_Adjoint_bx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx'))) - def test_Jvec_bxbz(self): self.assertTrue(testJvec(*getProb(rxTypes='bx,bz'))) - def test_Adjoint_bxbz(self): self.assertLess(*testAdjoint(*getProb(rxTypes='bx,bz'))) + def test_Jvec_bxbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz'))) + def test_Adjoint_bxbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz'))) - def test_Jvec_bxbz_2tx(self): self.assertTrue(testJvec(*getProb(rxTypes='bx,bz',nTx=2))) - def test_Adjoint_bxbz_2tx(self): self.assertLess(*testAdjoint(*getProb(rxTypes='bx,bz',nTx=2))) + def test_Jvec_bxbz_2tx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nTx=2))) + def test_Adjoint_bxbz_2tx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nTx=2))) - def test_Jvec_bxbzbz(self): self.assertTrue(testJvec(*getProb(rxTypes='bx,bz,bz'))) - def test_Adjoint_bxbzbz(self): self.assertLess(*testAdjoint(*getProb(rxTypes='bx,bz,bz'))) - - # def test_Jvec_ey(self): self.assertTrue(testJvec(*getProb(rxTypes='ey'))) - # def test_Adjoint_ey(self): self.assertLess(*testAdjoint(*getProb(rxTypes='ey'))) + def test_Jvec_bxbzbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz,bz'))) + def test_Adjoint_bxbzbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz,bz'))) + # def test_Jvec_ey(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey'))) + # def test_Adjoint_ey(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey'))) if __name__ == '__main__': From cb6781a13851dd4118f74b6201dd16fcf1167601 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 30 Apr 2014 21:53:59 -0700 Subject: [PATCH 7/7] turn off plotting in transect analytic test. --- simpegEM/Tests/test_FDEM_analytics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpegEM/Tests/test_FDEM_analytics.py b/simpegEM/Tests/test_FDEM_analytics.py index 5772be29..a1bb9fae 100644 --- a/simpegEM/Tests/test_FDEM_analytics.py +++ b/simpegEM/Tests/test_FDEM_analytics.py @@ -3,7 +3,7 @@ from SimPEG import * import simpegEM as EM from scipy.constants import mu_0 -plotIt = True +plotIt = False class FDEM_analyticTests(unittest.TestCase):