wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-09-12 12:51:28 +08:00
Code Issues Packages Projects Releases Wiki Activity

Futurize 1, futurize 2, pasteurize.

Browse Source
This commit is contained in:
Brendan Smithyman
2016-07-16 14:17:02 -05:00
parent 362975d2bd
commit ca8d8f8c2d
197 changed files with 2618 additions and 1235 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/Analytics/DC.py
+12 -4
View File
@@ -1,3 +1,11 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
import numpy as np
from scipy.constants import mu_0, pi
from scipy import special
@@ -24,8 +32,8 @@ def DCAnalyticHalf(txloc, rxlocs, sigma, earth_type="wholespace"):
rM = np.sqrt( (M[:,0]-txloc[0])**2 + (M[:,1]-txloc[1])**2 + (M[:,2]-txloc[1])**2 )
rN = np.sqrt( (N[:,0]-txloc[0])**2 + (N[:,1]-txloc[1])**2 + (N[:,2]-txloc[1])**2 )
phiM = 1./(4*np.pi*rM*sigma)
phiN = 1./(4*np.pi*rN*sigma)
phiM = old_div(1.,(4*np.pi*rM*sigma))
phiN = old_div(1.,(4*np.pi*rN*sigma))
phi = phiM - phiN
if earth_type == "halfspace":
@@ -69,8 +77,8 @@ def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
Pleg.append(special.legendre(i, monic=0))
rho = 1./sigma
rho1 = 1./sigma1
rho = old_div(1.,sigma)
rho1 = old_div(1.,sigma1)
# Center of the sphere should be aligned in txloc in y-direction
yc = txloc[1]
SimPEG/EM/Analytics/FDEM.py
+5
View File
@@ -1,4 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
from scipy.constants import mu_0, pi
from scipy.special import erf
SimPEG/EM/Analytics/FDEMDipolarfields.py
+5
View File
@@ -1,4 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
from scipy.constants import mu_0, pi, epsilon_0
from scipy.special import erf
SimPEG/EM/Analytics/FDEMcasing.py
+12 -5
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import Utils, np
from scipy.constants import mu_0, epsilon_0
from SimPEG.EM.Utils.EMUtils import k
@@ -6,7 +13,7 @@ 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))
return np.sqrt(old_div(b, a)) * np.exp(-1j*k(freq,sigma,mu,eps)*(b-a))
def _r2(xyz):
return np.sum(xyz**2,1)
@@ -34,7 +41,7 @@ def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones
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)
return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + old_div(1., sqrtr2z2))
def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
@@ -47,7 +54,7 @@ def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones
sqrtr2z2 = np.sqrt(r2z2)
k2 = k(freq,sigma[2],mu[2],eps)
return -HertzZ*dxyz[:,2] /sqrtr2z2 * (1j*k2 + 1./sqrtr2z2)
return -HertzZ*dxyz[:,2] /sqrtr2z2 * (1j*k2 + old_div(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)
@@ -62,7 +69,7 @@ def _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.o
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)
return dHertzZdr*(old_div(-z,sqrtr2z2))*(1j*k2+old_div(1.,sqrtr2z2)) + HertzZ*(z*r/sqrtr2z2**3)*(1j*k2 + old_div(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)
@@ -77,7 +84,7 @@ def _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.o
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)
return (dHertzZdz*z + HertzZ)/sqrtr2z2*(-1j*k2 - old_div(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)
SimPEG/EM/Analytics/TDEM.py
+11 -4
View File
@@ -1,12 +1,19 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
from scipy.constants import mu_0, pi
from scipy.special import erf
def hzAnalyticDipoleT(r, t, sigma):
theta = np.sqrt((sigma*mu_0)/(4*t))
theta = np.sqrt(old_div((sigma*mu_0),(4*t)))
tr = theta*r
etr = erf(tr)
t1 = (9/(2*tr**2) - 1)*etr
t2 = (1/np.sqrt(pi))*(9/tr + 4*tr)*np.exp(-tr**2)
hz = (t1 - t2)/(4*pi*r**3)
t1 = (old_div(9,(2*tr**2)) - 1)*etr
t2 = (old_div(1,np.sqrt(pi)))*(old_div(9,tr) + 4*tr)*np.exp(-tr**2)
hz = old_div((t1 - t2),(4*pi*r**3))
return hz
SimPEG/EM/Analytics/__init__.py
+11 -5
View File
@@ -1,5 +1,11 @@
from TDEM import hzAnalyticDipoleT
from FDEM import hzAnalyticDipoleF
from FDEMcasing import *
from DC import DCAnalyticHalf, DCAnalyticSphere
from FDEMDipolarfields import *
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .TDEM import hzAnalyticDipoleT
from .FDEM import hzAnalyticDipoleF
from .FDEMcasing import *
from .DC import DCAnalyticHalf, DCAnalyticSphere
from .FDEMDipolarfields import *
SimPEG/EM/Base.py
+8 -1
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
@@ -11,7 +18,7 @@ class EMPropMap(Maps.PropMap):
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))
mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = old_div(1.,mu_0), propertyLink=('mu', Maps.ReciprocalMap))
class BaseEMProblem(Problem.BaseProblem):
SimPEG/EM/FDEM/FieldsFDEM.py
+60 -52
View File
@@ -1,3 +1,11 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
import scipy.sparse as sp
import SimPEG
@@ -42,7 +50,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total electric field
"""
if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting e from %s is not implemented' %list(self.knownFields.keys())[0])
return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
@@ -56,7 +64,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total magnetic flux density
"""
if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting b from %s is not implemented' %list(self.knownFields.keys())[0])
return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
@@ -70,7 +78,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total magnetic field
"""
if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting h from %s is not implemented' %list(self.knownFields.keys())[0])
return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
@@ -84,7 +92,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total current density
"""
if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting j from %s is not implemented' %list(self.knownFields.keys())[0])
return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
@@ -100,7 +108,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -118,7 +126,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
@@ -136,7 +144,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
@@ -154,7 +162,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
@@ -285,7 +293,7 @@ class Fields3D_e(FieldsFDEM):
C = self._edgeCurl
b = (C * eSolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
b[:,i] *= old_div(- 1.,(1j*omega(src.freq)))
s_m, _ = src.eval(self.prob)
b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * s_m
return b
@@ -331,8 +339,8 @@ class Fields3D_e(FieldsFDEM):
:return: current density
"""
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(aveE2CCV.shape[0], self._nC)) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
return VI * (aveE2CCV * (self._MeSigma * self._e(eSolution,srcList) ) )
def _jDeriv_u(self, src, du_dm_v, adjoint = False):
@@ -345,8 +353,8 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
return self._eDeriv_u(src, self._MeSigma.T * (self._aveE2CCV.T * (VI.T * du_dm_v) ), adjoint = adjoint)
@@ -364,8 +372,8 @@ class Fields3D_e(FieldsFDEM):
:return: product of the current density derivative with respect to the inversion model with a vector
"""
e = self[src, 'e']
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
return self._MeSigmaDeriv(e).T * (self._aveE2CCV.T * (VI.T * v)) + self._eDeriv_m(src, self._aveE2CCV.T * (VI.T * v), adjoint=adjoint)
@@ -382,8 +390,8 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: magnetic field
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # Number of Components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
return VI * (self._aveF2CCV * (self._MfMui * self._b(eSolution, srcList)))
@@ -397,8 +405,8 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # Number of Components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * du_dm_v))
return self._bDeriv_u(src, v, adjoint=adjoint)
@@ -414,8 +422,8 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the magnetic field derivative with respect to the inversion model with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # Number of Components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * v))
return self._bDeriv_m(src, v, adjoint=adjoint)
@@ -607,8 +615,8 @@ class Fields3D_b(FieldsFDEM):
:return: primary current density
"""
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) )
@@ -624,8 +632,8 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
return self._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) )
return VI * (self._aveE2CCV * (self._edgeCurl.T * ( self._MfMui * du_dm_v ) ) )
@@ -652,8 +660,8 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: magnetic field
"""
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
return VI * (self._aveF2CCV * (self._MfMui * self._b(bSolution, srcList)))
def _hDeriv_u(self, src, du_dm_v, adjoint=False):
@@ -667,8 +675,8 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
return self._MfMui.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v) )
@@ -827,7 +835,7 @@ class Fields3D_j(FieldsFDEM):
h = (self._edgeCurl.T * (self._MfRho * jSolution) )
for i, src in enumerate(srcList):
h[:,i] *= -1./(1j*omega(src.freq))
h[:,i] *= old_div(-1.,(1j*omega(src.freq)))
s_m,_ = src.eval(self.prob)
h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (s_m)
return self._MeMuI * h
@@ -890,8 +898,8 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: electric field
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
return VI * (self._aveF2CCV * (self._MfRho * self._j(jSolution, srcList)))
def _eDeriv_u(self, src, du_dm_v, adjoint=False):
@@ -904,8 +912,8 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
return self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) )
return VI * (self._aveF2CCV * (self._MfRho * du_dm_v))
@@ -921,8 +929,8 @@ class Fields3D_j(FieldsFDEM):
:return: product of the derivative of the electric field with respect to the model with a vector
"""
jSolution = Utils.mkvc(self[src,'jSolution'])
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
return self._MfRhoDeriv(jSolution).T * ( self._aveF2CCV.T * ( VI.T * v ) )
return VI * (self._aveF2CCV * (self._MfRhoDeriv(jSolution) * v))
@@ -936,8 +944,8 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: secondary magnetic flux density
"""
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) )
@@ -951,8 +959,8 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
return -1./(1j*omega(src.freq)) * self._MfRho.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) )
@@ -969,8 +977,8 @@ class Fields3D_j(FieldsFDEM):
:return: product of the derivative of the magnetic flux density with respect to the model with a vector
"""
jSolution = self[src,'jSolution']
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) # number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
if adjoint:
@@ -1151,8 +1159,8 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: electric field
"""
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
return VI * (self._aveF2CCV * (self._MfRho * self._j(hSolution, srcList)))
def _eDeriv_u(self, src, du_dm_v, adjoint=False):
@@ -1165,8 +1173,8 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) )
return VI * (self._aveF2CCV * (self._MfRho * self._edgeCurl * du_dm_v ))
@@ -1182,8 +1190,8 @@ class Fields3D_h(FieldsFDEM):
:return: product of the electric field derivative with respect to the inversion model with a vector
"""
hSolution = Utils.mkvc(self[src,'hSolution'])
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveF2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) )
return VI * (self._aveF2CCV * (self._MfRhoDeriv(self._edgeCurl * hSolution) * v ))
@@ -1198,8 +1206,8 @@ class Fields3D_h(FieldsFDEM):
:return: magnetic flux density
"""
h = self._h(hSolution, srcList)
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
return VI * (self._aveE2CCV * (self._MeMu * h))
@@ -1213,8 +1221,8 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(old_div(self._aveE2CCV.shape[0], self._nC)) #number of components
VI = sdiag(np.kron(np.ones(n), old_div(1.,self.prob.mesh.vol)))
if adjoint:
return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v ))
return VI * (self._aveE2CCV * (self._MeMu * du_dm_v))
SimPEG/EM/FDEM/ProblemFDEM.py
+8 -2
View File
@@ -1,7 +1,13 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
from SurveyFDEM import Survey as SurveyFDEM
from FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
from .SurveyFDEM import Survey as SurveyFDEM
from .FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Utils import omega
SimPEG/EM/FDEM/RxFDEM.py
+7
View File
@@ -1,3 +1,10 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from builtins import super
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG import sp
SimPEG/EM/FDEM/SrcFDEM.py
+10 -3
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
from SimPEG import Survey, Problem, Utils, np, sp
from scipy.constants import mu_0
from SimPEG.EM.Utils import *
@@ -375,7 +382,7 @@ class MagDipole(BaseSrc):
formulation = prob._formulation
if formulation is 'EB':
mui_s = prob.curModel.mui - 1./self.mu
mui_s = prob.curModel.mui - old_div(1.,self.mu)
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
elif formulation is 'HJ':
@@ -489,7 +496,7 @@ class MagDipole_Bfield(BaseSrc):
formulation = prob._formulation
if formulation is 'EB':
mui_s = prob.curModel.mui - 1./self.mu
mui_s = prob.curModel.mui - old_div(1.,self.mu)
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
elif formulation is 'HJ':
@@ -601,7 +608,7 @@ class CircularLoop(BaseSrc):
formulation = prob._formulation
if formulation is 'EB':
mui_s = prob.curModel.mui - 1./self.mu
mui_s = prob.curModel.mui - old_div(1.,self.mu)
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
SimPEG/EM/FDEM/SurveyFDEM.py
+8 -2
View File
@@ -1,10 +1,16 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.EM.Utils import *
from SimPEG.EM.Base import BaseEMSurvey
from scipy.constants import mu_0
from SimPEG.Utils import Zero, Identity
import SrcFDEM as Src
import RxFDEM as Rx
from . import SrcFDEM as Src
from . import RxFDEM as Rx
from SimPEG import sp
class Survey(BaseEMSurvey):
SimPEG/EM/FDEM/__init__.py
+11 -5
View File
@@ -1,5 +1,11 @@
from SurveyFDEM import Survey
import SrcFDEM as Src
import RxFDEM as Rx
from ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
from FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .SurveyFDEM import Survey
from . import SrcFDEM as Src
from . import RxFDEM as Rx
from .ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
from .FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h
SimPEG/EM/Static/DC/BoundaryUtils.py
+37 -30
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
def getxBCyBC_CC(mesh, alpha, beta, gamma):
@@ -31,15 +38,15 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
# h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
h_xm, h_xp = mesh.hx[0], mesh.hx[-1]
a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
a_xm = old_div(gamma_xm,(0.5*alpha_xm-old_div(beta_xm,h_xm)))
b_xm = old_div((0.5*alpha_xm+old_div(beta_xm,h_xm)),(0.5*alpha_xm-old_div(beta_xm,h_xm)))
a_xp = old_div(gamma_xp,(0.5*alpha_xp-old_div(beta_xp,h_xp)))
b_xp = old_div((0.5*alpha_xp+old_div(beta_xp,h_xp)),(0.5*alpha_xp-old_div(beta_xp,h_xp)))
xBC_xm = 0.5*a_xm
xBC_xp = 0.5*a_xp/b_xp
yBC_xm = 0.5*(1.-b_xm)
yBC_xp = 0.5*(1.-1./b_xp)
yBC_xp = 0.5*(1.-old_div(1.,b_xp))
xBC = np.r_[xBC_xm, xBC_xp]
yBC = np.r_[yBC_xm, yBC_xp]
@@ -62,24 +69,24 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
a_xm = old_div(gamma_xm,(0.5*alpha_xm-old_div(beta_xm,h_xm)))
b_xm = old_div((0.5*alpha_xm+old_div(beta_xm,h_xm)),(0.5*alpha_xm-old_div(beta_xm,h_xm)))
a_xp = old_div(gamma_xp,(0.5*alpha_xp-old_div(beta_xp,h_xp)))
b_xp = old_div((0.5*alpha_xp+old_div(beta_xp,h_xp)),(0.5*alpha_xp-old_div(beta_xp,h_xp)))
a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
a_ym = old_div(gamma_ym,(0.5*alpha_ym-old_div(beta_ym,h_ym)))
b_ym = old_div((0.5*alpha_ym+old_div(beta_ym,h_ym)),(0.5*alpha_ym-old_div(beta_ym,h_ym)))
a_yp = old_div(gamma_yp,(0.5*alpha_yp-old_div(beta_yp,h_yp)))
b_yp = old_div((0.5*alpha_yp+old_div(beta_yp,h_yp)),(0.5*alpha_yp-old_div(beta_yp,h_yp)))
xBC_xm = 0.5*a_xm
xBC_xp = 0.5*a_xp/b_xp
yBC_xm = 0.5*(1.-b_xm)
yBC_xp = 0.5*(1.-1./b_xp)
yBC_xp = 0.5*(1.-old_div(1.,b_xp))
xBC_ym = 0.5*a_ym
xBC_yp = 0.5*a_yp/b_yp
yBC_ym = 0.5*(1.-b_ym)
yBC_yp = 0.5*(1.-1./b_yp)
yBC_yp = 0.5*(1.-old_div(1.,b_yp))
sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
@@ -114,33 +121,33 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp)
a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
a_xm = old_div(gamma_xm,(0.5*alpha_xm-old_div(beta_xm,h_xm)))
b_xm = old_div((0.5*alpha_xm+old_div(beta_xm,h_xm)),(0.5*alpha_xm-old_div(beta_xm,h_xm)))
a_xp = old_div(gamma_xp,(0.5*alpha_xp-old_div(beta_xp,h_xp)))
b_xp = old_div((0.5*alpha_xp+old_div(beta_xp,h_xp)),(0.5*alpha_xp-old_div(beta_xp,h_xp)))
a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
a_ym = old_div(gamma_ym,(0.5*alpha_ym-old_div(beta_ym,h_ym)))
b_ym = old_div((0.5*alpha_ym+old_div(beta_ym,h_ym)),(0.5*alpha_ym-old_div(beta_ym,h_ym)))
a_yp = old_div(gamma_yp,(0.5*alpha_yp-old_div(beta_yp,h_yp)))
b_yp = old_div((0.5*alpha_yp+old_div(beta_yp,h_yp)),(0.5*alpha_yp-old_div(beta_yp,h_yp)))
a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm)
b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm)
a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp)
b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp)
a_zm = old_div(gamma_zm,(0.5*alpha_zm-old_div(beta_zm,h_zm)))
b_zm = old_div((0.5*alpha_zm+old_div(beta_zm,h_zm)),(0.5*alpha_zm-old_div(beta_zm,h_zm)))
a_zp = old_div(gamma_zp,(0.5*alpha_zp-old_div(beta_zp,h_zp)))
b_zp = old_div((0.5*alpha_zp+old_div(beta_zp,h_zp)),(0.5*alpha_zp-old_div(beta_zp,h_zp)))
xBC_xm = 0.5*a_xm
xBC_xp = 0.5*a_xp/b_xp
yBC_xm = 0.5*(1.-b_xm)
yBC_xp = 0.5*(1.-1./b_xp)
yBC_xp = 0.5*(1.-old_div(1.,b_xp))
xBC_ym = 0.5*a_ym
xBC_yp = 0.5*a_yp/b_yp
yBC_ym = 0.5*(1.-b_ym)
yBC_yp = 0.5*(1.-1./b_yp)
yBC_yp = 0.5*(1.-old_div(1.,b_yp))
xBC_zm = 0.5*a_zm
xBC_zp = 0.5*a_zp/b_zp
yBC_zm = 0.5*(1.-b_zm)
yBC_zp = 0.5*(1.-1./b_zp)
yBC_zp = 0.5*(1.-old_div(1.,b_zp))
sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
SimPEG/EM/Static/DC/FieldsDC.py
+9 -3
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.Utils import Identity, Zero
import numpy as np
@@ -9,7 +15,7 @@ class Fields(SimPEG.Problem.Fields):
def _phiDeriv(self, src, du_dm_v, v, adjoint=False):
if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
@@ -18,7 +24,7 @@ class Fields(SimPEG.Problem.Fields):
def _eDeriv(self, src, du_dm_v, v, adjoint=False):
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -26,7 +32,7 @@ class Fields(SimPEG.Problem.Fields):
def _jDeriv(self, src, du_dm_v, v, adjoint=False):
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
SimPEG/EM/Static/DC/FieldsDC_2D.py
+9 -3
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.Utils import Identity, Zero
import numpy as np
@@ -32,7 +38,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _phiDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
@@ -41,7 +47,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _eDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
@@ -49,7 +55,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _jDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)
SimPEG/EM/Static/DC/ProblemDC.py
+9 -3
View File
@@ -1,11 +1,17 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Problem, Utils
from SimPEG.EM.Base import BaseEMProblem
from SurveyDC import Survey
from FieldsDC import Fields, Fields_CC, Fields_N
from .SurveyDC import Survey
from .FieldsDC import Fields, Fields_CC, Fields_N
from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from BoundaryUtils import getxBCyBC_CC
from .BoundaryUtils import getxBCyBC_CC
class BaseDCProblem(BaseEMProblem):
SimPEG/EM/Static/DC/ProblemDC_2D.py
+13 -5
View File
@@ -1,11 +1,19 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import Problem, Utils
from SimPEG.EM.Base import BaseEMProblem
from SurveyDC import Survey, Survey_ky
from FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
from .SurveyDC import Survey, Survey_ky
from .FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from BoundaryUtils import getxBCyBC_CC
from .BoundaryUtils import getxBCyBC_CC
class BaseDCProblem_2D(BaseEMProblem):
@@ -182,8 +190,8 @@ class Problem2D_CC(BaseDCProblem_2D):
MfRhoIDeriv = self.MfRhoIDeriv
rho = self.curModel.rho
if adjoint:
return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v
return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v
return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(old_div(-1.,rho**2)))*v
return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(old_div(-1.,rho**2)))*v
def getRHS(self, ky):
"""
SimPEG/EM/Static/DC/RxDC.py
+7
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
import SimPEG
import numpy as np
from SimPEG.Utils import Zero, closestPoints
SimPEG/EM/Static/DC/SrcDC.py
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
# from SimPEG.EM.Base import BaseEMSurvey
from SimPEG.Utils import Zero, closestPoints, mkvc
SimPEG/EM/Static/DC/SurveyDC.py
+8 -2
View File
@@ -1,9 +1,15 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.EM.Base import BaseEMSurvey
from SimPEG import sp, Survey
from SimPEG.Utils import Zero, Identity
from RxDC import BaseRx
from SrcDC import BaseSrc
from .RxDC import BaseRx
from .SrcDC import BaseSrc
class Survey(BaseEMSurvey):
rxPair = BaseRx
SimPEG/EM/Static/DC/Utils.py
+9 -1
View File
@@ -1,3 +1,11 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
import numpy as np
def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
@@ -5,7 +13,7 @@ def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
import SimPEG.EM.Static.DC as DC
elocs = np.arange(0,aSpacing*nElecs,aSpacing)
elocs -= (nElecs*aSpacing - aSpacing)/2
elocs -= old_div((nElecs*aSpacing - aSpacing),2)
space = 1
WENNER = np.zeros((0,),dtype=int)
for ii in range(nElecs):
SimPEG/EM/Static/DC/__init__.py
+14 -8
View File
@@ -1,8 +1,14 @@
from ProblemDC import Problem3D_CC, Problem3D_N
from ProblemDC_2D import Problem2D_CC, Problem2D_N
from SurveyDC import Survey, Survey_ky
import SrcDC as Src #Pole
import RxDC as Rx
from FieldsDC import Fields_CC
from BoundaryUtils import getxBCyBC_CC
import Utils
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .ProblemDC import Problem3D_CC, Problem3D_N
from .ProblemDC_2D import Problem2D_CC, Problem2D_N
from .SurveyDC import Survey, Survey_ky
from . import SrcDC as Src #Pole
from . import RxDC as Rx
from .FieldsDC import Fields_CC
from .BoundaryUtils import getxBCyBC_CC
from . import Utils
SimPEG/EM/Static/IP/ProblemIP.py
+7 -1
View File
@@ -1,3 +1,9 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Problem, Utils, Maps, Mesh
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
@@ -5,7 +11,7 @@ from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from SimPEG.EM.Static.DC import getxBCyBC_CC
from SurveyIP import Survey
from .SurveyIP import Survey
class IPPropMap(Maps.PropMap):
"""
SimPEG/EM/Static/IP/SurveyIP.py
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.EM.Base import BaseEMSurvey
from SimPEG import sp, Survey
SimPEG/EM/Static/IP/__init__.py
+8 -2
View File
@@ -1,2 +1,8 @@
from ProblemIP import Problem3D_CC, Problem3D_N
from SurveyIP import Survey
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .ProblemIP import Problem3D_CC, Problem3D_N
from .SurveyIP import Survey
SimPEG/EM/Static/SIP/ProblemSIP.py
+11 -2
View File
@@ -1,3 +1,12 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import Problem, Utils, Maps, Mesh
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
@@ -5,7 +14,7 @@ from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from SimPEG.EM.Static.DC import getxBCyBC_CC
from SurveySIP import Survey, Data
from .SurveySIP import Survey, Data
class ColeColePropMap(Maps.PropMap):
"""
@@ -105,7 +114,7 @@ class BaseSIPProblem(BaseEMProblem):
JvAll = []
#Assume only eta and tau (eta first then tau)
# v = [2*Mx1]
v = v.reshape((int(v.size/2), 2), order='F')
v = v.reshape((int(old_div(v.size,2)), 2), order='F')
for tind in range(len(self.survey.times)):
t = self.survey.times[tind]
SimPEG/EM/Static/SIP/Regularization.py
+7
View File
@@ -1,3 +1,10 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Utils, Maps, Mesh, sp, np
from SimPEG.Regularization import BaseRegularization, Simple
SimPEG/EM/Static/SIP/RxSIP.py
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
import numpy as np
from SimPEG.Utils import Zero, closestPoints
SimPEG/EM/Static/SIP/SrcSIP.py
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
# from SimPEG.EM.Base import BaseEMSurvey
from SimPEG.Utils import Zero, closestPoints, mkvc
SimPEG/EM/Static/SIP/SurveySIP.py
+7
View File
@@ -1,3 +1,10 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import str
import SimPEG
from SimPEG.EM.Base import BaseEMSurvey
from SimPEG import np, sp, Survey, Utils
SimPEG/EM/Static/SIP/__init__.py
+11 -5
View File
@@ -1,5 +1,11 @@
from ProblemSIP import Problem3D_CC, Problem3D_N
from SurveySIP import Survey, Data
import SrcSIP as Src #Pole
import RxSIP as Rx
from Regularization import MultiRegularization
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .ProblemSIP import Problem3D_CC, Problem3D_N
from .SurveySIP import Survey, Data
from . import SrcSIP as Src #Pole
from . import RxSIP as Rx
from .Regularization import MultiRegularization
SimPEG/EM/Static/Utils/StaticUtils.py
+37 -28
View File
@@ -1,3 +1,12 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import np
from SimPEG.EM.Static import DC, IP
@@ -51,7 +60,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
# Create mid-point location
Cmid = Tx[0]
Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
Pmid = old_div((Rx[0][:,0] + Rx[1][:,0]),2)
if DCsurvey.mesh.dim == 2:
zsrc = Tx[1]
elif DCsurvey.mesh.dim ==3:
@@ -64,12 +73,12 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
NB = np.abs(Tx[1][0] - Rx[1][:,0])
# Create mid-point location
Cmid = (Tx[0][0] + Tx[1][0])/2
Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
Cmid = old_div((Tx[0][0] + Tx[1][0]),2)
Pmid = old_div((Rx[0][:,0] + Rx[1][:,0]),2)
if DCsurvey.mesh.dim == 2:
zsrc = (Tx[0][1] + Tx[1][1])/2
zsrc = old_div((Tx[0][1] + Tx[1][1]),2)
elif DCsurvey.mesh.dim ==3:
zsrc = (Tx[0][2] + Tx[1][2])/2
zsrc = old_div((Tx[0][2] + Tx[1][2]),2)
# Change output for dtype
if dtype == 'volt':
@@ -85,16 +94,16 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
elif stype == 'dpdp':
leg = data * 2*np.pi / ( 1/MA - 1/MB + 1/NB - 1/NA )
LEG.append(1./(2*np.pi) *( 1/MA - 1/MB + 1/NB - 1/NA ))
leg = data * 2*np.pi / ( old_div(1,MA) - old_div(1,MB) + old_div(1,NB) - old_div(1,NA) )
LEG.append(1./(2*np.pi) *( old_div(1,MA) - old_div(1,MB) + old_div(1,NB) - old_div(1,NA) ))
else:
print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
print("""dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """)
break
if dtype == 'appc':
leg = np.log10(abs(1./leg))
leg = np.log10(abs(old_div(1.,leg)))
rho = np.hstack([rho,leg])
elif dtype == 'appr':
@@ -103,15 +112,15 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
rho = np.hstack([rho,leg])
else:
print """dtype must be 'appr' | 'appc' | 'volt' """
print("""dtype must be 'appr' | 'appc' | 'volt' """)
break
midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
midx = np.hstack([midx, old_div(( Cmid + Pmid ),2) ])
if DCsurvey.mesh.dim==3:
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ])
midz = np.hstack([midz, old_div(-np.abs(Cmid-Pmid),2) + zsrc ])
elif DCsurvey.mesh.dim==2:
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ])
midz = np.hstack([midz, old_div(-np.abs(Cmid-Pmid),2) + zsrc ])
ax = axs
# Grid points
@@ -184,14 +193,14 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
# Mesure survey length and direction
dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
dl_x = old_div(( endl[1,0] - endl[0,0] ), dl_len)
dl_y = old_div(( endl[1,1] - endl[0,1] ), dl_len)
nstn = np.floor( dl_len / a )
nstn = np.floor( old_div(dl_len, a) )
# Compute discrete pole location along line
stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
stn_x = endl[0,0] + np.array(list(range(int(nstn))))*dl_x*a
stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*a
if mesh.dim==2:
ztop = mesh.vectorNy[-1]
@@ -230,15 +239,15 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
# Number of receivers to fit
nstn = np.min([np.floor( (AB - b) / a ) , n])
nstn = np.min([np.floor( old_div((AB - b), a) ) , n])
# Check if there is enough space, else break the loop
if nstn <= 0:
continue
# Compute discrete pole location along line
stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
stn_x = N[ii,0] + dl_x*b + np.array(list(range(int(nstn))))*dl_x*a
stn_y = N[ii,1] + dl_y*b + np.array(list(range(int(nstn))))*dl_y*a
# Create receiver poles
@@ -275,17 +284,17 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
max_y = endl[1,1] - dl_y * b
box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
box_w = box_l/2.
box_w = old_div(box_l,2.)
nstn = np.floor( box_l / a )
nstn = np.floor( old_div(box_l, a) )
# Compute discrete pole location along line
stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
stn_x = min_x + np.array(list(range(int(nstn))))*dl_x*a
stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*a
# Define number of cross lines
nlin = int(np.floor( box_w / a ))
lind = range(-nlin,nlin+1)
nlin = int(np.floor( old_div(box_w, a) ))
lind = list(range(-nlin,nlin+1))
ngrad = nstn * len(lind)
@@ -310,7 +319,7 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
srcClass = DC.Src.Dipole([rxClass], M[0,:], N[-1,:])
SrcList.append(srcClass)
else:
print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
print("""stype must be either 'pdp', 'dpdp' or 'gradient'. """)
return SrcList
SimPEG/EM/Static/Utils/__init__.py
+7 -1
View File
@@ -1 +1,7 @@
from StaticUtils import *
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .StaticUtils import *
SimPEG/EM/Static/__init__.py
+9 -3
View File
@@ -1,3 +1,9 @@
import DC
import IP
import SIP
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from . import DC
from . import IP
from . import SIP
SimPEG/EM/TDEM/BaseTDEM.py
+21 -14
View File
@@ -1,3 +1,10 @@
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Solver, Problem
from SimPEG.Problem import BaseTimeProblem
from SimPEG.EM.Utils import *
@@ -47,7 +54,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
self.waveformType = "GENERAL"
def fields(self, m):
if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
if self.verbose: print('%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50))
self.curModel = m
# Create a fields storage object
F = self._FieldsForward_pair(self.mesh, self.survey)
@@ -55,7 +62,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
# Set the initial conditions
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)
if self.verbose: print('%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50))
return F
def forward(self, m, RHS, F=None):
@@ -70,13 +77,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
if Ainv is not None:
Ainv.clean()
A = self.getA(tInd)
if self.verbose: print 'Factoring... (dt = %e)'%dt
if self.verbose: print('Factoring... (dt = %e)'%dt)
Ainv = self.Solver(A, **self.solverOpts)
if self.verbose: print 'Done'
if self.verbose: print('Done')
rhs = RHS(tInd, F)
if self.verbose: print ' Solving... (tInd = %d)'%tInd
if self.verbose: print(' Solving... (tInd = %d)'%tInd)
sol = Ainv * rhs
if self.verbose: print ' Done...'
if self.verbose: print(' Done...')
if sol.ndim == 1:
sol.shape = (sol.size,1)
F[:,self.solType,tInd+1] = sol
@@ -95,13 +102,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
if Ainv is not None:
Ainv.clean()
A = self.getA(tInd)
if self.verbose: print 'Factoring (Adjoint)... (dt = %e)'%dt
if self.verbose: print('Factoring (Adjoint)... (dt = %e)'%dt)
Ainv = self.Solver(A, **self.solverOpts)
if self.verbose: print 'Done'
if self.verbose: print('Done')
rhs = RHS(tInd, F)
if self.verbose: print ' Solving (Adjoint)... (tInd = %d)'%tInd
if self.verbose: print(' Solving (Adjoint)... (tInd = %d)'%tInd)
sol = Ainv * rhs
if self.verbose: print ' Done...'
if self.verbose: print(' Done...')
if sol.ndim == 1:
sol.shape = (sol.size,1)
F[:,self.solType,tInd+1] = sol
@@ -123,14 +130,14 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
* Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
"""
if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
if self.verbose: print('%s\nCalculating J(v)\n%s'%('*'*50,'*'*50))
self.curModel = m
if f is None:
f = self.fields(m)
p = self.Gvec(m, v, f)
y = self.solveAh(m, p)
Jv = self.survey.evalDeriv(f, v=y)
if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
if self.verbose: print('%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50))
return - mkvc(Jv)
def Jtvec(self, m, v, f=None):
@@ -148,7 +155,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
* Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\)
"""
if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
if self.verbose: print('%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50))
self.curModel = m
if f is None:
f = self.fields(m)
@@ -159,6 +166,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
p = self.survey.evalDeriv(f, v=v, adjoint=True)
y = self.solveAht(m, p)
w = self.Gtvec(m, y, f)
if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
if self.verbose: print('%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50))
return - mkvc(w)
SimPEG/EM/TDEM/SurveyTDEM.py
+11 -5
View File
@@ -1,7 +1,13 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Utils, Survey, np
from SimPEG.Survey import BaseSurvey
from SimPEG.EM.Utils import *
from BaseTDEM import FieldsTDEM
from .BaseTDEM import FieldsTDEM
class RxTDEM(Survey.BaseTimeRx):
@@ -87,7 +93,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
def getInitialFields(self, mesh):
"""Vertical magnetic dipole, magnetic vector potential"""
if self.waveformType == "STEPOFF":
print ">> Step waveform: Non-zero initial condition"
print(">> Step waveform: Non-zero initial condition")
if mesh._meshType is 'CYL':
if mesh.isSymmetric:
MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
@@ -99,7 +105,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
raise Exception('Unknown mesh for VMD')
return {"b": mesh.edgeCurl*MVP}
elif self.waveformType == "GENERAL":
print ">> General waveform: Zero initial condition"
print(">> General waveform: Zero initial condition")
return {"b": np.zeros(mesh.nF)}
else:
raise NotImplementedError("Only use STEPOFF or GENERAL")
@@ -127,7 +133,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
def getInitialFields(self, mesh):
"""Circular Loop, magnetic vector potential"""
if self.waveformType == "STEPOFF":
print ">> Step waveform: Non-zero initial condition"
print(">> Step waveform: Non-zero initial condition")
if mesh._meshType is 'CYL':
if mesh.isSymmetric:
MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
@@ -139,7 +145,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
raise Exception('Unknown mesh for CircularLoop')
return {"b": mesh.edgeCurl*MVP}
elif self.waveformType == "GENERAL":
print ">> General waveform: Zero initial condition"
print(">> General waveform: Zero initial condition")
return {"b": np.zeros(mesh.nF)}
else:
raise NotImplementedError("Only use STEPOFF or GENERAL")
SimPEG/EM/TDEM/TDEM_b.py
+12 -4
View File
@@ -1,7 +1,15 @@
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from .BaseTDEM import BaseTDEMProblem, FieldsTDEM
from SimPEG.Utils import mkvc, sdiag
import numpy as np
from SurveyTDEM import SurveyTDEM
from .SurveyTDEM import SurveyTDEM
class FieldsTDEM_e_from_b(FieldsTDEM):
@@ -69,14 +77,14 @@ class ProblemTDEM_b(BaseTDEMProblem):
:return: A
"""
dt = self.timeSteps[tInd]
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + (1.0/dt)*self.MfMui
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + (old_div(1.0,dt))*self.MfMui
def getRHS(self, tInd, F):
dt = self.timeSteps[tInd]
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
RHS = (old_div(1.0,dt))*self.MfMui*B_n[0,:,:] #TODO: This is a hack
return RHS
####################################################
SimPEG/EM/TDEM/__init__.py
+9 -3
View File
@@ -1,3 +1,9 @@
from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
from TDEM_b import ProblemTDEM_b
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
from .BaseTDEM import BaseTDEMProblem, FieldsTDEM
from .TDEM_b import ProblemTDEM_b
SimPEG/EM/Utils/AnalyticUtils.py
+18 -10
View File
@@ -1,3 +1,11 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from past.utils import old_div
from SimPEG import *
from scipy.special import ellipk, ellipe
from scipy.constants import mu_0, pi
@@ -17,7 +25,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMo
#TODO: break this out!
if type(component) in [list, tuple]:
out = range(len(component))
out = list(range(len(component)))
for i, comp in enumerate(component):
out[i] = MagneticDipoleVectorPotential(srcLoc, obsLoc, comp, dipoleMoment=dipoleMoment)
return np.concatenate(out)
@@ -49,7 +57,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMo
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)
A[:, i] = +(old_div(mu,(4*pi))) * mCr[:,dimInd]/(r**3)
if nSrc == 1:
return A.flatten()
return A
@@ -90,11 +98,11 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, moment=1., mu = mu_0):
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)
B[:, i] = +(old_div(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)
B[:, i] = +(old_div(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)
B[:, i] = +(old_div(mu,(4*pi))) /(r**3) * (3*dR[:,2]**2/r**2-1)
else:
raise Exception("Not Implemented")
if nSrc == 1:
@@ -118,7 +126,7 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
"""
if type(component) in [list, tuple]:
out = range(len(component))
out = list(range(len(component)))
for i, comp in enumerate(component):
out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius, mu)
return np.concatenate(out)
@@ -148,7 +156,7 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_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 = old_div((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
@@ -158,11 +166,11 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
# % 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])
Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(old_div(radius, r[ind])) *((1. - old_div(m[ind], 2.)) * K[ind] - E[ind])
if component == 'x':
A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] )
A[ind, i] = Aphi[ind] * (old_div(-y[ind], r[ind]) )
elif component == 'y':
A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] )
A[ind, i] = Aphi[ind] * ( old_div(x[ind], r[ind]) )
else:
raise ValueError('Invalid component')
SimPEG/EM/Utils/EMUtils.py
+9 -2
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import numpy as np
from scipy.constants import mu_0, epsilon_0
@@ -9,8 +16,8 @@ def omega(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) )
alp = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (old_div(sigma, (eps*w)))**2 ) + 1) )
beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (old_div(sigma, (eps*w)))**2 ) - 1) )
return alp - 1j*beta
SimPEG/EM/Utils/__init__.py
+8 -2
View File
@@ -1,2 +1,8 @@
from EMUtils import omega, k
from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .EMUtils import omega, k
from .AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
SimPEG/EM/Utils/testingUtils.py
+15 -7
View File
@@ -1,3 +1,11 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from past.utils import old_div
import unittest
from SimPEG import *
from SimPEG import EM
@@ -24,7 +32,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
else:
mapping = Maps.ExpMap(mesh)
x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid
x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + old_div(cs,4.) #don't sample right by the source, slightly off alignment from either staggered grid
XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5))
Rx0 = getattr(EM.FDEM.Rx, 'Point_' + comp[0])
if comp[2] == 'r':
@@ -58,7 +66,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*S_m, S_e))
if verbose:
print ' Fetching %s problem' % (fdemType)
print(' Fetching %s problem' % (fdemType))
if fdemType == 'e':
survey = EM.FDEM.Survey(Src)
@@ -83,7 +91,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
try:
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
except ImportError as e:
prb.Solver = SolverLU
return prb
@@ -94,7 +102,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
mesh = prb1.mesh
print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
print('Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp))
logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
mu = np.ones(mesh.nC)*MU
@@ -112,7 +120,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
d1 = survey1.dpred(m)
if verbose:
print ' Problem 1 solved'
print(' Problem 1 solved')
prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose)
@@ -121,11 +129,11 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
d2 = survey2.dpred(m)
if verbose:
print ' Problem 2 solved'
print(' Problem 2 solved')
r = d2-d1
l2r = l2norm(r)
tol = np.max([TOL*(10**int(np.log10(0.5* (l2norm(d1) + l2norm(d2)) ))),FLR])
print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol
print(l2norm(d1), l2norm(d2), l2r , tol, l2r < tol)
return l2r < tol
SimPEG/EM/__init__.py
+12 -6
View File
@@ -1,7 +1,13 @@
import TDEM
import FDEM
import Static
import Base
import Analytics
import Utils
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from . import TDEM
from . import FDEM
from . import Static
from . import Base
from . import Analytics
from . import Utils
from scipy.constants import mu_0, epsilon_0