diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 5e465f36..dde77f15 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -19,7 +19,7 @@ class BaseFDEMProblem(BaseEMProblem): surveyPair = SurveyFDEM fieldsPair = FieldsFDEM - def fields(self, m): + def fields(self, m=None): self.curModel = m F = self.fieldsPair(self.mesh, self.survey) @@ -151,10 +151,6 @@ class BaseFDEMProblem(BaseEMProblem): return S_m, S_e - def getSourceTermDeriv(self,freq,m,v,u=None,adjoint=False): - raise NotImplementedError('getSourceTermDeriv not implemented yet') - return None, None - ########################################################################################## ################################ E-B Formulation ######################################### @@ -321,14 +317,14 @@ class ProblemFDEM_b(BaseFDEMProblem): def getRHSDeriv_m(self, src, v, adjoint=False): C = self.mesh.edgeCurl - S_m, S_e = self.getSourceTerm(src.freq) + S_m, S_e = src.eval(self) MfMui = self.MfMui if self._makeASymmetric and adjoint: v = self.MfMui * v if S_e is not None: - MeSigmaIDeriv = self.MeSigmaIDeriv(S_e) + MeSigmaIDeriv = self.MeSigmaIDeriv(Utils.mkvc(S_e)) if not adjoint: RHSderiv = C * (MeSigmaIDeriv * v) elif adjoint: @@ -577,7 +573,7 @@ class ProblemFDEM_h(BaseFDEMProblem): return RHS def getRHSDeriv_m(self, src, v, adjoint=False): - _, S_e = self.getSourceTerm(src.freq) + _, S_e = src.eval(self) C = self.mesh.edgeCurl MfRho = self.MfRho MfRhoDeriv = self.MfRhoDeriv(S_e) diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index 5838ddbe..341fd389 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -147,7 +147,7 @@ class FieldsFDEM_b(FieldsFDEM): for i,src in enumerate(srcList): _,S_e = src.eval(self.survey.prob) if S_e is not None: - e += -self._MeSigmaI*S_e + e[:,i] += -self._MeSigmaI*S_e return e def _eSecondaryDeriv_u(self, src, v, adjoint=False): @@ -162,7 +162,7 @@ class FieldsFDEM_b(FieldsFDEM): w = self._edgeCurl.T * (self._MfMui * bSolution) if S_e is not None: - w += -S_e + w += -Utils.mkvc(S_e,2) if not adjoint: de_dm = self._MeSigmaIDeriv(w) * v diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index 44dbfa58..ff60cb95 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -1,7 +1,7 @@ from SimPEG import Survey, Problem, Utils, np, sp from simpegEM.Utils import SrcUtils from simpegEM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b - +from scipy.constants import mu_0 #################################################### # Receivers @@ -189,11 +189,12 @@ class SrcFDEM_RawVec(SrcFDEM): class SrcFDEM_MagDipole(SrcFDEM): #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - def __init__(self, rxList, freq, loc, orientation='Z', moment=1.): + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): self.freq = float(freq) self.loc = loc self.orientation = orientation self.moment = moment + self.mu = mu SrcFDEM.__init__(self, rxList) def bPrimary(self,prob): @@ -216,13 +217,13 @@ class SrcFDEM_MagDipole(SrcFDEM): if not prob.mesh.isSymmetric: # TODO ? raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y') + a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu) else: srcfct = SrcUtils.MagneticDipoleVectorPotential - ax = srcfct(self.loc, gridX, 'x') - ay = srcfct(self.loc, gridY, 'y') - az = srcfct(self.loc, gridZ, 'z') + ax = srcfct(self.loc, gridX, 'x', mu=self.mu) + ay = srcfct(self.loc, gridY, 'y', mu=self.mu) + az = srcfct(self.loc, gridZ, 'z', mu=self.mu) a = np.concatenate((ax, ay, az)) return C*a @@ -235,17 +236,34 @@ class SrcFDEM_MagDipole(SrcFDEM): b_p = self.bPrimary(prob) return -1j*omega(self.freq)*b_p + def S_e(self,prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) class SrcFDEM_MagDipole_Bfield(SrcFDEM): #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that #TODO: neither does moment - def __init__(self, rxList, freq, loc, orientation='Z', moment=1.): + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): self.freq = float(freq) self.loc = loc self.orientation = orientation self.moment = moment + self.mu = mu SrcFDEM.__init__(self, rxList) def bPrimary(self,prob): @@ -268,13 +286,13 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM): if not prob.mesh.isSymmetric: # TODO ? raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - bx = srcfct(self.loc, gridX, 'x') - bz = srcfct(self.loc, gridZ, 'z') + bx = srcfct(self.loc, gridX, 'x', mu=self.mu) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu) b = np.concatenate((bx,bz)) else: - bx = srcfct(self.loc, gridX, 'x') - by = srcfct(self.loc, gridY, 'y') - bz = srcfct(self.loc, gridZ, 'z') + bx = srcfct(self.loc, gridX, 'x', mu=self.mu) + by = srcfct(self.loc, gridY, 'y', mu=self.mu) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu) b = np.concatenate((bx,by,bz)) return b diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index e56a92ba..bcddd91d 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -36,21 +36,42 @@ def getProblem(fdemType, comp): XYZ = Utils.ndgrid(x,x,np.r_[0.]) Rx0 = EM.FDEM.RxFDEM(XYZ, comp) Src0 = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.]) - Src1 = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq[1], loc=np.r_[0.,0.,0.]) - - survey = EM.FDEM.SurveyFDEM([Src0, Src1]) if verbose: print ' Fetching %s problem' % (fdemType) if fdemType == 'e': + S_m = np.zeros(mesh.nF) + S_e = np.zeros(mesh.nE) + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + Src1 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1]) prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) elif fdemType == 'b': + S_m = np.zeros(mesh.nF) + S_e = np.zeros(mesh.nE) + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + Src1 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1]) prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) elif fdemType == 'j': + S_m = np.zeros(mesh.nE) + S_e = np.zeros(mesh.nF) + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + Src1 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1]) prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) elif fdemType == 'h': + S_m = np.zeros(mesh.nE) + S_e = np.zeros(mesh.nF) + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + Src1 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1]) prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) else: raise NotImplementedError() @@ -69,15 +90,15 @@ def adjointTest(fdemType, comp): print 'Adjoint %s formulation - %s' % (fdemType, comp) m = np.log(np.ones(prb.mesh.nC)*CONDUCTIVITY) - mu = np.log(np.ones(prb.mesh.nC))*MU + mu = np.ones(prb.mesh.nC)*MU if addrandoms is True: - m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 + m = m + np.random.randn(prb.mesh.nC)*np.log(CONDUCTIVITY)*1e-1 mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 - prb.mu = mu survey = prb.survey - + prb.PropMap.PropModel.mu = mu + prb.PropMap.PropModel.mui = 1./mu u = prb.fields(m) v = np.random.rand(survey.nD) @@ -99,10 +120,12 @@ def derivTest(fdemType, comp): mu = np.log(np.ones(prb.mesh.nC)*MU) if addrandoms is True: - x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 + x0 = x0 + np.random.randn(prb.mesh.nC)*np.log(CONDUCTIVITY)*1e-1 mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 - prb.mu = mu + prb.PropMap.PropModel.mu = mu + prb.PropMap.PropModel.mui = 1./mu + survey = prb.survey def fun(x): return survey.dpred(x), lambda x: prb.Jvec(x0, x) @@ -120,10 +143,11 @@ def crossCheckTest(fdemType, comp): mu = np.log(np.ones(mesh.nC)*MU) if addrandoms is True: - m = m + np.random.randn(mesh.nC)*CONDUCTIVITY*1e-1 + m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1 mu = mu + np.random.randn(mesh.nC)*MU*1e-1 - prb1.mu = mu + prb1.PropMap.PropModel.mu = mu + prb1.PropMap.PropModel.mui = 1./mu survey1 = prb1.survey d1 = survey1.dpred(m)