In [1]:
from SimPEG import *
from MagAnalytics import spheremodel, MagSphereAnalFun, CongruousMagBC, MagSphereAnalFunA
%pylab inline


Warning: Python backend is being used for solver. Run setup.py from the command line.
Warning: mumps solver not available.
Warning: upgrade your scipy to 0.13.0
Populating the interactive namespace from numpy and matplotlib

Writing magnetostatic problem with total field formulation

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

Step:1 Generating mesh and operators


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

Step2: Compute Boundary indicies and set $\mathbf{B}_{bc}$


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]

Step3: Generating model $\mu$

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


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


C:\Users\SEOGI\AppData\Local\Enthought\Canopy\App\appdata\canopy-1.0.1.1189.win-x86_64\lib\site-packages\matplotlib\lines.py:483: RuntimeWarning: invalid value encountered in greater_equal
  return np.alltrue(x[1:]-x[0:-1]>=0)

In [9]:
M3.plotImage(np.log10(mu), imageType='CC')


Out[9]:
<matplotlib.collections.QuadMesh at 0xaa7f390>

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

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

$$\nabla \phi = \frac{1}{\mu}\mathbf{B}$$$$\nabla \cdot \mathbf{B} = q $$$$ (\phi)_{\partial\Omega} = \phi_{BC}$$

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

Things to remember

  • $\mathbf{P}_{BC}$ contains $\mathbf{diag}(\mathbf{area}_{BC})$ and a weired projection matrix, which has some +1 and -1 ....
  • Forward differential problem in weak form does not work well.., although inverse works really well
  • Mesh.getFaceInnerProduct or Mesh.getFaceMass should be used since we integrate $(\frac{1}{\mu}\mathbf{B}, \mathbf{F})$ in the cell, which means we need averaging from faces to cell centers
  • Weak formulation is a fancy representation of Finite Volume Approach. In terms of final system matrix \mathbf{A} is going to be exactly same, I guess..

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)


Wall time: 1.58 s

In [15]:
print info
print M3.nC


0
119646

In [16]:
B = -MfmuI*D.T*phi
rhsa = A*phi
print info
print np.linalg.norm(rhs-rhsa)/np.linalg.norm(rhs)


0
7.65372604673e-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 [17]:
figsize(16,5)
M3.plotImage(B, imageType='F')


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 plane


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


100.0 kang

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]:
[<matplotlib.lines.Line2D at 0xa11f550>]

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]:
<matplotlib.colorbar.Colorbar instance at 0x000000001244C608>

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]


[-285. -285. -285. -285. -285. -285. -285. -285. -285. -285. -285. -285.
 -285. -285. -285. -285. -285. -285. -285. -285. -285. -285. -285. -285.
 -285. -285. -285. -285. -285. -285. -285. -285. -285. -285. -285. -285.
 -285. -285. -285. -285. -285.]

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]:
[<matplotlib.lines.Line2D at 0x139014e0>]

Thoughts

  • It works well with non-uniform mesh!!
  • Actual accuray is ~10% relative error in secondary fields. As we pad more we can get better accuracy, since we did not consider secondary field at the boudnary ($\partial\Omega$).
  • Here, we can try primary secondary field approach to get better accuracy.
  • In addition, we can use the congruous sphere method to handle this secondary fields at boundaries