1
1
""" This demo implements a Poisson equations solver
2
2
based on the demo "dolfin/demo/pde/poisson/python/demo.py"
3
in Dolfin using Epetra matrices, the AztecOO CG solver and ML
3
in Dolfin using Epetra matrices, the AztecOO CG solver and ML
7
7
__author__ = "Kent-Andre Mardal (kent-and@simula.no)"
11
11
# Test for Trilinos:
13
from PyTrilinos import Epetra, AztecOO, TriUtils, ML
13
from PyTrilinos import Epetra, AztecOO, TriUtils, ML
15
15
print "You Need to have PyTrilinos with Epetra, AztecOO, TriUtils and ML installed for this demo to run",
19
19
from dolfin import *
21
if not has_linear_algebra_backend("Epetra"):
21
if not has_la_backend("Epetra"):
22
22
print "*** Warning: Dolfin is not compiled with Trilinos linear algebra backend"
26
dolfin_set("linear algebra backend", "Epetra")
26
parameters["linear_algebra_backend"] = "Epetra"
28
28
# Create mesh and finite element
29
29
mesh = UnitSquare(20,20)
37
37
# Define variational problem
38
38
v = TestFunction(V)
39
39
u = TrialFunction(V)
40
f = Function(V,"500.0 * exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
40
f = Expression("500.0 * exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", V= V)
43
42
a = dot(grad(v), grad(u))*dx
47
46
# Define boundary condition
49
48
bc = DirichletBC(V, u0, DirichletBoundary())
51
50
# Create linear system
52
A, b = assemble_system(a, L, bc)
51
A, b = assemble_system(a, L, bc)
57
# Fetch underlying epetra objects
58
A_epetra = down_cast(A).mat()
59
b_epetra = down_cast(b).vec()
60
x_epetra = down_cast(U.vector()).vec()
56
# Fetch underlying epetra objects
57
A_epetra = down_cast(A).mat()
58
b_epetra = down_cast(b).vec()
59
x_epetra = down_cast(U.vector()).vec()
62
61
# Sets up the parameters for ML using a python dictionary
63
MLList = {"max levels" : 3,
62
MLList = {"max levels" : 3,
65
64
"smoother: type" : "ML symmetric Gauss-Seidel",
66
65
"aggregation: type" : "Uncoupled",
67
66
"ML validate parameter list" : False
70
# Create the preconditioner
69
# Create the preconditioner
71
70
Prec = ML.MultiLevelPreconditioner(A_epetra, False)
72
71
Prec.SetParameterList(MLList)
73
72
Prec.ComputePreconditioner()
75
# Create solver and solve system
74
# Create solver and solve system
76
75
Solver = AztecOO.AztecOO(A_epetra, x_epetra, b_epetra)
77
76
Solver.SetPrecOperator(Prec)
78
77
Solver.SetAztecOption(AztecOO.AZ_solver, AztecOO.AZ_cg)
79
78
Solver.SetAztecOption(AztecOO.AZ_output, 16)
80
79
Solver.Iterate(1550, 1e-5)