

In [1]:

    
from SimPEG import *
from simpegPF.MagAnalytics import spheremodel, MagSphereAnalFun, CongruousMagBC
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib

Writing magnetostatic problem with secondary field formulation

- Here we focus on cell-centered system with FVM-Weak form

Step:1 Generating mesh and operators



In [10]:

    
cs = 25.
hxind = [(cs,5,-1.3), (cs, 41),(cs,5,1.3)]
hyind = [(cs,5,-1.3), (cs, 41),(cs,5,1.3)]
hzind = [(cs,5,-1.3), (cs, 40),(cs,5,1.3)]
M3 = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')

Step2: Generating model $\mu$

$\mu = \mu_0(1+\chi)$



In [11]:

    
mu0 = 4*np.pi*1e-7
chibkg = 0.
chiblk = 1.
chi = np.ones(M3.nC)*chibkg



In [12]:

    
sph_ind = spheremodel(M3, 0., 0., 0., 50)
chi[sph_ind] = chiblk
mu = (1.+chi)*mu0

Step3: Compute Boundary condition $\mathbf{B}_{sBC}$ using congruous sphere method



In [13]:

    
BC = [['neumann', 'neumann'], ['neumann', 'neumann'], ['neumann', 'neumann']]
Box = 1 # Primary field in x-direction (background)
Boy = 0 # Primary field in y-direction (background)
Boz = 0 # Primary field in z-direction (background)
B0 = np.r_[Box*np.ones(np.prod(M3.nFx)), Boy*np.ones(np.prod(M3.nFy)), Boz*np.ones(np.prod(M3.nFz))]



In [28]:

    
Bbc, Bbc_const = CongruousMagBC(M3, np.array([Box, Boy, Boz]), chi)



In [29]:

    
figsize(10,10)
M3.plotGrid()



In [30]:

    
M3.plotImage(np.log10(mu))









    Out[30]:





<matplotlib.collections.QuadMesh at 0x10166160>

Step4: get system matrix $\mathbf{A}$ and right handside

4.1 Backgrounds: Weak formulation for Poisson equation with cell-centered system

$$\frac{1}{\mu}\mathbf{B}_s = (\frac{1}{\mu}_0-\frac{1}{\mu})\mathbf{B}_0+\nabla\phi_s$$

$$\nabla \cdot \mathbf{B}_s = 0 $$

$$ (\mathbf{B}_s\cdot{\bar{\mathbf{n}}})_{\partial\Omega} = B_{sBC}$$

In discretized form we have

$$\mathbf{diag(v)}\mathbf{D}\mathbf{P}_{in}^T \mathbf{P}_{in}\mathbf{B}_s + \mathbf{diag(v)}\mathbf{D} \mathbf{P}_{out}^T B_{sBC}=0$$

$$\mathbf{M}^f_{\frac{1}{\mu}}\mathbf{B}_s = \mathbf{M}^f_{\frac{1}{\mu_0}}\mathbf{B}_0-\mathbf{M}^f_{\frac{1}{\mu}}\mathbf{B}_0 -\mathbf{P}_{in}\mathbf{P}_{in}^T\mathbf{D}^T\mathbf{diag(v)}\phi_s$$

$$ \mathbf{Div} = \mathbf{diag(v)}\mathbf{D}\mathbf{P}_{in}^T \mathbf{P}_{in}$$

$$ \mathbf{Div}\mathbf{B}_s + \mathbf{diag(v)}\mathbf{D} \mathbf{P}_{out}^T B_{sBC}=0$$

