box constrained (PETSc) krylov solvers in FEniCS

Question

Is it possible to use box constrained Krylov solvers in FEniCS project?

Can I put a non-negative restriction in z.vector() of the code below?

from dolfin import *

parameters["linear_algebra_backend"] = "PETSc"

n = 100
mesh = UnitSquareMesh(n, n)

U = FunctionSpace(mesh, "CG", 1)
V = FunctionSpace(mesh, "CG", 1)

Z = MixedFunctionSpace([U, V])

z = Function(Z)
(u, v) = split(z)

#variational problem
g = Expression("-10.0 * x[0] * x[0] * sin(25 * x[0])")

l = g * v * ds
F = inner(grad(u), grad(v))*dx - l

#functional
J = inner(u, u) * dx
# J = u * dx    # Is it possible to put non-negative constrains on u.vector() ? 

#optimaly conditions
L = J + F    # What about J = u * dx and subject to u.vector() >= 0 ?
kkt = derivative(L, z, TestFunction(Z))

#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
solver = KrylovSolver(kkt_lhs, "tfqmr", "icc")

solver.parameters["absolute_tolerance"] = 1e-6
solver.parameters["relative_tolerance"] = 1e-4
solver.parameters["maximum_iterations"] = 1000
solver.parameters["monitor_convergence"] = True

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

null_space = VectorSpaceBasis([null_vec])
solver.set_nullspace(null_space)

null_space.orthogonalize(kkt_rhs);

#solve
solver.solve(z.vector(), kkt_rhs)    # Can we do something like 
                                     # solver.solve(z.vector(), kkt_rhs, z.vector() >= 0.0) ?

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

Comments

I've done with with the PETScSNESSolver, but not the krylov linear solver. Would the TAOLinearBoundSolver be suitable?

http://fenicsproject.org/documentation/dolfin/dev/cpp/programmers-reference/nls/TAOLinearBoundSolver.html

---

Thank you very much for the feedback.

Unfortunately I really can not compile Dolfin with TAO and use this TAOLinearBoundSolver.

( I tried several suggestions of experts from some mailing lists.
Nothing worked on my Ubuntu 12.04 of Parallels 9 virtual machine
on osx 10.9. All the others solvers of FEniCS works very well! )

Mwelland, do you think that I can try with PETScSNESSolver puting some
(weak) non-linear regularization penality togheter the linear system?

Can you show some simple example of how to put the box constrains
restrictions on PETScSNESSolver (e.g., non-negativity) to solve a non-linear
system with FEniCS project?

Thank you very much for the help.

Andre Machado

---

A quick and dirty example below. The regular solve has the entire vector <0, the bounded limits it to the trivial / initial guess of 0 everywhere.

from dolfin import *


class nlp(NonlinearProblem):
    def __init__(self, a, L, bc):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
        self.bc = bc
    def F(self, b, x):
        assemble(self.L, tensor=b)
        self.bc.apply(b)
    def J(self, A, x):
        assemble(self.a, tensor=A) 
        self.bc.apply(A)


n = 100
mesh = UnitSquareMesh(n, n)

V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
du = TrialFunction(V)
u_test = TestFunction(V)

left = CompiledSubDomain("(std::abs(x[0]) < DOLFIN_EPS) && on_boundary")

bc = DirichletBC(V, 0, left)

F = inner(grad(u), grad(u_test))*dx + 1*u_test*ds

problem = nlp(derivative(F, u, du), F, bc)

solver = PETScSNESSolver()

# You need a vector for the lower and upper bound. Here are some simple ones.
lb = interpolate(Constant(0), V)
ub = interpolate(Constant(100), V)


solver.solve(problem, u.vector()) # regular solve
#solver.solve(problem, u.vector(), lb.vector(), ub.vector())  # Bounded solve


plot(u, title='u')
interactive()

---

Dear Mwelland, the exemple of this is script is perfect to me!

Thank you very much for all the help.

Andre Machado.

( p.s. I hope that in the future the TAOLinearBoundSolver is also available
as a standart method for box constrained linear problems... )