In [1]:
from SimPEG import *
from MagAnalytics import spheremodel, MagSphereAnalFun, CongruousMagBC, MagSphereAnalFunA
%pylab inline
In [2]:
hxind = ((5,25,1.3),(41, 12.5),(5,25,1.3))
hyind = ((5,25,1.3),(41, 12.5),(5,25,1.3))
hzind = ((5,25,1.3),(40, 12.5),(1,25,1.3))
hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)
M3 = Mesh.TensorMesh([hx, hy, hz], [-sum(hx)/2,-sum(hy)/2,-sum(hz)/2])
In [3]:
BC = [['neumann', 'neumann'], ['neumann', 'neumann'], ['neumann', 'neumann']]
indxd, indxu, indyd, indyu, indzd, indzu = M3.faceBoundaryInd
ind = np.r_[(indxd | indxu), (indyd | indyu), (indzd | indzu)]
In [4]:
Box = 1 # Primary field in x-direction (background)
Boy = 0 # Primary field in y-direction (background)
Boz = 0 # Primary field in z-direction (background)
Bbcx = np.zeros(np.prod(M3.nFx))
Bbcy = np.zeros(np.prod(M3.nFy))
Bbcz = np.zeros(np.prod(M3.nFz))
Bbcx[indxd] = Box
Bbcx[indxu] = Box
Bbcy[indyd] = Boy
Bbcy[indyu] = Boy
Bbcz[indzd] = Boz
Bbcz[indzu] = Boz
Bbc = np.r_[Bbcx, Bbcy, Bbcz]
In [5]:
mu0 = 4*np.pi*1e-7
chibkg = 0.
chiblk = 0.01
chi = np.ones(M3.nC)*chibkg
In [6]:
sph_ind = spheremodel(M3, 0, 0, 0, 100)
chi[sph_ind] = chiblk
mu = (1.+chi)*mu0
In [7]:
Bbc = Bbc[ind] + CongruousMagBC(M3, np.array([Box, Boy, Boz]), chi)
In [8]:
figsize(10,10)
M3.plotGrid()
In [9]:
M3.plotImage(np.log10(mu), imageType='CC')
Out[9]:
or $$ (\mathbf{B}\cdot{\bar{\mathbf{n}}})_{\partial\Omega} = B_{BC}$$
In discretized form we have
$$\mathbf{diag(v)}\mathbf{D}\mathbf{P}_{in}^T \mathbf{P}_{in}\mathbf{B} + \mathbf{diag(v)}\mathbf{D} \mathbf{P}_{out}^T B_{BC}= \mathbf{diag(v)}q $$$$\mathbf{M}^f_{\frac{1}{\mu}}\mathbf{B} = -\mathbf{P}_{in}\mathbf{P}_{in}^T\mathbf{D}^T\mathbf{diag(v)}\phi + \mathbf{P}_{BC}\phi_{BC}$$$$\mathbf{B} = (\mathbf{M}^f_{\frac{1}{\mu}})^{-1}(-\mathbf{P}_{in}\mathbf{P}_{in}^T\mathbf{D}\mathbf{diag(v)}\phi +\mathbf{P}_{BC}\phi_{BC})$$$$\mathbf{diag(v)}\mathbf{D}\mathbf{P}_{in}^T \mathbf{P}_{in} (\mathbf{M}^f_{\frac{1}{\mu}})^{-1}(-\mathbf{P}_{in}\mathbf{P}_{in}^T\mathbf{D}\mathbf{diag(v)}\phi +\mathbf{P}_{BC}\phi_{BC}) + \mathbf{diag(v)}\mathbf{D} \mathbf{P}_{out}^T B_{BC}= \mathbf{diag(v)}q $$$$-\mathbf{diag(v)}\mathbf{D}\mathbf{P}_{in}^T \mathbf{P}_{in} (\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{P}_{in}\mathbf{P}_{in}^T\mathbf{D}\mathbf{diag(v)}\phi = \mathbf{diag(v)}q -\mathbf{diag(v)}\mathbf{D}\mathbf{P}_{in}^T \mathbf{P}_{in} (\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{P}_{BC}\phi_{BC} - \mathbf{diag(v)}\mathbf{D} \mathbf{P}_{out}^T B_{BC}$$$$\mathbf{A} = -\mathbf{diag(v)}\mathbf{D}\mathbf{P}_{in}^T \mathbf{P}_{in} (\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{P}_{in}\mathbf{P}_{in}^T\mathbf{D}\mathbf{diag(v)}$$$$\mathbf{rhs} = \mathbf{diag(v)}q -\mathbf{diag(v)}\mathbf{D}\mathbf{P}_{in}^T \mathbf{P}_{in} (\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{P}_{BC}\phi_{BC} - \mathbf{diag(v)}\mathbf{D} \mathbf{P}_{out}^T B_{BC} $$In magnetostatic case we do not have $\mathbf{q}$ and $\phi_{BC}$ and with
$$ \mathbf{Div} = \mathbf{diag(v)}\mathbf{D}\mathbf{P}_{in}^T \mathbf{P}_{in}$$$$\mathbf{A} = - D(\mathbf{M}^f_{\frac{1}{\mu}})^{-1}D^{T}$$$$\mathbf{rhs} = - \mathbf{diag(v)}\mathbf{D} \mathbf{P}_{out}^T B_{BC} $$$$\mathbf{B} = (\mathbf{M}^f_{\frac{1}{\mu}})^{-1}(-\mathbf{Div}^{T}\phi)$$
In [10]:
Dface = M3.faceDiv
Pbc,Pin, Pout = M3.getBCProjWF(BC, discretization='CC')
Mc = Utils.sdiag(M3.vol)
D = Mc*Dface*Pin.T*Pin
MfmuI = Utils.sdiag(1/M3.getFaceInnerProduct(mu = 1/mu).diagonal())
A = -D*MfmuI*D.T
rhs = -Mc*Dface*Pout.T*Bbc
In [11]:
# import petsc4py
# import sys
# from sys import getrefcount
# petsc4py.init(sys.argv)
# from petsc4py import PETSc
# import PETScIO as IO
# Apetsc = PETSc.Mat().createAIJ(size=A.shape,csr=(A.indptr, A.indices, A.data))
# bpetsc = IO.arrayToVec(rhs)
# xpetsc = IO.arrayToVec(0*rhs)
In [12]:
# u1 = PETSc.Vec().createSeq(M3.nC)
# u1.set(1)
# u1.normalize()
# basis = [u1]
# nullsp = PETSc.NullSpace().create(False, basis, comm=PETSc.COMM_SELF)
In [13]:
# %%time
# ksp = PETSc.KSP().create()
# pc = PETSc.PC().create()
# ksp.setOperators(Apetsc)
# ksp.setType(ksp.Type.BCGS)
# pc = ksp.getPC()
# pc.setType(pc.Type.SOR)
# OptDB = PETSc.Options()
# OptDB["ksp_rtol"] = 1e-8
# OptDB["pc_factor_levels"] = 1
# ksp.setFromOptions()
# ksp.view()
# ksp.solve(bpetsc, xpetsc)
# print ksp.its# print ksp.its
# phi = IO.vecToArray(xpetsc)
# print np.linalg.norm(A*phi-rhs)/norm(rhs)
In [14]:
%%time
m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(-1/A.diagonal()))
phi, info = sp.linalg.bicgstab(A, rhs, tol=1e-6, maxiter = 1000, M =m1)
# m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(np.sqrt(1/(Utils.mkvc(-A.diagonal())))))
# m2 = m1
# phi, info = sp.linalg.qmr(A, rhs, tol = 1e-5, M1 = m1, M2=m2)
In [15]:
print info
print M3.nC
In [16]:
B = -MfmuI*D.T*phi
rhsa = A*phi
print info
print np.linalg.norm(rhs-rhsa)/np.linalg.norm(rhs)
In [17]:
figsize(16,5)
M3.plotImage(B, imageType='F')
Outside of the sphere $(r>R)$
$$\mathbf{H}_1 = H_0 \hat{x} + H_0\frac{R}{r^5}\frac{\mu_2-\mu_1}{\mu_2+2\mu_1}[(2x^2-y^2-z^2)\hat{x}+(3xy)\hat{y}+(3xz)\hat{z}]$$$$H_{x1} = H_0 + H_0\frac{R}{r^5}\frac{\mu_2-\mu_1}{\mu_2+2\mu_1}(2x^2-y^2-z^2)$$$$H_{y1} = H_0\frac{R}{r^5}\frac{\mu_2-\mu_1}{\mu_2+2\mu_1}(3xy)$$$$H_{z1} = H_0\frac{R}{r^5}\frac{\mu_2-\mu_1}{\mu_2+2\mu_1}(3xz)$$Inside of the sphere $(r\le R)$
$$\mathbf{H}_2 = H_0\frac{3\mu_1}{\mu_2+2\mu_1}\hat{x}$$$$H_{x2} = H_0\frac{3\mu_1}{\mu_2+2\mu_1}$$
In [18]:
xr = np.linspace(-300, 300, 41)
yr = np.linspace(-300, 300, 41)
X, Y = np.meshgrid(xr, yr)
Z = np.ones((size(xr), size(yr)))*0.
In [19]:
rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]
Qfx = M3.getInterpolationMat(rxLoc,'Fx')
Qfy = M3.getInterpolationMat(rxLoc,'Fy')
Qfz = M3.getInterpolationMat(rxLoc,'Fz')
In [20]:
Bxr = np.reshape(Qfx*B, (size(xr), size(yr)), order='F')
Byr = np.reshape(Qfy*B, (size(xr), size(yr)), order='F')
Bzr = np.reshape(Qfz*B, (size(xr), size(yr)), order='F')
H0 = Box/mu0
flag = 'secondary'
if flag=='secondary':
Bxr = Bxr-Box
# Bxra, Byra, Bzra = MagSphereAnalFun(X, Y, Z, 100, 0., 0., 0., mu0, mu0*(1+chiblk), H0, flag)
Bxra, Byra, Bzra = MagSphereAnalFunA(X, Y, Z, 100., 0., 0., 0., chiblk, np.array([1., 0., 0.]), flag)
Bxra = np.reshape(Bxra, (size(xr), size(yr)), order='F')
Byra = np.reshape(Byra, (size(xr), size(yr)), order='F')
Bzra = np.reshape(Bzra, (size(xr), size(yr)), order='F')
In [22]:
figsize(10, 4)
plot(Utils.mkvc(Bxra))
plot(Utils.mkvc(Bxr), 'k:')
plot(Utils.mkvc(Byra))
plot(Utils.mkvc(Byr), 'k:')
plot(Utils.mkvc(Bzra))
plot(Utils.mkvc(Bzr), 'k:')
Out[22]:
In [24]:
fig, ax = subplots(3,2, figsize = (10,15))
dat1 = ax[0,0].imshow(Bxr); fig.colorbar(dat1, ax=ax[0,0])
dat2 = ax[0,1].imshow(Bxra); fig.colorbar(dat2, ax=ax[0,1])
dat3 = ax[1,0].imshow(Byr); fig.colorbar(dat3, ax=ax[1,0])
dat4 = ax[1,1].imshow(Byra); fig.colorbar(dat4, ax=ax[1,1])
dat5 = ax[2,0].imshow(Bzr); fig.colorbar(dat5, ax=ax[2,0])
dat6 = ax[2,1].imshow(Bzra); fig.colorbar(dat6, ax=ax[2,1])
Out[24]:
In [25]:
id = 21
fig, axes = subplots(1,3, figsize=(14,5))
epsx = np.linalg.norm(Utils.mkvc(Bxr))*1e-6
epsy = np.linalg.norm(Utils.mkvc(Byr))*1e-6
epsz = np.linalg.norm(Utils.mkvc(Bzr))*1e-6
axes[0].plot(Y[:,id], Bxra[:,id], 'b', Y[:,id], Bxr[:,id], 'r.')
axes[1].plot(Y[:,id], Byra[:,id], 'b', Y[:,id], Byr[:,id], 'r.')
axes[2].plot(Y[:,id], Bzra[:,id], 'b', Y[:,id], Bzr[:,id], 'r.')
print X[:,1]
In [26]:
fig, axes = subplots(1,3, figsize=(14,5))
epsx = np.linalg.norm(Utils.mkvc(Bxr))*1e-6
epsy = np.linalg.norm(Utils.mkvc(Byr))*1e-6
epsz = np.linalg.norm(Utils.mkvc(Bzr))*1e-6
axes[0].plot(Y[:,1], abs((Bxr[:,1]-Bxra[:,1])/(Bxra[:,1]+epsx)), 'r.')
axes[1].plot(Y[:,1], abs((Byr[:,1]-Byra[:,1])/(Byra[:,1]+epsy)), 'r.')
axes[2].plot(Y[:,1], abs((Bzr[:,1]-Bzra[:,1])/(Bzra[:,1]+epsz)), 'r.')
Out[26]: