mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-28 13:38:02 +08:00
Brendan Smithyman
250 lines
8.2 KiB
Python
250 lines
8.2 KiB
Python
from __future__ import print_function
|
|
from __future__ import division
|
|
from __future__ import unicode_literals
|
|
from __future__ import absolute_import
|
|
from future import standard_library
|
|
standard_library.install_aliases()
|
|
import unittest
|
|
from SimPEG import *
|
|
from SimPEG import EM
|
|
from scipy.constants import mu_0
|
|
|
|
plotIt = False
|
|
tol_EBdipole = 1e-2
|
|
|
|
if plotIt:
|
|
import matplotlib.pylab
|
|
|
|
|
|
class FDEM_analyticTests(unittest.TestCase):
|
|
|
|
def setUp(self):
|
|
|
|
cs = 10.
|
|
ncx, ncy, ncz = 10, 10, 10
|
|
npad = 4
|
|
freq = 1e2
|
|
|
|
hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
|
|
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
|
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
|
|
mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC')
|
|
|
|
mapping = Maps.ExpMap(mesh)
|
|
|
|
x = np.linspace(-10,10,5)
|
|
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
|
|
rxList = EM.FDEM.Rx.Point_e(XYZ, orientation='x', component='imag')
|
|
Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
|
|
|
|
survey = EM.FDEM.Survey([Src0])
|
|
|
|
prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping)
|
|
prb.pair(survey)
|
|
|
|
try:
|
|
from pymatsolver import MumpsSolver
|
|
prb.Solver = MumpsSolver
|
|
except ImportError as e:
|
|
prb.Solver = SolverLU
|
|
|
|
sig = 1e-1
|
|
sigma = np.ones(mesh.nC)*sig
|
|
sigma[mesh.gridCC[:,2] > 0] = 1e-8
|
|
m = np.log(sigma)
|
|
|
|
self.prb = prb
|
|
self.mesh = mesh
|
|
self.m = m
|
|
self.Src0 = Src0
|
|
self.sig = sig
|
|
|
|
def test_Transect(self):
|
|
print('Testing Transect for analytic')
|
|
|
|
u = self.prb.fields(self.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.Src0.freq, self.sig)
|
|
|
|
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.Src0, 'b']))))
|
|
plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r')
|
|
plt.plot(x,diff,'g')
|
|
plt.show()
|
|
|
|
# We want the difference to be an orderMag less
|
|
# than the analytic solution. Note that right at
|
|
# the source, both the analytic and the numerical
|
|
# solution will be poor. Use plotIt up top to see that...
|
|
orderMag = 1.6
|
|
passed = np.abs(np.mean(diff - np.log10(np.abs(mu_0*np.imag(an))))) > orderMag
|
|
self.assertTrue(passed)
|
|
|
|
|
|
def test_CylMeshEBDipoles(self):
|
|
print('Testing CylMesh Electric and Magnetic Dipoles in a wholespace- Analytic: J-formulation')
|
|
sigmaback = 1.
|
|
mur = 2.
|
|
freq = 1.
|
|
skdpth = old_div(500.,np.sqrt(sigmaback*freq))
|
|
|
|
csx, ncx, npadx = 5, 50, 25
|
|
csz, ncz, npadz = 5, 50, 25
|
|
hx = Utils.meshTensor([(csx,ncx), (csx,npadx,1.3)])
|
|
hz = Utils.meshTensor([(csz,npadz,-1.3), (csz,ncz), (csz,npadz,1.3)])
|
|
mesh = Mesh.CylMesh([hx,1,hz], [0.,0.,old_div(-hz.sum(),2)]) # define the cylindrical mesh
|
|
|
|
if plotIt:
|
|
mesh.plotGrid()
|
|
|
|
# make sure mesh is big enough
|
|
self.assertTrue(mesh.hz.sum() > skdpth*2.)
|
|
self.assertTrue(mesh.hx.sum() > skdpth*2.)
|
|
|
|
SigmaBack = sigmaback*np.ones((mesh.nC))
|
|
MuBack = mur*mu_0*np.ones((mesh.nC))
|
|
|
|
# set up source
|
|
# test electric dipole
|
|
src_loc = np.r_[0.,0.,0.]
|
|
s_ind = Utils.closestPoints(mesh,src_loc,'Fz') + mesh.nFx
|
|
|
|
de = np.zeros(mesh.nF,dtype=complex)
|
|
de[s_ind] = old_div(1.,csz)
|
|
de_p = [EM.FDEM.Src.RawVec_e([],freq,old_div(de,mesh.area))]
|
|
|
|
dm_p = [EM.FDEM.Src.MagDipole([],freq,src_loc)]
|
|
|
|
|
|
# Pair the problem and survey
|
|
surveye = EM.FDEM.Survey(de_p)
|
|
surveym = EM.FDEM.Survey(dm_p)
|
|
|
|
mapping = [('sigma', Maps.IdentityMap(mesh)),('mu', Maps.IdentityMap(mesh))]
|
|
|
|
prbe = EM.FDEM.Problem3D_h(mesh, mapping=mapping)
|
|
prbm = EM.FDEM.Problem3D_e(mesh, mapping=mapping)
|
|
|
|
prbe.pair(surveye) # pair problem and survey
|
|
prbm.pair(surveym)
|
|
|
|
# solve
|
|
fieldsBackE = prbe.fields(np.r_[SigmaBack, MuBack]) # Done
|
|
fieldsBackM = prbm.fields(np.r_[SigmaBack, MuBack]) # Done
|
|
|
|
|
|
rlim = [20.,500.]
|
|
lookAtTx = de_p
|
|
r = mesh.vectorCCx[np.argmin(np.abs(mesh.vectorCCx-rlim[0])):np.argmin(np.abs(mesh.vectorCCx-rlim[1]))]
|
|
z = 100.
|
|
|
|
# where we choose to measure
|
|
XYZ = Utils.ndgrid(r, np.r_[0.], np.r_[z])
|
|
|
|
Pf = mesh.getInterpolationMat(XYZ, 'CC')
|
|
Zero = sp.csr_matrix(Pf.shape)
|
|
Pfx,Pfz = sp.hstack([Pf,Zero]),sp.hstack([Zero,Pf])
|
|
|
|
jn = fieldsBackE[de_p,'j']
|
|
bn = fieldsBackM[dm_p,'b']
|
|
|
|
Rho = Utils.sdiag(old_div(1.,SigmaBack))
|
|
Rho = sp.block_diag([Rho,Rho])
|
|
|
|
en = Rho*mesh.aveF2CCV*jn
|
|
bn = mesh.aveF2CCV*bn
|
|
|
|
ex,ez = Pfx*en, Pfz*en
|
|
bx,bz = Pfx*bn, Pfz*bn
|
|
|
|
# get analytic solution
|
|
exa, eya, eza = EM.Analytics.FDEM.ElectricDipoleWholeSpace(XYZ, src_loc, sigmaback, freq,orientation='Z',mu= mur*mu_0)
|
|
exa, eya, eza = Utils.mkvc(exa,2), Utils.mkvc(eya,2), Utils.mkvc(eza,2)
|
|
|
|
bxa, bya, bza = EM.Analytics.FDEM.MagneticDipoleWholeSpace(XYZ, src_loc, sigmaback, freq,orientation='Z',mu= mur*mu_0)
|
|
bxa, bya, bza = Utils.mkvc(bxa,2), Utils.mkvc(bya,2), Utils.mkvc(bza,2)
|
|
|
|
print(' comp, anayltic, numeric, num - ana, (num - ana)/ana')
|
|
print(' ex:', np.linalg.norm(exa), np.linalg.norm(ex), np.linalg.norm(exa-ex), old_div(np.linalg.norm(exa-ex),np.linalg.norm(exa)))
|
|
print(' ez:', np.linalg.norm(eza), np.linalg.norm(ez), np.linalg.norm(eza-ez), old_div(np.linalg.norm(eza-ez),np.linalg.norm(eza)))
|
|
|
|
print(' bx:', np.linalg.norm(bxa), np.linalg.norm(bx), np.linalg.norm(bxa-bx), old_div(np.linalg.norm(bxa-bx),np.linalg.norm(bxa)))
|
|
print(' bz:', np.linalg.norm(bza), np.linalg.norm(bz), np.linalg.norm(bza-bz), old_div(np.linalg.norm(bza-bz),np.linalg.norm(bza)))
|
|
|
|
if plotIt:
|
|
# Edipole
|
|
plt.subplot(221)
|
|
plt.plot(r,ex.real,'o',r,exa.real,linewidth=2)
|
|
plt.grid(which='both')
|
|
plt.title('Ex Real')
|
|
plt.xlabel('r (m)')
|
|
|
|
plt.subplot(222)
|
|
plt.plot(r,ex.imag,'o',r,exa.imag,linewidth=2)
|
|
plt.grid(which='both')
|
|
plt.title('Ex Imag')
|
|
plt.legend(['Num','Ana'],bbox_to_anchor=(1.5,0.5))
|
|
plt.xlabel('r (m)')
|
|
|
|
plt.subplot(223)
|
|
plt.plot(r,ez.real,'o',r,eza.real,linewidth=2)
|
|
plt.grid(which='both')
|
|
plt.title('Ez Real')
|
|
plt.xlabel('r (m)')
|
|
|
|
plt.subplot(224)
|
|
plt.plot(r,ez.imag,'o',r,eza.imag,linewidth=2)
|
|
plt.grid(which='both')
|
|
plt.title('Ez Imag')
|
|
plt.xlabel('r (m)')
|
|
|
|
plt.tight_layout()
|
|
|
|
# Bdipole
|
|
plt.subplot(221)
|
|
plt.plot(r,bx.real,'o',r,bxa.real,linewidth=2)
|
|
plt.grid(which='both')
|
|
plt.title('Bx Real')
|
|
plt.xlabel('r (m)')
|
|
|
|
plt.subplot(222)
|
|
plt.plot(r,bx.imag,'o',r,bxa.imag,linewidth=2)
|
|
plt.grid(which='both')
|
|
plt.title('Bx Imag')
|
|
plt.legend(['Num','Ana'],bbox_to_anchor=(1.5,0.5))
|
|
plt.xlabel('r (m)')
|
|
|
|
plt.subplot(223)
|
|
plt.plot(r,bz.real,'o',r,bza.real,linewidth=2)
|
|
plt.grid(which='both')
|
|
plt.title('Bz Real')
|
|
plt.xlabel('r (m)')
|
|
|
|
plt.subplot(224)
|
|
plt.plot(r,bz.imag,'o',r,bza.imag,linewidth=2)
|
|
plt.grid(which='both')
|
|
plt.title('Bz Imag')
|
|
plt.xlabel('r (m)')
|
|
|
|
plt.tight_layout()
|
|
|
|
self.assertTrue(old_div(np.linalg.norm(exa-ex),np.linalg.norm(exa)) < tol_EBdipole)
|
|
self.assertTrue(old_div(np.linalg.norm(eza-ez),np.linalg.norm(eza)) < tol_EBdipole)
|
|
|
|
self.assertTrue(old_div(np.linalg.norm(bxa-bx),np.linalg.norm(bxa)) < tol_EBdipole)
|
|
self.assertTrue(old_div(np.linalg.norm(bza-bz),np.linalg.norm(bza)) < tol_EBdipole)
|
|
|
|
|
|
|
|
if __name__ == '__main__':
|
|
unittest.main()
|