In [1]:
from dolfin import *
import numpy as np
import scipy.sparse as sp
import numpy
from scipy2Trilinos import scipy_csr_matrix2CrsMatrix
from PyTrilinos import Epetra, ML, AztecOO
nn = 64
mesh = RectangleMesh(0, 0, 1, 1, nn, nn,'crossed')
order = 1
Magnetic = FunctionSpace(mesh, "N1curl", order)
Lagrange = FunctionSpace(mesh, "CG", order)
L= FunctionSpace(mesh, "DG", order-1)
C = sp.csr_matrix((Magnetic.dim(),Lagrange.dim()))
parameters['linear_algebra_backend'] = 'uBLAS'
u = Function(Lagrange)
# for i in range(0, Lagrange.dim() ):
# zero = np.zeros((Lagrange.dim(),1))[:,0]
# zero[i] = 1
# u.vector()[:] = zero
# uu = grad(u)
# Pv = project(uu,Magnetic)
# P = Pv.vector().array()
# index = P.nonzero()
# index = index[0]
# for x in range(0,len(index)):
# if np.abs(P[index[x]]) < 1e-3:
# P[index[x]] = 0
# print P.shape, C.shape
# pn = P.nonzero()[0]
# for j in range(0,len(pn)):
# C[pn[j],i] = P[pn[j]]
# del P
# print C.todense()
(u) = TrialFunction(Magnetic)
(v) = TestFunction(Magnetic)
(p) = TrialFunction(Lagrange)
(q) = TestFunction(Lagrange)
a = dot(curl(u),curl(v))*dx + inner(u, v)*dx
l = (inner(grad(p),grad(q))*dx)
Acurl = assemble(a)
Anode = assemble(l)
Acurl = Acurl.sparray()
Anode = Anode.sparray()
D = C
x = numpy.random.rand(Acurl.shape[1],1)
b = Acurl*x
x0 = numpy.ones((Acurl.shape[1],1))
comm = Epetra.PyComm()
A = scipy_csr_matrix2CrsMatrix(Acurl, comm)
L = scipy_csr_matrix2CrsMatrix(Anode, comm)
D = scipy_csr_matrix2CrsMatrix(C, comm)
x0 = Epetra.Vector(x0)
b = Epetra.Vector(b)
# MLList = {
# "default values" : "maxwell",
# "max levels" : 10,
# "prec type" : "MGV",
# "increasing or decreasing" : "decreasing",
# "aggregation: type" : "Uncoupled-MIS",
# "aggregation: damping factor" : 1.333,
# "eigen-analysis: type" : "cg",
# "eigen-analysis: iterations" : 10,
# "aggregation: edge prolongator drop threshold" : 0.0,
# "smoother: sweeps" : 1,
# "smoother: damping factor" : 1.0,
# "smoother: pre or post" : "both",
# "smoother: type" : "Hiptmair",
# "smoother: Hiptmair efficient symmetric" : True,
# "subsmoother: type" : "Chebyshev",
# "subsmoother: Chebyshev alpha" : 20.0,
# "subsmoother: node sweeps" : 4,
# "subsmoother: edge sweeps" : 4,
# "coarse: type" : "Amesos-KLU",
# "coarse: max size" : 128
# }
ML_Hiptmair = ML.MultiLevelPreconditioner(A)
ML_Hiptmair.ComputePreconditioner()
# ML_Hiptmair.ComputePreconditioner()
solver = AztecOO.AztecOO(A, x0, b)
solver.SetPrecOperator(ML_Hiptmair)
solver.SetAztecOption(AztecOO.AZ_solver, AztecOO.AZ_gmres);
solver.SetAztecOption(AztecOO.AZ_output, 16);
err = solver.Iterate(1550, 1e-5)
In [2]:
x0 = x0.array
In [6]:
x0=x0[:,0]
In [8]:
norm(x-x0)