wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-27 19:48:52 +08:00
Code Issues Packages Projects Releases Wiki Activity

prim-sec for EB formulation on same meshes

Browse Source
This commit is contained in:
Lindsey Heagy
2016-03-29 09:59:13 -07:00
parent d6daa93d9a
commit d11f5c736b
4 changed files with 65 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/Base.py
+15 -15
View File
@@ -2,14 +2,14 @@ from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solve
from scipy.constants import mu_0
class EMPropMap(Maps.PropMap):
"""
"""
Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m)
"""
sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap))
mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap))
mu = Maps.Property("Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap))
rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap))
rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap))
mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap))
@@ -21,7 +21,7 @@ class BaseEMProblem(Problem.BaseProblem):
surveyPair = Survey.BaseSurvey
dataPair = Survey.Data
PropMap = EMPropMap
Solver = SimpegSolver
@@ -51,7 +51,7 @@ class BaseEMProblem(Problem.BaseProblem):
if self.mapping.muMap is not None or self.mapping.muiMap is not None:
toDelete += ['_MeMu', '_MeMuI','_MfMui','_MfMuiI']
return toDelete
@property
def Me(self):
"""
@@ -71,7 +71,7 @@ class BaseEMProblem(Problem.BaseProblem):
return self._Mf
# ----- Magnetic Permeability ----- #
# ----- Magnetic Permeability ----- #
@property
def MfMui(self):
"""
@@ -109,7 +109,7 @@ class BaseEMProblem(Problem.BaseProblem):
return self._MeMuI
# ----- Electrical Conductivity ----- #
# ----- Electrical Conductivity ----- #
#TODO: hardcoded to sigma as the model
@property
def MeSigma(self):
@@ -120,18 +120,18 @@ class BaseEMProblem(Problem.BaseProblem):
self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma)
return self._MeSigma
# TODO: This should take a vector
# TODO: This should take a vector
def MeSigmaDeriv(self, u):
"""
Derivative of MeSigma with respect to the model
"""
"""
return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
@property
def MeSigmaI(self):
"""
Inverse of the edge inner product matrix for \\(\\sigma\\).
Inverse of the edge inner product matrix for \\(\\sigma\\).
"""
if getattr(self, '_MeSigmaI', None) is None:
self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True)
@@ -140,8 +140,8 @@ class BaseEMProblem(Problem.BaseProblem):
# TODO: This should take a vector
def MeSigmaIDeriv(self, u):
"""
Derivative of :code:`MeSigma` with respect to the model
"""
Derivative of :code:`MeSigma` with respect to the model
"""
# TODO: only works for diagonal tensors. getEdgeInnerProductDeriv, invMat=True should be implemented in SimPEG
dMeSigmaI_dI = -self.MeSigmaI**2
@@ -163,7 +163,7 @@ class BaseEMProblem(Problem.BaseProblem):
# TODO: This should take a vector
def MfRhoDeriv(self,u):
"""
Derivative of :code:`MfRho` with respect to the model.
Derivative of :code:`MfRho` with respect to the model.
"""
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv)
# self.curModel.rhoDeriv
@@ -181,6 +181,6 @@ class BaseEMProblem(Problem.BaseProblem):
# TODO: This should take a vector
def MfRhoIDeriv(self,u):
"""
Derivative of :code:`MfRhoI` with respect to the model.
Derivative of :code:`MfRhoI` with respect to the model.
"""
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
SimPEG/EM/FDEM/SrcFDEM.py
+35 -23
View File
@@ -608,46 +608,58 @@ class CircularLoop(BaseSrc):
class PrimSecSigma(BaseSrc):
# TODO: This will only work for E-B formulation
def __init__(self, rxList, freq, sigmaPrimary, ePrimary=None, primaryProblem=None, primarySurvey=None):
def __init__(self, rxList, freq, mPrimary, primaryProblem, primarySurvey):
self.freq = float(freq)
self._ePrimary = ePrimary
self.sigmaPrimary = sigmaPrimary
self._primaryProblem = primaryProblem
self._primarySurvey = primarySurvey
self.mPrimary = mPrimary
self.primaryProblem = primaryProblem
self.primarySurvey = primarySurvey
self.primaryMesh = self.primaryProblem.mesh
self.integrate = False
BaseSrc.__init__(self, rxList)
def _solvePrimary(self):
# check if I have fields
# if not, solve the primary to get fields object
if self._primarySurvey.ispaired:
if self._primarySurvey.prob is not self._primaryProblem:
raise Exception "The survey object is already paired to a problem. Use survey.unpair()"
else:
self._primaryProblem.pair(self._primarySurvey)
@property
def MeSigmaPrimary(self, prob):
if getattr(self, '_MeSigmaPrimary', None) is None:
sigmaprimary = self.primaryProblem.mapping.sigmaMap * self.mPrimary
# if self.primaryMesh != prob.mesh:
# if self.mesh._meshType == 'CYL' and prob.mesh._meshType != 'CYL':
# P = self.mesh.getInterpolationMatCartMesh(prob.mesh,'CC')
# sigmaprimary =
self._MeSigmaPrimary = prob.mesh.getEdgeInnerProduct(sigmaprimary)
return self._MeSigmaPrimary
return
@property
def _primaryFields(self):
if getattr(self, '__primaryFields', None) is None:
# check if I have fields
# if not, solve the primary to get fields object
if self.primarySurvey.ispaired:
if self.primarySurvey.prob is not self.primaryProblem:
raise Exception('The survey object is already paired to a problem. Use survey.unpair()')
else:
self.primaryProblem.pair(self.primarySurvey)
self.__primaryFields = self.primaryProblem.fields(self.mPrimary)
return self.__primaryFields
def ePrimary(self,prob):
# check if a primary problem is defined
if getattr(self, '_ePrimary', None) is None:
if self._primaryProblem is None or self._primarySurvey is None:
raise Exception "if Not specifying a ePrimary, a primarySurvey and primaryProblem must be provided."
if self.primaryProblem is None or self.primarySurvey is None:
raise Exception('if Not specifying a ePrimary, a primarySurvey and primaryProblem must be provided.')
self._ePrimary = self._solvePrimary[:,'e']
return self._ePrimary
return Utils.mkvc(self._primaryFields[:,'e'])
def S_e(self,prob):
# todo --> see if the primary and secondary are on the same mesh or not
MeSigma = prob.MeSigma
MeSigmaPrimary = prob.mesh.getEdgeInnerProduct(self.sigmaPrimary)
return (MeSigma - MeSigmaPrimary) * self._ePrimary
print type(prob.mesh)
return Utils.mkvc((MeSigma - self.MeSigmaPrimary(prob)) * self.ePrimary(prob))
def S_eDeriv(self, prob, v, adjoint = False):
MeSigmaDeriv = prob.MeSigmaDeriv
if adjoint is not True:
return MeSigmaDeriv(Utils.mkvc(self._ePrimary)) * v
return MeSigmaDeriv(Utils.mkvc(self._ePrimary)) * v
return MeSigmaDeriv(self.ePrimary(prob)) * v
return MeSigmaDeriv(self.ePrimary(prob)) * v
SimPEG/EM/Utils/testingUtils.py
+11 -7
View File
@@ -20,7 +20,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
if useMu is True:
mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))]
mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))]
else:
mapping = Maps.ExpMap(mesh)
@@ -53,11 +53,15 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
elif SrcType is 'PrimSecSigma':
if fdemType is 'e' or fdemType is '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.
Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_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.
primSrc = EM.FDEM.Src.MagDipole([], freq, np.r_[0.,0.,0.])
primarySurvey = EM.FDEM.Survey([primSrc])
primaryProblem = EM.FDEM.Problem_e(mesh,mapping=mapping)
mPrimary = np.ones(mapping.nP)*np.log(CONDUCTIVITY)
Src.append(EM.FDEM.Src.PrimSecSigma([Rx0], freq, mPrimary, primaryProblem, primarySurvey))
if verbose:
print ' Fetching %s problem' % (fdemType)
@@ -97,7 +101,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
mesh = prb1.mesh
print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
mu = np.ones(mesh.nC)*MU
tests/em/fdem/inverse/derivs/test_FDEM_derivs.py
+4 -4
View File
@@ -8,8 +8,8 @@ from SimPEG.EM.Utils.testingUtils import getFDEMProblem
testE = True
testB = True
testH = True
testJ = True
testH = False
testJ = False
verbose = False
@@ -20,8 +20,8 @@ MU = mu_0
freq = 1e-1
addrandoms = True
SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
# SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
SrcType = ['PrimSecSigma']
def derivTest(fdemType, comp):