

In [2]:

    
from SimPEG import *
import scipy.sparse as sp
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib

Testing inhomogeneous boundary conditions

Here we test inhomogeneous boundary conditions with Finite volume approach.

$$\nabla \phi = \frac{1}{\mu}\mathbf{B}$$$$\nabla \cdot \mathbf{B} = q $$

In discretized form we have

$$\mathbf{M}^f\mathbf{G} \phi = \mathbf{M}^f_{\frac{1}{\mu}}\mathbf{B}$$$$\mathbf{D}\mathbf{M}^f\mathbf{B} = \mathbf{M}^cq $$$$\mathbf{D}\mathbf{M}^f(\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{M}^f\mathbf{G} \phi = \mathbf{M}^c q $$

1. Inhomogenous Dirichlet



In [3]:

    
n = 100
hx = np.ones(n)*pi/n
M = Mesh.TensorMesh([hx]) 
Phia = lambda x: -np.cos(x)
Ba = lambda x: np.sin(x)
qa = lambda x: np.cos(x)
M.setCellGradBC('dirichlet')
G = M.cellGrad
D = M.faceDiv
Bbc = M.cellGradBC
Mf = M.getFaceMass()



In [4]:

    
A = D*G
xsol = sp.linalg.minres(A, qa(M.gridCC)-D*Bbc*np.array([-1., 1.]))
q = A*Phia(M.gridCC)



In [5]:

    
plt.subplot(131)
plt.plot(M.gridCC, xsol[0], M.gridCC, Phia(M.gridCC), 'r.')
plt.subplot(132)
plt.plot(M.gridCC, qa(M.gridCC)+D*Bbc*np.array([1., -1.]), M.gridCC, q ,'r.')
plt.subplot(133)
plt.plot(M.gridN, G*xsol[0], 'b--', M.gridN, Ba(M.gridN), 'r-')









    Out[5]:





[<matplotlib.lines.Line2D at 0x4663e10>,
 <matplotlib.lines.Line2D at 0x4667050>]

2. Inhomogeneous neumann



In [6]:

    
n = 20
hx = np.ones(n)*pi/n
M = Mesh.TensorMesh([hx]) 
Phia = lambda x: np.sin(x)
Ba = lambda x: np.cos(x)
qa = lambda x: -np.sin(x)
M.setCellGradBC('neumann')
G = M.cellGrad
D = M.faceDiv
Mf = M.getFaceMass()



In [7]:

    
h = pi/n
bc = np.zeros(n)
bc[0] = -1.
bc[n-1] = -1.
bc = bc/h
A = -D*D.T



In [8]:

    
rhs = qa(M.gridCC)-bc
xsol = sp.linalg.minres(A, rhs)
q = A*Phia(M.gridCC)



In [9]:

    
figsize(12, 5)
plt.subplot(131)
plt.plot(M.gridCC, xsol[0], M.gridCC, Phia(M.gridCC), 'r.')
plt.subplot(132)
plt.plot(M.gridCC, rhs, 'b.', M.gridCC, q, 'r.')
plt.subplot(133)
plt.plot(M.gridN, G*xsol[0], 'b--', M.gridN, Ba(M.gridN), 'r-')









    Out[9]:





[<matplotlib.lines.Line2D at 0x4a75b90>,
 <matplotlib.lines.Line2D at 0x48d6c90>]

3. Questions?

For inhomogeneous neumann boundary condtion we need to specify bondary condtion for ou gradient operators?
Here I tried homogeneous dirichlet boundary condition.

Idea for implementing boundary condition for flux is

$$ \mathbf{D}\mathbf{B} = q $$$$ \mathbf{D}\mathbf{B}_{inner}+\mathbf{D}\mathbf{P}\mathbf{B}_{bc} = q $$

Here we need to use G for homogenous neuman case and $\mathbf{P}$ is a projection to boundary. Now, to make it work :), we need to generate

$$\mathbf{D}\mathbf{P}$$
$$\mathbf{B}_{bc}$$



In [10]:

    
from SimPEG import *
from test_boundary import ddxFaceDivBC, faceDivBC, faceBCind
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib

4. First generate $\mathbf{DP}$

Here I wrote the function "faceDivBC"



In [11]:

    
hxind = ((5,25,1.3),(20, 25),(5,25,1.3))
hx = Utils.meshTensors(hxind)
hx = hx/sum(hx)*np.pi
# hx = np.ones(n)*np.pi/n
M3 = Mesh.TensorMesh([hx, hx, hx])
M3.setCellGradBC('neumann')
G = M3.cellGrad
Bbc = M3.cellGradBC
D = M3.faceDiv
BC = [['neumann', 'neumann'], ['neumann', 'neumann'], ['neumann', 'neumann']]