$$\mathbf{B}_s = (\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{M}^f_{\frac{1}{\mu_0}}\mathbf{B}_0-\mathbf{B}_0 -(\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{Div}^T\phi_s$$

$$\mathbf{A} = - \mathbf{Div}(\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{Div}^{T}$$

$$\mathbf{rhs} = - \mathbf{Div}(\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{M}^f_{\frac{1}{\mu_0}}\mathbf{B}_0 + \mathbf{Div}\mathbf{B}_0-\mathbf{diag}(v)\mathbf{D} \mathbf{P}_{out}^T B_{sBC}$$



In [31]:

    
Dface = M3.faceDiv
Pbc,Pin, Pout = M3.getBCProjWF(BC, discretization='CC')
Mc = Utils.sdiag(M3.vol)
Div = Mc*Dface*Pin.T*Pin
MfmuI = Utils.sdiag(1/M3.getFaceInnerProduct(1/mu).diagonal())
Mfmu0 = M3.getFaceInnerProduct(1/mu0)
A = -Div*MfmuI*Div.T
rhs = -Div*MfmuI*Mfmu0*B0 + Div*B0 - Mc*Dface*Pout.T*Bbc



In [32]:

    
print Mc.shape
print Dface.shape
print Pout.shape
print Bbc.shape









    



(130050, 130050)
(130050, 397851)
(15402, 397851)
(15402L,)



In [33]:

    
print MfmuI.shape, Div.T.shape









    



(397851, 397851) (397851, 130050)



In [36]:

    
# 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 [35]:

    
# u1 = PETSc.Vec().createSeq(M3.nC)
# u1.set(1)
# u1.normalize()
# basis = [u1]
# nullsp = PETSc.NullSpace().create(False, basis, comm=PETSc.COMM_SELF)



In [37]:

    
# %%time
# ksp = PETSc.KSP().create()
# pc = PETSc.PC().create()
# ksp.setOperators(Apetsc)
# ksp.setType(ksp.Type.BCGS)
# pc = ksp.getPC()
# pc.setType(pc.Type.BJACOBI)
# 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 [40]:

    
%%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)









    



Wall time: 1.53 s



In [41]:

    
print info
print M3.nC



In [42]:

    
B = MfmuI*Mfmu0*B0-B0-MfmuI*Div.T*phi
rhsa = A*phi
print info
print np.linalg.norm(rhs-rhsa)/np.linalg.norm(rhs)









    



0
8.03408261327e-07

BICG, qmr or BICGstab with Jacobi preconditioner works well (scipy), minres does not work ... Need to play with some other iterative solvers. Remeber that our system has null-space, which is constant



In [53]:

    
figsize(5,5)
M3.plotSlice(phi)









    Out[53]:





(<matplotlib.collections.QuadMesh at 0xa55a1d0>,)

Step4: Compute analytic function (sphere in whole space)

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}$$

Step5: Projection to receiver locations



In [54]:

    
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)))*80



In [55]:

    
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 [56]:

    
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'

Bxra, Byra, Bzra = MagSphereAnalFun(X, Y, Z, 50., 0., 0., 0., mu0, mu0*(1+chiblk), H0, 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 [57]:

    
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[57]:





[<matplotlib.lines.Line2D at 0xaddff28>]



In [58]:

    
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[58]:





<matplotlib.colorbar.Colorbar instance at 0x0000000012124588>



In [59]:

    
fig, axes = subplots(1,3, figsize=(16,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
id = 21
axes[0].plot(Y[:,id], Bxra[:,id], 'b', Y[:,id], Bxr[:,id], 'r.'); axes[0].set_title('$B^s_x$', fontsize = 16)
axes[1].plot(Y[:,id], Byra[:,id], 'b', Y[:,id], Byr[:,id], 'r.'); axes[1].set_title('$B^s_y$', fontsize = 16)
axes[2].plot(Y[:,id], Bzra[:,id], 'b', Y[:,id], Bzr[:,id], 'r.'); axes[2].set_title('$B^s_z$', fontsize = 16)
print X[1,id]



In [60]:

    
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], abs((Bxr[:,id]-Bxra[:,id])/(Bxra[:,id]+epsx)), 'r.')
axes[1].plot(Y[:,id], abs((Byr[:,id]-Byra[:,id])/(Byra[:,id]+epsy)), 'r.')
axes[2].plot(Y[:,id], abs((Bzr[:,id]-Bzra[:,id])/(Bzra[:,id]+epsz)), 'r.')









    Out[60]:





[<matplotlib.lines.Line2D at 0x213e39e8>]



In [61]:

    
fig, axes = subplots(1,1, figsize=(16,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
id = 21
# axes.plot(X[id,:], Bzra[id,:], 'b', X[id,:], Bzr[id,:], 'r.'); axes.set_title('$B^s_z$', fontsize = 16)
axes.plot(X[id,:], Bxra[id,:], 'b', X[id,:], Bxr[id,:], 'r.'); axes.set_title('$B^s_z$', fontsize = 16)
print Y[id,0]



In [77]:

    
print 1/abs(min(Bxra[id,:])/max(Bxra[id,:]))









    



0.200975085599

Thoughts

- Congruous sphere method works well, which means we need right boundary condition.

- I am not sure for real situation when we have multiple bodies in our domain ... ?

- Secondary field shows better accuracy than total field approach as expected

- Anyway every piece works well for forward modeling!!



In [59]: