How to tell Fenics that the problem is linear?

Question

Hi everybody,

i'm kind of stuck with the following problem. I consider the simple unconstrained minimization problem (there are other terms, but they are linear in u, so they dont play a role)

min 0.5 \|y-yd\|^2 + 0.5*alpha*\|u\|^2
s.t. \Laplace y = u in Omega with Dirichlet boundary constraints

Since this problem is unconstrained I want to consider the Lagrange Function

F(u,y,p) = 0.5 \|y-yd\|^2 + 0.5*alpha*\|u\|^2 + (-\Delta y - u , p)

compute its derivate and solve the resulting linear system. I use C++, so here is the resulting ufl form:

state_element = FiniteElement("Lagrange", triangle, 1)


# Test function space
element = MixedElement([state_element, state_element, state_element]) 


z = Coefficient(element)
[u, y, p] = split( z )

# additional variables
yd          = Coefficient(state_element)
mult_z      = Coefficient(state_element)
mult_lmbda  = Coefficient(state_element)
lmbda       = Coefficient(state_element)

alpha = Constant(triangle)
sigma = Constant(triangle)


# Define Lagrange function
L = 0.5*(y-yd)**2*dx + 0.5*alpha*u**2*dx  - alpha*lmbda*u*dx + mult_lmbda*(u-mult_z)*dx + 0.5*sigma*(u-mult_z)**2*dx  +   inner( grad(y),grad(p) )*dx - u*p*dx 

a = derivative (L, z, TestFunction(element))

Up to here everything works fine. Now I want to solve this linear optimality system. By construction Fenics handles the form a as nonlinear, so i cannot use the simple solve command. Of course I can compute the second derivative and use a Nonlinear-Variational Solver, but this is kind of over the top.

Does anybody knows a nice trick how to solve a==0 in a simple way? I want to treat a as a matrix - or at least as a linear operator. The solve command accepts GenericLinearOperators and its subclasses. But i dont now how to convert a into this form.

Thanks in advance!

SM · Answer 1 · Jan 27, 2017

If your form is truly linear and FEniCS is just mistaking it as nonlinear, you can try rewriting the equation e.g.:

a_temp = a
a = lhs(a_temp)
L = rhs(a_temp)
eqn = Equation(a, L)

in your UFL form. Then, in your C++ code, you can use:

solve(eqn, z)

This might work. Some references (though to older documentation):

How to tell Fenics that the problem is linear?

1 Answer