From b750faa3e004eea3dcfdd57d37f8b2caf06d4dbe Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 28 Apr 2014 16:11:22 -0700 Subject: [PATCH] Multiple recs on one tx for the same field had a bug. These should add in the projection adjoint. --- simpegEM/TDEM/BaseTDEM.py | 11 +++++------ simpegEM/TDEM/SurveyTDEM.py | 10 ++++++++-- simpegEM/TDEM/TDEM_b.py | 4 ++-- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 5 +++-- 4 files changed, 18 insertions(+), 12 deletions(-) diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index cacf6eaa..8b3c681f 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -16,15 +16,14 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): surveyPair = SurveyTDEM - def calcFields(self, sol, solType, tInd): + def calcFields(self, sol, tInd): - if solType == 'b': + if self.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) + raise NotImplementedError('solType "%s" is not implemented in CalcFields.' % self.solType) return {'b':b, 'e':e} @@ -53,7 +52,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): sol = Asolve.solve(rhs) if sol.ndim == 1: sol.shape = (sol.size,1) - F[:,:,tInd+1] = CalcFields(sol, self.solType, tInd) + F[:,:,tInd+1] = CalcFields(sol, tInd) return F def adjoint(self, m, RHS, CalcFields, F=None): @@ -71,6 +70,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): sol = Asolve.solve(rhs) if sol.ndim == 1: sol.shape = (sol.size,1) - F[:,:,tInd+1] = CalcFields(sol, self.solType, tInd) + F[:,:,tInd+1] = CalcFields(sol, tInd) return F diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 0e04f10a..3df191cd 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -124,8 +124,14 @@ class SurveyTDEM(Survey.BaseSurvey): for tx in self.txList: for rx in tx.rxList: Ptv = rx.projectFieldsDeriv(tx, self.mesh, self.prob.timeMesh, u, v, adjoint=True) - Ptv = Ptv.reshape((-1, 1, self.prob.timeMesh.nN), order='F') - f[tx, rx.projField, :] = Ptv + if rx.projField not in f: # first time we are projecting + Ptv = Ptv.reshape((-1, 1, self.prob.timeMesh.nN), order='F') + f[tx, rx.projField, :] = Ptv + else: + Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F') + addedPtv = f[tx, rx.projField, :] + Ptv + addedPtv = addedPtv.reshape((-1, 1, self.prob.timeMesh.nN), order='F') + f[tx, rx.projField, :] = addedPtv return f diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 9622af14..207baf8a 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -128,7 +128,7 @@ class ProblemTDEM_b(BaseTDEMProblem): dt = self.timeSteps[tInd] return rhs + 1.0/dt*self.MfMui*y[:,'b',tInd] - def AhCalcFields(sol, solType, tInd): + def AhCalcFields(sol, tInd): y_b = sol if self.survey.nTx == 1: y_b = mkvc(y_b) @@ -168,7 +168,7 @@ class ProblemTDEM_b(BaseTDEMProblem): dt = self.timeSteps[tInd+1] return rhs + 1.0/dt*self.MfMui*y[:,'b',tInd+2] - def AhtCalcFields(sol, solType, tInd): + def AhtCalcFields(sol, tInd): y_b = sol if self.survey.nTx == 1: y_b = mkvc(y_b) diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 36176944..071e0de5 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -22,8 +22,9 @@ class TDEM_bDerivTests(unittest.TestCase): [Maps.ExpMap, Maps.Vertical1DMap, activeMap]) rxOffset = 40. - 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]) + 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) survey = EM.TDEM.SurveyTDEM([tx])