diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index c46272ec..44d97424 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -2,6 +2,8 @@ from SimPEG import Utils, Solver from SimPEG.Data import BaseData from SimPEG.Problem import BaseProblem from simpegEM.Utils import Sources +from scipy.constants import mu_0 +from SimPEG.Utils import sdiag, mkvc import numpy as np class DataTDEM1D(BaseData): @@ -19,9 +21,8 @@ class DataTDEM1D(BaseData): BaseData.__init__(self, **kwargs) Utils.setKwargs(self, **kwargs) - def dpred(self, sigma, F=None): - if F is None: F = self.prob.field(sigma) - return self.Qrx.dot(F.b[:,:,0].T) + def projectField(self, u): + return self.Qrx.dot(u.b[:,:,0].T) #################################################### # Interpolation Matrices @@ -110,10 +111,56 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem): """docstring for ProblemTDEM1D""" def __init__(self, mesh, model, **kwargs): BaseProblem.__init__(self, mesh, model, **kwargs) + + + #################################################### + # Physical Properties + #################################################### + + @property + def sigma(self): + return self._sigma + @sigma.setter + def sigma(self, value): + self._sigma = value + _sigma = None + + #################################################### + # Mass Matrices + #################################################### + + @property + def MfMui(self): return self._MfMui + + @property + def MeSigmaI(self): return self._MeSigmaI + + def makeMassMatrices(self, m): + MeSigma = self.mesh.getMass(m, loc='e') + self._MeSigmaI = sdiag(1/MeSigma.diagonal()) + self._MfMui = self.mesh.getMass(1/mu_0, loc='f') + + + def calcFields(self, sol, solType, tInd): + + if solType == 'b': + b = sol + e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b + # Todo: implement non-zero js + else: + errStr = 'solType: ' + solType + raise NotImplementedError(errStr) + + return {'b':b, 'e':e} solveOpts = {'factorize':True,'backend':'scipy'} - def field(self, m): + def field(self, m, useThisRhs=None, useThisCalcFields=None): + RHS = useThisRhs or self.getRHS + CalcFields = useThisCalcFields or self.calcFields + + self.makeMassMatrices(m) + F = self.getInitialFields() dtFact = None for tInd, t in enumerate(self.times): @@ -121,15 +168,18 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem): if dt!=dtFact: dtFact = dt A = self.getA(tInd) - print 'Factoring... (dt = ' + str(dt) + ')' + # print 'Factoring... (dt = ' + str(dt) + ')' Asolve = Solver(A,options=self.solveOpts) - print 'Done' - rhs = self.getRHS(tInd, F) + # print 'Done' + rhs = RHS(tInd, F) sol = Asolve.solve(rhs) if sol.ndim == 1: sol.shape = (sol.size,1) - F.update(sol, tInd, self.solType) + newFields = CalcFields(sol, self.solType, tInd) + F.update(newFields, tInd) return F + + class FieldsTDEM(object): """docstring for FieldsTDEM""" @@ -154,12 +204,9 @@ class FieldsTDEM(object): self.nTx = nTx #: Number of transmitters self.mesh = mesh - def update(self, sol, tInd, solType): - if solType == 'b': - self.set_b(sol, tInd) - else: - errStr = 'solType: ' + solType - raise NotImplementedError(errStr) + def update(self, newFields, tInd): + self.set_b(newFields['b'], tInd) + self.set_e(newFields['e'], tInd) #################################################### # Get Methods @@ -171,6 +218,12 @@ class FieldsTDEM(object): else: return self.b[ind,:,:] + def get_e(self, ind): + if ind == -1: + return self.e0 + else: + return self.e[ind,:,:] + #################################################### # Set Methods #################################################### @@ -179,4 +232,10 @@ class FieldsTDEM(object): if self.b is None: self.b = np.zeros((self.nTimes, np.sum(self.mesh.nF), self.nTx)) self.b[:] = np.nan - self.b[ind, :] = b + self.b[ind,:,:] = b + + def set_e(self, e, ind): + if self.e is None: + self.e = np.zeros((self.nTimes, np.sum(self.mesh.nE), self.nTx)) + self.e[:] = np.nan + self.e[ind,:,:] = e diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index fe560063..8230dfac 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -1,6 +1,6 @@ from BaseTDEM import ProblemBaseTDEM -from scipy.constants import mu_0 -from SimPEG.Utils import sdiag +from BaseTDEM import FieldsTDEM +import numpy as np class ProblemTDEM_b(ProblemBaseTDEM): """ @@ -10,43 +10,6 @@ class ProblemTDEM_b(ProblemBaseTDEM): ProblemBaseTDEM.__init__(self, mesh, model, **kwargs) solType = 'b' - - #################################################### - # Physical Properties - #################################################### - - @property - def sigma(self): - return self._sigma - @sigma.setter - def sigma(self, value): - self._sigma = value - _sigma = None - - #################################################### - # Mass Matrices - #################################################### - - @property - def MfMui(self): - if self._MfMui is None: - self._MfMui = self.mesh.getMass(1/mu_0, loc='f') - return self._MfMui - @MfMui.setter - def MfMui(self, value): - self._MfMui = value - _MfMui = None - - @property - def MeSigmaI(self): - if self._MeSigmaI is None: - MeSigma = self.mesh.getMass(self.sigma, loc='e') - self._MeSigmaI = sdiag(1/MeSigma.diagonal()) - return self._MeSigmaI - @MeSigmaI.setter - def MeSigmaI(self, value): - self._MeSigmaI = value - _MeSigmaI = None #################################################### # Internal Methods @@ -58,7 +21,48 @@ class ProblemTDEM_b(ProblemBaseTDEM): def getRHS(self, tInd, F): dt = self.getDt(tInd) - return (1/dt)*self.MfMui*F.get_b(tInd-1) + return (1/dt)*self.MfMui*F.get_b(tInd-1) + + + #################################################### + # Derivatives + #################################################### + + def J(self, m, v, u=None): + if u is None: + u = self.field(m) + p = self.G(m, v, u) + y = self.solveAh(m, p) + return self.data.projectField(y) + + def G(self, m, v, u=None): + if u is None: + u = self.field(m) + p = FieldsTDEM(self.mesh, 1, self.times.size, 'b') + c = self.mesh.getEdgeMassDeriv()*self.model.transformDeriv(m)*v + for i in range(self.times.size): + ei = u.get_e(i) + pVal = np.empty_like(ei) + for j in range(ei.shape[1]): + pVal[:,j] = -ei[:,j]*c + p.set_e(pVal,i) + return p + + def solveAh(self, m, p): + def AhRHS(tInd, u): + if tInd == 0: + return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + else: + dt = self.getDt(tInd) + return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + 1./dt*self.MfMui*u.get_b(tInd-1) + + def AhCalcFields(sol, solType, tInd): + b = sol + e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b - self.MeSigmaI*p.get_e(tInd) + return {'b':b, 'e':e} + + Y = self.field(m, useThisRhs=AhRHS, useThisCalcFields=AhCalcFields) + return Y if __name__ == '__main__': from SimPEG import * @@ -76,24 +80,37 @@ if __name__ == '__main__': mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2) model = Model.Vertical1DModel(mesh) - txLoc = 0. - txType = 'VMD_MVP' - rxLoc = np.r_[150., 0.] - rxType = 'bz' - timeCh = np.logspace(-4,-2,20) - dat = EM.TDEM.DataTDEM1D(txLoc=txLoc, txType=txType, rxLoc=rxLoc, rxType=rxType, timeCh=timeCh) + opts = {'txLoc':0., + 'txType':'VMD_MVP', + 'rxLoc':np.r_[150., 0.], + 'rxType':'bz', + 'timeCh':np.logspace(-4,-2,20), + } + dat = EM.TDEM.DataTDEM1D(**opts) prb = EM.TDEM.ProblemTDEM_b(mesh, model) prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150]) - prb.sigma = np.ones(mesh.nCz)*1e-8 - prb.sigma[mesh.vectorCCz<0] = 0.1 + sigma = np.ones(mesh.nCz)*1e-8 + sigma[mesh.vectorCCz<0] = 0.1 prb.pair(dat) + f = prb.field(sigma) + # prb.G(prb.sigma, prb.sigma) + # prb.solveAh(prb.sigma, f) + # prb.J(prb.sigma, prb.sigma, f) - bz_calc = dat.dpred(-999999) - bz_ana = mu_0*hzAnalyticDipoleT(dat.rxLoc[0], prb.times, prb.sigma[0]) + from SimPEG.Tests import checkDerivative + m0 = sigma + dx = np.zeros_like(sigma) + dx[prb.mesh.vectorCCz<0] = 1e-4 + derChk = lambda m: [dat.dpred(m), lambda mx: prb.J(m0, mx, u=f)] + passed = checkDerivative(derChk, m0, dx=dx, plotIt=False) + + # bz_calc = dat.dpred(sigma) + # bz_ana = mu_0*hzAnalyticDipoleT(dat.rxLoc[0], prb.times, sigma[0]) + + # plt.loglog(prb.times, np.abs(bz_calc.flatten()), label='TDEM_b') + # plt.loglog(prb.times, np.abs(bz_ana), 'r', label='Analytic') + # plt.legend() + # plt.show() - plt.loglog(prb.times, np.abs(bz_calc.flatten()), label='TDEM_b') - plt.loglog(prb.times, np.abs(bz_ana), 'r', label='Analytic') - plt.legend() - plt.show() diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py index 3a1d5919..1a856bdf 100644 --- a/simpegEM/Tests/test_forward_EMproblem.py +++ b/simpegEM/Tests/test_forward_EMproblem.py @@ -1,15 +1,52 @@ import unittest from SimPEG import * +import simpegEM as EM +from scipy.constants import mu_0 +from simpegEM.Utils.Ana import hzAnalyticDipoleT class EMProblemTests(unittest.TestCase): def setUp(self): - pass + + cs = 5. + ncx = 20 + ncy = 6 + npad = 20 + hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs))) + hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs))) + mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2) + model = Model.Vertical1DModel(mesh) + + opts = {'txLoc':0., + 'txType':'VMD_MVP', + 'rxLoc':np.r_[150., 0.], + 'rxType':'bz', + 'timeCh':np.logspace(-4,-2,20), + } + self.dat = EM.TDEM.DataTDEM1D(**opts) + + self.prb = EM.TDEM.ProblemTDEM_b(mesh, model) + self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150]) + self.sigma = np.ones(mesh.nCz)*1e-8 + self.sigma[mesh.vectorCCz<0] = 0.1 + self.prb.pair(self.dat) + + + def test_analitic_b(self): + bz_calc = self.dat.dpred(self.sigma) + bz_ana = mu_0*hzAnalyticDipoleT(self.dat.rxLoc[0], self.prb.times, self.sigma[0]) + + diff = np.linalg.norm(bz_calc.flatten() - bz_ana.flatten())/np.linalg.norm(bz_ana.flatten()) + print diff + self.assertTrue(diff<0.05) def test_awesome(self): self.assertTrue(True) + + + if __name__ == '__main__': unittest.main()