diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 47c3655e..507127f7 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -1,6 +1,7 @@ from SimPEG import Problem, Solver, Utils, np, sp from scipy.constants import mu_0 from SurveyFDEM import SurveyFDEM, DataFDEM, FieldsFDEM +from simpegEM.Utils import Sources def omega(freq): """Change frequency to angular frequency, omega""" @@ -60,16 +61,33 @@ class ProblemFDEM_e(Problem.BaseProblem): def getA(self, freq): """ - :param int fInd: Frequency index + :param float freq: Frequency :rtype: scipy.sparse.csr_matrix :return: A """ return self.mesh.edgeCurl.T*self.MfMui*self.mesh.edgeCurl + 1j*omega(freq)*self.MeSigma def getRHS(self, freq): - #TODO: this needs to also depend on your transmitter! + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nTx) + :return: RHS + """ + Txs = self.survey.getTransmitters(freq) + rhs = range(len(Txs)) + for i, tx in enumerate(Txs): + if tx.txType == 'VMD': + src = Sources.MagneticDipoleVectorPotential + else: + raise NotImplemented('%s txType is not implemented' % tx.txType) + SRCx = src(tx.loc, self.mesh.gridEx, 'x') + SRCy = src(tx.loc, self.mesh.gridEy, 'y') + SRCz = src(tx.loc, self.mesh.gridEz, 'z') + rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) - return -1j*omega(freq)*self.Me*self.j_s + #TODO: this is completely wrong. b0 = self.mesh.edgeCurl*rhs but we are doing an e formulation... + j_s = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') + return -1j*omega(freq)*self.Me*j_s def fields(self, m, useThisRhs=None): diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index a864b66b..23c53528 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -9,45 +9,31 @@ class FDEM_bDerivTests(unittest.TestCase): def setUp(self): cs = 5. - ncx = 2 - ncy = 2 - ncz = 2 + ncx, ncy, ncz = 2, 2, 2 npad = 3 hx = Utils.meshTensors(((npad,cs), (ncx,cs), (npad,cs))) hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs))) hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs))) mesh = Mesh.TensorMesh([hx,hy,hz]) - XY = Utils.ndgrid(np.linspace(20,50,3), np.linspace(20,50,3)) - rxLoc = np.c_[XY, np.ones(XY.shape[0])*40] - model = Model.LogModel(mesh) 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) + Tx0 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-2, rxList) x = np.linspace(5,10,3) XYZ = Utils.ndgrid(x,x,np.r_[0]) rxList = EM.FDEM.RxListFDEM(XYZ, 'Ey') - Tx1 = EM.FDEM.TxFDEM(None, 'VMD', 3, rxList) + Tx1 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-4, rxList) survey = EM.FDEM.SurveyFDEM([Tx0, Tx1]) prb = EM.FDEM.ProblemFDEM_e(model) prb.pair(survey) - sigma = np.log(np.ones(mesh.nC)*1e-3) - - j_sx = np.zeros(mesh.vnEx) - j_sx[4,4,4] = 1 - j_s = np.r_[Utils.mkvc(j_sx),np.zeros(mesh.nEy+mesh.nEz)] - - prb.j_s = np.c_[j_s, j_s] - f = prb.fields(sigma) - - self.sigma = sigma + self.sigma = np.log(np.ones(mesh.nC)*1e-3) self.prb = prb self.survey = survey