From ac3e7765fe370d9ca85cca2acebeb4f66fa0c0e1 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 23 Jun 2015 18:32:00 -0700 Subject: [PATCH] Mag dipole sources, defined using prim-sec with a zero-frequency primary actually have a non-zer S_e component if there is variable mu. This has been corrected in the vector potential definition of the sources, and an issue created for the remaining sources. Testing has also been made more robust by using a rawVec source with both S_m and S_e defined which showed a couple bugs in some of the getRHSDeriv_m (which have been corrected), along with a couple minor sizing things in fields derivs --- simpegEM/FDEM/FDEM.py | 12 ++++------ simpegEM/FDEM/FieldsFDEM.py | 4 ++-- simpegEM/FDEM/SurveyFDEM.py | 42 +++++++++++++++++++++++---------- simpegEM/Tests/test_FDEM.py | 46 ++++++++++++++++++++++++++++--------- 4 files changed, 71 insertions(+), 33 deletions(-) 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)