In [12]:

    
Bbc.shape









    Out[12]:





(83700, 5400)

So, the number of boundary faces elements is 600 (100$\times$6)



In [13]:

    
print M3.gridFx.shape
print M3.gridFy.shape
print M3.gridFz.shape
print M3.nF



In [14]:

    
indxd, indxu, indyd, indyu, indzd, indzu = faceBCind(M3)
ind = np.r_[(indxd | indxu), (indyd | indyu), (indzd | indzu)]



In [15]:

    
Dbc = faceDivBC(M3, BC, ind)



In [16]:

    
figsize(16,4)
M3.plotImage(ind, imageType='F')



In [17]:

    
figsize(6,6)
M3.plotGrid()









    



/usr/lib/pymodules/python2.7/matplotlib/lines.py:483: RuntimeWarning: invalid value encountered in greater_equal
  return np.alltrue(x[1:]-x[0:-1]>=0)

5. Second we generate $\mathbf{B}_{bc}$

$$\nabla\cdot\vec{B} = q$$$$\nabla\phi = \frac{1}{\mu}\vec{B}$$$$\phi_{true} = sin(x)-cos(y)-cos(z)$$$$\vec{B}_{true} = cos(x)\hat{i}+sin(y)\hat{j}+sin(z)\hat{k}$$$$\vec{q}_{true} = -sin(x)+cos(y)+cos(z)$$

In discrete form

