wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-29 05:18:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
Lindsey
2015-06-23 18:32:00 -07:00
parent 36f8eca25c
commit ac3e7765fe
4 changed files with 71 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/FDEM/FDEM.py
+4 -8
View File
@@ -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)
simpegEM/FDEM/FieldsFDEM.py
+2 -2
View File
@@ -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
simpegEM/FDEM/SurveyFDEM.py
+30 -12
View File
@@ -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
simpegEM/Tests/test_FDEM.py
+35 -11
View File
@@ -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)