From 93778d13f10a04cce2c4b5d1079fb93427433802 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Fri, 17 Apr 2015 10:53:25 -0700 Subject: [PATCH] ProblemFDEM_x.getSource now returns two matrices, S_m and S_e, which cleans up getRHS and the Fields --- simpegEM/FDEM/FDEM.py | 73 +++++++++++------------------------ simpegEM/FDEM/FieldsFDEM.py | 76 +++++++++++++++++-------------------- simpegEM/FDEM/SurveyFDEM.py | 3 -- simpegEM/Tests/test_FDEM.py | 6 +-- 4 files changed, 60 insertions(+), 98 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index fe89951f..870558b6 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -109,13 +109,21 @@ class BaseFDEMProblem(BaseEMProblem): :return: RHS """ Txs = self.survey.getTransmitters(freq) - j_m = range(len(Txs)) - j_e = range(len(Txs)) - for i, tx in enumerate(Txs): - j_m[i], j_e[i] = tx.getSource(self) + if self._eqLocs is 'FE': + S_m = 1j*np.zeros((self.mesh.nF,len(Txs))) + S_e = 1j*np.zeros((self.mesh.nE,len(Txs))) + elif self._eqLocs is 'EF': + S_m = 1j*np.zeros((self.mesh.nE,len(Txs))) + S_e = 1j*np.zeros((self.mesh.nF,len(Txs))) - return j_m, j_e - # return np.concatenate(rhs).reshape((-1, len(Txs)), order='F') #, np.concatenate(j_e).reshape((-1, len(Txs)), order='F') + for i, tx in enumerate(Txs): + smi, sei = tx.getSource(self) + if smi is not None: + S_m[:,i] = smi + if sei is not None: + S_e[:,i] = sei + + return S_m, S_e def getSourceDeriv(self,freq,adjoint=False): return None, None @@ -181,20 +189,11 @@ class ProblemFDEM_e(BaseFDEMProblem): :return: RHS """ - j_m, j_g = self.getSource(freq) - nTx_freq = self.survey.nTxByFreq[freq] - RHS = 1j*np.zeros([self.mesh.nE, nTx_freq]) - + S_m, S_e = self.getSource(freq) C = self.mesh.edgeCurl MfMui = self.MfMui - for ii in range(nTx_freq): - if j_m[ii] is not None: - - RHS[:, ii] += C.T * (MfMui * j_m[ii]) - - if j_g[ii] is not None: - RHS[:, ii] += -1j*omega(freq)*j_g[ii] + RHS = C.T * (MfMui * S_m) -1j*omega(freq)*S_e return RHS @@ -256,20 +255,11 @@ class ProblemFDEM_b(BaseFDEMProblem): :return: RHS """ - j_m, j_g = self.getSource(freq) - nTx_freq = self.survey.nTxByFreq[freq] - RHS = 1j*np.zeros([self.mesh.nF, nTx_freq]) - + S_m, S_e = self.getSource(freq) C = self.mesh.edgeCurl - MfSigmai = self.MfSigmai - - for ii in range(nTx_freq): - if j_m[ii] is not None: - RHS[:,ii] += j_m[ii] - - if j_g[ii] is not None: - RHS[:,ii] += C * ( MfSigmai * j_g[ii] ) + MeSigmaI = self.MeSigmaI + RHS = S_m + C * ( MeSigmaI * S_e ) if self._makeASymmetric is True: mui = self.MfMui @@ -373,21 +363,12 @@ class ProblemFDEM_j(BaseFDEMProblem): :return: RHS """ - j_m, j_g = self.getSource(freq) - nTx_freq = self.survey.nTxByFreq[freq] - RHS = 1j*np.zeros([self.mesh.nF, nTx_freq]) - + S_m, S_e = self.getSource(freq) C = self.mesh.edgeCurl MeMuI = self.MeMuI - for ii in range(nTx_freq): - if j_m[ii] is not None: - RHS[:,ii] += C * (MeMuI * j_m[ii]) - - if j_g[ii] is not None: - RHS[:,ii] += -1j * omega(freq) * j_g[ii] - + RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e if self._makeASymmetric is True: MfSigi = self.MfSigmai return MfSigi.T*RHS @@ -465,19 +446,11 @@ class ProblemFDEM_h(BaseFDEMProblem): :return: RHS """ - j_m, j_g = self.getSource(freq) - nTx_freq = self.survey.nTxByFreq[freq] - RHS = 1j*np.zeros([self.mesh.nE, nTx_freq]) - + S_m, S_e = self.getSource(freq) C = self.mesh.edgeCurl MfSigmai = self.MfSigmai - for ii in range(nTx_freq): - if j_m[ii] is not None: - RHS[:,ii] += j_m[ii] - - if j_g[ii] is not None: - RHS[:,ii] += C.T * ( MfSigmai * j_g[ii] ) + RHS = S_m + C.T * ( MfSigmai * S_e ) return RHS diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index 0f8c8596..bda8accf 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -4,9 +4,9 @@ from simpegEM.Utils.EMUtils import omega class FieldsFDEM(Problem.Fields): """Fancy Field Storage for a FDEM survey.""" - knownFields = None + knownFields = {} dtype = complex - + class FieldsFDEM_e(FieldsFDEM): knownFields = {'e':'E'} @@ -30,23 +30,21 @@ class FieldsFDEM_e(FieldsFDEM): return None def _b(self, e, tx): #adjoint=False - b = self._b_sec(e,tx) - j_m,_ = self.getSource(tx.freq) - if j_m[0] is not None: - b += 1./(1j*omega(tx.freq)) * np.array([j_m[0]]).T - return b + b_sec = self._b_sec(e,tx) + S_m,_ = self.getSource(tx.freq) + return b_sec + 1./(1j*omega(tx.freq)) * S_m def _bDeriv(self, e, tx, adjoint=False): - j_mDeriv,_ = self.getSourceDeriv(tx.freq, adjoint) + S_mDeriv,_ = self.getSourceDeriv(tx.freq, adjoint) b_secDeriv = self._b_secDeriv(e,tx.freq,adjoint) - if j_mDeriv is None & b_secDeriv is None: + if S_mDeriv is None & b_secDeriv is None: return None elif b_secDeriv is None: - return 1./(1j*omega(tx.freq)) * j_mDeriv - elif j_mDeriv is None: + return 1./(1j*omega(tx.freq)) * S_mDeriv + elif S_mDeriv is None: return b_secDeriv else: - return 1./(1j*omega(tx.freq)) * j_mDeriv + b_secDeriv + return 1./(1j*omega(tx.freq)) * S_mDeriv + b_secDeriv class FieldsFDEM_b(FieldsFDEM): @@ -73,24 +71,22 @@ class FieldsFDEM_b(FieldsFDEM): return None def _e(self, b, tx): - e = self._e_sec(b,tx) - _, j_g = self.getSource(tx.freq) - if j_g[0] is not None: - e += -np.array([j_g[0]]).T - return e + e_sec = self._e_sec(b,tx) + _, S_e = self.getSource(tx.freq) + return e_sec + S_e def _eDeriv(self, b, tx, adjoint=False): - _,j_gDeriv = self.getSourceDeriv(tx.freq, adjoint) + _,S_eDeriv = self.getSourceDeriv(tx.freq, adjoint) e_secDeriv = self._e_secDeriv(b,tx,adjoint) - if j_gDeriv is None & e_secDeriv is None: + if S_eDeriv is None & e_secDeriv is None: return None elif e_secDeriv is None: - return -j_gDeriv - elif j_gDeriv is None: + return -S_eDeriv + elif S_eDeriv is None: return e_secDeriv else: - return e_secDeriv - j_gDeriv + return e_secDeriv - S_eDeriv class FieldsFDEM_j(FieldsFDEM): @@ -129,23 +125,21 @@ class FieldsFDEM_j(FieldsFDEM): raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) def _h(self, j, tx): #adjoint=False - h = self._h_sec(j,tx) - j_m,_ = self.getSource(tx.freq) - if j_m[0] is not None: - h += 1./(1j*omega(tx.freq)) * self.MeMuI * np.array([j_m[0]]).T - return h + h_sec = self._h_sec(j,tx) + S_m,_ = self.getSource(tx.freq) + return h_sec + 1./(1j*omega(tx.freq)) * self.MeMuI * S_m def _hDeriv(self, j, tx, adjoint=False): - j_mDeriv,_ = self.getSourceDeriv(tx.freq, adjoint) + S_mDeriv,_ = self.getSourceDeriv(tx.freq, adjoint) h_secDeriv = self._h_secDeriv(j,tx.freq,adjoint) - if j_mDeriv is None & h_secDeriv is None: + if S_mDeriv is None & h_secDeriv is None: return None elif h_secDeriv is None: - return 1./(1j*omega(tx.freq)) * j_mDeriv - elif j_mDeriv is None: + return 1./(1j*omega(tx.freq)) * S_mDeriv + elif S_mDeriv is None: return h_secDeriv else: - return 1./(1j*omega(tx.freq)) * j_mDeriv + h_secDeriv + return 1./(1j*omega(tx.freq)) * S_mDeriv + h_secDeriv class FieldsFDEM_h(FieldsFDEM): knownFields = {'h':'E'} @@ -171,23 +165,21 @@ class FieldsFDEM_h(FieldsFDEM): return None def _j(self, h, tx): #adjoint=False - j = self._j_sec(h,tx) - _,j_g = self.getSource(tx.freq) - if j_g[0] is not None: - j += -np.array([j_g[0]]).T - return j + j_sec = self._j_sec(h,tx) + _,S_e = self.getSource(tx.freq) + return j_sec - S_e def _jDeriv(self, h, tx, adjoint=False): - _,j_gDeriv = self.getSourceDeriv(tx.freq, adjoint) + _,S_eDeriv = self.getSourceDeriv(tx.freq, adjoint) j_secDeriv = self._j_secDeriv(j,tx.freq,adjoint) - if j_gDeriv is None & j_secDeriv is None: + if S_eDeriv is None & j_secDeriv is None: return None elif j_secDeriv is None: - return - j_gDeriv - elif j_gDeriv is None: + return - S_eDeriv + elif S_eDeriv is None: return j_secDeriv else: - return - j_gDeriv + j_secDeriv + return - S_eDeriv + j_secDeriv # def calcFields(self, sol, freq, fieldType, adjoint=False): diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index e7553c43..b9aa03a0 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -182,10 +182,7 @@ class TxFDEM(Survey.BaseTx): a = SRC b_0 = C*a - # if solType == 'b' or solType == 'h': return -1j*omega(freq)*b_0, None - # elif solType == 'e' or solType == 'j': - # return -1j*omega(freq)*C.T*mui*b_0, None class SimpleTxFDEM_g(TxFDEM): diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index dc72f02c..9171caec 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -396,7 +396,7 @@ class FDEM_DerivTests(unittest.TestCase): def test_EB_CrossCheck_byr_Eform(self): self.assertTrue(crossCheckTest('e', 'byr')) def test_EB_CrossCheck_bzr_Eform(self): - self.assertTrue(crossCheckTest('e', 'bzr')) # Doesn't make sense to test this for p-s approach + self.assertTrue(crossCheckTest('e', 'bzr')) def test_EB_CrossCheck_bxi_Eform(self): self.assertTrue(crossCheckTest('e', 'bxi')) def test_EB_CrossCheck_byi_Eform(self): @@ -422,8 +422,8 @@ class FDEM_DerivTests(unittest.TestCase): self.assertTrue(crossCheckTest('j', 'hxr')) def test_HJ_CrossCheck_hyr_Jform(self): self.assertTrue(crossCheckTest('j', 'hyr')) - # def test_HJ_CrossCheck_hzr_Jform(self): - # self.assertTrue(crossCheckTest('j', 'hzr')) # Doesn't make sense to test this for p-s approach + def test_HJ_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest('j', 'hzr')) def test_HJ_CrossCheck_hxi_Jform(self): self.assertTrue(crossCheckTest('j', 'hxi')) def test_HJ_CrossCheck_hyi_Jform(self):