In [1]:
#!/usr/bin/python
import petsc4py
import sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
Print = PETSc.Sys.Print
# from MatrixOperations import *
from dolfin import *
import numpy as np
import matplotlib.pylab as plt
import os
import scipy.io
from PyTrilinos import Epetra, EpetraExt, AztecOO, ML, Amesos
from scipy2Trilinos import scipy_csr_matrix2CrsMatrix
import PETScIO as IO
m = 6
errL2u = np.zeros((m-1,1))
errL2p = np.zeros((m-1,1))
l2uorder = np.zeros((m-1,1))
l2porder = np.zeros((m-1,1))
NN = np.zeros((m-1,1))
DoF = np.zeros((m-1,1))
Vdim = np.zeros((m-1,1))
Qdim = np.zeros((m-1,1))
Wdim = np.zeros((m-1,1))
iterations = np.zeros((m-1,1))
SolTime = np.zeros((m-1,1))
udiv = np.zeros((m-1,1))
nn = 2
dim = 2
Solving = 'Iterative'
ShowResultPlots = 'no'
ShowErrorPlots = 'no'
EigenProblem = 'no'
SavePrecond = 'no'
case = 3
parameters['linear_algebra_backend'] = 'PETSc'
xx = 4
nn = 2**(xx)
# Create mesh and define function space
nn = int(nn)
NN[xx-1] = nn
mesh = RectangleMesh(0, 0, 1, 1, nn, nn,'right')
parameters['reorder_dofs_serial'] = False
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
parameters['reorder_dofs_serial'] = False
W = V*Q
Vdim[xx-1] = V.dim()
Qdim[xx-1] = Q.dim()
Wdim[xx-1] = W.dim()
print "\n\nV: ",Vdim[xx-1],"Q: ",Qdim[xx-1],"W: ",Wdim[xx-1],"\n\n"
def boundary(x, on_boundary):
return on_boundary
if case == 1:
u0 = Expression(("20*x[0]*pow(x[1],3)","5*pow(x[0],4)-5*pow(x[1],4)"))
p0 = Expression("60*pow(x[0],2)*x[1]-20*pow(x[1],3)")
elif case == 2:
u0 = Expression(("sin(pi*x[1])","sin(pi*x[0])"))
p0 = Expression("sin(x[1]*x[0])")
elif case == 3:
u0 =Expression(("sin(x[1])*exp(x[0])","cos(x[1])*exp(x[0])"))
p0 = Expression("sin(x[0])*cos(x[1])")
elif case == 4:
u0 = Expression(("sin(x[1])*exp(x[0])","cos(x[1])*exp(x[0])"))
p0 = Expression("sin(x[0])*cos(x[1])")
bc = DirichletBC(W.sub(0),u0, boundary)
bcs = [bc]
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
if case == 1:
f = Expression(("120*x[0]*x[1]*(1-mu)","60*(pow(x[0],2)-pow(x[1],2))*(1-mu)"), mu = 1e0)
elif case == 2:
f = Expression(("pi*pi*sin(pi*x[1])+x[1]*cos(x[1]*x[0])","pi*pi*sin(pi*x[0])+x[0]*cos(x[1]*x[0])"))
elif case == 3:
f = Expression(("cos(x[0])*cos(x[1])","-sin(x[0])*sin(x[1])"))
elif case == 4:
f = Expression(("cos(x[1])*cos(x[0])","-sin(x[1])*sin(x[0])"))
N = FacetNormal(mesh)
h = CellSize(mesh)
h_avg =avg(h)
alpha = 10.0
gamma =10.0
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg =avg(h)
d = 0
a11 = inner(grad(v), grad(u))*dx
a12 = div(v)*p*dx
a21 = div(u)*q*dx
L1 = inner(v,f)*dx
a = a11-a12-a21
i = p*q*dx
tic()
AA, bb = assemble_system(a, L1, bcs)
A = as_backend_type(AA).mat()
print toc()
b = bb.array()
zeros = 0*b
del bb
bb = IO.arrayToVec(b)
x = IO.arrayToVec(zeros)
PP, Pb = assemble_system(a11+i,L1,bcs)
P = as_backend_type(PP).mat()
In [14]:
(p) = TrialFunction(Q)
(q) = TestFunction(Q)
(u) = TrialFunction(V)
(v) = TestFunction(V)
p11 = inner(grad(v), grad(u))*dx
p22 = p*q*dx
In [36]:
P11 = assemble(p11)
bc.apply(P11)
P22 = assemble(p22)
A = as_backend_type(P11).mat()
n = A.size[1]
n
Out[36]:
In [16]:
class apply:
def __init__(self,P):
# self.A = A
# self.B = B
self.P = P
def mult(self,x,y):
"y <- M * x"
# A = as_backend_type(self.A).mat()
# B = as_backend_type(self.B).mat()
P = as_backend_type(self.P).mat()
ksp = PETSc.KSP().Create()
pc = ksp.getPC()
# n = A.size[0]
# X = x.array()[:n]
ksp.setType('preonly')
pc.setType('lu')
ksp.setFromOptions()
ksp.setOperators(P)
ksp.solve(y, x)
In [23]:
del P
P = PETSc.Mat().create()
P.setSizes([W.dim(), W.dim()])
P.setType('python')
shell = apply(PP)
P.setPythonContext(shell)
P.setUp()
apply.apply?
In [18]:
ksp = PETSc.KSP().create()
ksp.setOperators(A,P)
ksp.setTolerances(1e-10)
ksp.setType('minres')
pc = ksp.getPC()
pc.setType('shell')
ksp.view()
ksp.solve(bb, x)
In [ ]:
P.set