In [2]:
#!/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 [1]:
class ApplyPC:
def __init__(self,P):
self.P = P
def apply(self, _pc, x, y):
"y <- M * x"
ksp = PETSc.KSP().Create()
pc = ksp.getPC()
ksp.setType('preonly')
pc.setType('lu')
ksp.setFromOptions()
ksp.setOperators(self.P)
ksp.solve(y, x)
In [3]:
del P
# P = PETSc.Mat().create()
# P.setSizes([W.dim(), W.dim()])
# P.setType('python')
# shell = apply(PP)
# P.setPythonContext(shell)
# P.setUp()
In [10]:
P = PETSc.Mat().create()
P.setSizes([W.dim(), W.dim()])
P.setType('python')
# # shell = ApplyPC(PP)
# # P.setPythonContext(shell)
# # P.setUp()
# shell = ApplyPC(PP)
# P.setPythonContext(shell)
ksp = PETSc.KSP()
ksp.create(comm=PETSc.COMM_WORLD)
ksp.setTolerances(1e-10)
ksp.setType('minres')
pc = ksp.getPC()
pc.setType('python')
shell = ApplyPC(PP)
pc.setPythonContext(shell)
ksp.setOperators(A,P)
OptDB = PETSc.Options()
module = __name__
factory = 'ApplyPC'
# OptDB['pc_python_type'] = '%s.%s' % (module, factory)
ksp.setFromOptions()
ksp.solve(bb, x)
In [14]:
pc = ksp.getPC()
pc.setType('python')
shell = ApplyPC(PP)
pc.setPythonContext(shell)
ksp.setOperators(A,P)
OptDB = PETSc.Options()
module = __name__
factory = 'ApplyPC'
# OptDB['pc_python_type'] = '%s.%s' % (module, factory)
ksp.setFromOptions()
ksp.solve(bb, x)