From 200a1680cae5a768fbd729c563225bb8172ba9c5 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 24 Apr 2014 16:50:58 -0700 Subject: [PATCH] Switch to BaseTimeProblem --- simpegEM/TDEM/BaseTDEM.py | 66 ++++---------------- simpegEM/TDEM/FieldsTDEM.py | 10 +-- simpegEM/TDEM/SurveyTDEM.py | 6 +- simpegEM/TDEM/TDEM_b.py | 34 +++++----- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 60 +++++++++--------- simpegEM/Tests/test_TDEM_forward_Analytic.py | 10 +-- 6 files changed, 71 insertions(+), 115 deletions(-) diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index 05a084a1..c30c1137 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -1,9 +1,10 @@ from SimPEG import Solver -from SimPEG.Problem import BaseProblem +from SimPEG.Problem import BaseTimeProblem from simpegEM.Utils import Sources from FieldsTDEM import FieldsTDEM from scipy.constants import mu_0 from SimPEG.Utils import sdiag, mkvc +from SimPEG import Utils, Mesh import numpy as np @@ -37,63 +38,18 @@ class MixinInitialFieldCalc(object): raise Exception('Unknown mesh for VMD') # Initialize field object - F = FieldsTDEM(self.mesh, 1, self.times.size, store=self.storeTheseFields) + F = FieldsTDEM(self.mesh, 1, self.times[1:].size, store=self.storeTheseFields) # Set initial B F.b0 = self.mesh.edgeCurl*MVP return F -class MixinTimeStuff(object): - """docstring for MixinTimeStuff""" - def dt(): - doc = "Size of time steps" - def fget(self): - return self._dt - def fdel(self): - del self._dt - return locals() - dt = property(**dt()) - - def nsteps(): - doc = "Number of steps to take" - def fget(self): - return self._nsteps - def fdel(self): - del self._nsteps - return locals() - nsteps = property(**nsteps()) - - def times(): - doc = "Modeling times" - def fget(self): - t = np.r_[1:self.nsteps[0]+1]*self.dt[0] - for i in range(1,self.dt.size): - t = np.r_[t, np.r_[1:self.nsteps[i]+1]*self.dt[i]+t[-1]] - return t - return locals() - times = property(**times()) - - def getDt(self, tInd): - return np.concatenate([self.dt[i].repeat(self.nsteps[i]) for i in range(self.dt.size)])[tInd] - - def setTimes(self, dt, nsteps): - dt = np.array(dt) - nsteps = np.array(nsteps) - assert dt.size==nsteps.size, "dt, nsteps must be same length" - self._dt = dt - self._nsteps = nsteps - - @property - def nTimes(self): - return self.times.size - - -class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem): +class ProblemBaseTDEM(MixinInitialFieldCalc, BaseTimeProblem): """docstring for ProblemTDEM1D""" def __init__(self, mesh, mapping=None, **kwargs): - BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs) + BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs) #################################################### @@ -151,11 +107,11 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem): def forward(self, m, RHS, CalcFields, F=None): if F is None: - F = FieldsTDEM(self.mesh, self.survey.nTx, self.nTimes, store=self.storeTheseFields) + F = FieldsTDEM(self.mesh, self.survey.nTx, self.nT, store=self.storeTheseFields) dtFact = None - for tInd, t in enumerate(self.times): - dt = self.getDt(tInd) + for tInd, t in enumerate(self.times[1:]): + dt = self.timeSteps[tInd] if dt!=dtFact: dtFact = dt A = self.getA(tInd) @@ -172,11 +128,11 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem): def adjoint(self, m, RHS, CalcFields, F=None): if F is None: - F = FieldsTDEM(self.mesh, self.survey.nTx, self.nTimes, store=self.storeTheseFields) + F = FieldsTDEM(self.mesh, self.survey.nTx, self.nT, store=self.storeTheseFields) dtFact = None - for tInd, t in reversed(list(enumerate(self.times))): - dt = self.getDt(tInd) + for tInd, t in reversed(list(enumerate(self.times[1:]))): + dt = self.timeSteps[tInd] if dt!=dtFact: dtFact = dt A = self.getA(tInd) diff --git a/simpegEM/TDEM/FieldsTDEM.py b/simpegEM/TDEM/FieldsTDEM.py index 25234892..79e54ed2 100644 --- a/simpegEM/TDEM/FieldsTDEM.py +++ b/simpegEM/TDEM/FieldsTDEM.py @@ -18,9 +18,9 @@ class FieldsTDEM(object): j = None #: Current density h = None #: Magnetic field - def __init__(self, mesh, nTx, nTimes, store='b'): + def __init__(self, mesh, nTx, nT, store='b'): - self.nTimes = nTimes #: Number of times + self.nT = nT #: Number of times self.nTx = nTx #: Number of transmitters self.mesh = mesh @@ -30,7 +30,7 @@ class FieldsTDEM(object): def fieldVec(self): u = np.ndarray((0, self.nTx)) - for i in range(self.nTimes): + 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() @@ -58,7 +58,7 @@ class FieldsTDEM(object): def set_b(self, b, ind): if self.b is None: - self.b = np.zeros((self.nTimes, np.sum(self.mesh.nF), self.nTx)) + 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] @@ -66,7 +66,7 @@ class FieldsTDEM(object): 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.zeros((self.nT, np.sum(self.mesh.nE), self.nTx)) self.e[:] = np.nan if len(e.shape) == 1: e = e[:, np.newaxis] diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 4b1a91b3..4e712fd4 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -29,11 +29,11 @@ class SurveyTDEM1D(BaseSurvey): def projectFieldsAdjoint(self, d): # TODO: make the following self.nTimeCh - d = d.reshape((self.prob.nTimes, self.nTx), order='F') + 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.nTimes, 'b') - for ii in range(self.prob.nTimes): + 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) diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 5a901bf0..5a346965 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -34,11 +34,11 @@ class ProblemTDEM_b(ProblemBaseTDEM): :return: A """ - dt = self.getDt(tInd) + dt = self.timeSteps[tInd] return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + (1.0/dt)*self.MfMui def getRHS(self, tInd, F): - dt = self.getDt(tInd) + dt = self.timeSteps[tInd] return (1.0/dt)*self.MfMui*F.get_b(tInd-1) @@ -73,10 +73,10 @@ class ProblemTDEM_b(ProblemBaseTDEM): """ if u is None: u = self.fields(m) - p = FieldsTDEM(self.mesh, 1, self.times.size, 'b') + p = FieldsTDEM(self.mesh, 1, self.nT, 'b') curModel = self.mapping.transform(m) c = self.mesh.getEdgeInnerProductDeriv(curModel)*self.mapping.transformDeriv(m)*vec - for i in range(self.times.size): + for i in range(self.nT): ei = u.get_e(i) pVal = np.empty_like(ei) for j in range(ei.shape[1]): @@ -90,7 +90,7 @@ class ProblemTDEM_b(ProblemBaseTDEM): if u is None: u = self.fields(m) tmp = np.zeros((self.mesh.nE,self.survey.nTx)) - for i in range(self.nTimes): + for i in range(self.nT): tmp += v.get_e(i)*u.get_e(i) curModel = self.mapping.transform(m) @@ -102,7 +102,7 @@ class ProblemTDEM_b(ProblemBaseTDEM): rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + p.get_b(tInd) if tInd == 0: return rhs - dt = self.getDt(tInd) + dt = self.timeSteps[tInd] return rhs + 1.0/dt*self.MfMui*u.get_b(tInd-1) def AhCalcFields(sol, solType, tInd): @@ -117,9 +117,9 @@ class ProblemTDEM_b(ProblemBaseTDEM): def AhtRHS(tInd, u): rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + p.get_b(tInd) - if tInd == self.nTimes-1: + if tInd == self.nT-1: return rhs - dt = self.getDt(tInd+1) + dt = self.timeSteps[tInd+1] return rhs + 1.0/dt*self.MfMui*u.get_b(tInd+1) def AhtCalcFields(sol, solType, tInd): @@ -169,14 +169,14 @@ class ProblemTDEM_b(ProblemBaseTDEM): """ self.makeMassMatrices(sigma) - dt = self.getDt(0) + dt = self.timeSteps[0] b = 1.0/dt*self.MfMui*vec.get_b(0) + self.MfMui*self.mesh.edgeCurl*vec.get_e(0) e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(0) - self.MeSigma*vec.get_e(0) - f = FieldsTDEM(self.mesh, 1, self.times.size, 'b') + f = FieldsTDEM(self.mesh, 1, self.nT, 'b') f.set_b(b, 0) f.set_e(e, 0) - for i in range(1,self.nTimes): - dt = self.getDt(i) + for i in range(1,self.nT): + dt = self.timeSteps[i] b = 1.0/dt*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1.0/dt*self.MfMui*vec.get_b(i-1) e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(i) - self.MeSigma*vec.get_e(i) f.set_b(b, i) @@ -217,14 +217,14 @@ class ProblemTDEM_b(ProblemBaseTDEM): \\right] \\\\ """ self.makeMassMatrices(sigma) - f = FieldsTDEM(self.mesh, 1, self.times.size, 'b') - for i in range(self.nTimes-1): - b = 1/self.getDt(i)*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1/self.getDt(i+1)*self.MfMui*vec.get_b(i+1) + f = FieldsTDEM(self.mesh, 1, self.nT, 'b') + for i in range(self.nT-1): + b = 1.0/self.timeSteps[i]*self.MfMui*vec.get_b(i) + self.MfMui*self.mesh.edgeCurl*vec.get_e(i) - 1.0/self.timeSteps[i+1]*self.MfMui*vec.get_b(i+1) e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(i) - self.MeSigma*vec.get_e(i) f.set_b(b, i) f.set_e(e, i) - N = self.nTimes - 1 - b = 1/self.getDt(N)*self.MfMui*vec.get_b(N) + self.MfMui*self.mesh.edgeCurl*vec.get_e(N) + N = self.nT - 1 + b = 1.0/self.timeSteps[N]*self.MfMui*vec.get_b(N) + self.MfMui*self.mesh.edgeCurl*vec.get_e(N) e = self.mesh.edgeCurl.T*self.MfMui*vec.get_b(N) - self.MeSigma*vec.get_e(N) f.set_b(b, N) f.set_e(e, N) diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index fe5d98ae..7d136a29 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -29,7 +29,7 @@ class TDEM_bDerivTests(unittest.TestCase): self.dat = EM.TDEM.SurveyTDEM1D(**opts) self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) - self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10]) + self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (0.00025, 10)] self.sigma = np.ones(mesh.nCz)*1e-8 self.sigma[mesh.vectorCCz<0] = 1e-1 @@ -50,16 +50,16 @@ class TDEM_bDerivTests(unittest.TestCase): Ahu = prb.AhVec(sigma, u) V1 = Ahu.get_b(0) - V2 = 1./prb.getDt(0)*prb.MfMui*u.get_b(-1) + V2 = 1./prb.timeSteps[0]*prb.MfMui*u.get_b(-1) # print np.linalg.norm(V1-V2), np.linalg.norm(V2), np.linalg.norm(V1-V2)/np.linalg.norm(V2) # self.assertTrue(np.linalg.norm(V1-V2)/np.linalg.norm(V2) < 1.e-6) V1 = Ahu.get_e(0) self.assertTrue(np.linalg.norm(V1) < 1.e-6) - for i in range(1,u.nTimes): + for i in range(1,u.nT): - dt = prb.getDt(i) + dt = prb.timeSteps[i] V1 = Ahu.get_b(i) V2 = 1/dt*prb.MfMui*u.get_b(i-1) @@ -72,11 +72,11 @@ class TDEM_bDerivTests(unittest.TestCase): def test_AhVecVSMat_OneTS(self): prb = self.prb - prb.setTimes([1e-5], [1]) + prb.timeSteps = [(1e-05, 1)] sigma = self.sigma prb.makeMassMatrices(sigma) - dt = prb.getDt(0) + dt = prb.timeSteps[0] a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF) a12 = prb.MfMui*prb.mesh.edgeCurl a21 = prb.mesh.edgeCurl.T*prb.MfMui @@ -92,12 +92,12 @@ class TDEM_bDerivTests(unittest.TestCase): def test_solveAhVSMat_OneTS(self): prb = self.prb - prb.setTimes([1e-5], [1]) + prb.timeSteps = [(1e-05, 1)] sigma = self.sigma prb.makeMassMatrices(sigma) - dt = prb.getDt(0) + dt = prb.timeSteps[0] a11 = 1/dt*prb.MfMui*sp.eye(prb.mesh.nF) a12 = prb.MfMui*prb.mesh.edgeCurl a21 = prb.mesh.edgeCurl.T*prb.MfMui @@ -120,8 +120,8 @@ class TDEM_bDerivTests(unittest.TestCase): sigma = self.sigma self.prb.makeMassMatrices(sigma) - f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.times.size, 'b') - for i in range(f.nTimes): + f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') + for i in range(f.nT): f.set_b(np.zeros((mesh.nF, 1)), i) f.set_e(np.random.rand(mesh.nE, 1), i) @@ -152,7 +152,7 @@ class TDEM_bDerivTests(unittest.TestCase): def test_Deriv_dUdM(self): prb = self.prb - prb.setTimes([1e-5, 1e-4, 1e-3], [10, 10, 10]) + prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] mesh = self.mesh sigma = self.sigma @@ -168,7 +168,7 @@ class TDEM_bDerivTests(unittest.TestCase): def test_Deriv_J(self): prb = self.prb - prb.setTimes([1e-5, 1e-4, 1e-3], [10, 10, 10]) + prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] mesh = self.mesh sigma = self.sigma @@ -188,11 +188,11 @@ class TDEM_bDerivTests(unittest.TestCase): mesh = self.mesh # Generate random fields and data - f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.times.size, 'b') - for i in range(f.nTimes): + f = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') + for i in range(f.nT): f.set_b(np.random.rand(mesh.nF, 1), i) f.set_e(np.random.rand(mesh.nE, 1), i) - d = np.random.rand(dat.prob.nTimes, dat.nTx) + d = np.random.rand(dat.prob.nT, dat.nTx) # Check that d.T*Q*f = f.T*Q.T*d V1 = d.T.dot(dat.projectFields(f)) @@ -205,13 +205,13 @@ class TDEM_bDerivTests(unittest.TestCase): mesh = self.mesh sigma = self.sigma - f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b') - for i in range(f1.nTimes): + f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') + for i in range(f1.nT): f1.set_b(np.random.rand(mesh.nF, 1), i) f1.set_e(np.random.rand(mesh.nE, 1), i) - f2 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b') - for i in range(f2.nTimes): + f2 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') + for i in range(f2.nT): f2.set_b(np.random.rand(mesh.nF, 1), i) f2.set_e(np.random.rand(mesh.nE, 1), i) @@ -224,8 +224,8 @@ class TDEM_bDerivTests(unittest.TestCase): mesh = self.mesh sigma = np.random.rand(prb.mapping.nP) - f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nTimes, 'b') - for i in range(f1.nTimes): + f1 = EM.TDEM.FieldsTDEM(mesh, 1, prb.nT, 'b') + for i in range(f1.nT): f1.set_b(np.random.rand(mesh.nF, 1), i) f1.set_e(np.random.rand(mesh.nE, 1), i) @@ -241,13 +241,13 @@ class TDEM_bDerivTests(unittest.TestCase): mesh = self.mesh sigma = self.sigma - f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b') - for i in range(f1.nTimes): + f1 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') + for i in range(f1.nT): f1.set_b(np.random.rand(mesh.nF, 1), i) f1.set_e(np.random.rand(mesh.nE, 1), i) - f2 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b') - for i in range(f2.nTimes): + f2 = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') + for i in range(f2.nT): f2.set_b(np.random.rand(mesh.nF, 1), i) f2.set_e(np.random.rand(mesh.nE, 1), i) @@ -262,13 +262,13 @@ class TDEM_bDerivTests(unittest.TestCase): m = np.random.rand(prb.mapping.nP) sigma = np.random.rand(prb.mapping.nP) - u = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b') - for i in range(u.nTimes): + u = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') + for i in range(u.nT): u.set_b(np.random.rand(mesh.nF, 1), i) u.set_e(np.random.rand(mesh.nE, 1), i) - v = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nTimes, 'b') - for i in range(v.nTimes): + v = EM.TDEM.FieldsTDEM(prb.mesh, 1, prb.nT, 'b') + for i in range(v.nT): v.set_b(np.random.rand(mesh.nF, 1), i) v.set_e(np.random.rand(mesh.nE, 1), i) @@ -282,7 +282,7 @@ class TDEM_bDerivTests(unittest.TestCase): sigma = self.sigma m = np.random.rand(prb.mapping.nP) - d = np.random.rand(prb.nTimes) + d = np.random.rand(prb.nT) V1 = d.dot(prb.Jvec(sigma, m)) V2 = m.dot(prb.Jtvec(sigma, d)) diff --git a/simpegEM/Tests/test_TDEM_forward_Analytic.py b/simpegEM/Tests/test_TDEM_forward_Analytic.py index 63e81a5e..ceab2249 100644 --- a/simpegEM/Tests/test_TDEM_forward_Analytic.py +++ b/simpegEM/Tests/test_TDEM_forward_Analytic.py @@ -39,23 +39,23 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5, # except ImportError, e: # pass - prb.setTimes([1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4], [40, 40, 40, 40, 40, 40]) + prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40), (0.0001, 40), (0.0005, 40)] sigma = np.ones(mesh.nCz)*1e-8 sigma[active] = sig_half sigma = np.log(sigma[active]) prb.pair(survey) - bz_ana = mu_0*EM.Utils.Ana.hzAnalyticDipoleT(survey.rxLoc[0], prb.times, sig_half) + bz_ana = mu_0*EM.Utils.Ana.hzAnalyticDipoleT(survey.rxLoc[0]+1e-3, prb.times[1:], sig_half) bz_calc = survey.dpred(sigma) - ind = np.logical_and(prb.times > bounds[0],prb.times < bounds[1]) + ind = np.logical_and(prb.times[1:] > bounds[0],prb.times[1:] < bounds[1]) log10diff = np.linalg.norm(np.log10(np.abs(bz_calc[ind])) - np.log10(np.abs(bz_ana[ind])))/np.linalg.norm(np.log10(np.abs(bz_ana[ind]))) print 'Difference: ', log10diff if showIt == True: - plt.loglog(prb.times[bz_calc>0], bz_calc[bz_calc>0], 'r', prb.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--') - plt.loglog(prb.times, abs(bz_ana), 'b*') + plt.loglog(prb.times[1:][bz_calc>0], bz_calc[bz_calc>0], 'r', prb.times[1:][bz_calc<0], -bz_calc[bz_calc<0], 'r--') + plt.loglog(prb.times[1:], abs(bz_ana), 'b*') plt.title('sig_half = %e'%sig_half) plt.show()