← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to sandbox/la/trilinos/demo3.py

  • Committer: Garth N. Wells
  • Date: 2008-03-29 09:34:25 UTC
  • Revision ID: gnw20@cam.ac.uk-20080329093425-3ea2vhjvccq56zvi
Add some basic build & install instructions.

Show diffs side-by-side

added added

removed removed

Lines of Context:
sandbox/la/trilinos/demo3.py
1
 
from PyTrilinos import Epetra, AztecOO, TriUtils, ML 
2
 
from sys import path 
3
 
path.append("../poisson")
4
 
 
5
 
 
6
 
 
7
 
from dolfin import *
8
 
from BlockLinearAlgebra import *
9
 
from Krylov import *
10
 
from MinRes import *
11
 
 
12
 
class MLPreconditioner: 
13
 
    def __init__(self, A): 
14
 
        # Sets up the parameters for ML using a python dictionary
15
 
        MLList = {
16
 
              "max levels"        : 3, 
17
 
              "output"            : 0, # as little output as possible
18
 
              "smoother: type"    : "ML symmetric Gauss-Seidel",
19
 
              "aggregation: type" : "Uncoupled",
20
 
              "ML validate parameter list" : False
21
 
        }
22
 
        ml_prec = ML.MultiLevelPreconditioner(A.mat(), False)
23
 
        ml_prec.SetParameterList(MLList)
24
 
        ml_prec.ComputePreconditioner()
25
 
        self.ml_prec = ml_prec
26
 
    def __mul__(self, b):
27
 
        x = b.copy()
28
 
        x.zero()
29
 
        err = self.ml_prec.ApplyInverse(b.vec(),x.vec())
30
 
        if not err == 0: 
31
 
            print "err ", err
32
 
            return -1 
33
 
        return x
34
 
 
35
 
class SaddlePrec: 
36
 
    def __init__(self, M, N):  
37
 
        self.n = 2 
38
 
        self.M = M
39
 
        self.N = N
40
 
        self.M_prec = MLPreconditioner(M)
41
 
        self.N_prec = MLPreconditioner(N)
42
 
 
43
 
    def __mul__(self, b):
44
 
        if (not isinstance(b, BlockVector)):
45
 
            raise TypeError, "Can not multiply DiagBlockMatrix with %s" % (str(type(b)))
46
 
        if (not b.n == self.n):
47
 
            raise ValueError, "dim(BlockVector) != dim(Prec) "  
48
 
 
49
 
        b0 = b.data[0]
50
 
        b1 = b.data[1]
51
 
 
52
 
        x0 = self.M_prec*b0
53
 
        x1 = self.N_prec*b1
54
 
 
55
 
        xx = BlockVector(x0, x1)
56
 
        return xx 
57
 
 
58
 
 
59
 
# Function for no-slip boundary condition for velocity
60
 
class Noslip(Function):
61
 
    def __init__(self, mesh):
62
 
        Function.__init__(self, mesh)
63
 
 
64
 
    def eval(self, values, x):
65
 
        values[0] = 0.0
66
 
        values[1] = 0.0
67
 
 
68
 
    def rank(self):
69
 
        return 1
70
 
 
71
 
    def dim(self, i):
72
 
        return 2
73
 
 
74
 
# Function for inflow boundary condition for velocity
75
 
class Inflow(Function):
76
 
    def __init__(self, mesh):
77
 
        Function.__init__(self, mesh)
78
 
 
79
 
    def eval(self, values, x):
80
 
        values[0] = -1.0
81
 
        values[1] = 0.0
82
 
 
83
 
    def rank(self):
84
 
        return 1
85
 
 
86
 
    def dim(self, i):
87
 
        return 2
88
 
 
89
 
mesh = Mesh("../../../data/meshes/dolfin-2.xml.gz")
90
 
sub_domains = MeshFunction("uint", mesh, "../../../demo/pde/stokes/taylor-hood/subdomains.xml.gz")
91
 
 
92
 
# create Taylor-Hood elements
93
 
shape = "triangle"
94
 
V = VectorElement("CG", shape, 2)
95
 
