mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-28 07:34:52 +08:00
Rowan Cockett
+16
-16
@@ -1,9 +1,9 @@
|
||||
from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
|
||||
from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
|
||||
from scipy.constants import mu_0
|
||||
from SurveyFDEM import SurveyFDEM
|
||||
from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b, FieldsFDEM_h, FieldsFDEM_j
|
||||
from SurveyFDEM import Survey as SurveyFDEM
|
||||
from FieldsFDEM import Fields, Fields_e, Fields_b, Fields_h, Fields_j
|
||||
from SimPEG.EM.Base import BaseEMProblem
|
||||
from SimPEG.EM.Utils.EMUtils import omega
|
||||
from SimPEG.EM.Utils import omega
|
||||
|
||||
|
||||
class BaseFDEMProblem(BaseEMProblem):
|
||||
@@ -17,8 +17,8 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
\mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
|
||||
{\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{M^e} \mathbf{s_e}}
|
||||
|
||||
if using the E-B formulation (:code:`ProblemFDEM_e`
|
||||
or :code:`ProblemFDEM_b`) or the magnetic field
|
||||
if using the E-B formulation (:code:`Problem_e`
|
||||
or :code:`Problem_b`) or the magnetic field
|
||||
\\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
|
||||
|
||||
.. math ::
|
||||
@@ -26,13 +26,13 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\\\
|
||||
\mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
|
||||
|
||||
if using the H-J formulation (:code:`ProblemFDEM_j` or :code:`ProblemFDEM_h`).
|
||||
if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`).
|
||||
|
||||
The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
|
||||
"""
|
||||
|
||||
surveyPair = SurveyFDEM
|
||||
fieldsPair = FieldsFDEM
|
||||
fieldsPair = Fields
|
||||
|
||||
def fields(self, m=None):
|
||||
"""
|
||||
@@ -185,7 +185,7 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
################################ E-B Formulation #########################################
|
||||
##########################################################################################
|
||||
|
||||
class ProblemFDEM_e(BaseFDEMProblem):
|
||||
class Problem_e(BaseFDEMProblem):
|
||||
"""
|
||||
By eliminating the magnetic flux density using
|
||||
|
||||
@@ -205,7 +205,7 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
|
||||
_fieldType = 'e'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = FieldsFDEM_e
|
||||
fieldsPair = Fields_e
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, mesh, **kwargs)
|
||||
@@ -284,7 +284,7 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
return None
|
||||
|
||||
|
||||
class ProblemFDEM_b(BaseFDEMProblem):
|
||||
class Problem_b(BaseFDEMProblem):
|
||||
"""
|
||||
We eliminate \\\(\\\mathbf{e}\\\) using
|
||||
|
||||
@@ -304,7 +304,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
|
||||
|
||||
_fieldType = 'b'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = FieldsFDEM_b
|
||||
fieldsPair = Fields_b
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, mesh, **kwargs)
|
||||
@@ -425,7 +425,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
|
||||
##########################################################################################
|
||||
|
||||
|
||||
class ProblemFDEM_j(BaseFDEMProblem):
|
||||
class Problem_j(BaseFDEMProblem):
|
||||
"""
|
||||
We eliminate \\\(\\\mathbf{h}\\\) using
|
||||
|
||||
@@ -446,7 +446,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
|
||||
_fieldType = 'j'
|
||||
_eqLocs = 'EF'
|
||||
fieldsPair = FieldsFDEM_j
|
||||
fieldsPair = Fields_j
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, mesh, **kwargs)
|
||||
@@ -558,7 +558,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
|
||||
|
||||
|
||||
class ProblemFDEM_h(BaseFDEMProblem):
|
||||
class Problem_h(BaseFDEMProblem):
|
||||
"""
|
||||
We eliminate \\\(\\\mathbf{j}\\\) using
|
||||
|
||||
@@ -576,7 +576,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
|
||||
|
||||
_fieldType = 'h'
|
||||
_eqLocs = 'EF'
|
||||
fieldsPair = FieldsFDEM_h
|
||||
fieldsPair = Fields_h
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
@@ -1,13 +1,16 @@
|
||||
from SimPEG import Survey, Problem, Utils, np, sp
|
||||
from SimPEG.EM.Utils.EMUtils import omega
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import SimPEG
|
||||
from SimPEG import Utils
|
||||
from SimPEG.EM.Utils import omega
|
||||
|
||||
|
||||
class FieldsFDEM(Problem.Fields):
|
||||
class Fields(SimPEG.Problem.Fields):
|
||||
"""Fancy Field Storage for a FDEM survey."""
|
||||
knownFields = {}
|
||||
dtype = complex
|
||||
|
||||
class FieldsFDEM_e(FieldsFDEM):
|
||||
class Fields_e(Fields):
|
||||
knownFields = {'eSolution':'E'}
|
||||
aliasFields = {
|
||||
'e' : ['eSolution','E','_e'],
|
||||
@@ -19,7 +22,7 @@ class FieldsFDEM_e(FieldsFDEM):
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
|
||||
Fields.__init__(self,mesh,survey,**kwargs)
|
||||
|
||||
def startup(self):
|
||||
self.prob = self.survey.prob
|
||||
@@ -89,7 +92,7 @@ class FieldsFDEM_e(FieldsFDEM):
|
||||
return self._bSecondaryDeriv_m(src, v, adjoint)
|
||||
|
||||
|
||||
class FieldsFDEM_b(FieldsFDEM):
|
||||
class Fields_b(Fields):
|
||||
knownFields = {'bSolution':'F'}
|
||||
aliasFields = {
|
||||
'b' : ['bSolution','F','_b'],
|
||||
@@ -101,7 +104,7 @@ class FieldsFDEM_b(FieldsFDEM):
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
|
||||
Fields.__init__(self,mesh,survey,**kwargs)
|
||||
|
||||
def startup(self):
|
||||
self.prob = self.survey.prob
|
||||
@@ -190,7 +193,7 @@ class FieldsFDEM_b(FieldsFDEM):
|
||||
return self._eSecondaryDeriv_m(src, v, adjoint)
|
||||
|
||||
|
||||
class FieldsFDEM_j(FieldsFDEM):
|
||||
class Fields_j(Fields):
|
||||
knownFields = {'jSolution':'F'}
|
||||
aliasFields = {
|
||||
'j' : ['jSolution','F','_j'],
|
||||
@@ -202,7 +205,7 @@ class FieldsFDEM_j(FieldsFDEM):
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
|
||||
Fields.__init__(self,mesh,survey,**kwargs)
|
||||
|
||||
def startup(self):
|
||||
self.prob = self.survey.prob
|
||||
@@ -293,7 +296,7 @@ class FieldsFDEM_j(FieldsFDEM):
|
||||
return self._hSecondaryDeriv_m(src, v, adjoint)
|
||||
|
||||
|
||||
class FieldsFDEM_h(FieldsFDEM):
|
||||
class Fields_h(Fields):
|
||||
knownFields = {'hSolution':'E'}
|
||||
aliasFields = {
|
||||
'h' : ['hSolution','E','_h'],
|
||||
@@ -305,7 +308,7 @@ class FieldsFDEM_h(FieldsFDEM):
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
|
||||
Fields.__init__(self,mesh,survey,**kwargs)
|
||||
|
||||
def startup(self):
|
||||
self.prob = self.survey.prob
|
||||
|
||||
@@ -0,0 +1,347 @@
|
||||
from SimPEG import Survey, Problem, Utils, np, sp
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils import *
|
||||
# from SurveyFDEM import Rx
|
||||
|
||||
|
||||
class BaseSrc(Survey.BaseSrc):
|
||||
freq = None
|
||||
# rxPair = Rx
|
||||
integrate = True
|
||||
|
||||
def eval(self, prob):
|
||||
S_m = self.S_m(prob)
|
||||
S_e = self.S_e(prob)
|
||||
return S_m, S_e
|
||||
|
||||
def evalDeriv(self, prob, v, adjoint=False):
|
||||
return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
|
||||
|
||||
def bPrimary(self, prob):
|
||||
return None
|
||||
|
||||
def hPrimary(self, prob):
|
||||
return None
|
||||
|
||||
def ePrimary(self, prob):
|
||||
return None
|
||||
|
||||
def jPrimary(self, prob):
|
||||
return None
|
||||
|
||||
def S_m(self, prob):
|
||||
return None
|
||||
|
||||
def S_e(self, prob):
|
||||
return None
|
||||
|
||||
def S_mDeriv(self, prob, v, adjoint = False):
|
||||
return None
|
||||
|
||||
def S_eDeriv(self, prob, v, adjoint = False):
|
||||
return None
|
||||
|
||||
|
||||
class RawVec_e(BaseSrc):
|
||||
"""
|
||||
RawVec electric source. It is defined by the user provided vector S_e
|
||||
|
||||
:param numpy.array S_e: electric source term
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
|
||||
def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
|
||||
self._S_e = np.array(S_e,dtype=complex)
|
||||
self._ePrimary = ePrimary
|
||||
self._bPrimary = bPrimary
|
||||
self._hPrimary = hPrimary
|
||||
self._jPrimary = jPrimary
|
||||
self.freq = float(freq)
|
||||
BaseSrc.__init__(self, rxList)
|
||||
|
||||
def S_e(self, prob):
|
||||
return self._S_e
|
||||
|
||||
def ePrimary(self, prob):
|
||||
return self._ePrimary
|
||||
|
||||
def bPrimary(self, prob):
|
||||
return self._bPrimary
|
||||
|
||||
def hPrimary(self, prob):
|
||||
return self._hPrimary
|
||||
|
||||
def jPrimary(self, prob):
|
||||
return self._jPrimary
|
||||
|
||||
|
||||
class RawVec_m(BaseSrc):
|
||||
"""
|
||||
RawVec magnetic source. It is defined by the user provided vector S_m
|
||||
|
||||
:param numpy.array S_m: magnetic source term
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
|
||||
def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
|
||||
self._S_m = np.array(S_m,dtype=complex)
|
||||
self.freq = float(freq)
|
||||
self.integrate = integrate
|
||||
self._ePrimary = np.array(ePrimary,dtype=complex)
|
||||
self._bPrimary = np.array(bPrimary,dtype=complex)
|
||||
self._hPrimary = np.array(hPrimary,dtype=complex)
|
||||
self._jPrimary = np.array(jPrimary,dtype=complex)
|
||||
|
||||
BaseSrc.__init__(self, rxList)
|
||||
|
||||
def S_m(self, prob):
|
||||
return self._S_m
|
||||
|
||||
def ePrimary(self, prob):
|
||||
return self._ePrimary
|
||||
|
||||
def bPrimary(self, prob):
|
||||
return self._bPrimary
|
||||
|
||||
def hPrimary(self, prob):
|
||||
return self._hPrimary
|
||||
|
||||
def jPrimary(self, prob):
|
||||
return self._jPrimary
|
||||
|
||||
|
||||
class RawVec(BaseSrc):
|
||||
"""
|
||||
RawVec source. It is defined by the user provided vectors S_m, S_e
|
||||
|
||||
:param numpy.array S_m: magnetic source term
|
||||
:param numpy.array S_e: electric source term
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
def __init__(self, rxList, freq, S_m, S_e, integrate = True):
|
||||
self._S_m = np.array(S_m,dtype=complex)
|
||||
self._S_e = np.array(S_e,dtype=complex)
|
||||
self.freq = float(freq)
|
||||
self.integrate = integrate
|
||||
BaseSrc.__init__(self, rxList)
|
||||
|
||||
def S_m(self, prob):
|
||||
if prob._eqLocs is 'EF' and self.integrate is True:
|
||||
return prob.Me * self._S_m
|
||||
return self._S_m
|
||||
|
||||
def S_e(self, prob):
|
||||
if prob._eqLocs is 'FE' and self.integrate is True:
|
||||
return prob.Me * self._S_e
|
||||
return self._S_e
|
||||
|
||||
|
||||
class MagDipole(BaseSrc):
|
||||
|
||||
#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., mu = mu_0):
|
||||
self.freq = float(freq)
|
||||
self.loc = loc
|
||||
self.orientation = orientation
|
||||
self.moment = moment
|
||||
self.mu = mu
|
||||
self.integrate = False
|
||||
BaseSrc.__init__(self, rxList)
|
||||
|
||||
def bPrimary(self, prob):
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
if eqLocs is 'FE':
|
||||
gridX = prob.mesh.gridEx
|
||||
gridY = prob.mesh.gridEy
|
||||
gridZ = prob.mesh.gridEz
|
||||
C = prob.mesh.edgeCurl
|
||||
|
||||
elif eqLocs is 'EF':
|
||||
gridX = prob.mesh.gridFx
|
||||
gridY = prob.mesh.gridFy
|
||||
gridZ = prob.mesh.gridFz
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
|
||||
if prob.mesh._meshType is 'CYL':
|
||||
if not prob.mesh.isSymmetric:
|
||||
# TODO ?
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
|
||||
|
||||
else:
|
||||
srcfct = MagneticDipoleVectorPotential
|
||||
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
|
||||
|
||||
def hPrimary(self, prob):
|
||||
b = self.bPrimary(prob)
|
||||
return h_from_b(prob,b)
|
||||
|
||||
def S_m(self, prob):
|
||||
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 MagDipole_Bfield(BaseSrc):
|
||||
|
||||
#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., mu = mu_0):
|
||||
self.freq = float(freq)
|
||||
self.loc = loc
|
||||
self.orientation = orientation
|
||||
self.moment = moment
|
||||
self.mu = mu
|
||||
BaseSrc.__init__(self, rxList)
|
||||
|
||||
def bPrimary(self, prob):
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
if eqLocs is 'FE':
|
||||
gridX = prob.mesh.gridFx
|
||||
gridY = prob.mesh.gridFy
|
||||
gridZ = prob.mesh.gridFz
|
||||
C = prob.mesh.edgeCurl
|
||||
|
||||
elif eqLocs is 'EF':
|
||||
gridX = prob.mesh.gridEx
|
||||
gridY = prob.mesh.gridEy
|
||||
gridZ = prob.mesh.gridEz
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
srcfct = MagneticDipoleFields
|
||||
if prob.mesh._meshType is 'CYL':
|
||||
if not prob.mesh.isSymmetric:
|
||||
# TODO ?
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
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, 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
|
||||
|
||||
def hPrimary(self, prob):
|
||||
b = self.bPrimary(prob)
|
||||
return h_from_b(prob, 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))
|
||||
|
||||
|
||||
class CircularLoop(BaseSrc):
|
||||
|
||||
#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., mu=mu_0):
|
||||
self.freq = float(freq)
|
||||
self.orientation = orientation
|
||||
self.radius = radius
|
||||
self.mu = mu
|
||||
self.loc = loc
|
||||
self.integrate = False
|
||||
BaseSrc.__init__(self, rxList)
|
||||
|
||||
def bPrimary(self, prob):
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
if eqLocs is 'FE':
|
||||
gridX = prob.mesh.gridEx
|
||||
gridY = prob.mesh.gridEy
|
||||
gridZ = prob.mesh.gridEz
|
||||
C = prob.mesh.edgeCurl
|
||||
|
||||
elif eqLocs is 'EF':
|
||||
gridX = prob.mesh.gridFx
|
||||
gridY = prob.mesh.gridFy
|
||||
gridZ = prob.mesh.gridFz
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
if prob.mesh._meshType is 'CYL':
|
||||
if not prob.mesh.isSymmetric:
|
||||
# TODO ?
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
|
||||
|
||||
else:
|
||||
srcfct = MagneticDipoleVectorPotential
|
||||
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 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))
|
||||
|
||||
|
||||
@@ -1,13 +1,13 @@
|
||||
from SimPEG import Survey, Problem, Utils, np, sp
|
||||
from SimPEG.EM.Utils import SrcUtils
|
||||
from SimPEG.EM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
|
||||
import SimPEG
|
||||
from SimPEG.EM.Utils import *
|
||||
from scipy.constants import mu_0
|
||||
import SrcFDEM as Src
|
||||
|
||||
####################################################
|
||||
# Receivers
|
||||
####################################################
|
||||
|
||||
class RxFDEM(Survey.BaseRx):
|
||||
class Rx(SimPEG.Survey.BaseRx):
|
||||
|
||||
knownRxTypes = {
|
||||
'exr':['e', 'Ex', 'real'],
|
||||
@@ -41,7 +41,7 @@ class RxFDEM(Survey.BaseRx):
|
||||
radius = None
|
||||
|
||||
def __init__(self, locs, rxType):
|
||||
Survey.BaseRx.__init__(self, locs, rxType)
|
||||
SimPEG.Survey.BaseRx.__init__(self, locs, rxType)
|
||||
|
||||
@property
|
||||
def projField(self):
|
||||
@@ -87,366 +87,21 @@ class RxFDEM(Survey.BaseRx):
|
||||
return Pv
|
||||
|
||||
|
||||
####################################################
|
||||
# Sources
|
||||
####################################################
|
||||
|
||||
class SrcFDEM(Survey.BaseSrc):
|
||||
freq = None
|
||||
rxPair = RxFDEM
|
||||
integrate = True
|
||||
|
||||
def eval(self, prob):
|
||||
S_m = self.S_m(prob)
|
||||
S_e = self.S_e(prob)
|
||||
return S_m, S_e
|
||||
|
||||
def evalDeriv(self, prob, v, adjoint=False):
|
||||
return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
|
||||
|
||||
def bPrimary(self, prob):
|
||||
return None
|
||||
|
||||
def hPrimary(self, prob):
|
||||
return None
|
||||
|
||||
def ePrimary(self, prob):
|
||||
return None
|
||||
|
||||
def jPrimary(self, prob):
|
||||
return None
|
||||
|
||||
def S_m(self, prob):
|
||||
return None
|
||||
|
||||
def S_e(self, prob):
|
||||
return None
|
||||
|
||||
def S_mDeriv(self, prob, v, adjoint = False):
|
||||
return None
|
||||
|
||||
def S_eDeriv(self, prob, v, adjoint = False):
|
||||
return None
|
||||
|
||||
|
||||
class SrcFDEM_RawVec_e(SrcFDEM):
|
||||
"""
|
||||
RawVec electric source. It is defined by the user provided vector S_e
|
||||
|
||||
:param numpy.array S_e: electric source term
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
|
||||
def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
|
||||
self._S_e = np.array(S_e,dtype=complex)
|
||||
self._ePrimary = ePrimary
|
||||
self._bPrimary = bPrimary
|
||||
self._hPrimary = hPrimary
|
||||
self._jPrimary = jPrimary
|
||||
self.freq = float(freq)
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def S_e(self, prob):
|
||||
return self._S_e
|
||||
|
||||
def ePrimary(self, prob):
|
||||
return self._ePrimary
|
||||
|
||||
def bPrimary(self, prob):
|
||||
return self._bPrimary
|
||||
|
||||
def hPrimary(self, prob):
|
||||
return self._hPrimary
|
||||
|
||||
def jPrimary(self, prob):
|
||||
return self._jPrimary
|
||||
|
||||
|
||||
class SrcFDEM_RawVec_m(SrcFDEM):
|
||||
"""
|
||||
RawVec magnetic source. It is defined by the user provided vector S_m
|
||||
|
||||
:param numpy.array S_m: magnetic source term
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
|
||||
def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
|
||||
self._S_m = np.array(S_m,dtype=complex)
|
||||
self.freq = float(freq)
|
||||
self.integrate = integrate
|
||||
self._ePrimary = np.array(ePrimary,dtype=complex)
|
||||
self._bPrimary = np.array(bPrimary,dtype=complex)
|
||||
self._hPrimary = np.array(hPrimary,dtype=complex)
|
||||
self._jPrimary = np.array(jPrimary,dtype=complex)
|
||||
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def S_m(self, prob):
|
||||
return self._S_m
|
||||
|
||||
def ePrimary(self, prob):
|
||||
return self._ePrimary
|
||||
|
||||
def bPrimary(self, prob):
|
||||
return self._bPrimary
|
||||
|
||||
def hPrimary(self, prob):
|
||||
return self._hPrimary
|
||||
|
||||
def jPrimary(self, prob):
|
||||
return self._jPrimary
|
||||
|
||||
|
||||
class SrcFDEM_RawVec(SrcFDEM):
|
||||
"""
|
||||
RawVec source. It is defined by the user provided vectors S_m, S_e
|
||||
|
||||
:param numpy.array S_m: magnetic source term
|
||||
:param numpy.array S_e: electric source term
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
def __init__(self, rxList, freq, S_m, S_e, integrate = True):
|
||||
self._S_m = np.array(S_m,dtype=complex)
|
||||
self._S_e = np.array(S_e,dtype=complex)
|
||||
self.freq = float(freq)
|
||||
self.integrate = integrate
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def S_m(self, prob):
|
||||
if prob._eqLocs is 'EF' and self.integrate is True:
|
||||
return prob.Me * self._S_m
|
||||
return self._S_m
|
||||
|
||||
def S_e(self, prob):
|
||||
if prob._eqLocs is 'FE' and self.integrate is True:
|
||||
return prob.Me * self._S_e
|
||||
return self._S_e
|
||||
|
||||
|
||||
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., mu = mu_0):
|
||||
self.freq = float(freq)
|
||||
self.loc = loc
|
||||
self.orientation = orientation
|
||||
self.moment = moment
|
||||
self.mu = mu
|
||||
self.integrate = False
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def bPrimary(self, prob):
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
if eqLocs is 'FE':
|
||||
gridX = prob.mesh.gridEx
|
||||
gridY = prob.mesh.gridEy
|
||||
gridZ = prob.mesh.gridEz
|
||||
C = prob.mesh.edgeCurl
|
||||
|
||||
elif eqLocs is 'EF':
|
||||
gridX = prob.mesh.gridFx
|
||||
gridY = prob.mesh.gridFy
|
||||
gridZ = prob.mesh.gridFz
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
|
||||
if prob.mesh._meshType is 'CYL':
|
||||
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, moment=self.moment)
|
||||
|
||||
else:
|
||||
srcfct = SrcUtils.MagneticDipoleVectorPotential
|
||||
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
|
||||
|
||||
def hPrimary(self, prob):
|
||||
b = self.bPrimary(prob)
|
||||
return h_from_b(prob,b)
|
||||
|
||||
def S_m(self, prob):
|
||||
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., 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):
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
if eqLocs is 'FE':
|
||||
gridX = prob.mesh.gridFx
|
||||
gridY = prob.mesh.gridFy
|
||||
gridZ = prob.mesh.gridFz
|
||||
C = prob.mesh.edgeCurl
|
||||
|
||||
elif eqLocs is 'EF':
|
||||
gridX = prob.mesh.gridEx
|
||||
gridY = prob.mesh.gridEy
|
||||
gridZ = prob.mesh.gridEz
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
srcfct = SrcUtils.MagneticDipoleFields
|
||||
if prob.mesh._meshType is 'CYL':
|
||||
if not prob.mesh.isSymmetric:
|
||||
# TODO ?
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
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, 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
|
||||
|
||||
def hPrimary(self, prob):
|
||||
b = self.bPrimary(prob)
|
||||
return h_from_b(prob, 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))
|
||||
|
||||
|
||||
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., mu=mu_0):
|
||||
self.freq = float(freq)
|
||||
self.orientation = orientation
|
||||
self.radius = radius
|
||||
self.mu = mu
|
||||
self.loc = loc
|
||||
self.integrate = False
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def bPrimary(self, prob):
|
||||
eqLocs = prob._eqLocs
|
||||
|
||||
if eqLocs is 'FE':
|
||||
gridX = prob.mesh.gridEx
|
||||
gridY = prob.mesh.gridEy
|
||||
gridZ = prob.mesh.gridEz
|
||||
C = prob.mesh.edgeCurl
|
||||
|
||||
elif eqLocs is 'EF':
|
||||
gridX = prob.mesh.gridFx
|
||||
gridY = prob.mesh.gridFy
|
||||
gridZ = prob.mesh.gridFz
|
||||
C = prob.mesh.edgeCurl.T
|
||||
|
||||
if prob.mesh._meshType is 'CYL':
|
||||
if not prob.mesh.isSymmetric:
|
||||
# TODO ?
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
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, 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 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
|
||||
####################################################
|
||||
|
||||
class SurveyFDEM(Survey.BaseSurvey):
|
||||
class Survey(SimPEG.Survey.BaseSurvey):
|
||||
"""
|
||||
docstring for SurveyFDEM
|
||||
"""
|
||||
|
||||
srcPair = SrcFDEM
|
||||
srcPair = Src.BaseSrc
|
||||
|
||||
def __init__(self, srcList, **kwargs):
|
||||
# Sort these by frequency
|
||||
self.srcList = srcList
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
|
||||
_freqDict = {}
|
||||
for src in srcList:
|
||||
@@ -481,7 +136,7 @@ class SurveyFDEM(Survey.BaseSurvey):
|
||||
return self._freqDict[freq]
|
||||
|
||||
def projectFields(self, u):
|
||||
data = Survey.Data(self)
|
||||
data = SimPEG.Survey.Data(self)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.projectFields(src, self.mesh, u)
|
||||
|
||||
@@ -1,3 +1,3 @@
|
||||
from SurveyFDEM import *
|
||||
from FDEM import BaseFDEMProblem, ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h
|
||||
from SurveyFDEM import Rx, Src, Survey
|
||||
from FDEM import BaseFDEMProblem, Problem_e, Problem_b, Problem_j, Problem_h
|
||||
from FieldsFDEM import *
|
||||
@@ -1,6 +1,6 @@
|
||||
from SimPEG import Solver, Problem
|
||||
from SimPEG.Problem import BaseTimeProblem
|
||||
from SimPEG.EM.Utils import SrcUtils
|
||||
from SimPEG.EM.Utils import *
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.Utils import sdiag, mkvc
|
||||
from SimPEG import Utils, Mesh
|
||||
|
||||
@@ -1,6 +1,6 @@
|
||||
from SimPEG import Utils, Survey, np
|
||||
from SimPEG.Survey import BaseSurvey
|
||||
from SimPEG.EM.Utils import SrcUtils
|
||||
from SimPEG.EM.Utils import *
|
||||
from BaseTDEM import FieldsTDEM
|
||||
|
||||
|
||||
@@ -87,11 +87,11 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
|
||||
"""Vertical magnetic dipole, magnetic vector potential"""
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
|
||||
MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
|
||||
MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
|
||||
@@ -109,11 +109,11 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
|
||||
"""Circular Loop, magnetic vector potential"""
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
|
||||
MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
|
||||
MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
|
||||
else:
|
||||
raise Exception('Unknown mesh for CircularLoop')
|
||||
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
# import Sources
|
||||
# import Ana
|
||||
# import Solver
|
||||
import EMUtils
|
||||
import SrcUtils
|
||||
from EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
|
||||
from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
|
||||
@@ -17,50 +17,50 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
|
||||
|
||||
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)
|
||||
Rx0 = EM.FDEM.Rx(XYZ, comp)
|
||||
|
||||
Src = []
|
||||
|
||||
for SrcType in SrcList:
|
||||
if SrcType is 'MagDipole':
|
||||
Src.append(EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
|
||||
Src.append(EM.FDEM.Src.MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
|
||||
elif SrcType is 'MagDipole_Bfield':
|
||||
Src.append(EM.FDEM.SrcFDEM_MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
|
||||
Src.append(EM.FDEM.Src.MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
|
||||
elif SrcType is 'CircularLoop':
|
||||
Src.append(EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
|
||||
Src.append(EM.FDEM.Src.CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
|
||||
elif SrcType is 'RawVec':
|
||||
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.SrcFDEM_RawVec([Rx0], freq, S_m, S_e))
|
||||
Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
|
||||
|
||||
elif fdemType is 'h' or fdemType is '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.
|
||||
Src.append(EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e))
|
||||
Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
|
||||
|
||||
if verbose:
|
||||
print ' Fetching %s problem' % (fdemType)
|
||||
|
||||
if fdemType == 'e':
|
||||
survey = EM.FDEM.SurveyFDEM(Src)
|
||||
prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
|
||||
survey = EM.FDEM.Survey(Src)
|
||||
prb = EM.FDEM.Problem_e(mesh, mapping=mapping)
|
||||
|
||||
elif fdemType == 'b':
|
||||
survey = EM.FDEM.SurveyFDEM(Src)
|
||||
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
|
||||
survey = EM.FDEM.Survey(Src)
|
||||
prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
|
||||
|
||||
elif fdemType == 'j':
|
||||
survey = EM.FDEM.SurveyFDEM(Src)
|
||||
prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping)
|
||||
survey = EM.FDEM.Survey(Src)
|
||||
prb = EM.FDEM.Problem_j(mesh, mapping=mapping)
|
||||
|
||||
elif fdemType == 'h':
|
||||
survey = EM.FDEM.SurveyFDEM(Src)
|
||||
prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
|
||||
survey = EM.FDEM.Survey(Src)
|
||||
prb = EM.FDEM.Problem_h(mesh, mapping=mapping)
|
||||
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
@@ -28,12 +28,12 @@ class FDEM_analyticTests(unittest.TestCase):
|
||||
|
||||
x = np.linspace(-10,10,5)
|
||||
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
|
||||
rxList = EM.FDEM.RxFDEM(XYZ, 'exi')
|
||||
Src0 = EM.FDEM.SrcFDEM_MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
|
||||
rxList = EM.FDEM.Rx(XYZ, 'exi')
|
||||
Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Src0])
|
||||
survey = EM.FDEM.Survey([Src0])
|
||||
|
||||
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
|
||||
prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
|
||||
prb.pair(survey)
|
||||
|
||||
try:
|
||||
@@ -114,19 +114,19 @@ class FDEM_analyticTests(unittest.TestCase):
|
||||
|
||||
de = np.zeros(mesh.nF,dtype=complex)
|
||||
de[s_ind] = 1./csz
|
||||
de_p = [EM.FDEM.SrcFDEM_RawVec_e([],freq,de/mesh.area)]
|
||||
de_p = [EM.FDEM.Src.RawVec_e([],freq,de/mesh.area)]
|
||||
|
||||
dm_p = [EM.FDEM.SrcFDEM_MagDipole([],freq,src_loc)]
|
||||
dm_p = [EM.FDEM.Src.MagDipole([],freq,src_loc)]
|
||||
|
||||
|
||||
# Pair the problem and survey
|
||||
surveye = EM.FDEM.SurveyFDEM(de_p)
|
||||
surveym = EM.FDEM.SurveyFDEM(dm_p)
|
||||
surveye = EM.FDEM.Survey(de_p)
|
||||
surveym = EM.FDEM.Survey(dm_p)
|
||||
|
||||
mapping = [('sigma', Maps.IdentityMap(mesh)),('mu', Maps.IdentityMap(mesh))]
|
||||
|
||||
prbe = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
|
||||
prbm = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
|
||||
prbe = EM.FDEM.Problem_h(mesh, mapping=mapping)
|
||||
prbm = EM.FDEM.Problem_e(mesh, mapping=mapping)
|
||||
|
||||
prbe.pair(surveye) # pair problem and survey
|
||||
prbm.pair(surveym)
|
||||
|
||||
@@ -3,6 +3,7 @@ from SimPEG import *
|
||||
from SimPEG import EM
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils.testingUtils import getFDEMProblem
|
||||
|
||||
testDerivs = True
|
||||
testEB = True
|
||||
@@ -20,91 +21,9 @@ addrandoms = True
|
||||
SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
|
||||
|
||||
|
||||
def getProblem(fdemType, comp):
|
||||
cs = 5.
|
||||
ncx, ncy, ncz = 6, 6, 6
|
||||
npad = 3
|
||||
hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
|
||||
mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
|
||||
|
||||
mapping = Maps.ExpMap(mesh)
|
||||
|
||||
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)
|
||||
|
||||
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':
|
||||
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':
|
||||
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':
|
||||
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':
|
||||
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)
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
pass
|
||||
|
||||
return prb
|
||||
|
||||
|
||||
def derivTest(fdemType, comp):
|
||||
|
||||
prb = getProblem(fdemType, comp)
|
||||
prb = getFDEMProblem(fdemType, comp, [SrcType], freq)
|
||||
print '%s formulation - %s' % (fdemType, comp)
|
||||
x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
mu = np.log(np.ones(prb.mesh.nC)*MU)
|
||||
|
||||
Reference in New Issue
Block a user