mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-28 16:47:03 +08:00
Lindsey
@@ -19,9 +19,9 @@ model = Model.LogModel(mesh)
|
||||
x = np.linspace(-10,10,5)
|
||||
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
|
||||
rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex')
|
||||
Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, rxList)
|
||||
Src0 = EM.FDEM.SrcFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, rxList)
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Tx0])
|
||||
survey = EM.FDEM.SurveyFDEM([Src0])
|
||||
|
||||
prb = EM.FDEM.ProblemFDEM_b(model)
|
||||
prb.pair(survey)
|
||||
@@ -31,26 +31,26 @@ sigma = np.ones(mesh.nC)*sig
|
||||
sigma[mesh.gridCC[:,2] > 0] = 1e-8
|
||||
m = np.log(sigma)
|
||||
|
||||
skin = 500*np.sqrt(1/(sig*Tx0.freq))
|
||||
skin = 500*np.sqrt(1/(sig*Src0.freq))
|
||||
print 'The skin depth is: %4.2f m' % skin
|
||||
|
||||
prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.spsolve, factorize=False, checkAccuracy=True)
|
||||
|
||||
u = prb.fields(m)
|
||||
|
||||
plt.colorbar(mesh.plotImage(np.log10(np.abs(u[Tx0, 'b'].real)), 'Fz'))
|
||||
plt.colorbar(mesh.plotImage(np.log10(np.abs(u[Src0, 'b'].real)), 'Fz'))
|
||||
|
||||
bfz = mesh.r(u[Tx0, 'b'],'F','Fz','M')
|
||||
bfz = mesh.r(u[Src0, 'b'],'F','Fz','M')
|
||||
|
||||
x = np.linspace(-55,55,12)
|
||||
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
|
||||
|
||||
P = mesh.getInterpolationMat(XYZ, 'Fz')
|
||||
|
||||
an = EM.Utils.Ana.FEM.hzAnalyticDipoleF(x, Tx0.freq, sig)
|
||||
an = EM.Utils.Ana.FEM.hzAnalyticDipoleF(x, Src0.freq, sig)
|
||||
|
||||
plt.figure(2)
|
||||
plt.plot(x,np.log10(np.abs(P*np.imag(u[Tx0, 'b']))))
|
||||
plt.plot(x,np.log10(np.abs(P*np.imag(u[Src0, 'b']))))
|
||||
plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r')
|
||||
plt.xlabel('Distance, m')
|
||||
plt.ylabel('Log10 Response imag($B_z$)')
|
||||
|
||||
+65
-12
@@ -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) )
|
||||
@@ -34,11 +35,14 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True):
|
||||
|
||||
if secondary:
|
||||
hp =-1/(4*np.pi*r**3)
|
||||
return hz-hp
|
||||
hz = hz-hp
|
||||
|
||||
if hz.ndim == 1:
|
||||
hz = Utils.mkvc(hz,2)
|
||||
|
||||
return hz
|
||||
|
||||
def AnalyticMagDipoleWholeSpace(XYZ, txLoc, sig, f, m=1., orientation='X'):
|
||||
def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, moment=1., orientation='X', mu = mu_0):
|
||||
"""
|
||||
Analytical solution for a dipole in a whole-space.
|
||||
|
||||
@@ -67,15 +71,15 @@ def AnalyticMagDipoleWholeSpace(XYZ, txLoc, sig, f, m=1., orientation='X'):
|
||||
|
||||
XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
|
||||
|
||||
dx = XYZ[:,0]-txLoc[0]
|
||||
dy = XYZ[:,1]-txLoc[1]
|
||||
dz = XYZ[:,2]-txLoc[2]
|
||||
dx = XYZ[:,0]-srcLoc[0]
|
||||
dy = XYZ[:,1]-srcLoc[1]
|
||||
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)
|
||||
front = moment / (4.*pi * r**3.) * np.exp(-1j*kr)
|
||||
mid = -kr**2. + 3.*1j*kr + 3.
|
||||
|
||||
if orientation.upper() == 'X':
|
||||
@@ -93,8 +97,57 @@ def AnalyticMagDipoleWholeSpace(XYZ, txLoc, 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)
|
||||
|
||||
if By.ndim is 1:
|
||||
By = Utils.mkvc(By,2)
|
||||
|
||||
if Bz.ndim is 1:
|
||||
Bz = Utils.mkvc(Bz,2)
|
||||
|
||||
return Bx, By, Bz
|
||||
|
||||
|
||||
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]
|
||||
dy = XYZ[:,1]-srcLoc[1]
|
||||
dz = XYZ[:,2]-srcLoc[2]
|
||||
|
||||
r = np.sqrt( dx**2. + dy**2. + dz**2.)
|
||||
k = np.sqrt( -1j*2.*np.pi*f*mu*sig )
|
||||
kr = 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
|
||||
|
||||
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':
|
||||
# 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':
|
||||
# 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
|
||||
@@ -0,0 +1,98 @@
|
||||
from SimPEG import Utils, np
|
||||
from scipy.constants import mu_0, epsilon_0
|
||||
from simpegEM.Utils.EMUtils import k
|
||||
|
||||
def getKc(freq,sigma,a,b,mu=mu_0,eps=epsilon_0):
|
||||
a = float(a)
|
||||
b = float(b)
|
||||
# return 1./(2*np.pi) * np.sqrt(b / a) * np.exp(-1j*k(freq,sigma,mu,eps)*(b-a))
|
||||
return np.sqrt(b / a) * np.exp(-1j*k(freq,sigma,mu,eps)*(b-a))
|
||||
|
||||
def _r2(xyz):
|
||||
return np.sum(xyz**2,1)
|
||||
|
||||
def _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
|
||||
Kc1 = getKc(freq,sigma[1],a,b,mu[1],eps)
|
||||
|
||||
nobs = obsloc.shape[0]
|
||||
dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
|
||||
|
||||
r2 = _r2(dxyz[:,:2])
|
||||
sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2)
|
||||
k2 = k(freq,sigma[2],mu[2],eps)
|
||||
|
||||
return Kc1 * moment / (4.*np.pi) *np.exp(-1j*k2*sqrtr2z2) / sqrtr2z2
|
||||
|
||||
|
||||
def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
|
||||
HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
|
||||
nobs = obsloc.shape[0]
|
||||
dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
|
||||
|
||||
r2 = _r2(dxyz[:,:2])
|
||||
sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2)
|
||||
k2 = k(freq,sigma[2],mu[2],eps)
|
||||
|
||||
return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1./ sqrtr2z2)
|
||||
|
||||
|
||||
def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
|
||||
HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
|
||||
nobs = obsloc.shape[0]
|
||||
dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
|
||||
|
||||
r2z2 = _r2(dxyz)
|
||||
sqrtr2z2 = np.sqrt(r2z2)
|
||||
k2 = k(freq,sigma[2],mu[2],eps)
|
||||
|
||||
return -HertzZ*dxyz[:,2] /sqrtr2z2 * (1j*k2 + 1./sqrtr2z2)
|
||||
|
||||
def _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
|
||||
HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
dHertzZdr = _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
|
||||
nobs = obsloc.shape[0]
|
||||
dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
|
||||
|
||||
r2 = _r2(dxyz[:,:2])
|
||||
r = np.sqrt(r2)
|
||||
z = dxyz[:,2]
|
||||
sqrtr2z2 = np.sqrt(r2 + z**2)
|
||||
k2 = k(freq,sigma[2],mu[2],eps)
|
||||
|
||||
return dHertzZdr*(-z/sqrtr2z2)*(1j*k2+1./sqrtr2z2) + HertzZ*(z*r/sqrtr2z2**3)*(1j*k2 + 2./sqrtr2z2)
|
||||
|
||||
def _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
|
||||
HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
dHertzZdz = _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
|
||||
nobs = obsloc.shape[0]
|
||||
dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
|
||||
|
||||
r2 = _r2(dxyz[:,:2])
|
||||
r = np.sqrt(r2)
|
||||
z = dxyz[:,2]
|
||||
sqrtr2z2 = np.sqrt(r2 + z**2)
|
||||
k2 = k(freq,sigma[2],mu[2],eps)
|
||||
|
||||
return (dHertzZdz*z + HertzZ)/sqrtr2z2*(-1j*k2 - 1./sqrtr2z2) + HertzZ*z/sqrtr2z2**3*(1j*k2*z + 2.*z/sqrtr2z2)
|
||||
|
||||
def getCasingEphiMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
|
||||
return 1j * omega(freq) * mu * _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
|
||||
def getCasingHrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
|
||||
return _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
|
||||
def getCasingHzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
|
||||
d2HertzZdz2 = _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
k2 = k(freq,sigma[2],mu[2],eps)
|
||||
HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
return d2HertzZdz2 + k2**2 * HertzZ
|
||||
|
||||
def getCasingBrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
|
||||
return mu_0 * getCasingHrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
|
||||
def getCasingBzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
|
||||
return mu_0 * getCasingHzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
|
||||
@@ -1,2 +1,3 @@
|
||||
from TDEM import hzAnalyticDipoleT
|
||||
from FDEM import hzAnalyticDipoleF
|
||||
from FDEMcasing import *
|
||||
|
||||
+95
-59
@@ -1,16 +1,24 @@
|
||||
from SimPEG import Survey, Problem, Utils, Models, np, sp, Solver as SimpegSolver
|
||||
from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver
|
||||
from scipy.constants import mu_0
|
||||
|
||||
class EMPropMap(Maps.PropMap):
|
||||
sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap))
|
||||
mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',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))
|
||||
|
||||
|
||||
class BaseEMProblem(Problem.BaseProblem):
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
solType = None
|
||||
storeTheseFields = ['e', 'b']
|
||||
|
||||
surveyPair = Survey.BaseSurvey
|
||||
dataPair = Survey.Data
|
||||
|
||||
PropMap = EMPropMap
|
||||
|
||||
Solver = SimpegSolver
|
||||
solverOpts = {}
|
||||
@@ -26,84 +34,112 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
self.__makeASymmetric = True
|
||||
return self.__makeASymmetric
|
||||
|
||||
####################################################
|
||||
# Mu Model
|
||||
####################################################
|
||||
@property
|
||||
def mu(self):
|
||||
if getattr(self, '_mu', None) is None:
|
||||
self._mu = mu_0
|
||||
return self._mu
|
||||
@mu.setter
|
||||
def mu(self, value):
|
||||
if getattr(self, '_MfMui', None) is not None:
|
||||
del self._MfMui
|
||||
if getattr(self, '_MeMu', None) is not None:
|
||||
del delf._MeMu
|
||||
if getattr(self, '_MeMuI', None) is not None:
|
||||
del self._MeMuI
|
||||
self._mu = value
|
||||
|
||||
|
||||
|
||||
####################################################
|
||||
# Mass Matrices
|
||||
####################################################
|
||||
|
||||
@property
|
||||
def MfMui(self):
|
||||
if getattr(self, '_MfMui', None) is None:
|
||||
self._MfMui = self.mesh.getFaceInnerProduct(1/self.mu)
|
||||
return self._MfMui
|
||||
|
||||
@property
|
||||
def MeMuI(self):
|
||||
# TODO: Assuming isotropic mu
|
||||
if getattr(self, '_MeMuI', None) is None:
|
||||
self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True)
|
||||
return self._MeMuI
|
||||
|
||||
@property
|
||||
def MeMu(self):
|
||||
#TODO: Assuming isotropic mu
|
||||
if getattr(self, '_MeMu', None) is None:
|
||||
self._MeMu = self.mesh.getEdgeInnerProduct(self.mu)
|
||||
return self._MeMu
|
||||
|
||||
def deleteTheseOnModelUpdate(self):
|
||||
toDelete = []
|
||||
if self.mapping.sigmaMap is not None or self.mapping.rhoMap is not None:
|
||||
toDelete += ['_MeSigma', '_MeSigmaI','_MfRho','_MfRhoI']
|
||||
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):
|
||||
if getattr(self, '_Me', None) is None:
|
||||
self._Me = self.mesh.getEdgeInnerProduct()
|
||||
return self._Me
|
||||
|
||||
@property
|
||||
def Mf(self):
|
||||
if getattr(self, '_Mf', None) is None:
|
||||
self._Mf = self.mesh.getFaceInnerProduct()
|
||||
return self._Mf
|
||||
|
||||
|
||||
# ----- Magnetic Permeability ----- #
|
||||
@property
|
||||
def MfMui(self):
|
||||
if getattr(self, '_MfMui', None) is None:
|
||||
self._MfMui = self.mesh.getFaceInnerProduct(self.curModel.mui)
|
||||
return self._MfMui
|
||||
|
||||
@property
|
||||
def MfMuiI(self):
|
||||
if getattr(self, '_MfMuiI', None) is None:
|
||||
self._MfMuiI = self.mesh.getFaceInnerProduct(self.curModel.mui, invMat=True)
|
||||
return self._MfMuiI
|
||||
|
||||
@property
|
||||
def MeMu(self):
|
||||
if getattr(self, '_MeMu', None) is None:
|
||||
self._MeMu = self.mesh.getEdgeInnerProduct(self.curModel.mu)
|
||||
return self._MeMu
|
||||
|
||||
@property
|
||||
def MeMuI(self):
|
||||
if getattr(self, '_MeMuI', None) is None:
|
||||
self._MeMuI = self.mesh.getEdgeInnerProduct(self.curModel.mu, invMat=True)
|
||||
return self._MeMuI
|
||||
|
||||
# ----- Electrical Conductivity ----- #
|
||||
#TODO: hardcoded to sigma as the model
|
||||
@property
|
||||
def MeSigma(self):
|
||||
#TODO: hardcoded to sigma as the model
|
||||
if getattr(self, '_MeSigma', None) is None:
|
||||
sigma = self.curModel.transform
|
||||
self._MeSigma = self.mesh.getEdgeInnerProduct(sigma)
|
||||
self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma)
|
||||
return self._MeSigma
|
||||
|
||||
# TODO: This should take a vector
|
||||
def MeSigmaDeriv(self, u):
|
||||
"""
|
||||
Deriv of MeSigma wrt sigma
|
||||
"""
|
||||
return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
|
||||
|
||||
|
||||
@property
|
||||
def MeSigmaI(self):
|
||||
#TODO: hardcoded to sigma as the model
|
||||
if getattr(self, '_MeSigmaI', None) is None:
|
||||
sigma = self.curModel.transform
|
||||
self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True)
|
||||
self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True)
|
||||
return self._MeSigmaI
|
||||
|
||||
# TODO: This should take a vector
|
||||
def MeSigmaIDeriv(self, u):
|
||||
"""
|
||||
Deriv of MeSigma wrt sigma
|
||||
"""
|
||||
# TODO: only works for diagonal tensors. getEdgeInnerProductDeriv, invMat=True should be implemented in SimPEG
|
||||
|
||||
dMeSigmaI_dI = -self.MeSigmaI**2
|
||||
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u)
|
||||
dsig_dm = self.curModel.sigmaDeriv
|
||||
return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm))
|
||||
# return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u)
|
||||
|
||||
|
||||
@property
|
||||
def MfSigmai(self):
|
||||
#TODO: hardcoded to sigma as the model
|
||||
#TODO: hardcoded to sigma diagonal
|
||||
if getattr(self, '_MfSigmai', None) is None:
|
||||
sigma = self.curModel.transform
|
||||
self._MfSigmai = self.mesh.getFaceInnerProduct(1/sigma)
|
||||
return self._MfSigmai
|
||||
def MfRho(self):
|
||||
if getattr(self, '_MfRho', None) is None:
|
||||
self._MfRho = self.mesh.getFaceInnerProduct(self.curModel.rho)
|
||||
return self._MfRho
|
||||
|
||||
deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai']
|
||||
# TODO: This should take a vector
|
||||
def MfRhoDeriv(self,u):
|
||||
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv)
|
||||
# self.curModel.rhoDeriv
|
||||
|
||||
def fields(self, m):
|
||||
self.curModel = m
|
||||
F = self.forward(m, self.getRHS, self.calcFields)
|
||||
return F
|
||||
@property
|
||||
def MfRhoI(self):
|
||||
if getattr(self, '_MfRhoI', None) is None:
|
||||
self._MfRhoI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True)
|
||||
return self._MfRhoI
|
||||
|
||||
# TODO: This isn't going to work yet
|
||||
# TODO: This should take a vector
|
||||
def MfRhoIDeriv(self,u):
|
||||
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
|
||||
|
||||
@@ -36,8 +36,8 @@ if plotIt:
|
||||
|
||||
rxOffset=1e-3
|
||||
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
|
||||
tx = EM.TDEM.TxTDEM(np.array([0., 0., 80]), 'VMD_MVP', [rx])
|
||||
survey = EM.TDEM.SurveyTDEM([tx])
|
||||
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], np.array([0., 0., 80]))
|
||||
survey = EM.TDEM.SurveyTDEM([src])
|
||||
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
|
||||
prb.Solver = SolverLU
|
||||
|
||||
+328
-327
@@ -1,108 +1,10 @@
|
||||
from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
|
||||
from scipy.constants import mu_0
|
||||
from SurveyFDEM import SurveyFDEM, FieldsFDEM
|
||||
from simpegEM import Sources
|
||||
from SurveyFDEM import SurveyFDEM
|
||||
from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b, FieldsFDEM_h, FieldsFDEM_j
|
||||
from simpegEM.Base import BaseEMProblem
|
||||
from simpegEM.Utils.EMUtils import omega
|
||||
|
||||
def omega(freq):
|
||||
"""Change frequency to angular frequency, omega"""
|
||||
return 2.*np.pi*freq
|
||||
|
||||
def getSource(self,freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
:rtype: numpy.ndarray (nE, nTx)
|
||||
:return: RHS
|
||||
"""
|
||||
Txs = self.survey.getTransmitters(freq)
|
||||
rhs = range(len(Txs))
|
||||
|
||||
solType = self.solType
|
||||
|
||||
if solType == 'e' or solType == 'b':
|
||||
gridEJx = self.mesh.gridEx
|
||||
gridEJy = self.mesh.gridEy
|
||||
gridEJz = self.mesh.gridEz
|
||||
nEJ = self.mesh.nE
|
||||
|
||||
gridBHx = self.mesh.gridFx
|
||||
gridBHy = self.mesh.gridFy
|
||||
gridBHz = self.mesh.gridFz
|
||||
nBH = self.mesh.nF
|
||||
|
||||
|
||||
C = self.mesh.edgeCurl
|
||||
mui = self.MfMui
|
||||
|
||||
elif solType == 'h' or solType == 'j':
|
||||
gridEJx = self.mesh.gridFx
|
||||
gridEJy = self.mesh.gridFy
|
||||
gridEJz = self.mesh.gridFz
|
||||
nEJ = self.mesh.nF
|
||||
|
||||
gridBHx = self.mesh.gridEx
|
||||
gridBHy = self.mesh.gridEy
|
||||
gridBHz = self.mesh.gridEz
|
||||
nBH = self.mesh.nE
|
||||
|
||||
C = self.mesh.edgeCurl.T
|
||||
mui = self.MeMuI
|
||||
|
||||
else:
|
||||
NotImplementedError('Only E or F sources')
|
||||
|
||||
for i, tx in enumerate(Txs):
|
||||
if self.mesh._meshType is 'CYL':
|
||||
if self.mesh.isSymmetric:
|
||||
if tx.txType == 'VMD':
|
||||
SRC = Sources.MagneticDipoleVectorPotential(tx.loc, gridEJy, 'y')
|
||||
elif tx.txType =='CircularLoop':
|
||||
SRC = Sources.MagneticLoopVectorPotential(tx.loc, gridEJy, 'y', tx.radius)
|
||||
else:
|
||||
raise NotImplementedError('Only VMD and CircularLoop')
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
|
||||
elif self.mesh._meshType is 'TENSOR':
|
||||
|
||||
if tx.txType == 'VMD':
|
||||
src = Sources.MagneticDipoleVectorPotential
|
||||
SRCx = src(tx.loc, gridEJx, 'x')
|
||||
SRCy = src(tx.loc, gridEJy, 'y')
|
||||
SRCz = src(tx.loc, gridEJz, 'z')
|
||||
|
||||
elif tx.txType == 'VMD_B':
|
||||
src = Sources.MagneticDipoleFields
|
||||
SRCx = src(tx.loc, gridBHx, 'x')
|
||||
SRCy = src(tx.loc, gridBHy, 'y')
|
||||
SRCz = src(tx.loc, gridBHz, 'z')
|
||||
|
||||
elif tx.txType == 'CircularLoop':
|
||||
src = Sources.MagneticLoopVectorPotential
|
||||
SRCx = src(tx.loc, gridEJx, 'x', tx.radius)
|
||||
SRCy = src(tx.loc, gridEJy, 'y', tx.radius)
|
||||
SRCz = src(tx.loc, gridEJz, 'z', tx.radius)
|
||||
else:
|
||||
|
||||
raise NotImplemented('%s txType is not implemented' % tx.txType)
|
||||
SRC = np.concatenate((SRCx, SRCy, SRCz))
|
||||
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
|
||||
rhs[i] = SRC
|
||||
|
||||
# b-forumlation
|
||||
if tx.txType == 'VMD_B':
|
||||
b_0 = np.concatenate(rhs).reshape((nBH, len(Txs)), order='E')
|
||||
else:
|
||||
a = np.concatenate(rhs).reshape((nEJ, len(Txs)), order='F')
|
||||
b_0 = C*a
|
||||
|
||||
if solType == 'b' or solType == 'h':
|
||||
return b_0
|
||||
elif solType == 'e' or solType == 'j':
|
||||
return C.T*mui*b_0
|
||||
|
||||
class BaseFDEMProblem(BaseEMProblem):
|
||||
"""
|
||||
@@ -110,58 +12,71 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
|
||||
.. math::
|
||||
|
||||
\\nabla \\times \\vec{E} + i \\omega \\vec{B} = 0 \\\\
|
||||
\\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{J_s}
|
||||
\\nabla \\times \\vec{E} + i \\omega \\vec{B} = \\vec{S_m} \\\\
|
||||
\\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{S_e}
|
||||
|
||||
"""
|
||||
surveyPair = SurveyFDEM
|
||||
fieldsPair = FieldsFDEM
|
||||
|
||||
def forward(self, m, RHS, CalcFields):
|
||||
|
||||
F = FieldsFDEM(self.mesh, self.survey)
|
||||
def fields(self, m=None):
|
||||
self.curModel = m
|
||||
F = self.fieldsPair(self.mesh, self.survey)
|
||||
|
||||
for freq in self.survey.freqs:
|
||||
A = self.getA(freq)
|
||||
rhs = RHS(freq)
|
||||
rhs = self.getRHS(freq)
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
sol = Ainv * rhs
|
||||
for fieldType in self.storeTheseFields:
|
||||
Txs = self.survey.getTransmitters(freq)
|
||||
F[Txs, fieldType] = CalcFields(sol, freq, fieldType)
|
||||
Srcs = self.survey.getSrcByFreq(freq)
|
||||
ftype = self._fieldType + 'Solution'
|
||||
F[Srcs, ftype] = sol
|
||||
|
||||
return F
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
def Jvec(self, m, v, f=None):
|
||||
if f is None:
|
||||
f = self.fields(m)
|
||||
|
||||
self.curModel = m
|
||||
|
||||
Jv = self.dataPair(self.survey)
|
||||
|
||||
for freq in self.survey.freqs:
|
||||
A = self.getA(freq)
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
dA_du = self.getA(freq) #
|
||||
dA_duI = self.Solver(dA_du, **self.solverOpts)
|
||||
|
||||
for tx in self.survey.getTransmitters(freq):
|
||||
u_tx = u[tx, self.solType]
|
||||
w = self.getADeriv(freq, u_tx, v)
|
||||
Ainvw = Ainv * w
|
||||
for rx in tx.rxList:
|
||||
fAinvw = self.calcFields(Ainvw, freq, rx.projField)
|
||||
P = lambda v: rx.projectFieldsDeriv(tx, self.mesh, u, v)
|
||||
for src in self.survey.getSrcByFreq(freq):
|
||||
ftype = self._fieldType + 'Solution'
|
||||
u_src = f[src, ftype]
|
||||
dA_dm = self.getADeriv_m(freq, u_src, v)
|
||||
dRHS_dm = self.getRHSDeriv_m(src, v)
|
||||
if dRHS_dm is None:
|
||||
du_dm = dA_duI * ( - dA_dm )
|
||||
else:
|
||||
du_dm = dA_duI * ( - dA_dm + dRHS_dm )
|
||||
for rx in src.rxList:
|
||||
# df_duFun = u.deriv_u(rx.fieldsUsed, m)
|
||||
df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
|
||||
df_du = df_duFun(src, du_dm, adjoint=False)
|
||||
if df_du is not None:
|
||||
du_dm = df_du
|
||||
|
||||
df_dm = self.calcFieldsDeriv(u_tx, freq, rx.projField, v)
|
||||
if df_dm is None:
|
||||
Jv[tx, rx] = - P(fAinvw)
|
||||
else:
|
||||
Jv[tx, rx] = - P(fAinvw) + P(df_dm)
|
||||
df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
|
||||
df_dm = df_dmFun(src, v, adjoint=False)
|
||||
if df_dm is not None:
|
||||
du_dm += df_dm
|
||||
|
||||
P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m
|
||||
|
||||
|
||||
Jv[src, rx] = P(du_dm)
|
||||
|
||||
return Utils.mkvc(Jv)
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
def Jtvec(self, m, v, f=None):
|
||||
if f is None:
|
||||
f = self.fields(m)
|
||||
|
||||
self.curModel = m
|
||||
|
||||
@@ -169,39 +84,73 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
if not isinstance(v, self.dataPair):
|
||||
v = self.dataPair(self.survey, v)
|
||||
|
||||
Jtv = np.zeros(self.mapping.nP)
|
||||
Jtv = np.zeros(m.size)
|
||||
|
||||
for freq in self.survey.freqs:
|
||||
AT = self.getA(freq).T
|
||||
ATinv = self.Solver(AT, **self.solverOpts)
|
||||
|
||||
for tx in self.survey.getTransmitters(freq):
|
||||
u_tx = u[tx, self.solType]
|
||||
for src in self.survey.getSrcByFreq(freq):
|
||||
ftype = self._fieldType + 'Solution'
|
||||
u_src = f[src, ftype]
|
||||
|
||||
for rx in tx.rxList:
|
||||
PTv = rx.projectFieldsDeriv(tx, self.mesh, u, v[tx, rx], adjoint=True)
|
||||
fPTv = self.calcFields(PTv, freq, rx.projField, adjoint=True)
|
||||
for rx in src.rxList:
|
||||
PTv = rx.projectFieldsDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
|
||||
|
||||
w = ATinv * fPTv
|
||||
Jtv_rx = - self.getADeriv(freq, u_tx, w, adjoint=True)
|
||||
df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
|
||||
df_duT = df_duTFun(src, PTv, adjoint=True)
|
||||
if df_duT is not None:
|
||||
dA_duIT = ATinv * df_duT
|
||||
else:
|
||||
dA_duIT = ATinv * PTv
|
||||
|
||||
df_dm = self.calcFieldsDeriv(u_tx, freq, rx.projField, PTv, adjoint=True)
|
||||
dA_dmT = self.getADeriv_m(freq, u_src, dA_duIT, adjoint=True)
|
||||
|
||||
if df_dm is not None:
|
||||
Jtv_rx += df_dm
|
||||
dRHS_dmT = self.getRHSDeriv_m(src, dA_duIT, adjoint=True)
|
||||
|
||||
if dRHS_dmT is None:
|
||||
du_dmT = - dA_dmT
|
||||
else:
|
||||
du_dmT = -dA_dmT + dRHS_dmT
|
||||
|
||||
df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
|
||||
dfT_dm = df_dmFun(src, PTv, adjoint=True)
|
||||
if dfT_dm is not None:
|
||||
du_dmT += dfT_dm
|
||||
|
||||
real_or_imag = rx.projComp
|
||||
if real_or_imag == 'real':
|
||||
Jtv += Jtv_rx.real
|
||||
Jtv += du_dmT.real
|
||||
elif real_or_imag == 'imag':
|
||||
Jtv += - Jtv_rx.real
|
||||
Jtv += - du_dmT.real
|
||||
else:
|
||||
raise Exception('Must be real or imag')
|
||||
|
||||
return Jtv
|
||||
|
||||
def getSource(self,freq):
|
||||
return self.getSource(freq)
|
||||
def getSourceTerm(self, freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
:rtype: numpy.ndarray (nE or nF, nSrc)
|
||||
:return: RHS
|
||||
"""
|
||||
Srcs = self.survey.getSrcByFreq(freq)
|
||||
if self._eqLocs is 'FE':
|
||||
S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
|
||||
S_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
|
||||
elif self._eqLocs is 'EF':
|
||||
S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
|
||||
S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
|
||||
|
||||
for i, src in enumerate(Srcs):
|
||||
smi, sei = src.eval(self)
|
||||
if smi is not None:
|
||||
S_m[:,i] = Utils.mkvc(smi)
|
||||
if sei is not None:
|
||||
S_e[:,i] = Utils.mkvc(sei)
|
||||
|
||||
return S_m, S_e
|
||||
|
||||
|
||||
##########################################################################################
|
||||
################################ E-B Formulation #########################################
|
||||
@@ -225,10 +174,13 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
|
||||
|
||||
"""
|
||||
solType = 'e'
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
_fieldType = 'e'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = FieldsFDEM_e
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
def getA(self, freq):
|
||||
"""
|
||||
@@ -236,62 +188,77 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: A
|
||||
"""
|
||||
mui = self.MfMui
|
||||
sig = self.MeSigma
|
||||
MfMui = self.MfMui
|
||||
MeSigma = self.MeSigma
|
||||
C = self.mesh.edgeCurl
|
||||
|
||||
return C.T*mui*C + 1j*omega(freq)*sig
|
||||
return C.T*MfMui*C + 1j*omega(freq)*MeSigma
|
||||
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
sig = self.curModel.transform
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(u)
|
||||
def getADeriv_m(self, freq, u, v, adjoint=False):
|
||||
dsig_dm = self.curModel.sigmaDeriv
|
||||
dMe_dsig = self.MeSigmaDeriv(u)
|
||||
|
||||
if adjoint:
|
||||
return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
|
||||
return 1j * omega(freq) * ( dMe_dsig.T * v )
|
||||
|
||||
return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
|
||||
return 1j * omega(freq) * ( dMe_dsig * v )
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
:rtype: numpy.ndarray (nE, nTx)
|
||||
:rtype: numpy.ndarray (nE, nSrc)
|
||||
:return: RHS
|
||||
"""
|
||||
|
||||
j_s = getSource(self,freq)
|
||||
return -1j*omega(freq)*j_s
|
||||
S_m, S_e = self.getSourceTerm(freq)
|
||||
C = self.mesh.edgeCurl
|
||||
MfMui = self.MfMui
|
||||
|
||||
def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
e = sol
|
||||
if fieldType == 'e':
|
||||
return e
|
||||
elif fieldType == 'b':
|
||||
if not adjoint:
|
||||
b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl * e )
|
||||
RHS = C.T * (MfMui * S_m) -1j*omega(freq)*S_e
|
||||
|
||||
return RHS
|
||||
|
||||
def getRHSDeriv_m(self, src, v, adjoint=False):
|
||||
C = self.mesh.edgeCurl
|
||||
MfMui = self.MfMui
|
||||
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
|
||||
|
||||
if adjoint:
|
||||
dRHS = MfMui * (C * v)
|
||||
S_mDerivv = S_mDeriv(dRHS)
|
||||
S_eDerivv = S_eDeriv(v)
|
||||
if S_mDerivv is not None and S_eDerivv is not None:
|
||||
return S_mDerivv - 1j*omega(freq)*S_eDerivv
|
||||
elif S_mDerivv is not None:
|
||||
return S_mDerivv
|
||||
elif S_eDerivv is not None:
|
||||
return - 1j*omega(freq)*S_eDerivv
|
||||
else:
|
||||
b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl.T * e )
|
||||
return b
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
return None
|
||||
else:
|
||||
S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v)
|
||||
|
||||
def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
e = sol
|
||||
if fieldType == 'e':
|
||||
return None
|
||||
elif fieldType == 'b':
|
||||
return None
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
if S_mDerivv is not None and S_eDerivv is not None:
|
||||
return C.T * (MfMui * S_mDerivv) -1j*omega(freq)*S_eDerivv
|
||||
elif S_mDerivv is not None:
|
||||
return C.T * (MfMui * S_mDerivv)
|
||||
elif S_eDerivv is not None:
|
||||
return -1j*omega(freq)*S_eDerivv
|
||||
else:
|
||||
return None
|
||||
|
||||
|
||||
class ProblemFDEM_b(BaseFDEMProblem):
|
||||
"""
|
||||
Solving for b!
|
||||
"""
|
||||
solType = 'b'
|
||||
_fieldType = 'b'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = FieldsFDEM_b
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
def getA(self, freq):
|
||||
"""
|
||||
@@ -299,87 +266,100 @@ class ProblemFDEM_b(BaseFDEMProblem):
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: A
|
||||
"""
|
||||
mui = self.MfMui
|
||||
sigI = self.MeSigmaI
|
||||
MfMui = self.MfMui
|
||||
MeSigmaI = self.MeSigmaI
|
||||
C = self.mesh.edgeCurl
|
||||
iomega = 1j * omega(freq) * sp.eye(self.mesh.nF)
|
||||
|
||||
A = C*sigI*C.T*mui + iomega
|
||||
A = C*MeSigmaI*C.T*MfMui + iomega
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return mui.T*A
|
||||
return A
|
||||
return MfMui.T*A
|
||||
return A
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
def getADeriv_m(self, freq, u, v, adjoint=False):
|
||||
|
||||
mui = self.MfMui
|
||||
MfMui = self.MfMui
|
||||
C = self.mesh.edgeCurl
|
||||
sig = self.curModel.transform
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
#TODO: This only works if diagonal (no tensors)...
|
||||
dMeSigmaI_dI = - self.MeSigmaI**2
|
||||
MeSigmaIDeriv = self.MeSigmaIDeriv
|
||||
vec = C.T*(MfMui*u)
|
||||
|
||||
vec = (C.T*(mui*u))
|
||||
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec)
|
||||
MeSigmaIDeriv = MeSigmaIDeriv(vec)
|
||||
|
||||
if adjoint:
|
||||
if self._makeASymmetric is True:
|
||||
v = mui * v
|
||||
return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * v ) ) )
|
||||
v = MfMui * v
|
||||
return MeSigmaIDeriv.T * (C.T * v)
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return mui.T * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) )
|
||||
return C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) )
|
||||
return MfMui.T * ( C * ( MeSigmaIDeriv * v ) )
|
||||
return C * ( MeSigmaIDeriv * v )
|
||||
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
:rtype: numpy.ndarray (nE, nTx)
|
||||
:rtype: numpy.ndarray (nE, nSrc)
|
||||
:return: RHS
|
||||
"""
|
||||
|
||||
b_0 = getSource(self,freq)
|
||||
|
||||
rhs = -1j*omega(freq)*b_0
|
||||
S_m, S_e = self.getSourceTerm(freq)
|
||||
C = self.mesh.edgeCurl
|
||||
MeSigmaI = self.MeSigmaI
|
||||
|
||||
RHS = S_m + C * ( MeSigmaI * S_e )
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
mui = self.MfMui
|
||||
return mui.T*rhs
|
||||
return rhs
|
||||
MfMui = self.MfMui
|
||||
return MfMui.T*RHS
|
||||
|
||||
def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
b = sol
|
||||
if fieldType == 'e':
|
||||
return RHS
|
||||
|
||||
def getRHSDeriv_m(self, src, v, adjoint=False):
|
||||
C = self.mesh.edgeCurl
|
||||
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(Utils.mkvc(S_e))
|
||||
if not adjoint:
|
||||
e = self.MeSigmaI * ( self.mesh.edgeCurl.T * ( self.MfMui * b ) )
|
||||
else:
|
||||
e = self.MfMui.T * ( self.mesh.edgeCurl * ( self.MeSigmaI.T * b ) )
|
||||
return e
|
||||
elif fieldType == 'b':
|
||||
return b
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
RHSderiv = C * (MeSigmaIDeriv * v)
|
||||
elif adjoint:
|
||||
RHSderiv = MeSigmaIDeriv.T * (C.T * v)
|
||||
else:
|
||||
RHSderiv = None
|
||||
|
||||
def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
b = sol
|
||||
if fieldType == 'e':
|
||||
sig = self.curModel.transform
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
|
||||
C = self.mesh.edgeCurl
|
||||
mui = self.MfMui
|
||||
|
||||
#TODO: This only works if diagonal (no tensors)...
|
||||
dMeSigmaI_dI = - self.MeSigmaI**2
|
||||
|
||||
vec = C.T * ( mui * b )
|
||||
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec)
|
||||
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
|
||||
S_mDeriv, S_eDeriv = S_mDeriv(v), S_eDeriv(v)
|
||||
if S_mDeriv is not None and S_eDeriv is not None:
|
||||
if not adjoint:
|
||||
return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) )
|
||||
else:
|
||||
return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * v ) )
|
||||
elif fieldType == 'b':
|
||||
return None
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
SrcDeriv = S_mDeriv + C * (self.MeSigmaI * S_eDeriv)
|
||||
elif adjoint:
|
||||
SrcDeriv = S_mDeriv + Self.MeSigmaI.T * ( C.T * S_eDeriv)
|
||||
elif S_mDeriv is not None:
|
||||
SrcDeriv = S_mDeriv
|
||||
elif S_eDeriv is not None:
|
||||
if not adjoint:
|
||||
SrcDeriv = C * (self.MeSigmaI * S_eDeriv)
|
||||
elif adjoint:
|
||||
SrcDeriv = self.MeSigmaI.T * ( C.T * S_eDeriv)
|
||||
else:
|
||||
SrcDeriv = None
|
||||
|
||||
if RHSderiv is not None and SrcDeriv is not None:
|
||||
RHSderiv += SrcDeriv
|
||||
elif SrcDeriv is not None:
|
||||
RHSderiv = SrcDeriv
|
||||
|
||||
if RHSderiv is not None:
|
||||
if self._makeASymmetric is True and not adjoint:
|
||||
return MfMui.T * RHSderiv
|
||||
|
||||
return RHSderiv
|
||||
|
||||
|
||||
|
||||
##########################################################################################
|
||||
@@ -401,7 +381,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
|
||||
.. math::
|
||||
\\nabla \\times ( \\mu^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} ) + i\\omega \\vec{J} = - i\\omega\\vec{J_s}
|
||||
|
||||
|
||||
We discretize this to:
|
||||
|
||||
.. math::
|
||||
@@ -412,15 +392,16 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
|
||||
"""
|
||||
|
||||
solType = 'j'
|
||||
storeTheseFields = ['j','h']
|
||||
_fieldType = 'j'
|
||||
_eqLocs = 'EF'
|
||||
fieldsPair = FieldsFDEM_j
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
def getA(self, freq):
|
||||
"""
|
||||
Here, we form the operator \(\\mathbf{A}\) to solce
|
||||
Here, we form the operator \(\\mathbf{A}\) to solce
|
||||
.. math::
|
||||
\\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega
|
||||
|
||||
@@ -430,96 +411,98 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
"""
|
||||
|
||||
MeMuI = self.MeMuI
|
||||
MfSigi = self.MfSigmai
|
||||
MfRho = self.MfRho
|
||||
C = self.mesh.edgeCurl
|
||||
iomega = 1j * omega(freq) * sp.eye(self.mesh.nF)
|
||||
|
||||
A = C * MeMuI * C.T * MfSigi + iomega
|
||||
A = C * MeMuI * C.T * MfRho + iomega
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return MfSigi.T*A
|
||||
return A
|
||||
return MfRho.T*A
|
||||
return A
|
||||
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
def getADeriv_m(self, freq, u, v, adjoint=False):
|
||||
"""
|
||||
In this case, we assume that electrical conductivity, \(\\sigma\) is the physical property of interest (i.e. \(\sigma\) = model.transform). Then we want
|
||||
.. math::
|
||||
.. math::
|
||||
\\frac{\mathbf{A(\\sigma)} \mathbf{v}}{d \\mathbf{m}} &= \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{m}}
|
||||
&= \\mathbf{C} \\mathbf{M^e_{mu}^{-1}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{\\sigma^{-1}}} \\frac{d \\mathbf{\\sigma^{-1}}}{d \\mathbf{\\sigma}} \\frac{d \\mathbf{\\sigma}}{d \\mathbf{m}}
|
||||
"""
|
||||
|
||||
MeMuI = self.MeMuI
|
||||
MfSigi = self.MfSigmai
|
||||
MfRho = self.MfRho
|
||||
C = self.mesh.edgeCurl
|
||||
sig = self.curModel.transform
|
||||
sigi = 1/sig
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u)
|
||||
MfRhoDeriv_m = self.MfRhoDeriv(u)
|
||||
|
||||
if adjoint:
|
||||
if self._makeASymmetric is True:
|
||||
v = MfSigi * v
|
||||
return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) )
|
||||
v = MfRho * v
|
||||
return MfRhoDeriv_m.T * (C * (MeMuI.T * (C.T * v)))
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return MfSigi.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) )
|
||||
return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) )
|
||||
if self._makeASymmetric is True:
|
||||
return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv_m * v) )))
|
||||
return C * (MeMuI * (C.T * (MfRhoDeriv_m * v)))
|
||||
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
:rtype: numpy.ndarray (nE, nTx)
|
||||
:rtype: numpy.ndarray (nE, nSrc)
|
||||
:return: RHS
|
||||
"""
|
||||
j_s = getSource(self,freq)
|
||||
rhs = -1j*omega(freq)*j_s
|
||||
|
||||
S_m, S_e = self.getSourceTerm(freq)
|
||||
C = self.mesh.edgeCurl
|
||||
MeMuI = self.MeMuI
|
||||
|
||||
|
||||
RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e
|
||||
if self._makeASymmetric is True:
|
||||
MfSigi = self.MfSigmai
|
||||
return MfSigi.T*rhs
|
||||
return rhs
|
||||
MfRho = self.MfRho
|
||||
return MfRho.T*RHS
|
||||
|
||||
def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
j = sol
|
||||
if fieldType == 'j':
|
||||
return j
|
||||
elif fieldType == 'h':
|
||||
MeMuI = self.MeMuI
|
||||
C = self.mesh.edgeCurl
|
||||
MfSigi = self.MfSigmai
|
||||
if not adjoint:
|
||||
h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfSigi * j ) )
|
||||
return RHS
|
||||
|
||||
def getRHSDeriv_m(self, src, v, adjoint=False):
|
||||
C = self.mesh.edgeCurl
|
||||
MeMuI = self.MeMuI
|
||||
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
|
||||
|
||||
if adjoint:
|
||||
if self._makeASymmetric:
|
||||
MfRho = self.MfRho
|
||||
v = MfRho*v
|
||||
S_mDerivv = S_mDeriv(MeMuI.T * (C.T * v))
|
||||
S_eDerivv = S_eDeriv(v)
|
||||
if S_mDerivv is not None and S_eDerivv is not None:
|
||||
return S_mDerivv - 1j*omega(freq)*S_eDerivv
|
||||
elif S_mDerivv is not None:
|
||||
return S_mDerivv
|
||||
elif S_eDerivv is not None:
|
||||
return - 1j*omega(freq)*S_eDerivv
|
||||
else:
|
||||
h = -(1./(1j*omega(freq))) * MfSigi.T * ( C * ( MeMuI.T * j ) )
|
||||
return h
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
return None
|
||||
else:
|
||||
S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v)
|
||||
|
||||
def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
j = sol
|
||||
if fieldType == 'j':
|
||||
return None
|
||||
elif fieldType == 'h':
|
||||
MeMuI = self.MeMuI
|
||||
C = self.mesh.edgeCurl
|
||||
sig = self.curModel.transform
|
||||
sigi = 1/sig
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j)
|
||||
sigi = self.MfSigmai
|
||||
if not adjoint:
|
||||
return -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) )
|
||||
if S_mDerivv is not None and S_eDerivv is not None:
|
||||
RHSDeriv = C * (MeMuI * S_mDerivv) - 1j * omega(freq) * S_eDerivv
|
||||
elif S_mDerivv is not None:
|
||||
RHSDeriv = C * (MeMuI * S_mDerivv)
|
||||
elif S_eDerivv is not None:
|
||||
RHSDeriv = - 1j * omega(freq) * S_eDerivv
|
||||
else:
|
||||
return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMuI.T * v ) ) ) )
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
return None
|
||||
|
||||
if self._makeASymmetric:
|
||||
MfRho = self.MfRho
|
||||
return MfRho.T * RHSDeriv
|
||||
return RHSDeriv
|
||||
|
||||
|
||||
|
||||
|
||||
# Solving for h! - using primary- secondary approach
|
||||
class ProblemFDEM_h(BaseFDEMProblem):
|
||||
"""
|
||||
Using the H-J formulation of Maxwell's equations
|
||||
@@ -545,11 +528,12 @@ class ProblemFDEM_h(BaseFDEMProblem):
|
||||
|
||||
"""
|
||||
|
||||
solType = 'h'
|
||||
storeTheseFields = ['j','h']
|
||||
_fieldType = 'h'
|
||||
_eqLocs = 'EF'
|
||||
fieldsPair = FieldsFDEM_h
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
def getA(self, freq):
|
||||
"""
|
||||
@@ -559,47 +543,64 @@ class ProblemFDEM_h(BaseFDEMProblem):
|
||||
"""
|
||||
|
||||
MeMu = self.MeMu
|
||||
MfSigi = self.MfSigmai
|
||||
MfRho = self.MfRho
|
||||
C = self.mesh.edgeCurl
|
||||
|
||||
return C.T * MfSigi * C + 1j*omega(freq)*MeMu
|
||||
return C.T * MfRho * C + 1j*omega(freq)*MeMu
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
def getADeriv_m(self, freq, u, v, adjoint=False):
|
||||
|
||||
MeMu = self.MeMu
|
||||
C = self.mesh.edgeCurl
|
||||
sig = self.curModel.transform
|
||||
sigi = 1/sig
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
|
||||
dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(C*u)
|
||||
MfRhoDeriv_m = self.MfRhoDeriv(C*u)
|
||||
|
||||
if adjoint:
|
||||
return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v))))
|
||||
return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v))))
|
||||
|
||||
|
||||
return MfRhoDeriv_m.T * (C * v)
|
||||
return C.T * (MfRhoDeriv_m * v)
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
:param float freq: Frequency
|
||||
:rtype: numpy.ndarray (nE, nTx)
|
||||
:rtype: numpy.ndarray (nE, nSrc)
|
||||
:return: RHS
|
||||
"""
|
||||
b_0 = getSource(self,freq)
|
||||
return -1j*omega(freq)*b_0
|
||||
|
||||
def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
h = sol
|
||||
if fieldType == 'j':
|
||||
C = self.mesh.edgeCurl
|
||||
if adjoint:
|
||||
return C.T*h
|
||||
return C*h
|
||||
elif fieldType == 'h':
|
||||
return h
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
return None
|
||||
S_m, S_e = self.getSourceTerm(freq)
|
||||
C = self.mesh.edgeCurl
|
||||
MfRho = self.MfRho
|
||||
|
||||
RHS = S_m + C.T * ( MfRho * S_e )
|
||||
|
||||
return RHS
|
||||
|
||||
def getRHSDeriv_m(self, src, v, adjoint=False):
|
||||
_, S_e = src.eval(self)
|
||||
C = self.mesh.edgeCurl
|
||||
MfRho = self.MfRho
|
||||
|
||||
RHSDeriv = None
|
||||
|
||||
if S_e is not None:
|
||||
MfRhoDeriv = self.MfRhoDeriv(S_e)
|
||||
if not adjoint:
|
||||
RHSDeriv = C.T * (MfRhoDeriv * v)
|
||||
elif adjoint:
|
||||
RHSDeriv = MfRhoDeriv.T * (C * v)
|
||||
|
||||
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
|
||||
|
||||
S_mDeriv = S_mDeriv(v)
|
||||
S_eDeriv = S_eDeriv(v)
|
||||
if S_mDeriv is not None:
|
||||
if RHSDeriv is not None:
|
||||
RHSDeriv += S_mDeriv(v)
|
||||
else:
|
||||
RHSDeriv = S_mDeriv(v)
|
||||
if S_eDeriv is not None:
|
||||
if RHSDeriv is not None:
|
||||
RHSDeriv += C.T * (MfRho * S_e)
|
||||
else:
|
||||
RHSDeriv = C.T * (MfRho * S_e)
|
||||
|
||||
return RHSDeriv
|
||||
|
||||
|
||||
@@ -0,0 +1,373 @@
|
||||
from SimPEG import Survey, Problem, Utils, np, sp
|
||||
from simpegEM.Utils.EMUtils import omega
|
||||
|
||||
|
||||
class FieldsFDEM(Problem.Fields):
|
||||
"""Fancy Field Storage for a FDEM survey."""
|
||||
knownFields = {}
|
||||
dtype = complex
|
||||
|
||||
class FieldsFDEM_e(FieldsFDEM):
|
||||
knownFields = {'eSolution':'E'}
|
||||
aliasFields = {
|
||||
'e' : ['eSolution','E','_e'],
|
||||
'ePrimary' : ['eSolution','E','_ePrimary'],
|
||||
'eSecondary' : ['eSolution','E','_eSecondary'],
|
||||
'b' : ['eSolution','F','_b'],
|
||||
'bPrimary' : ['eSolution','F','_bPrimary'],
|
||||
'bSecondary' : ['eSolution','F','_bSecondary']
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
|
||||
|
||||
def startup(self):
|
||||
self.prob = self.survey.prob
|
||||
self._edgeCurl = self.survey.prob.mesh.edgeCurl
|
||||
|
||||
def _ePrimary(self, eSolution, srcList):
|
||||
ePrimary = np.zeros_like(eSolution)
|
||||
for i, src in enumerate(srcList):
|
||||
ep = src.ePrimary(self.prob)
|
||||
if ep is not None:
|
||||
ePrimary[:,i] = ep
|
||||
return ePrimary
|
||||
|
||||
def _eSecondary(self, eSolution, srcList):
|
||||
return eSolution
|
||||
|
||||
def _e(self, eSolution, srcList):
|
||||
return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
|
||||
|
||||
def _eDeriv_u(self, src, v, adjoint = False):
|
||||
return None
|
||||
|
||||
def _eDeriv_m(self, src, v, adjoint = False):
|
||||
# assuming primary does not depend on the model
|
||||
return None
|
||||
|
||||
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.prob)
|
||||
if bp is not None:
|
||||
bPrimary[:,i] += bp
|
||||
return bPrimary
|
||||
|
||||
def _bSecondary(self, eSolution, srcList):
|
||||
C = self._edgeCurl
|
||||
b = (C * eSolution)
|
||||
for i, src in enumerate(srcList):
|
||||
b[:,i] *= - 1./(1j*omega(src.freq))
|
||||
S_m, _ = src.eval(self.prob)
|
||||
if S_m is not None:
|
||||
b[:,i] += 1./(1j*omega(src.freq)) * S_m
|
||||
return b
|
||||
|
||||
def _bSecondaryDeriv_u(self, src, v, adjoint = False):
|
||||
C = self._edgeCurl
|
||||
if adjoint:
|
||||
return - 1./(1j*omega(src.freq)) * (C.T * v)
|
||||
return - 1./(1j*omega(src.freq)) * (C * v)
|
||||
|
||||
def _bSecondaryDeriv_m(self, src, v, adjoint = False):
|
||||
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
|
||||
return None
|
||||
|
||||
def _b(self, eSolution, srcList):
|
||||
return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
|
||||
|
||||
def _bDeriv_u(self, src, v, adjoint=False):
|
||||
# Primary does not depend on u
|
||||
return self._bSecondaryDeriv_u(src, v, adjoint)
|
||||
|
||||
def _bDeriv_m(self, src, v, adjoint=False):
|
||||
# Assuming the primary does not depend on the model
|
||||
return self._bSecondaryDeriv_m(src, v, adjoint)
|
||||
|
||||
|
||||
class FieldsFDEM_b(FieldsFDEM):
|
||||
knownFields = {'bSolution':'F'}
|
||||
aliasFields = {
|
||||
'b' : ['bSolution','F','_b'],
|
||||
'bPrimary' : ['bSolution','F','_bPrimary'],
|
||||
'bSecondary' : ['bSolution','F','_bSecondary'],
|
||||
'e' : ['bSolution','E','_e'],
|
||||
'ePrimary' : ['bSolution','E','_ePrimary'],
|
||||
'eSecondary' : ['bSolution','E','_eSecondary'],
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
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
|
||||
self._MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
|
||||
|
||||
def _bPrimary(self, bSolution, srcList):
|
||||
bPrimary = np.zeros_like(bSolution)
|
||||
for i, src in enumerate(srcList):
|
||||
bp = src.bPrimary(self.prob)
|
||||
if bp is not None:
|
||||
bPrimary[:,i] = bp
|
||||
return bPrimary
|
||||
|
||||
def _bSecondary(self, bSolution, srcList):
|
||||
return bSolution
|
||||
|
||||
def _b(self, bSolution, srcList):
|
||||
return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList)
|
||||
|
||||
def _bDeriv_u(self, src, v, adjoint=False):
|
||||
return None
|
||||
|
||||
def _bDeriv_m(self, src, v, adjoint=False):
|
||||
# assuming primary does not depend on the model
|
||||
return None
|
||||
|
||||
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.prob)
|
||||
if ep is not None:
|
||||
ePrimary[:,i] = ep
|
||||
return ePrimary
|
||||
|
||||
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.prob)
|
||||
if S_e is not None:
|
||||
e[:,i] += -self._MeSigmaI*S_e
|
||||
return e
|
||||
|
||||
def _eSecondaryDeriv_u(self, src, v, adjoint=False):
|
||||
if not adjoint:
|
||||
return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * v) )
|
||||
else:
|
||||
return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * v))
|
||||
|
||||
def _eSecondaryDeriv_m(self, src, v, adjoint=False):
|
||||
bSolution = self[[src],'bSolution']
|
||||
_,S_e = src.eval(self.prob)
|
||||
|
||||
w = self._edgeCurl.T * (self._MfMui * bSolution)
|
||||
if S_e is not None:
|
||||
w += -Utils.mkvc(S_e,2)
|
||||
|
||||
if not adjoint:
|
||||
de_dm = self._MeSigmaIDeriv(w) * v
|
||||
elif adjoint:
|
||||
de_dm = self._MeSigmaIDeriv(w).T * v
|
||||
|
||||
_, S_eDeriv = src.evalDeriv(self.prob, adjoint)
|
||||
Se_Deriv = S_eDeriv(v)
|
||||
|
||||
if Se_Deriv is not None:
|
||||
de_dm += -self._MeSigmaI * Se_Deriv
|
||||
|
||||
return de_dm
|
||||
|
||||
def _e(self, bSolution, srcList):
|
||||
return self._ePrimary(bSolution, srcList) + self._eSecondary(bSolution, srcList)
|
||||
|
||||
def _eDeriv_u(self, src, v, adjoint=False):
|
||||
return self._eSecondaryDeriv_u(src, v, adjoint)
|
||||
|
||||
def _eDeriv_m(self, src, v, adjoint=False):
|
||||
# assuming primary doesn't depend on model
|
||||
return self._eSecondaryDeriv_m(src, v, adjoint)
|
||||
|
||||
|
||||
class FieldsFDEM_j(FieldsFDEM):
|
||||
knownFields = {'jSolution':'F'}
|
||||
aliasFields = {
|
||||
'j' : ['jSolution','F','_j'],
|
||||
'jPrimary' : ['jSolution','F','_jPrimary'],
|
||||
'jSecondary' : ['jSolution','F','_jSecondary'],
|
||||
'h' : ['jSolution','E','_h'],
|
||||
'hPrimary' : ['jSolution','E','_hPrimary'],
|
||||
'hSecondary' : ['jSolution','E','_hSecondary'],
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
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
|
||||
self._MfRhoDeriv = self.survey.prob.MfRhoDeriv
|
||||
|
||||
def _jPrimary(self, jSolution, srcList):
|
||||
jPrimary = np.zeros_like(jSolution,dtype = complex)
|
||||
for i, src in enumerate(srcList):
|
||||
jp = src.jPrimary(self.prob)
|
||||
if jp is not None:
|
||||
jPrimary[:,i] += jp
|
||||
return jPrimary
|
||||
|
||||
def _jSecondary(self, jSolution, srcList):
|
||||
return jSolution
|
||||
|
||||
def _j(self, jSolution, srcList):
|
||||
return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList)
|
||||
|
||||
def _jDeriv_u(self, src, v, adjoint=False):
|
||||
return None
|
||||
|
||||
def _jDeriv_m(self, src, v, adjoint=False):
|
||||
# assuming primary does not depend on the model
|
||||
return None
|
||||
|
||||
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.prob)
|
||||
if hp is not None:
|
||||
hPrimary[:,i] = hp
|
||||
return hPrimary
|
||||
|
||||
def _hSecondary(self, jSolution, srcList):
|
||||
MeMuI = self._MeMuI
|
||||
C = self._edgeCurl
|
||||
MfRho = self._MfRho
|
||||
h = MeMuI * (C.T * (MfRho * jSolution) )
|
||||
for i, src in enumerate(srcList):
|
||||
h[:,i] *= -1./(1j*omega(src.freq))
|
||||
S_m,_ = src.eval(self.prob)
|
||||
if S_m is not None:
|
||||
h[:,i] += 1./(1j*omega(src.freq)) * MeMuI * S_m
|
||||
return h
|
||||
|
||||
def _hSecondaryDeriv_u(self, src, v, adjoint=False):
|
||||
MeMuI = self._MeMuI
|
||||
C = self._edgeCurl
|
||||
MfRho = self._MfRho
|
||||
if not adjoint:
|
||||
return -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRho * v) )
|
||||
elif adjoint:
|
||||
return -1./(1j*omega(src.freq)) * MfRho.T * (C * ( MeMuI.T * v))
|
||||
|
||||
def _hSecondaryDeriv_m(self, src, v, adjoint=False):
|
||||
jSolution = self[[src],'jSolution']
|
||||
MeMuI = self._MeMuI
|
||||
C = self._edgeCurl
|
||||
MfRho = self._MfRho
|
||||
MfRhoDeriv = self._MfRhoDeriv
|
||||
|
||||
if not adjoint:
|
||||
hDeriv_m = -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) )
|
||||
elif adjoint:
|
||||
hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) )
|
||||
|
||||
S_mDeriv,_ = src.evalDeriv(self.prob, adjoint)
|
||||
|
||||
if not adjoint:
|
||||
S_mDeriv = S_mDeriv(v)
|
||||
if S_mDeriv is not None:
|
||||
hDeriv_m += 1./(1j*omega(src.freq)) * MeMuI * S_mDeriv
|
||||
elif adjoint:
|
||||
S_mDeriv = S_mDeriv(MeMuI.T * v)
|
||||
if S_mDeriv is not None:
|
||||
hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv
|
||||
return hDeriv_m
|
||||
|
||||
|
||||
def _h(self, jSolution, srcList):
|
||||
return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList)
|
||||
|
||||
def _hDeriv_u(self, src, v, adjoint=False):
|
||||
return self._hSecondaryDeriv_u(src, v, adjoint)
|
||||
|
||||
def _hDeriv_m(self, src, v, adjoint=False):
|
||||
# assuming the primary doesn't depend on the model
|
||||
return self._hSecondaryDeriv_m(src, v, adjoint)
|
||||
|
||||
|
||||
class FieldsFDEM_h(FieldsFDEM):
|
||||
knownFields = {'hSolution':'E'}
|
||||
aliasFields = {
|
||||
'h' : ['hSolution','E','_h'],
|
||||
'hPrimary' : ['hSolution','E','_hPrimary'],
|
||||
'hSecondary' : ['hSolution','E','_hSecondary'],
|
||||
'j' : ['hSolution','F','_j'],
|
||||
'jPrimary' : ['hSolution','F','_jPrimary'],
|
||||
'jSecondary' : ['hSolution','F','_jSecondary']
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
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
|
||||
|
||||
def _hPrimary(self, hSolution, srcList):
|
||||
hPrimary = np.zeros_like(hSolution,dtype = complex)
|
||||
for i, src in enumerate(srcList):
|
||||
hp = src.hPrimary(self.prob)
|
||||
if hp is not None:
|
||||
hPrimary[:,i] += hp
|
||||
return hPrimary
|
||||
|
||||
def _hSecondary(self, hSolution, srcList):
|
||||
return hSolution
|
||||
|
||||
def _h(self, hSolution, srcList):
|
||||
return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList)
|
||||
|
||||
def _hDeriv_u(self, src, v, adjoint=False):
|
||||
return None
|
||||
|
||||
def _hDeriv_m(self, src, v, adjoint=False):
|
||||
# assuming primary does not depend on the model
|
||||
return None
|
||||
|
||||
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.prob)
|
||||
if jp is not None:
|
||||
jPrimary[:,i] = jp
|
||||
return jPrimary
|
||||
|
||||
def _jSecondary(self, hSolution, srcList):
|
||||
j = self._edgeCurl*hSolution
|
||||
for i, src in enumerate(srcList):
|
||||
_,S_e = src.eval(self.prob)
|
||||
if S_e is not None:
|
||||
j[:,i] += -S_e
|
||||
return j
|
||||
|
||||
def _jSecondaryDeriv_u(self, src, v, adjoint=False):
|
||||
if not adjoint:
|
||||
return self._edgeCurl*v
|
||||
elif adjoint:
|
||||
return self._edgeCurl.T*v
|
||||
|
||||
def _jSecondaryDeriv_m(self, src, v, adjoint=False):
|
||||
_,S_eDeriv = src.evalDeriv(self.prob, adjoint)
|
||||
S_eDeriv = S_eDeriv(v)
|
||||
if S_eDeriv is not None:
|
||||
return -S_eDeriv
|
||||
return None
|
||||
|
||||
def _j(self, hSolution, srcList):
|
||||
return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList)
|
||||
|
||||
def _jDeriv_u(self, src, v, adjoint=False):
|
||||
return self._jSecondaryDeriv_u(src,v,adjoint)
|
||||
|
||||
def _jDeriv_m(self, src, v, adjoint=False):
|
||||
# assuming the primary does not depend on the model
|
||||
return self._jSecondaryDeriv_m(src,v,adjoint)
|
||||
+324
-32
@@ -1,4 +1,11 @@
|
||||
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
|
||||
####################################################
|
||||
|
||||
class RxFDEM(Survey.BaseRx):
|
||||
|
||||
@@ -51,15 +58,15 @@ class RxFDEM(Survey.BaseRx):
|
||||
"""Component projection (real/imag)"""
|
||||
return self.knownRxTypes[self.rxType][2]
|
||||
|
||||
def projectFields(self, tx, mesh, u):
|
||||
def projectFields(self, src, mesh, u):
|
||||
P = self.getP(mesh)
|
||||
u_part_complex = u[tx, self.projField]
|
||||
u_part_complex = u[src, self.projField]
|
||||
# get the real or imag component
|
||||
real_or_imag = self.projComp
|
||||
u_part = getattr(u_part_complex, real_or_imag)
|
||||
return P*u_part
|
||||
|
||||
def projectFieldsDeriv(self, tx, mesh, u, v, adjoint=False):
|
||||
def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False):
|
||||
P = self.getP(mesh)
|
||||
|
||||
if not adjoint:
|
||||
@@ -80,45 +87,330 @@ class RxFDEM(Survey.BaseRx):
|
||||
return Pv
|
||||
|
||||
|
||||
class TxFDEM(Survey.BaseTx):
|
||||
|
||||
freq = None #: Frequency (float)
|
||||
####################################################
|
||||
# Sources
|
||||
####################################################
|
||||
|
||||
class SrcFDEM(Survey.BaseSrc):
|
||||
freq = None
|
||||
rxPair = RxFDEM
|
||||
|
||||
knownTxTypes = ['VMD', 'VMD_B', 'CircularLoop']
|
||||
def eval(self, prob):
|
||||
S_m = self.S_m(prob)
|
||||
S_e = self.S_e(prob)
|
||||
return S_m, S_e
|
||||
|
||||
radius = None
|
||||
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 __init__(self, loc, txType, freq, rxList):
|
||||
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):
|
||||
self._S_e = np.array(S_e,dtype=complex)
|
||||
self.freq = float(freq)
|
||||
Survey.BaseTx.__init__(self, loc, txType, rxList)
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def S_e(self, prob):
|
||||
return self._S_e
|
||||
|
||||
|
||||
class SrcFDEM_RawVec_m(SrcFDEM):
|
||||
"""
|
||||
RawVec magnetic source. It is defined by the user provided vector S_m
|
||||
|
||||
class FieldsFDEM(Problem.Fields):
|
||||
"""Fancy Field Storage for a FDEM survey."""
|
||||
knownFields = {'b': 'F', 'e': 'E', 'j': 'F', 'h': 'E'}
|
||||
dtype = complex
|
||||
:param numpy.array S_m: magnetic source term
|
||||
:param float freq: frequency
|
||||
:param rxList: receiver list
|
||||
"""
|
||||
|
||||
def __init__(self, rxList, freq, S_m):
|
||||
self._S_m = np.array(S_m,dtype=complex)
|
||||
self.freq = float(freq)
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def S_m(self, prob):
|
||||
return self._S_m
|
||||
|
||||
|
||||
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):
|
||||
self._S_m = np.array(S_m,dtype=complex)
|
||||
self._S_e = np.array(S_e,dtype=complex)
|
||||
self.freq = float(freq)
|
||||
SrcFDEM.__init__(self, rxList)
|
||||
|
||||
def S_m(self, prob):
|
||||
return self._S_m
|
||||
|
||||
def S_e(self, prob):
|
||||
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
|
||||
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
|
||||
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):
|
||||
"""
|
||||
docstring for SurveyFDEM
|
||||
"""
|
||||
|
||||
txPair = TxFDEM
|
||||
srcPair = SrcFDEM
|
||||
|
||||
def __init__(self, txList, **kwargs):
|
||||
def __init__(self, srcList, **kwargs):
|
||||
# Sort these by frequency
|
||||
self.txList = txList
|
||||
self.srcList = srcList
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
|
||||
_freqDict = {}
|
||||
for tx in txList:
|
||||
if tx.freq not in _freqDict:
|
||||
_freqDict[tx.freq] = []
|
||||
_freqDict[tx.freq] += [tx]
|
||||
for src in srcList:
|
||||
if src.freq not in _freqDict:
|
||||
_freqDict[src.freq] = []
|
||||
_freqDict[src.freq] += [src]
|
||||
|
||||
self._freqDict = _freqDict
|
||||
self._freqs = sorted([f for f in self._freqDict])
|
||||
@@ -134,24 +426,24 @@ class SurveyFDEM(Survey.BaseSurvey):
|
||||
return len(self._freqDict)
|
||||
|
||||
@property
|
||||
def nTxByFreq(self):
|
||||
if getattr(self, '_nTxByFreq', None) is None:
|
||||
self._nTxByFreq = {}
|
||||
def nSrcByFreq(self):
|
||||
if getattr(self, '_nSrcByFreq', None) is None:
|
||||
self._nSrcByFreq = {}
|
||||
for freq in self.freqs:
|
||||
self._nTxByFreq[freq] = len(self.getTransmitters(freq))
|
||||
return self._nTxByFreq
|
||||
self._nSrcByFreq[freq] = len(self.getSrcByFreq(freq))
|
||||
return self._nSrcByFreq
|
||||
|
||||
def getTransmitters(self, freq):
|
||||
"""Returns the transmitters associated with a specific frequency."""
|
||||
def getSrcByFreq(self, freq):
|
||||
"""Returns the sources associated with a specific frequency."""
|
||||
assert freq in self._freqDict, "The requested frequency is not in this survey."
|
||||
return self._freqDict[freq]
|
||||
|
||||
def projectFields(self, u):
|
||||
data = Survey.Data(self)
|
||||
for tx in self.txList:
|
||||
for rx in tx.rxList:
|
||||
data[tx, rx] = rx.projectFields(tx, self.mesh, u)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.projectFields(src, self.mesh, u)
|
||||
return data
|
||||
|
||||
def projectFieldsDeriv(self, u):
|
||||
raise Exception('Use Transmitters to project fields deriv.')
|
||||
raise Exception('Use Sources to project fields deriv.')
|
||||
|
||||
@@ -1,2 +1,3 @@
|
||||
from SurveyFDEM import *
|
||||
from FDEM import ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h
|
||||
from FDEM import BaseFDEMProblem, ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h
|
||||
from FieldsFDEM import *
|
||||
@@ -1,101 +0,0 @@
|
||||
from SimPEG import *
|
||||
from scipy.special import ellipk, ellipe
|
||||
|
||||
def MagneticLoopVectorPotential(txLoc, obsLoc, component, radius):
|
||||
"""
|
||||
Calculate the vector potential of horizontal circular loop
|
||||
at given locations
|
||||
|
||||
:param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z)
|
||||
:param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh
|
||||
:param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list
|
||||
:param numpy.ndarray I: Input current of the loop
|
||||
:param numpy.ndarray radius: radius of the loop
|
||||
:rtype: numpy.ndarray
|
||||
:return: The vector potential each dipole at each observation location
|
||||
"""
|
||||
|
||||
if type(component) in [list, tuple]:
|
||||
out = range(len(component))
|
||||
for i, comp in enumerate(component):
|
||||
out[i] = MagneticLoopVectorPotential(txLoc, obsLoc, comp, radius)
|
||||
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(txLoc, getattr(mesh,'grid'+component), component[1], radius)
|
||||
|
||||
txLoc = np.atleast_2d(txLoc)
|
||||
obsLoc = np.atleast_2d(obsLoc)
|
||||
|
||||
n = obsLoc.shape[0]
|
||||
nTx = txLoc.shape[0]
|
||||
|
||||
if component=='z':
|
||||
A = np.zeros((n, nTx))
|
||||
if nTx ==1:
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
else:
|
||||
|
||||
A = np.zeros((n, nTx))
|
||||
for i in range (nTx):
|
||||
x = obsLoc[:, 0] - txLoc[i, 0]
|
||||
y = obsLoc[:, 1] - txLoc[i, 1]
|
||||
z = obsLoc[:, 2] - txLoc[i, 2]
|
||||
r = np.sqrt(x**2 + y**2)
|
||||
m = (4 * radius * r) / ((radius + r)**2 + z**2)
|
||||
m[m > 1.] = 1.
|
||||
# m might be slightly larger than 1 due to rounding errors
|
||||
# but ellipke requires 0 <= m <= 1
|
||||
K = ellipk(m)
|
||||
E = ellipe(m)
|
||||
ind = (r > 0) & (m < 1)
|
||||
# % 1/r singular at r = 0 and K(m) singular at m = 1
|
||||
Aphi = np.zeros(n)
|
||||
# % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi.
|
||||
Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
|
||||
if component == 'x':
|
||||
A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] )
|
||||
elif component == 'y':
|
||||
A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] )
|
||||
else:
|
||||
raise ValueError('Invalid component')
|
||||
|
||||
if nTx == 1:
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG import Mesh
|
||||
import matplotlib.pyplot as plt
|
||||
cs = 20
|
||||
ncx, ncy, ncz = 41, 41, 40
|
||||
hx = np.ones(ncx)*cs
|
||||
hy = np.ones(ncy)*cs
|
||||
hz = np.ones(ncz)*cs
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz], 'CCC')
|
||||
txLoc = np.r_[0., 0., 0.]
|
||||
Ax = MagneticLoopVectorPotential(txLoc, mesh.gridEx, 'x', 200)
|
||||
Ay = MagneticLoopVectorPotential(txLoc, mesh.gridEy, 'y', 200)
|
||||
Az = MagneticLoopVectorPotential(txLoc, mesh.gridEz, 'z', 200)
|
||||
A = np.r_[Ax, Ay, Az]
|
||||
B0 = mesh.edgeCurl*A
|
||||
J0 = mesh.edgeCurl.T*B0
|
||||
|
||||
# mesh.plotImage(A, vType = 'Ex')
|
||||
# mesh.plotImage(A, vType = 'Ey')
|
||||
|
||||
mesh.plotImage(B0, vType = 'Fx')
|
||||
mesh.plotImage(B0, vType = 'Fy')
|
||||
mesh.plotImage(B0, vType = 'Fz')
|
||||
|
||||
# # mesh.plotImage(J0, vType = 'Ex')
|
||||
# mesh.plotImage(J0, vType = 'Ey')
|
||||
# mesh.plotImage(J0, vType = 'Ez')
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
@@ -1,3 +0,0 @@
|
||||
from magneticDipole import MagneticDipoleVectorPotential
|
||||
from CircularLoop import MagneticLoopVectorPotential
|
||||
from magneticDipole import MagneticDipoleFields
|
||||
@@ -1,100 +0,0 @@
|
||||
import numpy as np
|
||||
from scipy.constants import mu_0, pi
|
||||
from SimPEG import Mesh
|
||||
|
||||
def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)):
|
||||
"""
|
||||
Calculate the vector potential of a set of magnetic dipoles
|
||||
at given locations 'ref. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
|
||||
|
||||
:param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z)
|
||||
:param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh
|
||||
:param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list
|
||||
:param numpy.ndarray dipoleMoment: The vector dipole moment
|
||||
:rtype: numpy.ndarray
|
||||
:return: The vector potential each dipole at each observation location
|
||||
"""
|
||||
|
||||
if type(component) in [list, tuple]:
|
||||
out = range(len(component))
|
||||
for i, comp in enumerate(component):
|
||||
out[i] = MagneticDipoleVectorPotential(txLoc, obsLoc, comp, dipoleMoment=dipoleMoment)
|
||||
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 MagneticDipoleVectorPotential(txLoc, getattr(mesh,'grid'+component), component[1], dipoleMoment=dipoleMoment)
|
||||
|
||||
if component == 'x':
|
||||
dimInd = 0
|
||||
elif component == 'y':
|
||||
dimInd = 1
|
||||
elif component == 'z':
|
||||
dimInd = 2
|
||||
else:
|
||||
raise ValueError('Invalid component')
|
||||
|
||||
txLoc = np.atleast_2d(txLoc)
|
||||
obsLoc = np.atleast_2d(obsLoc)
|
||||
dipoleMoment = np.atleast_2d(dipoleMoment)
|
||||
|
||||
nEdges = obsLoc.shape[0]
|
||||
nTx = txLoc.shape[0]
|
||||
|
||||
m = np.array(dipoleMoment).repeat(nEdges, axis=0)
|
||||
A = np.empty((nEdges, nTx))
|
||||
for i in range(nTx):
|
||||
dR = obsLoc - txLoc[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)
|
||||
if nTx == 1:
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
def MagneticDipoleFields(txLoc, obsLoc, component, dipoleMoment=1.):
|
||||
"""
|
||||
Calculate the vector potential of a set of magnetic dipoles
|
||||
at given locations 'ref. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
|
||||
|
||||
:param numpy.ndarray txLoc: Location of the transmitter(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)
|
||||
:rtype: numpy.ndarray
|
||||
:return: The vector potential each dipole at each observation location
|
||||
"""
|
||||
|
||||
if component=='x':
|
||||
dimInd = 0
|
||||
elif component=='y':
|
||||
dimInd = 1
|
||||
elif component=='z':
|
||||
dimInd = 2
|
||||
else:
|
||||
raise ValueError('Invalid component')
|
||||
|
||||
txLoc = np.atleast_2d(txLoc)
|
||||
obsLoc = np.atleast_2d(obsLoc)
|
||||
dipoleMoment = np.atleast_2d(dipoleMoment)
|
||||
|
||||
nFaces = obsLoc.shape[0]
|
||||
nTx = txLoc.shape[0]
|
||||
|
||||
m = np.array(dipoleMoment).repeat(nFaces, axis=0)
|
||||
B = np.empty((nFaces, nTx))
|
||||
for i in range(nTx):
|
||||
dR = obsLoc - txLoc[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)
|
||||
elif dimInd == 1:
|
||||
B[:, i] = +(mu_0/(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)
|
||||
else:
|
||||
raise Exception("Not Implemented")
|
||||
if nTx == 1:
|
||||
return B.flatten()
|
||||
return B
|
||||
@@ -1,6 +1,6 @@
|
||||
from SimPEG import Solver, Problem
|
||||
from SimPEG.Problem import BaseTimeProblem
|
||||
from simpegEM import Sources
|
||||
from simpegEM.Utils import SrcUtils
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.Utils import sdiag, mkvc
|
||||
from SimPEG import Utils, Mesh
|
||||
@@ -13,21 +13,21 @@ class FieldsTDEM(Problem.TimeFields):
|
||||
knownFields = {'b': 'F', 'e': 'E'}
|
||||
|
||||
def tovec(self):
|
||||
nTx, nF, nE = self.survey.nTx, self.mesh.nF, self.mesh.nE
|
||||
u = np.empty(0 if nTx == 1 else (0, nTx))
|
||||
nSrc, nF, nE = self.survey.nSrc, self.mesh.nF, self.mesh.nE
|
||||
u = np.empty((0,nSrc)) #((0,1) if nSrc == 1 else (0, nSrc))
|
||||
|
||||
for i in range(self.survey.prob.nT):
|
||||
if 'b' in self:
|
||||
b = self[:,'b',i+1]
|
||||
else:
|
||||
b = np.zeros(nF if nTx == 1 else (nF, nTx))
|
||||
b = np.zeros((nF,nSrc)) # if nSrc == 1 else (nF, nSrc))
|
||||
|
||||
if 'e' in self:
|
||||
e = self[:,'e',i+1]
|
||||
else:
|
||||
e = np.zeros(nE if nTx == 1 else (nE, nTx))
|
||||
e = np.zeros((nE,nSrc)) # if nSrc == 1 else (nE, nSrc))
|
||||
u = np.concatenate((u, b, e))
|
||||
return Utils.mkvc(u)
|
||||
return Utils.mkvc(u,nSrc)
|
||||
|
||||
|
||||
class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
@@ -42,9 +42,9 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
self.curModel = m
|
||||
# Create a fields storage object
|
||||
F = self._FieldsForward_pair(self.mesh, self.survey)
|
||||
for tx in self.survey.txList:
|
||||
for src in self.survey.srcList:
|
||||
# Set the initial conditions
|
||||
F[tx,:,0] = tx.getInitialFields(self.mesh)
|
||||
F[src,:,0] = src.getInitialFields(self.mesh)
|
||||
F = self.forward(m, self.getRHS, F=F)
|
||||
if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
|
||||
return F
|
||||
|
||||
+45
-31
@@ -1,6 +1,6 @@
|
||||
from SimPEG import Utils, Survey, np
|
||||
from SimPEG.Survey import BaseSurvey
|
||||
from simpegEM import Sources
|
||||
from simpegEM.Utils import SrcUtils
|
||||
from BaseTDEM import FieldsTDEM
|
||||
|
||||
|
||||
@@ -51,76 +51,90 @@ class RxTDEM(Survey.BaseTimeRx):
|
||||
else:
|
||||
return timeMesh.getInterpolationMat(self.times, self.projTLoc)
|
||||
|
||||
def projectFields(self, tx, mesh, timeMesh, u):
|
||||
def projectFields(self, src, mesh, timeMesh, u):
|
||||
P = self.getP(mesh, timeMesh)
|
||||
u_part = Utils.mkvc(u[tx, self.projField, :])
|
||||
u_part = Utils.mkvc(u[src, self.projField, :])
|
||||
return P*u_part
|
||||
|
||||
def projectFieldsDeriv(self, tx, mesh, timeMesh, u, v, adjoint=False):
|
||||
def projectFieldsDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
|
||||
P = self.getP(mesh, timeMesh)
|
||||
|
||||
if not adjoint:
|
||||
return P * Utils.mkvc(v[tx, self.projField, :])
|
||||
return P * Utils.mkvc(v[src, self.projField, :])
|
||||
elif adjoint:
|
||||
return P.T * v[tx, self]
|
||||
return P.T * v[src, self]
|
||||
|
||||
|
||||
class TxTDEM(Survey.BaseTx):
|
||||
class SrcTDEM(Survey.BaseSrc):
|
||||
rxPair = RxTDEM
|
||||
radius = None
|
||||
knownTxTypes = ['VMD_MVP', 'CircularLoop_MVP']
|
||||
|
||||
def getInitialFields(self, mesh):
|
||||
F0 = getattr(self, '_getInitialFields_' + self.txType)(mesh)
|
||||
F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
|
||||
return F0
|
||||
|
||||
def _getInitialFields_VMD_MVP(self, mesh):
|
||||
def getJs(self, mesh, time):
|
||||
return None
|
||||
|
||||
|
||||
class SrcTDEM_VMD_MVP(SrcTDEM):
|
||||
|
||||
def __init__(self,rxList,loc):
|
||||
self.loc = loc
|
||||
SrcTDEM.__init__(self,rxList)
|
||||
|
||||
def getInitialFields(self, mesh):
|
||||
"""Vertical magnetic dipole, magnetic vector potential"""
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = Sources.MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
|
||||
MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = Sources.MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
|
||||
MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
|
||||
else:
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
|
||||
return {"b": mesh.edgeCurl*MVP}
|
||||
|
||||
def _getInitialFields_CircularLoop_MVP(self, mesh):
|
||||
|
||||
class SrcTDEM_CircularLoop_MVP(SrcTDEM):
|
||||
|
||||
def __init__(self,rxList,loc):
|
||||
self.loc = loc
|
||||
SrcTDEM.__init__(self,rxList)
|
||||
|
||||
def getInitialFields_(self, mesh):
|
||||
"""Circular Loop, magnetic vector potential"""
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = Sources.MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
|
||||
MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
|
||||
else:
|
||||
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
|
||||
elif mesh._meshType is 'TENSOR':
|
||||
MVP = Sources.MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
|
||||
MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
|
||||
else:
|
||||
raise Exception('Unknown mesh for CircularLoop')
|
||||
|
||||
return {"b": mesh.edgeCurl*MVP}
|
||||
|
||||
def getJs(self, mesh, time):
|
||||
return None
|
||||
|
||||
class SurveyTDEM(Survey.BaseSurvey):
|
||||
"""
|
||||
docstring for SurveyTDEM
|
||||
"""
|
||||
txPair = TxTDEM
|
||||
srcPair = SrcTDEM
|
||||
|
||||
def __init__(self, txList, **kwargs):
|
||||
def __init__(self, srcList, **kwargs):
|
||||
# Sort these by frequency
|
||||
self.txList = txList
|
||||
self.srcList = srcList
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
|
||||
def projectFields(self, u):
|
||||
data = Survey.Data(self)
|
||||
for tx in self.txList:
|
||||
for rx in tx.rxList:
|
||||
data[tx, rx] = rx.projectFields(tx, self.mesh, self.prob.timeMesh, u)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.projectFields(src, self.mesh, self.prob.timeMesh, u)
|
||||
return data
|
||||
|
||||
def projectFieldsDeriv(self, u, v=None, adjoint=False):
|
||||
@@ -128,20 +142,20 @@ class SurveyTDEM(Survey.BaseSurvey):
|
||||
|
||||
if not adjoint:
|
||||
data = Survey.Data(self)
|
||||
for tx in self.txList:
|
||||
for rx in tx.rxList:
|
||||
data[tx, rx] = rx.projectFieldsDeriv(tx, self.mesh, self.prob.timeMesh, u, v)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v)
|
||||
return data
|
||||
else:
|
||||
f = FieldsTDEM(self.mesh, self)
|
||||
for tx in self.txList:
|
||||
for rx in tx.rxList:
|
||||
Ptv = rx.projectFieldsDeriv(tx, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
Ptv = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
|
||||
Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
|
||||
if rx.projField not in f: # first time we are projecting
|
||||
f[tx, rx.projField, :] = Ptv
|
||||
f[src, rx.projField, :] = Ptv
|
||||
else: # there are already fields, so let's add to them!
|
||||
f[tx, rx.projField, :] += Ptv
|
||||
f[src, rx.projField, :] += Ptv
|
||||
return f
|
||||
|
||||
|
||||
|
||||
+26
-21
@@ -14,7 +14,7 @@ class FieldsTDEM_e_from_b(FieldsTDEM):
|
||||
self.edgeCurlT = self.survey.prob.mesh.edgeCurl.T
|
||||
self.MfMui = self.survey.prob.MfMui
|
||||
|
||||
def e_from_b(self, b, txInd, timeInd):
|
||||
def e_from_b(self, b, srcInd, timeInd):
|
||||
# TODO: implement non-zero js
|
||||
return self.MeSigmaI*(self.edgeCurlT*(self.MfMui*b))
|
||||
|
||||
@@ -32,10 +32,10 @@ class FieldsTDEM_e_from_b_Ah(FieldsTDEM):
|
||||
self.edgeCurlT = self.survey.prob.mesh.edgeCurl.T
|
||||
self.MfMui = self.survey.prob.MfMui
|
||||
|
||||
def e_from_b(self, y_b, txInd, tInd):
|
||||
def e_from_b(self, y_b, srcInd, tInd):
|
||||
y_e = self.MeSigmaI*(self.edgeCurlT*(self.MfMui*y_b))
|
||||
if 'e' in self.p:
|
||||
y_e = y_e - self.MeSigmaI*self.p[txInd,'e',tInd]
|
||||
y_e = y_e - self.MeSigmaI*self.p[srcInd,'e',tInd]
|
||||
return y_e
|
||||
|
||||
class ProblemTDEM_b(BaseTDEMProblem):
|
||||
@@ -73,8 +73,10 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
|
||||
def getRHS(self, tInd, F):
|
||||
dt = self.timeSteps[tInd]
|
||||
B_n = np.c_[[F[tx,'b',tInd] for tx in self.survey.txList]].T
|
||||
RHS = (1.0/dt)*self.MfMui*B_n
|
||||
B_n = np.c_[[F[src,'b',tInd] for src in self.survey.srcList]].T
|
||||
if B_n.shape[0] is not 1:
|
||||
raise NotImplementedError('getRHS not implemented for this shape of B_n')
|
||||
RHS = (1.0/dt)*self.MfMui*B_n[0,:,:] #TODO: This is a hack
|
||||
return RHS
|
||||
|
||||
####################################################
|
||||
@@ -95,7 +97,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
u = self.fields(m)
|
||||
self.curModel = m
|
||||
|
||||
# Note: Fields has shape (nF/E, nTx, nT+1)
|
||||
# Note: Fields has shape (nF/E, nSrc, nT+1)
|
||||
# However, p will only really fill (:,:,1:nT+1)
|
||||
# meaning the 'initial fields' are zero (:,:,0)
|
||||
p = FieldsTDEM(self.mesh, self.survey)
|
||||
@@ -106,15 +108,17 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
|
||||
# fake initial 'e' fields
|
||||
p[:, 'e', 0] = 0.0
|
||||
dMdsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform)
|
||||
dsigdm_x_v = self.curModel.transformDeriv*vec
|
||||
dMdsig = self.MeSigmaDeriv
|
||||
# self.mesh.getEdgeInnerProductDeriv(self.curModel.transform)
|
||||
# dsigdm_x_v = self.curModel.sigmaDeriv*vec
|
||||
# dsigdm_x_v = self.curModel.transformDeriv*vec
|
||||
for i in range(1,self.nT+1):
|
||||
# TODO: G[1] may be dependent on the model
|
||||
# for a galvanic source (deriv of the dc problem)
|
||||
#
|
||||
# Do multiplication for all tx in self.survey.txList
|
||||
for tx in self.survey.txList:
|
||||
p[tx, 'e', i] = - dMdsig(u[tx,'e',i]) * dsigdm_x_v
|
||||
# Do multiplication for all src in self.survey.srcList
|
||||
for src in self.survey.srcList:
|
||||
p[src, 'e', i] = - dMdsig(u[src,'e',i]) * vec
|
||||
return p
|
||||
|
||||
def Gtvec(self, m, vec, u=None):
|
||||
@@ -130,20 +134,21 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
self.curModel = m
|
||||
dMdsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform)
|
||||
dsigdm = self.curModel.transformDeriv
|
||||
# dMdsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform)
|
||||
# dsigdm = self.curModel.transformDeriv
|
||||
MeSigmaDeriv = self.MeSigmaDeriv
|
||||
|
||||
nTx = self.survey.nTx
|
||||
nSrc = self.survey.nSrc
|
||||
VUs = None
|
||||
# Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go.
|
||||
for i in range(1,self.nT+1):
|
||||
vu = None
|
||||
for tx in self.survey.txList:
|
||||
vutx = dMdsig(u[tx,'e',i]).T * vec[tx,'e',i]
|
||||
vu = vutx if vu is None else vu + vutx
|
||||
for src in self.survey.srcList:
|
||||
vusrc = MeSigmaDeriv(u[src,'e',i]).T * vec[src,'e',i]
|
||||
vu = vusrc if vu is None else vu + vusrc
|
||||
VUs = vu if VUs is None else VUs + vu
|
||||
p = -dsigdm.T*VUs
|
||||
return p
|
||||
# p = -dsigdm.T*VUs
|
||||
return -VUs
|
||||
|
||||
def solveAh(self, m, p):
|
||||
"""
|
||||
@@ -241,8 +246,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
|
||||
# 1 (tInd=1 uses fields 2 and 3)
|
||||
|
||||
def AhtRHS(tInd, y):
|
||||
nTx, nF = self.survey.nTx, self.mesh.nF
|
||||
rhs = np.zeros(nF if nTx == 1 else (nF, nTx))
|
||||
nSrc, nF = self.survey.nSrc, self.mesh.nF
|
||||
rhs = np.zeros((nF,1) if nSrc == 1 else (nF, nSrc))
|
||||
|
||||
if 'e' in p:
|
||||
rhs += self.MfMui*(self.mesh.edgeCurl*(self.MeSigmaI*p[:,'e',tInd+1]))
|
||||
|
||||
@@ -1,3 +1,3 @@
|
||||
from SurveyTDEM import SurveyTDEM, RxTDEM, TxTDEM
|
||||
from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
|
||||
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
|
||||
from TDEM_b import ProblemTDEM_b
|
||||
|
||||
+339
-269
@@ -1,9 +1,17 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import simpegEM as EM
|
||||
import sys
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
|
||||
testDerivs = True
|
||||
testCrossCheck = True
|
||||
testAdjoint = True
|
||||
testEB = True
|
||||
testHJ = True
|
||||
|
||||
verbose = False
|
||||
|
||||
TOL = 1e-4
|
||||
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
|
||||
CONDUCTIVITY = 1e1
|
||||
@@ -11,6 +19,9 @@ MU = mu_0
|
||||
freq = 1e-1
|
||||
addrandoms = True
|
||||
|
||||
SrcType = 'MagDipole' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
|
||||
|
||||
|
||||
def getProblem(fdemType, comp):
|
||||
cs = 5.
|
||||
ncx, ncy, ncz = 6, 6, 6
|
||||
@@ -22,22 +33,64 @@ def getProblem(fdemType, comp):
|
||||
|
||||
mapping = Maps.ExpMap(mesh)
|
||||
|
||||
x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the transmitter
|
||||
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)
|
||||
Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', freq, [Rx0])
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Tx0])
|
||||
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)
|
||||
@@ -54,38 +107,42 @@ def adjointTest(fdemType, comp):
|
||||
prb = getProblem(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)
|
||||
m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
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.mapping.nP)*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)
|
||||
w = np.random.rand(prb.mapping.nP)
|
||||
w = np.random.rand(prb.mesh.nC)
|
||||
|
||||
vJw = v.dot(prb.Jvec(m, w, u=u))
|
||||
wJtv = w.dot(prb.Jtvec(m, v, u=u))
|
||||
vJw = v.dot(prb.Jvec(m, w, u))
|
||||
wJtv = w.dot(prb.Jtvec(m, v, u))
|
||||
tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
|
||||
print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
|
||||
return np.abs(vJw - wJtv) < tol
|
||||
|
||||
|
||||
def derivTest(fdemType, comp):
|
||||
|
||||
prb = getProblem(fdemType, comp)
|
||||
print '%s formulation - %s' % (fdemType, comp)
|
||||
x0 = np.log(np.ones(prb.mesh.nC)*CONDUCTIVITY)
|
||||
x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
mu = np.log(np.ones(prb.mesh.nC)*MU)
|
||||
|
||||
if addrandoms is True:
|
||||
x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1
|
||||
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
|
||||
x0 = x0 + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
|
||||
mu = mu + np.random.randn(prb.mapping.nP)*MU*1e-1
|
||||
|
||||
# prb.PropMap.PropModel.mu = mu
|
||||
# prb.PropMap.PropModel.mui = 1./mu
|
||||
|
||||
prb.mu = mu
|
||||
survey = prb.survey
|
||||
def fun(x):
|
||||
return survey.dpred(x), lambda x: prb.Jvec(x0, x)
|
||||
@@ -103,16 +160,16 @@ 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)
|
||||
|
||||
u1 = prb1.fields(m)
|
||||
d1 = Utils.mkvc(survey1.projectFields(u1))
|
||||
|
||||
prb1.unpair
|
||||
if verbose:
|
||||
print ' Problem 1 solved'
|
||||
|
||||
if fdemType == 'e':
|
||||
prb2 = getProblem('b', comp)
|
||||
@@ -125,11 +182,12 @@ def crossCheckTest(fdemType, comp):
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
prb2.mu = mu
|
||||
# prb2.mu = mu
|
||||
survey2 = prb2.survey
|
||||
d2 = survey2.dpred(m)
|
||||
|
||||
u2 = prb2.fields(m)
|
||||
d2 = Utils.mkvc(survey2.projectFields(u2))
|
||||
if verbose:
|
||||
print ' Problem 2 solved'
|
||||
|
||||
r = d2-d1
|
||||
l2r = l2norm(r)
|
||||
@@ -144,265 +202,277 @@ def crossCheckTest(fdemType, comp):
|
||||
|
||||
class FDEM_DerivTests(unittest.TestCase):
|
||||
|
||||
def test_Jvec_exr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'exr'))
|
||||
def test_Jvec_exr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'exr'))
|
||||
def test_Jvec_eyr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'eyr'))
|
||||
def test_Jvec_eyr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'eyr'))
|
||||
def test_Jvec_ezr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'ezr'))
|
||||
def test_Jvec_ezr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'ezr'))
|
||||
def test_Jvec_exi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'exi'))
|
||||
def test_Jvec_exi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'exi'))
|
||||
def test_Jvec_eyi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'eyi'))
|
||||
def test_Jvec_eyi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'eyi'))
|
||||
def test_Jvec_ezi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'ezi'))
|
||||
def test_Jvec_ezi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'ezi'))
|
||||
if testDerivs:
|
||||
if testEB:
|
||||
def test_Jvec_exr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'exr'))
|
||||
def test_Jvec_eyr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'eyr'))
|
||||
def test_Jvec_ezr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'ezr'))
|
||||
def test_Jvec_exi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'exi'))
|
||||
def test_Jvec_eyi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'eyi'))
|
||||
def test_Jvec_ezi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'ezi'))
|
||||
|
||||
def test_Jvec_bxr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bxr'))
|
||||
def test_Jvec_bxr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bxr'))
|
||||
def test_Jvec_byr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'byr'))
|
||||
def test_Jvec_byr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'byr'))
|
||||
def test_Jvec_bzr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bzr'))
|
||||
def test_Jvec_bzr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bzr'))
|
||||
def test_Jvec_bxi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bxi'))
|
||||
def test_Jvec_bxi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bxi'))
|
||||
def test_Jvec_byi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'byi'))
|
||||
def test_Jvec_byi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'byi'))
|
||||
def test_Jvec_bzi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bzi'))
|
||||
def test_Jvec_bzi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bzi'))
|
||||
def test_Jvec_bxr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bxr'))
|
||||
def test_Jvec_byr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'byr'))
|
||||
def test_Jvec_bzr_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bzr'))
|
||||
def test_Jvec_bxi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bxi'))
|
||||
def test_Jvec_byi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'byi'))
|
||||
def test_Jvec_bzi_Eform(self):
|
||||
self.assertTrue(derivTest('e', 'bzi'))
|
||||
|
||||
def test_Jvec_exr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'exr'))
|
||||
def test_Jvec_eyr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'eyr'))
|
||||
def test_Jvec_ezr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'ezr'))
|
||||
def test_Jvec_exi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'exi'))
|
||||
def test_Jvec_eyi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'eyi'))
|
||||
def test_Jvec_ezi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'ezi'))
|
||||
|
||||
def test_Jvec_bxr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bxr'))
|
||||
def test_Jvec_byr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'byr'))
|
||||
def test_Jvec_bzr_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bzr'))
|
||||
def test_Jvec_bxi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bxi'))
|
||||
def test_Jvec_byi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'byi'))
|
||||
def test_Jvec_bzi_Bform(self):
|
||||
self.assertTrue(derivTest('b', 'bzi'))
|
||||
|
||||
if testHJ:
|
||||
def test_Jvec_jxr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jxr'))
|
||||
def test_Jvec_jyr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jyr'))
|
||||
def test_Jvec_jzr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jzr'))
|
||||
def test_Jvec_jxi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jxi'))
|
||||
def test_Jvec_jyi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jyi'))
|
||||
def test_Jvec_jzi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jzi'))
|
||||
|
||||
def test_Jvec_hxr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hxr'))
|
||||
def test_Jvec_hyr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hyr'))
|
||||
def test_Jvec_hzr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hzr'))
|
||||
def test_Jvec_hxi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hxi'))
|
||||
def test_Jvec_hyi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hyi'))
|
||||
def test_Jvec_hzi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hzi'))
|
||||
|
||||
def test_Jvec_hxr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hxr'))
|
||||
def test_Jvec_hyr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hyr'))
|
||||
def test_Jvec_hzr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hzr'))
|
||||
def test_Jvec_hxi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hxi'))
|
||||
def test_Jvec_hyi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hyi'))
|
||||
def test_Jvec_hzi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hzi'))
|
||||
|
||||
def test_Jvec_hxr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jxr'))
|
||||
def test_Jvec_hyr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jyr'))
|
||||
def test_Jvec_hzr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jzr'))
|
||||
def test_Jvec_hxi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jxi'))
|
||||
def test_Jvec_hyi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jyi'))
|
||||
def test_Jvec_hzi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jzi'))
|
||||
|
||||
|
||||
if testAdjoint:
|
||||
if testEB:
|
||||
def test_Jtvec_adjointTest_exr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'exr'))
|
||||
def test_Jtvec_adjointTest_eyr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'eyr'))
|
||||
def test_Jtvec_adjointTest_ezr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'ezr'))
|
||||
def test_Jtvec_adjointTest_exi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'exi'))
|
||||
def test_Jtvec_adjointTest_eyi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'eyi'))
|
||||
def test_Jtvec_adjointTest_ezi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'ezi'))
|
||||
|
||||
def test_Jtvec_adjointTest_exr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'exr'))
|
||||
def test_Jtvec_adjointTest_exr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'exr'))
|
||||
def test_Jtvec_adjointTest_eyr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'eyr'))
|
||||
def test_Jtvec_adjointTest_eyr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'eyr'))
|
||||
def test_Jtvec_adjointTest_ezr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'ezr'))
|
||||
def test_Jtvec_adjointTest_ezr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'ezr'))
|
||||
def test_Jtvec_adjointTest_exi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'exi'))
|
||||
def test_Jtvec_adjointTest_exi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'exi'))
|
||||
def test_Jtvec_adjointTest_eyi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'eyi'))
|
||||
def test_Jtvec_adjointTest_eyi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'eyi'))
|
||||
def test_Jtvec_adjointTest_ezi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'ezi'))
|
||||
def test_Jtvec_adjointTest_ezi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'ezi'))
|
||||
def test_Jtvec_adjointTest_bxr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bxr'))
|
||||
def test_Jtvec_adjointTest_byr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'byr'))
|
||||
def test_Jtvec_adjointTest_bzr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bzr'))
|
||||
def test_Jtvec_adjointTest_bxi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bxi'))
|
||||
def test_Jtvec_adjointTest_byi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'byi'))
|
||||
def test_Jtvec_adjointTest_bzi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_bxr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bxr'))
|
||||
def test_Jtvec_adjointTest_bxr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bxr'))
|
||||
def test_Jtvec_adjointTest_byr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'byr'))
|
||||
def test_Jtvec_adjointTest_byr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'byr'))
|
||||
def test_Jtvec_adjointTest_bzr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bzr'))
|
||||
def test_Jtvec_adjointTest_bzr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bzr'))
|
||||
def test_Jtvec_adjointTest_bxi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bxi'))
|
||||
def test_Jtvec_adjointTest_bxi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bxi'))
|
||||
def test_Jtvec_adjointTest_byi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'byi'))
|
||||
def test_Jtvec_adjointTest_byi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'byi'))
|
||||
def test_Jtvec_adjointTest_bzi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bzi'))
|
||||
def test_Jtvec_adjointTest_bzi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bzi'))
|
||||
def test_Jtvec_adjointTest_exr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'exr'))
|
||||
def test_Jtvec_adjointTest_eyr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'eyr'))
|
||||
def test_Jtvec_adjointTest_ezr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'ezr'))
|
||||
def test_Jtvec_adjointTest_exi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'exi'))
|
||||
def test_Jtvec_adjointTest_eyi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'eyi'))
|
||||
def test_Jtvec_adjointTest_ezi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'ezi'))
|
||||
def test_Jtvec_adjointTest_bxr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bxr'))
|
||||
def test_Jtvec_adjointTest_byr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'byr'))
|
||||
def test_Jtvec_adjointTest_bzr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bzr'))
|
||||
def test_Jtvec_adjointTest_bxi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bxi'))
|
||||
def test_Jtvec_adjointTest_byi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'byi'))
|
||||
def test_Jtvec_adjointTest_bzi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bzi'))
|
||||
|
||||
|
||||
def test_Jvec_jxr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jxr'))
|
||||
def test_Jvec_jyr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jyr'))
|
||||
def test_Jvec_jzr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jzr'))
|
||||
def test_Jvec_jxi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jxi'))
|
||||
def test_Jvec_jyi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jyi'))
|
||||
def test_Jvec_jzi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'jzi'))
|
||||
if testHJ:
|
||||
def test_Jtvec_adjointTest_jxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxr'))
|
||||
def test_Jtvec_adjointTest_jyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyr'))
|
||||
def test_Jtvec_adjointTest_jzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzr'))
|
||||
def test_Jtvec_adjointTest_jxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxi'))
|
||||
def test_Jtvec_adjointTest_jyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyi'))
|
||||
def test_Jtvec_adjointTest_jzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzi'))
|
||||
|
||||
def test_Jvec_hxr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hxr'))
|
||||
def test_Jvec_hyr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hyr'))
|
||||
def test_Jvec_hzr_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hzr'))
|
||||
def test_Jvec_hxi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hxi'))
|
||||
def test_Jvec_hyi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hyi'))
|
||||
def test_Jvec_hzi_Jform(self):
|
||||
self.assertTrue(derivTest('j', 'hzi'))
|
||||
def test_Jtvec_adjointTest_hxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzi'))
|
||||
|
||||
def test_Jvec_hxr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hxr'))
|
||||
def test_Jvec_hyr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hyr'))
|
||||
def test_Jvec_hzr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hzr'))
|
||||
def test_Jvec_hxi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hxi'))
|
||||
def test_Jvec_hyi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hyi'))
|
||||
def test_Jvec_hzi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'hzi'))
|
||||
def test_Jtvec_adjointTest_hxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzi'))
|
||||
|
||||
def test_Jvec_hxr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jxr'))
|
||||
def test_Jvec_hyr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jyr'))
|
||||
def test_Jvec_hzr_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jzr'))
|
||||
def test_Jvec_hxi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jxi'))
|
||||
def test_Jvec_hyi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jyi'))
|
||||
def test_Jvec_hzi_Hform(self):
|
||||
self.assertTrue(derivTest('h', 'jzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_jxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxr'))
|
||||
def test_Jtvec_adjointTest_jyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyr'))
|
||||
def test_Jtvec_adjointTest_jzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzr'))
|
||||
def test_Jtvec_adjointTest_jxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxi'))
|
||||
def test_Jtvec_adjointTest_jyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyi'))
|
||||
def test_Jtvec_adjointTest_jzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzi'))
|
||||
def test_Jtvec_adjointTest_hxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzi'))
|
||||
|
||||
|
||||
def test_EB_CrossCheck_exr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exr'))
|
||||
def test_EB_CrossCheck_eyr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyr'))
|
||||
def test_EB_CrossCheck_ezr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezr'))
|
||||
def test_EB_CrossCheck_exi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exi'))
|
||||
def test_EB_CrossCheck_eyi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyi'))
|
||||
def test_EB_CrossCheck_ezi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezi'))
|
||||
if testCrossCheck:
|
||||
if testEB:
|
||||
def test_EB_CrossCheck_exr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exr'))
|
||||
def test_EB_CrossCheck_eyr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyr'))
|
||||
def test_EB_CrossCheck_ezr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezr'))
|
||||
def test_EB_CrossCheck_exi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exi'))
|
||||
def test_EB_CrossCheck_eyi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyi'))
|
||||
def test_EB_CrossCheck_ezi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezi'))
|
||||
|
||||
def test_EB_CrossCheck_bxr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxr'))
|
||||
def test_EB_CrossCheck_byr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byr'))
|
||||
# def test_EB_CrossCheck_bzr_Eform(self):
|
||||
# self.assertTrue(crossCheckTest('e', 'bzr')) # Doesn't make sense to test this for p-s approach
|
||||
def test_EB_CrossCheck_bxi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxi'))
|
||||
def test_EB_CrossCheck_byi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byi'))
|
||||
def test_EB_CrossCheck_bzi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzi'))
|
||||
def test_EB_CrossCheck_bxr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxr'))
|
||||
def test_EB_CrossCheck_byr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byr'))
|
||||
def test_EB_CrossCheck_bzr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzr'))
|
||||
def test_EB_CrossCheck_bxi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxi'))
|
||||
def test_EB_CrossCheck_byi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byi'))
|
||||
def test_EB_CrossCheck_bzi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzi'))
|
||||
|
||||
def test_HJ_CrossCheck_jxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxr'))
|
||||
def test_HJ_CrossCheck_jyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyr'))
|
||||
def test_HJ_CrossCheck_jzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzr'))
|
||||
def test_HJ_CrossCheck_jxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxi'))
|
||||
def test_HJ_CrossCheck_jyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyi'))
|
||||
def test_HJ_CrossCheck_jzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzi'))
|
||||
if testHJ:
|
||||
def test_HJ_CrossCheck_jxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxr'))
|
||||
def test_HJ_CrossCheck_jyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyr'))
|
||||
def test_HJ_CrossCheck_jzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzr'))
|
||||
def test_HJ_CrossCheck_jxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxi'))
|
||||
def test_HJ_CrossCheck_jyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyi'))
|
||||
def test_HJ_CrossCheck_jzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzi'))
|
||||
|
||||
def test_HJ_CrossCheck_hxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxr'))
|
||||
def test_HJ_CrossCheck_hyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyr'))
|
||||
def test_HJ_CrossCheck_hzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzr'))
|
||||
def test_HJ_CrossCheck_hxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxi'))
|
||||
def test_HJ_CrossCheck_hyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyi'))
|
||||
def test_HJ_CrossCheck_hzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzi'))
|
||||
|
||||
def test_HJ_CrossCheck_hxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxr'))
|
||||
def test_HJ_CrossCheck_hyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyr'))
|
||||
# def test_HJ_CrossCheck_hzr_Jform(self):
|
||||
# self.assertTrue(crossCheckTest('j', 'hzr')) # Doesn't make sense to test this for p-s approach
|
||||
def test_HJ_CrossCheck_hxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxi'))
|
||||
def test_HJ_CrossCheck_hyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyi'))
|
||||
def test_HJ_CrossCheck_hzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzi'))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -0,0 +1,61 @@
|
||||
from SimPEG import Tests, Utils, np
|
||||
import simpegEM.Analytics.FDEMcasing as Casing
|
||||
import unittest
|
||||
from scipy.constants import mu_0
|
||||
|
||||
|
||||
n = 50.
|
||||
freq = 1.
|
||||
a = 5e-2
|
||||
b = a + 1e-2
|
||||
sigma = np.r_[10., 5.5e6, 1e-1]
|
||||
mu = mu_0*np.r_[1.,100.,1.]
|
||||
srcloc = np.r_[0., 0., 0.]
|
||||
xobs = np.random.rand(n)+10.
|
||||
yobs = np.zeros(n)
|
||||
zobs = np.random.randn(n)
|
||||
plotit = False
|
||||
|
||||
def CasingMagDipoleDeriv_r(x):
|
||||
obsloc = np.vstack([x, yobs, zobs]).T
|
||||
|
||||
f = Casing._getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu)
|
||||
g = Utils.sdiag(Casing._getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu))
|
||||
|
||||
return f,g
|
||||
|
||||
def CasingMagDipoleDeriv_z(z):
|
||||
obsloc = np.vstack([xobs, yobs, z]).T
|
||||
|
||||
f = Casing._getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu)
|
||||
g = Utils.sdiag(Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu))
|
||||
|
||||
return f,g
|
||||
|
||||
def CasingMagDipole2Deriv_z_r(x):
|
||||
obsloc = np.vstack([x, yobs, zobs]).T
|
||||
|
||||
f = Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu)
|
||||
g = Utils.sdiag(Casing._getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu))
|
||||
|
||||
return f,g
|
||||
|
||||
def CasingMagDipole2Deriv_z_z(z):
|
||||
obsloc = np.vstack([xobs, yobs, z]).T
|
||||
|
||||
f = Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu)
|
||||
g = Utils.sdiag(Casing._getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu))
|
||||
|
||||
return f,g
|
||||
|
||||
|
||||
|
||||
class Casing_DerivTest(unittest.TestCase):
|
||||
Tests.checkDerivative(CasingMagDipoleDeriv_r,np.ones(n)*10+np.random.randn(n),plotIt=False)
|
||||
Tests.checkDerivative(CasingMagDipoleDeriv_z,np.random.randn(n),plotIt=False)
|
||||
Tests.checkDerivative(CasingMagDipole2Deriv_z_r,np.ones(n)*10+np.random.randn(n),plotIt=False)
|
||||
Tests.checkDerivative(CasingMagDipole2Deriv_z_z,np.random.randn(n),plotIt=False)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -4,6 +4,7 @@ import simpegEM as EM
|
||||
from scipy.constants import mu_0
|
||||
|
||||
plotIt = False
|
||||
freq = 1e2
|
||||
|
||||
class FDEM_analyticTests(unittest.TestCase):
|
||||
|
||||
@@ -22,9 +23,9 @@ 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')
|
||||
Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, [rxList])
|
||||
Src0 = EM.FDEM.SrcFDEM_MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Tx0])
|
||||
survey = EM.FDEM.SurveyFDEM([Src0])
|
||||
|
||||
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
|
||||
prb.pair(survey)
|
||||
@@ -43,7 +44,7 @@ class FDEM_analyticTests(unittest.TestCase):
|
||||
self.prb = prb
|
||||
self.mesh = mesh
|
||||
self.m = m
|
||||
self.Tx0 = Tx0
|
||||
self.Src0 = Src0
|
||||
self.sig = sig
|
||||
|
||||
def test_Transect(self):
|
||||
@@ -51,20 +52,19 @@ class FDEM_analyticTests(unittest.TestCase):
|
||||
|
||||
u = self.prb.fields(self.m)
|
||||
|
||||
bfz = self.mesh.r(u[self.Tx0, 'b'],'F','Fz','M')
|
||||
|
||||
bfz = self.mesh.r(u[self.Src0, 'b'],'F','Fz','M')
|
||||
x = np.linspace(-55,55,12)
|
||||
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
|
||||
|
||||
P = self.mesh.getInterpolationMat(XYZ, 'Fz')
|
||||
|
||||
an = EM.Analytics.FDEM.hzAnalyticDipoleF(x, self.Tx0.freq, self.sig)
|
||||
an = EM.Analytics.FDEM.hzAnalyticDipoleF(x, self.Src0.freq, self.sig)
|
||||
|
||||
diff = np.log10(np.abs(P*np.imag(u[self.Tx0, 'b']) - mu_0*np.imag(an)))
|
||||
diff = np.log10(np.abs(P*np.imag(u[self.Src0, 'b']) - mu_0*np.imag(an)))
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
plt.plot(x,np.log10(np.abs(P*np.imag(u[self.Tx0, 'b']))))
|
||||
plt.plot(x,np.log10(np.abs(P*np.imag(u[self.Src0, 'b']))))
|
||||
plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r')
|
||||
plt.plot(x,diff,'g')
|
||||
plt.show()
|
||||
|
||||
@@ -1,103 +0,0 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import simpegEM as EM
|
||||
|
||||
class FieldsTest(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
|
||||
x = np.linspace(5,10,3)
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0.])
|
||||
txLoc = np.r_[0,0,0.]
|
||||
rxList0 = EM.FDEM.RxFDEM(XYZ, 'exi')
|
||||
Tx0 = EM.FDEM.TxFDEM(txLoc, 'VMD', 3., [rxList0])
|
||||
rxList1 = EM.FDEM.RxFDEM(XYZ, 'bxi')
|
||||
Tx1 = EM.FDEM.TxFDEM(txLoc, 'VMD', 3., [rxList1])
|
||||
rxList2 = EM.FDEM.RxFDEM(XYZ, 'bxi')
|
||||
Tx2 = EM.FDEM.TxFDEM(txLoc, 'VMD', 2., [rxList2])
|
||||
rxList3 = EM.FDEM.RxFDEM(XYZ, 'bxi')
|
||||
Tx3 = EM.FDEM.TxFDEM(txLoc, 'VMD', 2., [rxList3])
|
||||
Tx4 = EM.FDEM.TxFDEM(txLoc, 'VMD', 1., [rxList0, rxList1, rxList2, rxList3])
|
||||
txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
|
||||
survey = EM.FDEM.SurveyFDEM(txList)
|
||||
self.F = EM.FDEM.FieldsFDEM(mesh, survey)
|
||||
self.Tx0 = Tx0
|
||||
self.Tx1 = Tx1
|
||||
self.mesh = mesh
|
||||
self.XYZ = XYZ
|
||||
|
||||
def test_SetGet(self):
|
||||
F = self.F
|
||||
for freq in F.survey.freqs:
|
||||
nFreq = F.survey.nTxByFreq[freq]
|
||||
Txs = F.survey.getTransmitters(freq)
|
||||
e = np.random.rand(F.mesh.nE, nFreq)
|
||||
F[Txs, 'e'] = e
|
||||
b = np.random.rand(F.mesh.nF, nFreq)
|
||||
F[Txs, 'b'] = b
|
||||
if nFreq == 1:
|
||||
F[Txs, 'b'] = Utils.mkvc(b)
|
||||
if e.shape[1] == 1:
|
||||
e, b = Utils.mkvc(e), Utils.mkvc(b)
|
||||
self.assertTrue(np.all(F[Txs, 'e'] == e))
|
||||
self.assertTrue(np.all(F[Txs, 'b'] == b))
|
||||
F[Txs] = {'b':b,'e':e}
|
||||
self.assertTrue(np.all(F[Txs, 'e'] == e))
|
||||
self.assertTrue(np.all(F[Txs, 'b'] == b))
|
||||
|
||||
lastFreq = F[Txs]
|
||||
self.assertTrue(type(lastFreq) is dict)
|
||||
self.assertTrue(sorted([k for k in lastFreq]) == ['b','e'])
|
||||
self.assertTrue(np.all(lastFreq['b'] == b))
|
||||
self.assertTrue(np.all(lastFreq['e'] == e))
|
||||
|
||||
Tx_f3 = F.survey.getTransmitters(3.)
|
||||
self.assertTrue(F[Tx_f3,'b'].shape == (F.mesh.nF, 2))
|
||||
|
||||
b = np.random.rand(F.mesh.nF, 2)
|
||||
Tx_f0 = F.survey.getTransmitters(self.Tx0.freq)
|
||||
F[Tx_f0,'b'] = b
|
||||
self.assertTrue(F[self.Tx0]['b'].shape == (F.mesh.nF,))
|
||||
self.assertTrue(F[self.Tx0,'b'].shape == (F.mesh.nF,))
|
||||
self.assertTrue(np.all(F[self.Tx0,'b'] == b[:,0]))
|
||||
self.assertTrue(np.all(F[self.Tx1,'b'] == b[:,1]))
|
||||
|
||||
def test_assertions(self):
|
||||
freq = self.F.survey.freqs[0]
|
||||
Txs = self.F.survey.getTransmitters(freq)
|
||||
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nTxByFreq[freq])
|
||||
def fun(): self.F[Txs, 'b'] = bWrongSize
|
||||
self.assertRaises(ValueError, fun)
|
||||
def fun(): self.F[-999.]
|
||||
self.assertRaises(KeyError, fun)
|
||||
def fun(): self.F['notRight']
|
||||
self.assertRaises(KeyError, fun)
|
||||
def fun(): self.F[Txs,'notThere']
|
||||
self.assertRaises(KeyError, fun)
|
||||
|
||||
def test_FieldProjections(self):
|
||||
F = self.F
|
||||
for freq in F.survey.freqs:
|
||||
nFreq = F.survey.nTxByFreq[freq]
|
||||
Txs = F.survey.getTransmitters(freq)
|
||||
e = np.random.rand(F.mesh.nE, nFreq)
|
||||
b = np.random.rand(F.mesh.nF, nFreq)
|
||||
F[Txs] = {'b':b,'e':e}
|
||||
|
||||
Txs = F.survey.getTransmitters(freq)
|
||||
for ii, tx in enumerate(Txs):
|
||||
for jj, rx in enumerate(tx.rxList):
|
||||
dat = rx.projectFields(tx, self.mesh, F)
|
||||
self.assertTrue(dat.dtype == float)
|
||||
fieldType = rx.projField
|
||||
u = {'b':b[:,ii], 'e': e[:,ii]}[fieldType]
|
||||
real_or_imag = rx.projComp
|
||||
u = getattr(u, real_or_imag)
|
||||
gloc = rx.projGLoc
|
||||
d = self.mesh.getInterpolationMat(self.XYZ, gloc)*u
|
||||
self.assertTrue(np.all(dat == d))
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -3,6 +3,7 @@ from SimPEG import *
|
||||
import simpegEM as EM
|
||||
|
||||
plotIt = False
|
||||
tol = 1e-6
|
||||
|
||||
class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
@@ -22,9 +23,9 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
rxOffset = 40.
|
||||
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
|
||||
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
|
||||
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM([tx])
|
||||
survey = EM.TDEM.SurveyTDEM([src])
|
||||
|
||||
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
# self.prb.timeSteps = [1e-5]
|
||||
@@ -60,7 +61,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6)
|
||||
|
||||
V1 = Ahu[:,'e',1]
|
||||
self.assertLess(np.linalg.norm(V1), 1.e-6)
|
||||
return np.linalg.norm(V1) < 1.e-6
|
||||
|
||||
for i in range(2,prb.nT):
|
||||
|
||||
@@ -74,7 +75,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
V1 = Ahu[:,'e',i]
|
||||
V2 = prb.MeSigma*u[:,'e',i]
|
||||
# print np.linalg.norm(V1), np.linalg.norm(V2)
|
||||
self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
|
||||
return np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6
|
||||
|
||||
def test_AhVecVSMat_OneTS(self):
|
||||
|
||||
@@ -120,7 +121,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
u1 = prb.solveAh(sigma,f).tovec().flatten()
|
||||
u2 = sp.linalg.spsolve(A.tocsr(),f.tovec())
|
||||
|
||||
self.assertLess(np.linalg.norm(u1-u2),1e-8)
|
||||
self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
|
||||
|
||||
def test_solveAhVsAhVec(self):
|
||||
|
||||
@@ -156,7 +157,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
|
||||
print '\ntest_DerivG'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
return passed
|
||||
|
||||
def test_Deriv_dUdM(self):
|
||||
|
||||
@@ -171,8 +172,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
|
||||
print '\n'
|
||||
print 'test_Deriv_dUdM'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
|
||||
def test_Deriv_J(self):
|
||||
|
||||
@@ -188,8 +188,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
|
||||
print '\n'
|
||||
print 'test_Deriv_J'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
|
||||
|
||||
def test_projectAdjoint(self):
|
||||
prb = self.prb
|
||||
@@ -208,7 +207,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec())
|
||||
V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec())
|
||||
|
||||
self.assertLess((V1-V2)/np.abs(V1), 1e-6)
|
||||
self.assertTrue((V1-V2)/np.abs(V1) < tol)
|
||||
|
||||
def test_adjointAhVsAht(self):
|
||||
prb = self.prb
|
||||
@@ -227,7 +226,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec())
|
||||
V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec())
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
|
||||
|
||||
# def test_solveAhtVsAhtVec(self):
|
||||
# prb = self.prb
|
||||
@@ -292,7 +291,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
V1 = m.dot(prb.Gtvec(sigma, v, u))
|
||||
V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
|
||||
|
||||
def test_adjointJvecVsJtvec(self):
|
||||
mesh = self.mesh
|
||||
@@ -304,8 +303,10 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
V1 = d.dot(prb.Jvec(sigma, m))
|
||||
V2 = m.dot(prb.Jtvec(sigma, d))
|
||||
print 'AdjointTest', V1, V2
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
passed = np.abs(V1-V2)/np.abs(V1) < tol
|
||||
print 'AdjointTest', V1, V2, passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
+14
-17
@@ -22,11 +22,11 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
rxOffset = 40.
|
||||
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
|
||||
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
|
||||
src = EM.TDEM.SrcTDEM_VMD_MVP( [rx], loc=np.array([0., 0., 0.]))
|
||||
rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz')
|
||||
tx2 = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx2])
|
||||
src2 = EM.TDEM.SrcTDEM_VMD_MVP( [rx2], loc=np.array([0., 0., 0.]))
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM([tx,tx2])
|
||||
survey = EM.TDEM.SurveyTDEM([src,src2])
|
||||
|
||||
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
# self.prb.timeSteps = [1e-5]
|
||||
@@ -60,8 +60,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
|
||||
derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
|
||||
print '\ntest_DerivG'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
|
||||
def test_Deriv_dUdM(self):
|
||||
|
||||
@@ -76,8 +75,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
|
||||
print '\n'
|
||||
print 'test_Deriv_dUdM'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
|
||||
|
||||
def test_Deriv_J(self):
|
||||
|
||||
@@ -93,28 +91,27 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
|
||||
print '\n'
|
||||
print 'test_Deriv_J'
|
||||
passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
|
||||
self.assertTrue(passed)
|
||||
Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
|
||||
|
||||
def test_projectAdjoint(self):
|
||||
prb = self.prb
|
||||
survey = prb.survey
|
||||
nTx = survey.nTx
|
||||
nSrc = survey.nSrc
|
||||
mesh = self.mesh
|
||||
|
||||
# Generate random fields and data
|
||||
f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
|
||||
for i in range(prb.nT):
|
||||
f[:,'b',i] = np.random.rand(mesh.nF, nTx)
|
||||
f[:,'e',i] = np.random.rand(mesh.nE, nTx)
|
||||
f[:,'b',i] = np.random.rand(mesh.nF, nSrc)
|
||||
f[:,'e',i] = np.random.rand(mesh.nE, nSrc)
|
||||
d_vec = np.random.rand(survey.nD)
|
||||
d = Survey.Data(survey,v=d_vec)
|
||||
|
||||
# Check that d.T*Q*f = f.T*Q.T*d
|
||||
V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec())
|
||||
V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec())
|
||||
V2 = np.sum((f.tovec())*(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec()))
|
||||
|
||||
self.assertLess((V1-V2)/np.abs(V1), 1e-6)
|
||||
self.assertTrue((V1-V2)/np.abs(V1) < 1e-6)
|
||||
|
||||
def test_adjointGvecVsGtvec(self):
|
||||
mesh = self.mesh
|
||||
@@ -134,8 +131,8 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
v[:,'e',i] = np.random.rand(mesh.nE, 2)
|
||||
|
||||
V1 = m.dot(prb.Gtvec(sigma, v, u))
|
||||
V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
V2 = np.sum(v.tovec()*prb.Gvec(sigma, m, u).tovec())
|
||||
self.assertTrue(np.abs(V1-V2)/np.abs(V1) <1e-6)
|
||||
|
||||
def test_adjointJvecVsJtvec(self):
|
||||
mesh = self.mesh
|
||||
@@ -148,7 +145,7 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
V1 = d.dot(prb.Jvec(sigma, m))
|
||||
V2 = m.dot(prb.Jtvec(sigma, d))
|
||||
print 'AdjointTest', V1, V2
|
||||
self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
|
||||
self.assertTrue(np.abs(V1-V2)/np.abs(V1) < 1e-6)
|
||||
|
||||
|
||||
|
||||
@@ -4,7 +4,7 @@ import simpegEM as EM
|
||||
|
||||
plotIt = False
|
||||
|
||||
def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1):
|
||||
def getProb(meshType='CYL',rxTypes='bx,bz',nSrc=1):
|
||||
cs = 5.
|
||||
ncx = 20
|
||||
ncy = 6
|
||||
@@ -19,12 +19,12 @@ def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1):
|
||||
|
||||
rxOffset = 40.
|
||||
|
||||
txs = []
|
||||
for ii in range(nTx):
|
||||
srcs = []
|
||||
for ii in range(nSrc):
|
||||
rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')]
|
||||
txs += [EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', rxs)]
|
||||
srcs += [EM.TDEM.SrcTDEM_VMD_MVP(rxs,np.array([0., 0., 0.]))]
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM(txs)
|
||||
survey = EM.TDEM.SurveyTDEM(srcs)
|
||||
|
||||
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
# prb.timeSteps = [1e-5]
|
||||
@@ -68,8 +68,8 @@ class TDEM_bDerivTests(unittest.TestCase):
|
||||
def test_Jvec_bxbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz')))
|
||||
def test_Adjoint_bxbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz')))
|
||||
|
||||
def test_Jvec_bxbz_2tx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nTx=2)))
|
||||
def test_Adjoint_bxbz_2tx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nTx=2)))
|
||||
def test_Jvec_bxbz_2src(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nSrc=2)))
|
||||
def test_Adjoint_bxbz_2src(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nSrc=2)))
|
||||
|
||||
def test_Jvec_bxbzbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz,bz')))
|
||||
def test_Adjoint_bxbzbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz,bz')))
|
||||
|
||||
@@ -28,9 +28,10 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
|
||||
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
|
||||
|
||||
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
|
||||
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
|
||||
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
|
||||
# src = EM.TDEM.SrcTDEM([rx], loc=np.array([0., 0., 0.]))
|
||||
|
||||
survey = EM.TDEM.SurveyTDEM([tx])
|
||||
survey = EM.TDEM.SurveyTDEM([src])
|
||||
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
|
||||
prb.Solver = MumpsSolver
|
||||
|
||||
@@ -60,26 +61,26 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
|
||||
|
||||
class TDEM_bTests(unittest.TestCase):
|
||||
|
||||
def test_analitic_p2_CYL_50m(self):
|
||||
def test_analytic_p2_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+2) < 0.01)
|
||||
def test_analitic_p1_CYL_50m(self):
|
||||
def test_analytic_p1_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+1) < 0.01)
|
||||
def test_analitic_p0_CYL_50m(self):
|
||||
def test_analytic_p0_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+0) < 0.01)
|
||||
def test_analitic_m1_CYL_50m(self):
|
||||
def test_analytic_m1_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-1) < 0.01)
|
||||
def test_analitic_m2_CYL_50m(self):
|
||||
def test_analytic_m2_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-2) < 0.01)
|
||||
def test_analitic_m3_CYL_50m(self):
|
||||
def test_analytic_m3_CYL_50m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-3) < 0.02)
|
||||
|
||||
def test_analitic_p0_CYL_1m(self):
|
||||
def test_analytic_p0_CYL_1m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e+0) < 0.01)
|
||||
def test_analitic_m1_CYL_1m(self):
|
||||
def test_analytic_m1_CYL_1m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-1) < 0.01)
|
||||
def test_analitic_m2_CYL_1m(self):
|
||||
def test_analytic_m2_CYL_1m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-2) < 0.01)
|
||||
def test_analitic_m3_CYL_1m(self):
|
||||
def test_analytic_m3_CYL_1m(self):
|
||||
self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-3) < 0.02)
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,49 @@
|
||||
import numpy as np
|
||||
from scipy.constants import mu_0, epsilon_0
|
||||
|
||||
# useful params
|
||||
def omega(freq):
|
||||
"""Angular frequency, omega"""
|
||||
return 2.*np.pi*freq
|
||||
|
||||
def k(freq, sigma, mu=mu_0, eps=epsilon_0):
|
||||
""" Eq 1.47 - 1.49 in Ward and Hohmann """
|
||||
w = omega(freq)
|
||||
alp = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) + 1) )
|
||||
beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) )
|
||||
return alp - 1j*beta
|
||||
|
||||
# Constitutive relations
|
||||
def e_from_j(prob,j):
|
||||
eqLocs = prob._eqLocs
|
||||
if eqLocs is 'FE':
|
||||
MSigmaI = prob.MeSigmaI
|
||||
elif eqLocs is 'EF':
|
||||
MSigmaI = prob.MfRho
|
||||
return MSigmaI*j
|
||||
|
||||
def j_from_e(prob,e):
|
||||
eqLocs = prob._eqLocs
|
||||
if eqLocs is 'FE':
|
||||
MSigma = prob.MeSigma
|
||||
elif eqLocs is 'EF':
|
||||
MSigma = prob.MfRhoI
|
||||
return MSigma*e
|
||||
|
||||
def b_from_h(prob,h):
|
||||
eqLocs = prob._eqLocs
|
||||
if eqLocs is 'FE':
|
||||
MMu = prob.MfMuiI
|
||||
elif eqLocs is 'EF':
|
||||
MMu = prob.MeMu
|
||||
return MMu*h
|
||||
|
||||
def h_from_b(prob,b):
|
||||
eqLocs = prob._eqLocs
|
||||
if eqLocs is 'FE':
|
||||
MMuI = prob.MfMui
|
||||
elif eqLocs is 'EF':
|
||||
MMuI = prob.MeMuI
|
||||
return MMuI*b
|
||||
|
||||
|
||||
@@ -0,0 +1,203 @@
|
||||
from SimPEG import *
|
||||
from scipy.special import ellipk, ellipe
|
||||
from scipy.constants import mu_0, pi
|
||||
|
||||
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. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
|
||||
|
||||
:param numpy.ndarray srcLoc: Location of the source(s) (x, y, z)
|
||||
:param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh
|
||||
:param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list
|
||||
:param numpy.ndarray dipoleMoment: The vector dipole moment
|
||||
:rtype: numpy.ndarray
|
||||
: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):
|
||||
out[i] = MagneticDipoleVectorPotential(srcLoc, obsLoc, comp, dipoleMoment=dipoleMoment)
|
||||
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 MagneticDipoleVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], dipoleMoment=dipoleMoment)
|
||||
|
||||
if component == 'x':
|
||||
dimInd = 0
|
||||
elif component == 'y':
|
||||
dimInd = 1
|
||||
elif component == 'z':
|
||||
dimInd = 2
|
||||
else:
|
||||
raise ValueError('Invalid component')
|
||||
|
||||
srcLoc = np.atleast_2d(srcLoc)
|
||||
obsLoc = np.atleast_2d(obsLoc)
|
||||
dipoleMoment = np.atleast_2d(dipoleMoment)
|
||||
|
||||
nEdges = obsLoc.shape[0]
|
||||
nSrc = srcLoc.shape[0]
|
||||
|
||||
m = np.array(dipoleMoment).repeat(nEdges, axis=0)
|
||||
A = np.empty((nEdges, nSrc))
|
||||
for i in range(nSrc):
|
||||
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/(4*pi)) * mCr[:,dimInd]/(r**3)
|
||||
if nSrc == 1:
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
|
||||
def MagneticDipoleFields(srcLoc, obsLoc, component, moment=1., mu = mu_0):
|
||||
"""
|
||||
Calculate the vector potential of a set of magnetic dipoles
|
||||
at given locations 'ref. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
|
||||
|
||||
: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 moment: The vector dipole moment (vertical)
|
||||
:rtype: numpy.ndarray
|
||||
:return: The vector potential each dipole at each observation location
|
||||
"""
|
||||
|
||||
if component=='x':
|
||||
dimInd = 0
|
||||
elif component=='y':
|
||||
dimInd = 1
|
||||
elif component=='z':
|
||||
dimInd = 2
|
||||
else:
|
||||
raise ValueError('Invalid component')
|
||||
|
||||
srcLoc = np.atleast_2d(srcLoc)
|
||||
obsLoc = np.atleast_2d(obsLoc)
|
||||
moment = np.atleast_2d(moment)
|
||||
|
||||
nFaces = obsLoc.shape[0]
|
||||
nSrc = srcLoc.shape[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)
|
||||
r = np.sqrt((dR**2).sum(axis=1))
|
||||
if dimInd == 0:
|
||||
B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2)
|
||||
elif dimInd == 1:
|
||||
B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2)
|
||||
elif dimInd == 2:
|
||||
B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1)
|
||||
else:
|
||||
raise Exception("Not Implemented")
|
||||
if nSrc == 1:
|
||||
return B.flatten()
|
||||
return B
|
||||
|
||||
|
||||
|
||||
def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
|
||||
"""
|
||||
Calculate the vector potential of horizontal circular loop
|
||||
at given locations
|
||||
|
||||
:param numpy.ndarray srcLoc: Location of the source(s) (x, y, z)
|
||||
:param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh
|
||||
:param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list
|
||||
:param numpy.ndarray I: Input current of the loop
|
||||
:param numpy.ndarray radius: radius of the loop
|
||||
:rtype: numpy.ndarray
|
||||
:return: The vector potential each dipole at each observation location
|
||||
"""
|
||||
|
||||
if type(component) in [list, tuple]:
|
||||
out = range(len(component))
|
||||
for i, comp in enumerate(component):
|
||||
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, mu)
|
||||
|
||||
srcLoc = np.atleast_2d(srcLoc)
|
||||
obsLoc = np.atleast_2d(obsLoc)
|
||||
|
||||
n = obsLoc.shape[0]
|
||||
nSrc = srcLoc.shape[0]
|
||||
|
||||
if component=='z':
|
||||
A = np.zeros((n, nSrc))
|
||||
if nSrc ==1:
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
else:
|
||||
|
||||
A = np.zeros((n, nSrc))
|
||||
for i in range (nSrc):
|
||||
x = obsLoc[:, 0] - srcLoc[i, 0]
|
||||
y = obsLoc[:, 1] - srcLoc[i, 1]
|
||||
z = obsLoc[:, 2] - srcLoc[i, 2]
|
||||
r = np.sqrt(x**2 + y**2)
|
||||
m = (4 * radius * r) / ((radius + r)**2 + z**2)
|
||||
m[m > 1.] = 1.
|
||||
# m might be slightly larger than 1 due to rounding errors
|
||||
# but ellipke requires 0 <= m <= 1
|
||||
K = ellipk(m)
|
||||
E = ellipe(m)
|
||||
ind = (r > 0) & (m < 1)
|
||||
# % 1/r singular at r = 0 and K(m) singular at m = 1
|
||||
Aphi = np.zeros(n)
|
||||
# % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi.
|
||||
Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
|
||||
if component == 'x':
|
||||
A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] )
|
||||
elif component == 'y':
|
||||
A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] )
|
||||
else:
|
||||
raise ValueError('Invalid component')
|
||||
|
||||
if nSrc == 1:
|
||||
return A.flatten()
|
||||
return A
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG import Mesh
|
||||
import matplotlib.pyplot as plt
|
||||
cs = 20
|
||||
ncx, ncy, ncz = 41, 41, 40
|
||||
hx = np.ones(ncx)*cs
|
||||
hy = np.ones(ncy)*cs
|
||||
hz = np.ones(ncz)*cs
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz], 'CCC')
|
||||
srcLoc = np.r_[0., 0., 0.]
|
||||
Ax = MagneticLoopVectorPotential(srcLoc, mesh.gridEx, 'x', 200)
|
||||
Ay = MagneticLoopVectorPotential(srcLoc, mesh.gridEy, 'y', 200)
|
||||
Az = MagneticLoopVectorPotential(srcLoc, mesh.gridEz, 'z', 200)
|
||||
A = np.r_[Ax, Ay, Az]
|
||||
B0 = mesh.edgeCurl*A
|
||||
J0 = mesh.edgeCurl.T*B0
|
||||
|
||||
# mesh.plotImage(A, vType = 'Ex')
|
||||
# mesh.plotImage(A, vType = 'Ey')
|
||||
|
||||
mesh.plotImage(B0, vType = 'Fx')
|
||||
mesh.plotImage(B0, vType = 'Fy')
|
||||
mesh.plotImage(B0, vType = 'Fz')
|
||||
|
||||
# # mesh.plotImage(J0, vType = 'Ex')
|
||||
# mesh.plotImage(J0, vType = 'Ey')
|
||||
# mesh.plotImage(J0, vType = 'Ez')
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
@@ -1,3 +1,5 @@
|
||||
# import Sources
|
||||
# import Ana
|
||||
# import Solver
|
||||
import EMUtils
|
||||
import SrcUtils
|
||||
@@ -2,5 +2,5 @@
|
||||
import TDEM
|
||||
import FDEM
|
||||
import Base
|
||||
import Sources
|
||||
import Analytics
|
||||
import Utils
|
||||
|
||||
Reference in New Issue
Block a user