KKT system assembler with derivative UFL form operator

Question

Hi,
I have an elementar question. I would like to assembler KKT systems
using the FEniCS project. The code below illustrates my attempt.

from dolfin import *

mesh = UnitSquareMesh(100, 100)

U = FunctionSpace(mesh, "CG", 1) #foward solution space
V = FunctionSpace(mesh, "CG", 1) #ajoint solution space
M = FunctionSpace(mesh, "DG", 0) #control space

Z = MixedFunctionSpace([U, V, M])

z = Function(Z)
(u, lmbd, m) = split(z)

F = inner(grad(u), grad(lmbd))*dx - m*lmbd*dx

x = triangle.x
ud = exp(-1.0/(1.0-x[0]*x[0])-1.0/(1.0-x[1]*x[1]))
J = inner(u - ud, u - ud)*dx + Constant(1e-6)*inner(m, m)*dx

bc_u = DirichletBC(Z.sub(0), 0.0, "on_boundary")
bc_lmbd = DirichletBC(Z.sub(1), 0.0, "on_boundary")
bcs = [bc_u, bc_lmbd]

L = J + F
kkt = derivative(L, z, TestFunction(Z))    # !

KKT = assemble(kkt)    # ?

The derivative(L, z, TestFunction(Z)) returns a UFL form.
But the assembler(kkt) is a vector. How can I get the KKT system matrix
and the righ hand side vector? Thank you for your time.

( The documentation for derivative is rather sparse
http://fenicsproject.org/documentation/dolfin/dev/python/programmers-reference/fem/formmanipulations/derivative.html )

Comments

Please format your code properly (indent four spaces to get code markup).

Answer 1

The problem is that you do not use a dolfin.TrialFunction in your form definition. Hence your assembled kkt form becomes only a rank-1 tensor.

You have 2 options:

Either keep your current code and use the non-linear solver routine
in FEniCS to solve the KKT system:

   solve(kkt == 0, z)

Or replace "z" by a TrialFunction and then assemble the lhs and rhs:

kkt_lin = replace(kkt, {z: TrialFunction(Z)}) 
kkt_lhs = assemble(lhs(kkt_lin)) 
kkt_rhs = assemble(rhs(kkt_lin))

Comments

Simon,

The second option is perfect! Thank you very much.

I want to bypass the singular kkt systems with petsc krylov solvers orthogonal to null space:

http://fenicsproject.org/documentation/dolfin/1.3.0/python/demo/documented/singular-poisson/python/documentation.html

Here you have a simple example to ilustrate:

from dolfin import *

parameters["linear_algebra_backend"] = "PETSc"

#discrete (mixed) functions spaces
n = 100
mesh = UnitSquareMesh(n, n)

U = FunctionSpace(mesh, "CG", 1) #foward solution space
V = FunctionSpace(mesh, "CG", 1) #ajoint solution space

Z = MixedFunctionSpace([U, V])

#functions
z = Function(Z)
(u, lmbd) = split(z)

#variational problem
g = Expression("-10.0 * x[0] * x[0] * sin(25 * x[0])")
l = g * lmbd * ds
F = inner(grad(u), grad(lmbd))*dx - l

#functional
J = inner(u, u) * dx

#optimaly conditions
L = J + F
kkt = derivative(L, z, TestFunction(Z))

## WRONG SOLVER for singular problems
#solve(kkt == 0, z)

#kkt system assembler (without boudary conditions)
bcs = []
kkt_lin = replace(kkt, {z: TrialFunction(Z)})
kkt_lhs, kkt_rhs = assemble_system(lhs(kkt_lin), rhs(kkt_lin), bcs)

#Krylov solver (for singular kkt systems)
#solver = KrylovSolver(kkt_lhs, "gmres", "ilu")
solver = KrylovSolver(kkt_lhs, "tfqmr", "icc")

solver.parameters["absolute_tolerance"] = 1e-6
solver.parameters["relative_tolerance"] = 1e-5
solver.parameters["maximum_iterations"] = 1000
solver.parameters["monitor_convergence"] = True
#print solver.parameters.keys()

#vector that spans the null space
z = Function(Z)
null_vec = Vector(z.vector())
#V.dofmap().set(null_vec, 1.0)
Z.dofmap().set(null_vec, 1.0)
#U.dofmap().set(null_vec, 1.0)
null_vec *= 1.0/null_vec.norm("l2")

#null space basis object and attach to Krylov solver
null_space = VectorSpaceBasis([null_vec])
solver.set_nullspace(null_space)

#orthogonalize b with respect to the null space
null_space.orthogonalize(kkt_rhs);

#solve
solver.solve(z.vector(), kkt_rhs)

#plot
(u, lmbd) = split(z)
plot(u, title='u')
interactive()