Q = FiniteElement("CG", shape, 1)
96
 
 
97
 
# Test and trial functions
98
 
v = TestFunction(V)
99
 
u = TrialFunction(V)
100
 
q = TestFunction(Q)
101
 
p = TrialFunction(Q)
102
 
tau = Function(Q, mesh, 0.0)
103
 
f = Function(V, mesh, 0.0)
104
 
 
105
 
 
106
 
# velocity terms 
107
 
a00 = dot(grad(v), grad(u))*dx
108
 
 
109
 
# Divergence constraint 
110
 
a10 = -div(u)*q*dx
111
 
 
112
 
# gradient of p 
113
 
a01 = -div(v)*p*dx 
114
 
 
115
 
# stabilization of p 
116
 
a11 = tau*dot(grad(p), grad(q))*dx
117
 
 
118
 
 
119
 
# right hand side 
120
 
L0 = dot(v, f)*dx 
121
 
#FIXME: does not remember stabilization terms on rhs right now
122
 
L1 = tau*q*dx  
123
 
 
124
 
 
125
 
# Create functions for boundary conditions
126
 
noslip = Noslip(mesh)
127
 
inflow = Inflow(mesh)
128
 
 
129
 
# No-slip boundary condition for velocity
130
 
bc0 = DirichletBC(noslip, sub_domains, 0)
131
 
 
132
 
# Inflow boundary condition for velocity
133
 
bc1 = DirichletBC(inflow, sub_domains, 1)
134
 
 
135
 
# Collect boundary conditions
136
 
bcs = [bc0, bc1]
137
 
 
138
 
# fetch the EpetraFactory backend
139
 
backend = EpetraFactory.instance()
140
 
 
141
 
# create matrices
142
 
A00 = assemble(a00, mesh, backend=backend)
143
 
A10 = assemble(a10, mesh, backend=backend)
144
 
A01 = assemble(a01, mesh, backend=backend)
145
 
A11 = assemble(a11, mesh, backend=backend)
146
 
 
147
 
# create right hand sides
148
 
b0 = assemble(L0, mesh, backend=backend)
149
 
b1 = assemble(L1, mesh, backend=backend)
150
 
 
151
 
# apply bc 
152
 
for bc in bcs: 
153
 
    bc.apply(A00, b0, a00)
154
 
    bc.zero(A01, a00)
155
 
 
156
 
 
157
 
# ensure that bc are in start vector
158
 
x0 = b0.copy()
159
 
x1 = b1.copy()
160
 
 
161
 
# Functions
162
 
U = Function(V, mesh, x0)
163
 
P = Function(Q, mesh, x1)
164
 
 
165
 
 
166
 
# create block matrix
167
 
AA = BlockMatrix(2,2)
168
 
AA[0,0] = A00 
169
 
AA[1,0] = A10 
170
 
AA[0,1] = A01 
171
 
AA[1,1] = A11 
172
 
 
173
 
#file = File("A00.m"); file << A00
174
 
#file = File("A10.m"); file << A10
175
 
#file = File("A01.m"); file << A01
176
 
#file = File("A11.m"); file << A11
177
 
 
178
 
# create block vector 
179
 
bb = BlockVector(2)  
180
 
bb[0] = b0
181
 
bb[1] = b1
182
 
 
183
 
 
184
 
# create solution vector
185
 
xx = BlockVector(2)  
186
 
xx[0] = x0
187
 
xx[1] = x1
188
 
 
189
 
 
190
 
# create matrices needed in the preconditioner 
191
 
c11 = p*q*dx
192
 
C11 = assemble(c11, mesh, backend=backend)
193
 
#file = File("C11.m"); file << C11
194
 
 
195
 
 
196
 
# create block preconditioner
197
 
BB = SaddlePrec(AA[0,0], C11)
198
 
 
199
 
xx, i, rho  = precondBiCGStab(BB, AA, xx, bb, 10e-8, True, 200)
200
 
 
201
 
 
202
 
 
203
 
# Plot solution
204
 
#plot(mesh)
205
 
plot(U)
206
 
plot(P)
207
 
interactive()
208
 
 
209
 
 
210
 
 
211