$$ \mathbf{G} \mathbf{M}^c\phi + \mathbf{G}_{bc}\phi_{bc}= \mathbf{M}^f_{\frac{1}{\mu}}\mathbf{B}_{inner}$$$$ \mathbf{B}_{inner} = (\mathbf{M}^f_{\frac{1}{\mu}})^{-1}(\mathbf{G} \mathbf{M}^c\phi + \mathbf{G}_{bc}\phi_{bc})$$$$ \mathbf{M}^c\mathbf{D}\mathbf{B}_{inner}+\mathbf{M}^c\mathbf{D}_{bc}\mathbf{B}_{bc} = \mathbf{M}^cq $$$$ \mathbf{M}^c\mathbf{D}(\mathbf{M}^f_{\frac{1}{\mu}})^{-1}(\mathbf{G} \mathbf{M}^c\phi +\mathbf{G}_{bc}\phi_{bc})+\mathbf{M}^c\mathbf{D}_{bc}\mathbf{B}_{bc} = \mathbf{M}^cq $$$$ \mathbf{M}^c\mathbf{D}(\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{G} \mathbf{M}^c\phi = \mathbf{M}^cq -\mathbf{M}^c\mathbf{D}(\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{G}_{bc}\phi_{bc}-\mathbf{M}^c\mathbf{D}_{bc}\mathbf{B}_{bc}$$

6. Now we test inhomogeneous neumann for 3D case

Since $\mathbf{G}_{bc}$ is going to be zero matrix in this case, we have

$$ \mathbf{M}^c\mathbf{D}(\mathbf{M}^f_{\frac{1}{\mu}})^{-1}\mathbf{G} \mathbf{M}^c\phi = \mathbf{M}^cq -\mathbf{M}^c\mathbf{D}_{bc}\mathbf{B}_{bc}$$

And by assuming uniform mesh and constant $\mu$

$$ \mathbf{D}\mathbf{G}\phi = q - \mathbf{D}_{bc}\mathbf{B}_{bc}$$



In [18]:

    
def phitrue(x,y,z):
    phi = np.sin(x)-np.cos(y)+np.cos(z)
    return phi

def Btrue(x,y,z):
    Bx = np.cos(x)
    By = np.sin(y)
    Bz = np.sin(z)
    B = np.r_[Bx, By, Bz]
    return B
def qtrue(x,y,z):
    q = -np.sin(x)+np.cos(y)+np.cos(z)
    return q



In [19]:

    
# def phitrue(x,y,z):
#     phi = np.sin(x)+np.sin(y)+np.sin(z)
#     return phi

# def Btrue(x,y,z):
#     Bx = np.cos(x)
#     By = np.cos(y)
#     Bz = np.cos(z)
#     B = np.r_[Bx, By, Bz]
#     return B
# def qtrue(x,y,z):
#     q = -np.sin(x)-np.sin(y)-np.sin(z)
#     return q



In [20]:

    
Mc = Utils.sdiag(M3.vol) 
Mfmu = Utils.sdiag(1/M3.getFaceMass().diagonal())
# A = Mc*D*Mfmu*G*Mc

A = D*G



In [21]:

    
qa = qtrue(M3.gridCC[:,0], M3.gridCC[:,1], M3.gridCC[:,2])
Ba = Btrue(M3.gridFx[:,0], M3.gridFy[:,1], M3.gridFz[:,2])
Phia = phitrue(M3.gridCC[:,0], M3.gridCC[:,1], M3.gridCC[:,2])



In [22]:

    
# M3.plotImage(Ba, imageType='F')
# M3.plotImage(qa, imageType='CC')



In [23]:

    
Bbcx = np.zeros(np.prod(M3.nFz))
Bbcy = np.zeros(np.prod(M3.nFz))
Bbcz = np.zeros(np.prod(M3.nFz))



In [24]:

    
# Bbcx[indxd] = 1
# Bbcx[indxu] = -1
# Bbcy[indyd] = 1
# Bbcy[indyu] = -1
# Bbcz[indzd] = 1
# Bbcz[indzu] = -1
# Bbc = np.r_[Bbcx, Bbcy, Bbcz]



In [25]:

    
Bbcx[indxd] = 1
Bbcx[indxu] = -1
Bbcy[indyd] = 0
Bbcy[indyu] = 0
Bbcz[indzd] = 0
Bbcz[indzu] = 0
Bbc = np.r_[Bbcx, Bbcy, Bbcz]



In [26]:

    
figsize(16,4)
M3.plotImage(Bbc, imageType='F')



In [27]:

    
# m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(-1/A.diagonal()))

If I get rid of minus sign there ... that does not work



In [28]:

    
# %%time
# rhs = Mc*qa-Mc*Dbc*Bbc[ind]
rhs = qa-Dbc*Bbc[ind]
# Try to use Jacobi preconditioner, but does not work well...
# x, info = sp.linalg.minres(A, rhs, M = m1,tol=1e-6, maxiter = 2000, show=False)
# x, info = sp.linalg.minres(A, rhs, tol=1e-6, maxiter = 1000)
# x, info = sp.linalg.minres(A, rhs, tol=1e-6, maxiter = 2000, show=False)



In [29]:

    
# rhs = Mc*qa-Mc*Dbc*Bbc[ind]
rhs = qa-Dbc*Bbc[ind]

import petsc4py
import sys
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 [30]:

    
%%time
ksp = PETSc.KSP().create()
pc = PETSc.PC().create()
ksp.setOperators(Apetsc)
ksp.setType(ksp.Type.GMRES)
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)
x = IO.vecToArray(xpetsc)









    



CPU times: user 11.1 s, sys: 409 µs, total: 11.1 s
Wall time: 11.1 s

Compute error for right handside



In [31]:

    
rhsa = A*x
print np.linalg.norm(rhs-rhsa)/np.linalg.norm(rhs)
print info









    



0.000208241391358
<function info at 0x2c3fd70>



In [32]:

    
print A.shape, Phia.shape, qa.shape, Dbc.shape, Bbc[ind].shape









    



(27000, 27000) (27000,) (27000,) (27000, 5400) (5400,)



In [33]:

    
fig, ax = plt.subplots(1,2, figsize=(16,4))
M3.plotImage(x, imageType='CC', ax = ax[0])
M3.plotImage(Phia, imageType='CC', ax = ax[1])









    Out[33]:





<matplotlib.collections.QuadMesh at 0x3f27810>



In [34]:

    
M3.plotImage(Ba, imageType='F')



In [44]:

    
Bnum = G*x
# Bnum = G*x

M3.plotImage(Bnum, imageType='F')



In [45]:

    
print Mc.diagonal()









    



[ 0.01307576  0.01005828  0.00773714 ...,  0.00773714  0.01005828
  0.01307576]

Compute error for B



In [46]:

    
print np.linalg.norm(Ba[~ind]-Bnum[~ind])/np.linalg.norm(Ba[~ind])









    



0.00104365335183



In [47]:

    
clf;
plot(Ba[~ind]); hold
plot(Bnum[~ind], 'r')









    Out[47]:





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

Compute error for $\phi$



In [38]:

    
print np.linalg.norm(x-Phia)/np.linalg.norm(Phia)









    



1.31230938236



In [39]:

    
plot(Phia); hold
plot(x, 'r')









    Out[39]:





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

Here we can recognize that there are some static... since we do not treat null-space of $A$ properly