From 7c01c84514bea0d650dd3d6071c0e23e592a6d31 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 23 Jun 2015 18:22:44 -0700 Subject: [PATCH 1/7] Analytics should take other values for mu --- simpegEM/Analytics/FDEM.py | 43 +++++++++++++++++++++++++------------- 1 file changed, 28 insertions(+), 15 deletions(-) diff --git a/simpegEM/Analytics/FDEM.py b/simpegEM/Analytics/FDEM.py index 36a2de63..a6c051cd 100644 --- a/simpegEM/Analytics/FDEM.py +++ b/simpegEM/Analytics/FDEM.py @@ -5,7 +5,8 @@ from scipy.special import erf import matplotlib.pyplot as plt from SimPEG import Utils -def hzAnalyticDipoleF(r, freq, sigma, secondary=True): + +def hzAnalyticDipoleF(r, freq, sigma, secondary=True, mu=mu_0): """ 4.56 in Ward and Hohmann @@ -25,7 +26,7 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True): """ r = np.abs(r) - k = np.sqrt(-1j*2.*np.pi*freq*mu_0*sigma) + k = np.sqrt(-1j*2.*np.pi*freq*mu*sigma) m = 1 front = m / (2. * np.pi * (k**2) * (r**5) ) @@ -41,7 +42,7 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True): return hz -def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'): +def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X', mu = mu_0): """ Analytical solution for a dipole in a whole-space. @@ -75,7 +76,7 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'): dz = XYZ[:,2]-srcLoc[2] r = np.sqrt( dx**2. + dy**2. + dz**2.) - k = np.sqrt( -1j*2.*np.pi*f*mu_0*sig ) + k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) kr = k*r front = m / (4.*pi * r**3.) * np.exp(-1j*kr) @@ -96,9 +97,9 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'): Hy = front*( (dy*dz/r**2.) * mid ) Hz = front*( (dz/r)**2. * mid + (kr**2. - 1j*kr - 1.) ) - Bx = mu_0*Hx - By = mu_0*Hy - Bz = mu_0*Hz + Bx = mu*Hx + By = mu*Hy + Bz = mu*Hz if Bx.ndim is 1: Bx = Utils.mkvc(Bx,2) @@ -112,7 +113,7 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'): return Bx, By, Bz -def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'): +def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): XYZ = Utils.asArray_N_x_Dim(XYZ, 3) dx = XYZ[:,0]-srcLoc[0] @@ -120,21 +121,33 @@ def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'): dz = XYZ[:,2]-srcLoc[2] r = np.sqrt( dx**2. + dy**2. + dz**2.) - k = np.sqrt( -1j*2.*np.pi*f*mu_0*sig ) + k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) kr = k*r - front = moment / (4. * np.pi * sig * r**3) * exp(-1j*k*r) + front = current * length / (4. * np.pi * sig * r**3) * np.exp(-1j*k*r) mid = -k**2 * r**2 + 3*1j*k*r + 3 - Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r)) - Ey = front*(dx*dy / r**2)*mid - Ez = front*(dx*dz / r**2)*mid + # Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r)) + # Ey = front*(dx*dy / r**2)*mid + # Ez = front*(dx*dz / r**2)*mid if orientation.upper() == 'X': + Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ey = front*(dx*dy / r**2)*mid + Ez = front*(dx*dz / r**2)*mid return Ex, Ey, Ez elif orientation.upper() == 'Y': - return Ez, Ex, Ey + # x--> y, y--> z, z-->x + Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ez = front*(dy*dz / r**2)*mid + Ex = front*(dy*dx / r**2)*mid + return Ex, Ey, Ez elif orientation.upper() == 'Z': - return Ey, Ez, Ex \ No newline at end of file + # x --> z, y --> x, z --> y + Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ex = front*(dz*dx / r**2)*mid + Ey = front*(dz*dy / r**2)*mid + return Ex, Ey, Ez + # return Ey, Ez, Ex \ No newline at end of file From 5732e85a2f7ff40788756c9952cd111d3ef03ba4 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 23 Jun 2015 18:24:40 -0700 Subject: [PATCH 2/7] fixed typo in MfRhoDeriv. note that it still won't work (but we don't call it anywhere... yet --- simpegEM/Base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 4205c518..5a17c69b 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -141,5 +141,5 @@ class BaseEMProblem(Problem.BaseProblem): # TODO: This isn't going to work yet # TODO: This should take a vector - def dMfRhoIDeriv(self,u): + def MfRhoIDeriv(self,u): return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv From 36f8eca25c04e443b1ec32059904013b1414b71f Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 23 Jun 2015 18:25:29 -0700 Subject: [PATCH 3/7] Sources can take mu values other than mu_0 --- simpegEM/Utils/SrcUtils.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/simpegEM/Utils/SrcUtils.py b/simpegEM/Utils/SrcUtils.py index ec638a0f..7c672bf4 100644 --- a/simpegEM/Utils/SrcUtils.py +++ b/simpegEM/Utils/SrcUtils.py @@ -2,7 +2,7 @@ from SimPEG import * from scipy.special import ellipk, ellipe from scipy.constants import mu_0, pi -def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)): +def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0., 1.), mu = mu_0): """ Calculate the vector potential of a set of magnetic dipoles at given locations 'ref. ' @@ -48,13 +48,13 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0 dR = obsLoc - srcLoc[i, np.newaxis].repeat(nEdges, axis=0) mCr = np.cross(m, dR) r = np.sqrt((dR**2).sum(axis=1)) - A[:, i] = +(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3) + A[:, i] = +(mu/(4*pi)) * mCr[:,dimInd]/(r**3) if nSrc == 1: return A.flatten() return A -def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.): +def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1., mu = mu_0): """ Calculate the vector potential of a set of magnetic dipoles at given locations 'ref. ' @@ -89,11 +89,11 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.): dR = obsLoc - srcLoc[i, np.newaxis].repeat(nFaces, axis=0) r = np.sqrt((dR**2).sum(axis=1)) if dimInd == 0: - B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2) + B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2) elif dimInd == 1: - B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2) + B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2) elif dimInd == 2: - B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1) + B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1) else: raise Exception("Not Implemented") if nSrc == 1: From ac3e7765fe370d9ca85cca2acebeb4f66fa0c0e1 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Tue, 23 Jun 2015 18:32:00 -0700 Subject: [PATCH 4/7] 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) From d9e899633661e98cd68b0f349368b7abcb398cb2 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 23 Jun 2015 23:02:33 -0700 Subject: [PATCH 5/7] fixed loop source and mag dipole fields for variable mu, added both to testing --- simpegEM/Analytics/FDEM.py | 4 +-- simpegEM/FDEM/SurveyFDEM.py | 66 ++++++++++++++++++++++++++++--------- simpegEM/Tests/test_FDEM.py | 19 ++++++----- simpegEM/Utils/SrcUtils.py | 17 +++++----- 4 files changed, 72 insertions(+), 34 deletions(-) diff --git a/simpegEM/Analytics/FDEM.py b/simpegEM/Analytics/FDEM.py index a6c051cd..9abb0a15 100644 --- a/simpegEM/Analytics/FDEM.py +++ b/simpegEM/Analytics/FDEM.py @@ -42,7 +42,7 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True, mu=mu_0): return hz -def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X', mu = mu_0): +def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, moment=1., orientation='X', mu = mu_0): """ Analytical solution for a dipole in a whole-space. @@ -79,7 +79,7 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X', mu = k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) kr = k*r - front = m / (4.*pi * r**3.) * np.exp(-1j*kr) + front = moment / (4.*pi * r**3.) * np.exp(-1j*kr) mid = -kr**2. + 3.*1j*kr + 3. if orientation.upper() == 'X': diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index ff60cb95..8a6937e1 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -217,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', mu=self.mu) + a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) else: srcfct = SrcUtils.MagneticDipoleVectorPotential - 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) + ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) a = np.concatenate((ax, ay, az)) return C*a @@ -286,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', mu=self.mu) - bz = srcfct(self.loc, gridZ, 'z', mu=self.mu) + bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) b = np.concatenate((bx,bz)) else: - 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) + bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) b = np.concatenate((bx,by,bz)) return b @@ -305,14 +305,33 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM): b = self.bPrimary(prob) return -1j*omega(self.freq)*b + 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_CircularLoop(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', radius = 1.): + def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): self.freq = float(freq) self.orientation = orientation self.radius = radius + self.mu = mu + self.loc = loc SrcFDEM.__init__(self, rxList) def bPrimary(self,prob): @@ -334,25 +353,42 @@ class SrcFDEM_CircularLoop(SrcFDEM): if not prob.mesh.isSymmetric: # TODO ? raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - a = SrcUtils.MagneticDipoleVectorPotential(src.loc, gridY, 'y', self.radius) + a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) else: srcfct = SrcUtils.MagneticDipoleVectorPotential - ax = srcfct(self.loc, gridX, 'x', self.radius) - ay = srcfct(self.loc, gridY, 'y', self.radius) - az = srcfct(self.loc, gridZ, 'z', self.radius) + ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu) + ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu) + az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu) a = np.concatenate((ax, ay, az)) return C*a def hPrimary(self,prob): b = self.bPrimary(prob) - return h_from_b + return 1./self.mu*b def S_m(self, prob): b = self.bPrimary(prob) return -1j*omega(self.freq)*b + 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)) + #################################################### # Survey diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index bcddd91d..e3cbfc7b 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -36,7 +36,8 @@ 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_Bfield([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.]) + Src2 = EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.]) if verbose: print ' Fetching %s problem' % (fdemType) @@ -46,32 +47,32 @@ def getProblem(fdemType, comp): 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]) + Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) 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]) + Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) 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]) + Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) 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]) + Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) else: raise NotImplementedError() diff --git a/simpegEM/Utils/SrcUtils.py b/simpegEM/Utils/SrcUtils.py index 7c672bf4..1827f6b2 100644 --- a/simpegEM/Utils/SrcUtils.py +++ b/simpegEM/Utils/SrcUtils.py @@ -2,7 +2,7 @@ from SimPEG import * from scipy.special import ellipk, ellipe from scipy.constants import mu_0, pi -def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0., 1.), mu = mu_0): +def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMoment=(0., 0., 1.), mu = mu_0): """ Calculate the vector potential of a set of magnetic dipoles at given locations 'ref. ' @@ -15,6 +15,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0 :return: The vector potential each dipole at each observation location """ #TODO: break this out! + if type(component) in [list, tuple]: out = range(len(component)) for i, comp in enumerate(component): @@ -54,7 +55,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0 return A -def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1., mu = mu_0): +def MagneticDipoleFields(srcLoc, obsLoc, component, moment=1., mu = mu_0): """ Calculate the vector potential of a set of magnetic dipoles at given locations 'ref. ' @@ -62,7 +63,7 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1., mu = mu_0): :param numpy.ndarray srcLoc: Location of the source(s) (x, y, z) :param numpy.ndarray obsLoc: Where the potentials will be calculated (x, y, z) :param str component: The component to calculate - 'x', 'y', or 'z' - :param numpy.ndarray dipoleMoment: The vector dipole moment (vertical) + :param numpy.ndarray moment: The vector dipole moment (vertical) :rtype: numpy.ndarray :return: The vector potential each dipole at each observation location """ @@ -78,12 +79,12 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1., mu = mu_0): srcLoc = np.atleast_2d(srcLoc) obsLoc = np.atleast_2d(obsLoc) - dipoleMoment = np.atleast_2d(dipoleMoment) + moment = np.atleast_2d(moment) nFaces = obsLoc.shape[0] nSrc = srcLoc.shape[0] - m = np.array(dipoleMoment).repeat(nFaces, axis=0) + m = np.array(moment).repeat(nFaces, axis=0) B = np.empty((nFaces, nSrc)) for i in range(nSrc): dR = obsLoc - srcLoc[i, np.newaxis].repeat(nFaces, axis=0) @@ -102,7 +103,7 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1., mu = mu_0): -def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius): +def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0): """ Calculate the vector potential of horizontal circular loop at given locations @@ -119,13 +120,13 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius): if type(component) in [list, tuple]: out = range(len(component)) for i, comp in enumerate(component): - out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius) + out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius, mu) return np.concatenate(out) if isinstance(obsLoc, Mesh.BaseMesh): mesh = obsLoc assert component in ['Ex','Ey','Ez','Fx','Fy','Fz'], "Components must be in: ['Ex','Ey','Ez','Fx','Fy','Fz']" - return MagneticLoopVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], radius) + return MagneticLoopVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], radius, mu) srcLoc = np.atleast_2d(srcLoc) obsLoc = np.atleast_2d(obsLoc) From 1fbf40568d45278a5eaf7714c0f960cf0349cbeb Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 24 Jun 2015 16:57:15 -0700 Subject: [PATCH 6/7] shortened the test for travis --- simpegEM/Tests/test_FDEM.py | 74 +++++++++++++++++++++++-------------- 1 file changed, 46 insertions(+), 28 deletions(-) diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index e3cbfc7b..540ff24f 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -17,9 +17,11 @@ TOL = 1e-4 FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order CONDUCTIVITY = 1e1 MU = mu_0 -freq = [1e-1, 2e-1] +freq = 1e-1 addrandoms = True +SrcType = 'MagDipole' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' + def getProblem(fdemType, comp): cs = 5. @@ -35,45 +37,61 @@ def getProblem(fdemType, comp): x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source 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_Bfield([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.]) - Src2 = EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.]) + + if SrcType is 'MagDipole': + Src = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) + elif SrcType is 'MagDipole_Bfield': + Src = EM.FDEM.SrcFDEM_MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) + elif SrcType is 'CircularLoop': + Src2 = EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) 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. - Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) - survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) + if SrcType is 'RawVec': + 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. + Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) + + survey = EM.FDEM.SurveyFDEM([Src]) 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. - Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) - survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) + if SrcType is 'RawVec': + 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. + Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) + + survey = EM.FDEM.SurveyFDEM([Src]) 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. - Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) - survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) + if SrcType is 'RawVec': + 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. + Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) + + survey = EM.FDEM.SurveyFDEM([Src]) 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. - Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) - survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) + if SrcType is 'RawVec': + 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. + Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) + + survey = EM.FDEM.SurveyFDEM([Src]) prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) + else: raise NotImplementedError() prb.pair(survey) From 5ec86300f6b4cc941fa16111a19d9ac60bc4a892 Mon Sep 17 00:00:00 2001 From: Lindsey Date: Mon, 29 Jun 2015 18:05:32 -0700 Subject: [PATCH 7/7] have fields keep track of the problem so that calls to src.eval use the proper mass matrices if a model is switched out when running multiple forwards. This way fields also stores curModel --- simpegEM/FDEM/FieldsFDEM.py | 39 ++++++++++++++++++++----------------- 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index 341fd389..46b67fa7 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -7,7 +7,6 @@ class FieldsFDEM(Problem.Fields): knownFields = {} dtype = complex - class FieldsFDEM_e(FieldsFDEM): knownFields = {'eSolution':'E'} aliasFields = { @@ -23,6 +22,7 @@ class FieldsFDEM_e(FieldsFDEM): FieldsFDEM.__init__(self,mesh,survey,**kwargs) def startup(self): + self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl # def getDeriv_u(self, fieldsList, src, v, adjoint=False): @@ -32,7 +32,7 @@ class FieldsFDEM_e(FieldsFDEM): def _ePrimary(self, eSolution, srcList): ePrimary = np.zeros_like(eSolution) for i, src in enumerate(srcList): - ep = src.ePrimary(self.survey.prob) + ep = src.ePrimary(self.prob) if ep is not None: ePrimary[:,i] = ep return ePrimary @@ -53,7 +53,7 @@ class FieldsFDEM_e(FieldsFDEM): def _bPrimary(self, eSolution, srcList): bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex) for i, src in enumerate(srcList): - bp = src.bPrimary(self.survey.prob) + bp = src.bPrimary(self.prob) if bp is not None: bPrimary[:,i] += bp return bPrimary @@ -63,7 +63,7 @@ class FieldsFDEM_e(FieldsFDEM): b = (C * eSolution) for i, src in enumerate(srcList): b[:,i] *= - 1./(1j*omega(src.freq)) - S_m, _ = src.eval(self.survey.prob) + S_m, _ = src.eval(self.prob) if S_m is not None: b[:,i] += 1./(1j*omega(src.freq)) * S_m return b @@ -75,7 +75,7 @@ class FieldsFDEM_e(FieldsFDEM): return - 1./(1j*omega(src.freq)) * (C * v) def _bSecondaryDeriv_m(self, src, v, adjoint = False): - S_mDeriv, _ = src.evalDeriv(self.survey.prob, adjoint) + S_mDeriv, _ = src.evalDeriv(self.prob, adjoint) S_mDeriv = S_mDeriv(v) if S_mDeriv is not None: return 1./(1j * omega(src.freq)) * S_mDeriv @@ -108,6 +108,7 @@ class FieldsFDEM_b(FieldsFDEM): FieldsFDEM.__init__(self,mesh,survey,**kwargs) def startup(self): + self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeSigmaI = self.survey.prob.MeSigmaI self._MfMui = self.survey.prob.MfMui @@ -116,7 +117,7 @@ class FieldsFDEM_b(FieldsFDEM): def _bPrimary(self, bSolution, srcList): bPrimary = np.zeros_like(bSolution) for i, src in enumerate(srcList): - bp = src.bPrimary(self.survey.prob) + bp = src.bPrimary(self.prob) if bp is not None: bPrimary[:,i] = bp return bPrimary @@ -137,7 +138,7 @@ class FieldsFDEM_b(FieldsFDEM): def _ePrimary(self, bSolution, srcList): ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex) for i,src in enumerate(srcList): - ep = src.ePrimary(self.survey.prob) + ep = src.ePrimary(self.prob) if ep is not None: ePrimary[:,i] = ep return ePrimary @@ -145,7 +146,7 @@ class FieldsFDEM_b(FieldsFDEM): def _eSecondary(self, bSolution, srcList): e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution)) for i,src in enumerate(srcList): - _,S_e = src.eval(self.survey.prob) + _,S_e = src.eval(self.prob) if S_e is not None: e[:,i] += -self._MeSigmaI*S_e return e @@ -158,7 +159,7 @@ class FieldsFDEM_b(FieldsFDEM): def _eSecondaryDeriv_m(self, src, v, adjoint=False): bSolution = self[[src],'bSolution'] - _,S_e = src.eval(self.survey.prob) + _,S_e = src.eval(self.prob) w = self._edgeCurl.T * (self._MfMui * bSolution) if S_e is not None: @@ -169,7 +170,7 @@ class FieldsFDEM_b(FieldsFDEM): elif adjoint: de_dm = self._MeSigmaIDeriv(w).T * v - _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint) + _, S_eDeriv = src.evalDeriv(self.prob, adjoint) Se_Deriv = S_eDeriv(v) if Se_Deriv is not None: @@ -203,6 +204,7 @@ class FieldsFDEM_j(FieldsFDEM): FieldsFDEM.__init__(self,mesh,survey,**kwargs) def startup(self): + self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeMuI = self.survey.prob.MeMuI self._MfRho = self.survey.prob.MfRho @@ -211,7 +213,7 @@ class FieldsFDEM_j(FieldsFDEM): def _jPrimary(self, jSolution, srcList): jPrimary = np.zeros_like(jSolution,dtype = complex) for i, src in enumerate(srcList): - jp = src.jPrimary(self.survey.prob) + jp = src.jPrimary(self.prob) if jp is not None: jPrimary[:,i] += jp return jPrimary @@ -232,7 +234,7 @@ class FieldsFDEM_j(FieldsFDEM): def _hPrimary(self, jSolution, srcList): hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex) for i, src in enumerate(srcList): - hp = src.hPrimary(self.survey.prob) + hp = src.hPrimary(self.prob) if hp is not None: hPrimary[:,i] = hp return hPrimary @@ -244,7 +246,7 @@ class FieldsFDEM_j(FieldsFDEM): h = MeMuI * (C.T * (MfRho * jSolution) ) for i, src in enumerate(srcList): h[:,i] *= -1./(1j*omega(src.freq)) - S_m,_ = src.eval(self.survey.prob) + S_m,_ = src.eval(self.prob) if S_m is not None: h[:,i] += 1./(1j*omega(src.freq)) * MeMuI * S_m return h @@ -270,7 +272,7 @@ class FieldsFDEM_j(FieldsFDEM): elif adjoint: hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) ) - S_mDeriv,_ = src.evalDeriv(self.survey.prob, adjoint) + S_mDeriv,_ = src.evalDeriv(self.prob, adjoint) if not adjoint: S_mDeriv = S_mDeriv(v) @@ -309,6 +311,7 @@ class FieldsFDEM_h(FieldsFDEM): FieldsFDEM.__init__(self,mesh,survey,**kwargs) def startup(self): + self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeMuI = self.survey.prob.MeMuI self._MfRho = self.survey.prob.MfRho @@ -316,7 +319,7 @@ class FieldsFDEM_h(FieldsFDEM): def _hPrimary(self, hSolution, srcList): hPrimary = np.zeros_like(hSolution,dtype = complex) for i, src in enumerate(srcList): - hp = src.hPrimary(self.survey.prob) + hp = src.hPrimary(self.prob) if hp is not None: hPrimary[:,i] += hp return hPrimary @@ -337,7 +340,7 @@ class FieldsFDEM_h(FieldsFDEM): def _jPrimary(self, hSolution, srcList): jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]]) for i, src in enumerate(srcList): - jp = src.jPrimary(self.survey.prob) + jp = src.jPrimary(self.prob) if jp is not None: jPrimary[:,i] = jp return jPrimary @@ -345,7 +348,7 @@ class FieldsFDEM_h(FieldsFDEM): def _jSecondary(self, hSolution, srcList): j = self._edgeCurl*hSolution for i, src in enumerate(srcList): - _,S_e = src.eval(self.survey.prob) + _,S_e = src.eval(self.prob) if S_e is not None: j[:,i] += -S_e return j @@ -357,7 +360,7 @@ class FieldsFDEM_h(FieldsFDEM): return self._edgeCurl.T*v def _jSecondaryDeriv_m(self, src, v, adjoint=False): - _,S_eDeriv = src.evalDeriv(self.survey.prob, adjoint) + _,S_eDeriv = src.evalDeriv(self.prob, adjoint) S_eDeriv = S_eDeriv(v) if S_eDeriv is not None: return -S_eDeriv