From d5967d20b98d5be3c89b03d5c39bf7990851c916 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 22 Feb 2016 17:42:22 -0800 Subject: [PATCH] forward is running, but not passing --- SimPEG/EM/TDEM/FieldsTDEM.py | 10 +- SimPEG/EM/TDEM/SurveyTDEM.py | 119 +++++++++++++++++--- SimPEG/EM/TDEM/TDEM.py | 64 +++++++++-- tests/em/tdem/test_TDEM_forward_Analytic.py | 54 ++++----- 4 files changed, 193 insertions(+), 54 deletions(-) diff --git a/SimPEG/EM/TDEM/FieldsTDEM.py b/SimPEG/EM/TDEM/FieldsTDEM.py index 83197cd6..79601850 100644 --- a/SimPEG/EM/TDEM/FieldsTDEM.py +++ b/SimPEG/EM/TDEM/FieldsTDEM.py @@ -5,7 +5,7 @@ from SimPEG import Utils from SimPEG.EM.Utils import omega from SimPEG.Utils import Zero, Identity -class Fields(SimPEG.Problem.Fields): +class Fields(SimPEG.Problem.TimeFields): """ Fancy Field Storage for a TDEM survey. Only one field type is stored for @@ -32,8 +32,6 @@ class Fields(SimPEG.Problem.Fields): dtype = float - - class Fields_b(Fields): """Fancy Field Storage for a TDEM survey.""" knownFields = {'bSolution': 'F'} @@ -57,7 +55,11 @@ class Fields_b(Fields): return Zero() def _e(self, bSolution, srcList, tInd): - return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) ) + e = self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) ) + for i, src in enumerate(srcList): + _, S_e = src.eval(self.prob, tInd) + e[:,i,tInd] = e[:,i,tInd] - self.MeSigmaI * S_e + return e def _eDeriv_u(self, src, du_dm_v, adjoint = False): raise NotImplementedError diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index b25d0888..05dae215 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/SimPEG/EM/TDEM/SurveyTDEM.py @@ -1,13 +1,14 @@ -from SimPEG import Survey, np -from SimPEG.Survey import BaseSurvey +import SimPEG +from SimPEG import np, Utils from SimPEG.Utils import Zero, Identity from scipy.constants import mu_0 +from SimPEG.EM.Utils import * #################################################### # Receivers #################################################### -class Rx(Survey.BaseTimeRx): +class Rx(SimPEG.Survey.BaseTimeRx): knownRxTypes = { 'ex':['e', 'Ex', 'N'], @@ -24,7 +25,7 @@ class Rx(Survey.BaseTimeRx): } def __init__(self, locs, times, rxType): - Survey.BaseTimeRx.__init__(self, locs, times, rxType) + SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType) @property def projField(self): @@ -77,6 +78,10 @@ class BaseWaveform(object): self.offTime = offTime self.hasInitialFields = hasInitialFields + def _assertMatchesPair(self, pair): + assert (isinstance(self, pair) + ), "Waveform object must be an instance of a %s BaseWaveform class."%(pair.__name__) + def eval(self, time): raise NotImplementedError @@ -91,11 +96,11 @@ class StepOffWaveform(BaseWaveform): -class BaseSrc(Survey.BaseSrc): +class BaseSrc(SimPEG.Survey.BaseSrc): rxPair = Rx integrate = True - waveformPair = StepOffWaveform + waveformPair = BaseWaveform @property def waveform(self): @@ -103,17 +108,17 @@ class BaseSrc(Survey.BaseSrc): return getattr(self, '_waveform', None) @waveform.setter def waveform(self, val): - if self.waveform is None: - val._assertMatchesPair(self.waveformPair) - self._waveform = val + # if self.waveform is None: + val._assertMatchesPair(self.waveformPair) + self._waveform = val def __init__(self, rxList, waveform = None): self.waveform = waveform - Survey.BaseSrc.__init__(self, rxList) + SimPEG.Survey.BaseSrc.__init__(self, rxList) - def bInitial(self, mesh): + def bInitial(self): return Zero() def eval(self, prob, time): @@ -135,7 +140,7 @@ class BaseSrc(Survey.BaseSrc): class MagDipole(BaseSrc): - def __init__(self, rxList, waveform, loc, orientation='Z', moment=1., mu=mu_0, **kwargs): + def __init__(self, rxList, waveform=None, loc=None, orientation='Z', moment=1., mu=mu_0): self.loc = loc self.orientation = orientation @@ -143,10 +148,98 @@ class MagDipole(BaseSrc): self.moment = moment self.mu = mu self.integrate = False - BaseSrc.__init__(self, rxList) + BaseSrc.__init__(self, rxList, waveform) + + def _bfromVectorPotential(self, prob): + if prob._eqLocs is 'FE': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl + + elif prob._eqLocs is 'EF': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl.T + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + + else: + srcfct = MagneticDipoleVectorPotential + ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + a = np.concatenate((ax, ay, az)) + + return C*a + def bInitial(self, prob): + eqLocs = prob._eqLocs + + if self.waveform.hasInitialFields is False: + return Zero() + + return self._bfromVectorPotential(prob) + + def S_m(self, prob, time): + if self.waveform.hasInitialFields is False: + raise NotImplementedError + return Zero() + + def S_e(self, prob, time): + if self.waveform.hasInitialFields is False: + raise NotImplementedError + return Zero() + +#################################################### +# Survey +#################################################### + +class Survey(SimPEG.Survey.BaseSurvey): + """ + Time domain electromagnetic survey + """ + + srcPair = BaseSrc + rxPair = Rx + + def __init__(self, srcList, **kwargs): + # Sort these by frequency + self.srcList = srcList + SimPEG.Survey.BaseSurvey.__init__(self, **kwargs) + + def eval(self, u): + data = SimPEG.Survey.Data(self) + for src in self.srcList: + for rx in src.rxList: + data[src, rx] = rx.eval(src, self.mesh, self.prob.timeMesh, u) + return data + + def evalDeriv(self, u, v=None, adjoint=False): + assert v is not None, 'v to multiply must be provided.' + + if not adjoint: + data = SimPEG.Survey.Data(self) + for src in self.srcList: + for rx in src.rxList: + data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v) + return data + else: + f = FieldsTDEM(self.mesh, self) + for src in self.srcList: + for rx in src.rxList: + Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True) + Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F') + if rx.projField not in f: # first time we are projecting + f[src, rx.projField, :] = Ptv + else: # there are already fields, so let's add to them! + f[src, rx.projField, :] += Ptv + return f diff --git a/SimPEG/EM/TDEM/TDEM.py b/SimPEG/EM/TDEM/TDEM.py index 6fb8398e..e6d80009 100644 --- a/SimPEG/EM/TDEM/TDEM.py +++ b/SimPEG/EM/TDEM/TDEM.py @@ -29,20 +29,59 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): tic = time.time() self.curModel = m + F = self.fieldsPair(self.mesh, self.survey) - # for + # set initial fields + for i, src in enumerate(self.survey.srcList): + F[src,'bSolution',0] = src.bInitial(self) + # timestep to solve forward + Ainv = None + for tInd, dt in enumerate(self.timeSteps): + print dt, self.timeSteps[tInd] + if Ainv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same + Ainv.clean() + + if Ainv is None: + A = self.getA(tInd) + if self.verbose: print 'Factoring... (dt = %e)'%dt + Ainv = self.Solver(A, **self.solverOpts) + if self.verbose: print 'Done' + + rhs = self.getRHS(tInd, F) + if self.verbose: print ' Solving... (tInd = %d)'%tInd + sol = Ainv * rhs + if self.verbose: print ' Done...' + if sol.ndim == 1: + sol.shape = (sol.size,1) + F[:,self._fieldType+'Solution',tInd+1] = sol + Ainv.clean() + return F def Jvec(self, m, v, u=None): - return None + raise NotImplementedError def Jtvec(self, m, v, u=None): - return None + raise NotImplementedError def getSourceTerm(self, tInd): - return None - # return S_m, S_e + + Srcs = self.survey.srcList + + if self._eqLocs is 'FE': + S_m = np.zeros((self.mesh.nF,len(Srcs))) + S_e = np.zeros((self.mesh.nE,len(Srcs))) + elif self._eqLocs is 'EF': + S_m = np.zeros((self.mesh.nE,len(Srcs))) + S_e = np.zeros((self.mesh.nF,len(Srcs))) + + for i, src in enumerate(Srcs): + smi, sei = src.eval(self, self.times[tInd]) + S_m[:,i] = S_m[:,i] + smi + S_e[:,i] = S_e[:,i] + sei + + return S_m, S_e @@ -114,7 +153,7 @@ class Problem_b(BaseTDEMProblem): A = I + dt * ( C * ( MeSigmaI * (C.T * MfMui ) ) ) if self._makeASymmetric is True: - return MeMui.T * A + return MfMui.T * A return A def getADeriv(self, freq, u, v, adjoint=False): @@ -135,21 +174,24 @@ class Problem_b(BaseTDEMProblem): return ADeriv - def getRHS(self, tInd): + def getRHS(self, tInd, F): dt = self.timeSteps[tInd] C = self.mesh.edgeCurl MeSigmaI = self.MeSigmaI MfMui = self.MfMui - S_m, S_e = self.getSourceTerm(tInd+1) # I think this is tInd+1 ? + S_m, S_e = self.getSourceTerm(tInd+1) - B_n = np.c_[[F[src,'b',tInd] for src in self.survey.srcList]].T + B_n = np.c_[[F[src,'bSolution',tInd] for src in self.survey.srcList]].T if B_n.shape[0] is not 1: raise NotImplementedError('getRHS not implemented for this shape of B_n') - return B_n + dt * (C * (MeSigmaIDeriv * S_e) + S_m) + rhs = B_n[:,:,0].T + dt * (C * (MeSigmaI * S_e) + S_m) + if self._makeASymmetric: + return MfMui.T * rhs + return rhs - def getRHSDeriv(self, tInd, src, v, adjoint=False): + def getRHSDeriv(self, tInd, F, src, v, adjoint=False): raise NotImplementedError diff --git a/tests/em/tdem/test_TDEM_forward_Analytic.py b/tests/em/tdem/test_TDEM_forward_Analytic.py index 3d85140b..e90dbd3d 100644 --- a/tests/em/tdem/test_TDEM_forward_Analytic.py +++ b/tests/em/tdem/test_TDEM_forward_Analytic.py @@ -10,7 +10,7 @@ except ImportError, e: MumpsSolver = SolverLU -def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,1e-3], showIt=False): +def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,1e-3], showIt=True): if meshType == 'CYL': cs, ncx, ncz, npad = 5., 30, 10, 15 hx = [(cs,ncx), (cs,npad,1.3)] @@ -28,11 +28,11 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5, mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz') - src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.])) + src = EM.TDEM.SurveyTDEM.MagDipole([rx], waveform= EM.TDEM.SurveyTDEM.StepOffWaveform(), loc=np.array([0., 0., 0.])) # src = EM.TDEM.SrcTDEM([rx], loc=np.array([0., 0., 0.])) - survey = EM.TDEM.SurveyTDEM([src]) - prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) + survey = EM.TDEM.Survey([src]) + prb = EM.TDEM.Problem_b(mesh, mapping=mapping) prb.Solver = MumpsSolver prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40), (0.0001, 40), (0.0005, 40)] @@ -48,6 +48,8 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5, ind = np.logical_and(rx.times > bounds[0],rx.times < 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 ' |bz_ana| = ',np.linalg.norm(bz_ana), ' |bz_num| = ', np.linalg.norm(bz_calc), ' |bz_ana - bz_num| =', np.linalg.norm(bz_ana-bz_calc) print 'Difference: ', log10diff if showIt == True: @@ -61,35 +63,35 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5, class TDEM_SimpleSrcTests(unittest.TestCase): def test_source(self): - waveform = EM.TDEM.SurveyTDEM.StepOffWaveform([]) + waveform = EM.TDEM.SurveyTDEM.StepOffWaveform() assert waveform.eval(0.) == 0. -# class TDEM_bTests(unittest.TestCase): +class TDEM_bTests(unittest.TestCase): -# def test_analytic_p2_CYL_50m(self): -# self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+2) < 0.01) -# def test_analytic_p1_CYL_50m(self): -# self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+1) < 0.01) -# def test_analytic_p0_CYL_50m(self): -# self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+0) < 0.01) -# def test_analytic_m1_CYL_50m(self): -# self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-1) < 0.01) -# def test_analytic_m2_CYL_50m(self): -# self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-2) < 0.01) -# def test_analytic_m3_CYL_50m(self): -# self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-3) < 0.02) + def test_analytic_p2_CYL_50m(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+2) < 0.01) + def test_analytic_p1_CYL_50m(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+1) < 0.01) + def test_analytic_p0_CYL_50m(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+0) < 0.01) + def test_analytic_m1_CYL_50m(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-1) < 0.01) + def test_analytic_m2_CYL_50m(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-2) < 0.01) + def test_analytic_m3_CYL_50m(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-3) < 0.02) -# def test_analytic_p0_CYL_1m(self): -# self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e+0) < 0.01) -# def test_analytic_m1_CYL_1m(self): -# self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-1) < 0.01) -# def test_analytic_m2_CYL_1m(self): -# self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-2) < 0.01) -# def test_analytic_m3_CYL_1m(self): -# self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-3) < 0.02) + def test_analytic_p0_CYL_1m(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e+0) < 0.01) + def test_analytic_m1_CYL_1m(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-1) < 0.01) + def test_analytic_m2_CYL_1m(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-2) < 0.01) + def test_analytic_m3_CYL_1m(self): + self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-3) < 0.02)