From e1010f9ecaca827627841086fe3fc3b002eb7a8e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 10 Mar 2014 22:47:43 -0700 Subject: [PATCH] Initial fields and data implementation. --- simpegEM/FDEM/FDEM.py | 10 ++---- simpegEM/FDEM/SurveyFDEM.py | 50 +++++++++++------------------ simpegEM/Tests/test_FDEM.py | 13 ++++---- simpegEM/Tests/test_FieldsObject.py | 10 ++++-- 4 files changed, 35 insertions(+), 48 deletions(-) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 260d0e06..435990cf 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -93,8 +93,6 @@ class ProblemFDEM_e(Problem.BaseProblem): def Jvec(self, m, v, u=None): - # TODO: only 1 transmitter for now - # TODO: all P's the same if u is None: u = self.fields(m) @@ -108,17 +106,15 @@ class ProblemFDEM_e(Problem.BaseProblem): for tx in self.survey.getTransmitters(freq): dMe_dsig = self.mesh.getEdgeInnerProductDeriv(m, v=e) dsig_dm = self.model.transformDeriv(m) + b = 1j*omega(freq) * ( dMe_dsig * ( dsig_dm * v ) ) Ab = solver.solve(b) - #TODO: look at Rx for this... - P = self.survey.projectFieldsDeriv(u) + P = tx.projectFieldsDeriv(self.mesh, u) Jv[tx] = -P*Ab - Jv = np.concatenate(Jvs) - - return Jv + return Utils.mkvc(Jv) def Jtvec(self, m, v, u=None): diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index 0fadb93f..a0e5a33b 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -35,15 +35,18 @@ class TxFDEM(Survey.BaseTx): def projectFields(self, mesh, u): + if self.rxList.rxType in ['Ex', 'Ey', 'Ez']: + u_part = u[self, 'e'] + else: + raise NotImplemented('Unknown receiver type.') + P = self.rxList.getP(mesh) - - u_part = u[self] - Pu = P*u_part return Pu def projectFieldsDeriv(self, mesh, u): - pass + P = self.rxList.getP(mesh) + return P class FieldsFDEM(object): @@ -119,8 +122,10 @@ class FieldsFDEM(object): elif isinstance(key, TxFDEM): raise Exception('Cannot set one transmitter at a time.') - for field in newFields: + for name in newFields: field = self._initStore(name) + if field[freq].shape[1] == 1: + newFields[name] = Utils.mkvc(newFields[name],2) assert field[freq].shape == newFields[name].shape, 'Must be correct shape (n%s x nTx[freq])' % self.knownFields[name] field[freq] = newFields[name] @@ -173,13 +178,13 @@ class DataFDEM(object): self._ensureCorrectKey(key) return self._dataDict[key] - def toarray(self): + def tovec(self): D = self._dataDict return np.concatenate([D[k] for k in D]) -class SurveyFDEM(Survey.MixinFancyProjection, Survey.BaseSurvey): +class SurveyFDEM(Survey.BaseSurvey): """ docstring for SurveyFDEM """ @@ -228,35 +233,16 @@ class SurveyFDEM(Survey.MixinFancyProjection, Survey.BaseSurvey): self._nTx[freq] = len(self.getTransmitters(freq)) return self._nTx - def getTransmitters(self, freq): """Returns the transmitters associated with a specific frequency.""" assert freq in self._freqDict, "The requested frequency is not in this survey." return self._freqDict[freq] def projectFields(self, u): - P = sp.identity(self.prob.mesh.nE) - Pes = range(self.nFreq) - #TODO: this is hardcoded to 1Tx - for i, freqInd in enumerate(range(self.nFreq)): - e = u.get_e(freqInd) - Pes[i] = P*e - Pe = np.concatenate(Pes) - return Pe + data = DataFDEM(self) + for tx in self.txList: + data[tx] = tx.projectFields(self.mesh, u) + return data - - def projectFieldsDerivVec(self, u, v=None): - raise NotImplemented('projectFieldsDerivVec is not yet implemented') - - def projectAdjointFieldsDerivVec(self, u, v=None): - raise NotImplemented('projectAdjointFieldsDerivVec is not yet implemented') - - - #################################################### - # Interpolation Matrices - #################################################### - - @Utils.requires('prob') - def getP(self, Tx): - # TODO: store these in a mesh lookup table by transmitter? - pass + def projectFieldsDeriv(self, u): + raise Exception('Use Transmitters to project fields deriv.') diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index f40e5dd6..a2378d75 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -21,13 +21,12 @@ class FDEM_bDerivTests(unittest.TestCase): model = Model.LogModel(mesh) - opts = {'txLoc':0., - 'txType':'VMD_MVP', - 'rxLoc': rxLoc, - 'rxType':'bz', - 'freq': np.logspace(0,3,4), - } - dat = EM.FDEM.SurveyFDEM(**opts) + x = np.linspace(5,10,3) + XYZ = Utils.ndgrid(x,x,np.r_[0]) + rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex') + Tx0 = EM.FDEM.TxFDEM(None, 'VMD', 3, rxList) + + dat = EM.FDEM.SurveyFDEM([Tx0]) prb = EM.FDEM.ProblemFDEM_e(model) prb.pair(dat) diff --git a/simpegEM/Tests/test_FieldsObject.py b/simpegEM/Tests/test_FieldsObject.py index 302c4882..27137c75 100644 --- a/simpegEM/Tests/test_FieldsObject.py +++ b/simpegEM/Tests/test_FieldsObject.py @@ -22,10 +22,16 @@ class FieldsTest(unittest.TestCase): def test_SetGet(self): F = self.F for freq in F.survey.freqs: - e = np.random.rand(F.mesh.nE, F.survey.nTx[freq]) + nFreq = F.survey.nTx[freq] + e = np.random.rand(F.mesh.nE, nFreq) F[freq, 'e'] = e - b = np.random.rand(F.mesh.nF, F.survey.nTx[freq]) + b = np.random.rand(F.mesh.nF, nFreq) F[freq, 'b'] = b + if nFreq == 1: + F[freq, 'b'] = Utils.mkvc(b) + self.assertTrue(np.all(F[freq, 'e'] == e)) + self.assertTrue(np.all(F[freq, 'b'] == b)) + F[freq] = {'b':b,'e':e} self.assertTrue(np.all(F[freq, 'e'] == e)) self.assertTrue(np.all(F[freq, 'b'] == b))