From 0be942730a7408e81223b5ca1f021002bcbdc123 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 18 Mar 2016 15:47:31 -0700 Subject: [PATCH] start of Problem_e (e Jvec working) --- SimPEG/EM/TDEM/FieldsTDEM.py | 59 ++++++++++++++--- SimPEG/EM/TDEM/SurveyTDEM.py | 37 ++++++++++- SimPEG/EM/TDEM/TDEM.py | 77 +++++++++++++++++++++-- SimPEG/EM/TDEM/__init__.py | 4 +- tests/em/tdem/test_TDEM_b_DerivAdjoint.py | 42 +++++++------ tests/em/tdem/test_TDEM_crosscheck.py | 76 ++++++++++++++++++++++ 6 files changed, 261 insertions(+), 34 deletions(-) create mode 100644 tests/em/tdem/test_TDEM_crosscheck.py diff --git a/SimPEG/EM/TDEM/FieldsTDEM.py b/SimPEG/EM/TDEM/FieldsTDEM.py index 033dd68f..611b929f 100644 --- a/SimPEG/EM/TDEM/FieldsTDEM.py +++ b/SimPEG/EM/TDEM/FieldsTDEM.py @@ -31,7 +31,15 @@ class Fields(SimPEG.Problem.TimeFields): knownFields = {} dtype = float + def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): + if adjoint is True: + return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint) + return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v) + def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): + if adjoint is True: + return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint) + return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v) class Fields_Derivs(Fields): knownFields = { @@ -65,10 +73,10 @@ class Fields_b(Fields): def _bDeriv_m(self, tInd, src, v, adjoint=False): return Zero() - def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): - if adjoint is True: - return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint) - return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v) + # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): + # if adjoint is True: + # return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint) + # return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v) def _e(self, bSolution, srcList, tInd): e = self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) ) @@ -93,7 +101,42 @@ class Fields_b(Fields): return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution)) * v - self.MeSigmaI * S_eDeriv(v) - def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): - if adjoint is True: - return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint) - return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v) + + +class Fields_e(Fields): + """Fancy Field Storage for a TDEM survey.""" + knownFields = {'eSolution': 'E'} + aliasFields = { + 'e': ['eSolution', 'E', '_e'], + 'b': ['eSolution', 'F', '_b'], + } + + def startup(self): + self.MeSigmaI = self.survey.prob.MeSigmaI + self.MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv + self.edgeCurl = self.survey.prob.mesh.edgeCurl + self.MfMui = self.survey.prob.MfMui + + + def _e(self, eSolution, srcList, tInd): + return eSolution + + def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False): + return dun_dm_v + + def _eDeriv_m(self, tInd, src, v, adjoint = False): + return Zero() + + def _b(self, eSolution, srcList, tInd): + raise NotImplementedError + + def _bDeriv_u(self, tInd, src, dun_dm_v, adjoint=False): + raise NotImplementedError + + def _bDeriv_m(self, tInd, src, v, adjoint=False): + raise NotImplementedError + + # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): + # if adjoint is True: + # return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint) + # return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v) diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index acd154ba..2e8d3ce7 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/SimPEG/EM/TDEM/SurveyTDEM.py @@ -144,13 +144,18 @@ class BaseSrc(SimPEG.Survey.BaseSrc): SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs) - def bInitial(self, prob): return Zero() def bInitialDeriv(self, prob, v=None, adjoint=False): return Zero() + def eInitial(self, prob): + return Zero() + + def eInitialDeriv(self, prob, v=None, adjoint=False): + return Zero() + def eval(self, prob, time): S_m = self.S_m(prob, time) S_e = self.S_e(prob, time) @@ -218,13 +223,41 @@ class MagDipole(BaseSrc): def bInitial(self, prob): - eqLocs = prob._eqLocs if self.waveform.hasInitialFields is False: return Zero() return self._bfromVectorPotential(prob) + def eInitial(self, prob): + + if self.waveform.hasInitialFields is False: + return Zero() + + b = self.bInitial(prob) + MeSigmaI = prob.MeSigmaI + MfMui = prob.MfMui + C = prob.mesh.edgeCurl + + return MeSigmaI * (C.T * (MfMui * b)) + + def eInitialDeriv(self, prob, v=None, adjoint=False): + + if self.waveform.hasInitialFields is False: + return Zero() + + b = self.bInitial(prob) + MeSigmaIDeriv = prob.MeSigmaIDeriv + MfMui = prob.MfMui + C = prob.mesh.edgeCurl + S_e = self.S_e(prob, prob.t0) + + if adjoint: + raise NotImplementedError + + return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ) * v + + def S_m(self, prob, time): if self.waveform.hasInitialFields is False: raise NotImplementedError diff --git a/SimPEG/EM/TDEM/TDEM.py b/SimPEG/EM/TDEM/TDEM.py index b087231e..bc7c6b4c 100644 --- a/SimPEG/EM/TDEM/TDEM.py +++ b/SimPEG/EM/TDEM/TDEM.py @@ -116,9 +116,13 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): dA_dm_v = self.getAdiagDeriv(tInd, un_src, v) # cell centered on time mesh dRHS_dm_v = self.getRHSDeriv(tInd+1, src, v) # on nodes of time mesh - # dAsubdiag_dm_v = 0 - JRHS = dRHS_dm_v - dA_dm_v # - dAsubdiag_dm_v (which is zero) + # if tInd >= 1: + dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, u[src,ftype,tInd], v) + # else: + # dAsubdiag_dm_v = Zero() + + JRHS = dRHS_dm_v - dAsubdiag_dm_v - dA_dm_v # step in time and overwrite if tInd != len(self.timeSteps+1): @@ -180,8 +184,6 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): df_duT_v[src, '%sDeriv'%self._fieldType, tInd] = df_duT_v[src, '%sDeriv'%self._fieldType, tInd] + Utils.mkvc(cur[0],2) JTv = JTv + cur[1] - - del PT_v # no longer need this @@ -309,6 +311,8 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem): ################################ E-B Formulation ######################################### ########################################################################################## +# ------------------------------- Problem_b -------------------------------------------- # + class Problem_b(BaseTDEMProblem): """ Starting from the quasi-static E-B formulation of Maxwell's equations (semi-discretized) @@ -404,6 +408,9 @@ class Problem_b(BaseTDEMProblem): return Asubdiag + def getAsubdiagDeriv(self, tInd, u, v, adjoint=False): + return Zero() * v + def getRHS(self, tInd): @@ -450,6 +457,68 @@ class Problem_b(BaseTDEMProblem): return RHSDeriv +# ------------------------------- Problem_e -------------------------------------------- # + +class Problem_e(BaseTDEMProblem): + + _fieldType = 'e' + _eqLocs = 'FE' + fieldsPair = Fields_e + surveyPair = SurveyTDEM + + def __init__(self, mesh, mapping=None, **kwargs): + BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs) + + def getAdiag(self, tInd): + """ + System matrix at a given time index + + """ + assert tInd >= 0 and tInd < self.nT + + dt = self.timeSteps[tInd] + C = self.mesh.edgeCurl + MfMui = self.MfMui + MeSigma = self.MeSigma + + return C.T * ( MfMui * C ) + 1./dt * MeSigma + + + def getAdiagDeriv(self, tInd, u, v, adjoint=False): + + dt = self.timeSteps[tInd] + C = self.mesh.edgeCurl + MfMui = self.MfMui + MeSigmaDeriv = self.MeSigmaDeriv(u) + + if adjoint: + raise 1./dt * MeSigmaDeriv.T * v + + return 1./dt * MeSigmaDeriv * v + + + def getAsubdiag(self, tInd): + assert tInd >= 0 and tInd < self.nT + + dt = self.timeSteps[tInd] + + return - 1./dt * self.MeSigma + + def getAsubdiagDeriv(self, tInd, u, v, adjoint=False): + dt = self.timeSteps[tInd] + + if adjoint: + return - 1./dt * self.MeSigmaDeriv(u).T * v + + return - 1./dt * self.MeSigmaDeriv(u) * v + + def getRHS(self, tInd): + return Zero() + + def getRHSDeriv(self, tInd, src, v, adjoint=False): + return Zero() + + diff --git a/SimPEG/EM/TDEM/__init__.py b/SimPEG/EM/TDEM/__init__.py index be3b9d8a..13010109 100644 --- a/SimPEG/EM/TDEM/__init__.py +++ b/SimPEG/EM/TDEM/__init__.py @@ -1,3 +1,3 @@ -from TDEM import BaseTDEMProblem, Problem_b +from TDEM import BaseTDEMProblem, Problem_b, Problem_e from FieldsTDEM import Fields, Fields_b -from SurveyTDEM import Survey, BaseSrc, Rx \ No newline at end of file +from SurveyTDEM import Survey, BaseSrc, Rx diff --git a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py index c5d9ec5d..a297e6cb 100644 --- a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py @@ -5,11 +5,11 @@ from SimPEG import EM plotIt = False testDeriv = True -testAdjoint = True +testAdjoint = False TOL = 1e-5 -def setUp(rxcomp='bz'): +def setUp(prbtype='b', rxcomp='bz'): cs = 5. ncx = 20 ncy = 15 @@ -28,7 +28,11 @@ def setUp(rxcomp='bz'): survey = EM.TDEM.Survey([src]) - prb = EM.TDEM.Problem_b(mesh, mapping=mapping) + if prbtype == 'b': + prb = EM.TDEM.Problem_b(mesh, mapping=mapping) + elif prbtype == 'e': + prb = EM.TDEM.Problem_e(mesh, mapping=mapping) + prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] # prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)] @@ -129,40 +133,42 @@ class TDEM_DerivTests(unittest.TestCase): if testDeriv: - def JvecTest(self, rxcomp): - prb, m, mesh = setUp(rxcomp) + def JvecTest(self, prbtype, rxcomp): + prb, m, mesh = setUp(prbtype, rxcomp) derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)] print '\n' - print 'test_Jvec_%s' %(rxcomp) + print 'test_Jvec_%s_%s' %(prbtype, rxcomp) Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20) def test_Jvec_b_bx(self): - self.JvecTest('bx') + self.JvecTest('b','bx') def test_Jvec_b_bz(self): - self.JvecTest('bz') + self.JvecTest('b','bz') def test_Jvec_b_dbxdt(self): - self.JvecTest('dbxdt') + self.JvecTest('b','dbxdt') def test_Jvec_b_dbzdt(self): - self.JvecTest('dbzdt') + self.JvecTest('b','dbzdt') def test_Jvec_b_ey(self): - self.JvecTest('ey') + self.JvecTest('b','ey') + def test_Jvec_e_ey(self): + self.JvecTest('e','ey') # ====== TEST Jtvec ========== # if testAdjoint: - def JvecVsJtvecTest(self, rxcomp='bz'): + def JvecVsJtvecTest(self, prbtype='b', rxcomp='bz'): print '\nAdjoint Testing Jvec, Jtvec %s' %(rxcomp) - prb, m0, mesh = setUp(rxcomp) + prb, m0, mesh = setUp(prbtype, rxcomp) m = np.random.rand(prb.mapping.nP) d = np.random.randn(prb.survey.nD) V1 = d.dot(prb.Jvec(m0, m)) @@ -174,19 +180,19 @@ class TDEM_DerivTests(unittest.TestCase): self.assertTrue(passed) def test_Jvec_adjoint_b_bx(self): - self.JvecVsJtvecTest('bx') + self.JvecVsJtvecTest('b', 'bx') def test_Jvec_adjoint_b_bz(self): - self.JvecVsJtvecTest('bz') + self.JvecVsJtvecTest('b', 'bz') def test_Jvec_adjoint_b_dbxdt(self): - self.JvecVsJtvecTest('bx') + self.JvecVsJtvecTest('b', 'bx') def test_Jvec_adjoint_b_dbzdt(self): - self.JvecVsJtvecTest('bz') + self.JvecVsJtvecTest('b', 'bz') def test_Jvec_adjoint_b_ey(self): - self.JvecVsJtvecTest('ey') + self.JvecVsJtvecTest('b', 'ey') diff --git a/tests/em/tdem/test_TDEM_crosscheck.py b/tests/em/tdem/test_TDEM_crosscheck.py new file mode 100644 index 00000000..ef2b8a83 --- /dev/null +++ b/tests/em/tdem/test_TDEM_crosscheck.py @@ -0,0 +1,76 @@ +import unittest +from SimPEG import * +from SimPEG import EM + +TOL = 1e-5 +FLR = 1e-20 + +np.random.seed(seed=25) # set a seed so that the same conductivity model is used for all runs + +def setUp(prbtype = 'b', rxcomp='bz'): + cs = 5. + ncx = 20 + ncy = 15 + 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.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) + mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap + + rxOffset = 10. + rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,] + src = EM.TDEM.SurveyTDEM.MagDipole([rx], loc=np.array([0., 0., 0.])) + + survey = EM.TDEM.Survey([src]) + + if prbtype == 'b': + prb = EM.TDEM.Problem_b(mesh, mapping=mapping) + elif prbtype == 'e': + prb = EM.TDEM.Problem_e(mesh, mapping=mapping) + + prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] + # prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)] + + try: + from pymatsolver import MumpsSolver + prb.Solver = MumpsSolver + except ImportError, e: + prb.Solver = SolverLU + + m = np.log(1e-1)*np.ones(prb.mapping.nP) #+ 1e-2*np.random.randn(prb.mapping.nP) + + prb.pair(survey) + mesh = mesh + + return prb, m, mesh + +def CrossCheck(prbtype1='b', prbtype2='e', rxcomp='bz'): + + prb1,m1,mesh1 = setUp(prbtype1, rxcomp) + prb2,m2,mesh2 = setUp(prbtype2, rxcomp) + + assert (m1 == m2).all(), 'Models for two formulations are different' + + d1 = prb1.survey.dpred(m1) + d2 = prb2.survey.dpred(m2) + + + check = np.linalg.norm(d1 - d2) + tol = 0.5 * (np.linalg.norm(d1) + np.linalg.norm(d2)) * TOL + passed = check < tol + + print 'Checking %s, %s for %s data'%(prbtype1, prbtype2, rxcomp) + print ' ', np.linalg.norm(d1), np.linalg.norm(d2), np.linalg.norm(check), tol, passed + + assert passed + +class TDEM_cross_check_EB(unittest.TestCase): + def test_EB_ey(self): + CrossCheck('b','e','ey') + +if __name__ == '__main__': + unittest.main() +