notebook.community

Edit and run 





In [1]:



    


from dolfin import *
import petsc4py, sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
import matplotlib.pylab as plt
import PETScIO as IO
import numpy as np
import scipy.sparse as sparse
import CheckPetsc4py as CP
import scipy.sparse.linalg as sparselin
import scipy as sp
import time



    











In [2]:



    


nn = 4
mesh = RectangleMesh(0, 0, 1, 1, nn, nn,'left')

order  = 2
parameters['reorder_dofs_serial'] = False
V = FunctionSpace(mesh, "N2curl", order)
Q = FunctionSpace(mesh, "CG", order)
W = MixedFunctionSpace([V,Q])



    


















    






DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.

















In [3]:



    


parameters['linear_algebra_backend'] = 'uBLAS'



    











In [4]:



    


b0 = Expression(("x[1]*x[1]*(x[1]-1)","x[0]*x[0]*(x[0]-1)"))
print V.dim(), Q.dim()



    


















    






112 25

















In [5]:



    


def boundary(x, on_boundary):
    return on_boundary
bcb = DirichletBC(V,b0, boundary)
bc = [bcb]



    











In [6]:



    


(v) = TrialFunction(V)
(u) = TestFunction(V)
(uMix,pMix) = TrialFunctions(W)
(vMix,qMix) = TestFunctions(W)
CurlCurl = Expression(("-6*x[1]+2","-6*x[0]+2"))+b0
f = CurlCurl



    











In [7]:



    


a = inner(curl(v),curl(u))*dx
m = inner(u,v)*dx
b = inner(vMix,grad(pMix))*dx



    











In [8]:



    


A = assemble(a)
M = assemble(m)
Ms = M.sparray()
A = A.sparray()



    


















    






DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.

















In [9]:



    


B = assemble(b)
B = B.sparray()[W.dim()-V.dim():,W.dim()-Q.dim():]



    


















    






DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.

















In [10]:



    


ksp = PETSc.KSP().create()
parameters['linear_algebra_backend'] = 'PETSc'
M = assemble(m)
M = CP.Assemble(M)
ksp.setOperators(M)
x = M.getVecLeft()
ksp.setFromOptions()
ksp.setType(ksp.Type.CG)
ksp.pc.setType(ksp.pc.Type.BJACOBI)



    


















    






DEBUG:FFC:Reusing form from cache.

















In [11]:



    


OptDB = PETSc.Options()
# OptDB["pc_factor_mat_ordering_type"] = "rcm"
# OptDB["pc_factor_mat_solver_package"] = "cholmod"
ksp.setFromOptions()
C = sparse.csr_matrix((V.dim(),Q.dim()))



    











In [38]:



    


C = sparse.csr_matrix((V.dim(),Q.dim()))
(v) = TrialFunction(V)
(u) = TestFunction(V)
tic()
for i in range(0,Q.dim()):
    uOut = Function(V)
    uu = Function(Q)
    x = M.getVecRight()
    zero = np.zeros((Q.dim(),1))[:,0]
    zero[i] = 1
    uu.vector()[:] = zero
    L = assemble(inner(u, grad(uu))*dx)
    rhs = IO.arrayToVec(L.array())
    ksp.solve(rhs,x)
#     x = project(grad(uu),V)
    P = x.array
    uOut.vector()[:] = P
    low_values_indices = np.abs(P) < 1e-3
    P[low_values_indices] = 0  
    pn = P.nonzero()[0]
    for j in range(0,len(pn)):
        C[pn[j],i] = P[pn[j]]
    del uu
print toc()



    


















    






DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.










    






0.13164305687

















In [35]:



    


B = A*C
print np.min(np.abs(B.toarray()))
print np.max(np.abs(B.toarray()))



    


















    






1.3850177567e-18
0.000363568441818

















In [57]:



    


rhs = PETSc.Vec().create()
rhs.createWithArray(B[:,0].toarray())
V.thisown



    


















    
Out[57]:








True

















In [23]:



    


print C > 1e-6



    


















    






  (111, 24)	True

















In [49]:



    


import scipy.sparse.linalg



    











In [51]:



    


M.norm()



    


















    
Out[51]:








90.68455705845307

















In [14]:



    


print np.min(np.abs(B.toarray()))
print np.max(np.abs(B.toarray()))



    


















    






0.0
16.0000106667

















In [ ]:



    






    











In [58]:



    


ksp.Type.BICG



    


















    
Out[58]:








'bicg'

















In [ ]:



    






    











In [59]:



    


plt.spy(C.todense())
plt.show()



    











In [ ]:



    






    











In [29]:



    


low_values_indices = x < 1e-4  # Where values are low
x[low_values_indices] = 0  # All



    











In [30]:



    


x



    


















    
Out[30]:








array([ 0.,  0.,  0., ...,  0.,  0.,  0.])

Content source: wathen/PhD

Similar notebooks:

notebook.community | gallery | about