Lower and upper bounds of a solution

I am trying to solve the problem \nabla u = f where the solution is u has lower and upper bounds i.e 0 \leq u \leq 1 and \int u\;dx=c_0. I follow steps and example posted in fenics demos. So I used SNES solver for upper and lower bounds and Lengrange Function Space for integral constraint. Integral constraint is working fine, I am having problems with the lower and upper bounds. Any help will be much appreciated

import numpy
from dolfin import *

mesh = UnitSquareMesh(64, 64)

V = FiniteElement("CG", mesh.ufl_cell(), 1)
R = FiniteElement("R", mesh.ufl_cell(), 0)
W_elem = MixedElement([V, R])
W = FunctionSpace(mesh, W_elem)

c0 = Constant(5)

u_lower = Constant(0.0) # lower bound 
u_upper = Constant(1.0) # upper bound 

# Define variational problem
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("-sin(5*x[0])", degree=2)
a = inner(grad(u), grad(v)) * dx + c * v * dx + u * d * dx
L = f * v * dx + c0 * d * dx
F = a - L

# Compute solution
w = Function(W)
bc = []

# J = derivative(F, w, TrialFunction(W))
J = derivative(L, w, TestFunction(W))
H = derivative(L, w, TrialFunction(W))


snes_solver_parameters = {"nonlinear_solver": "snes",
                          "snes_solver": {"linear_solver": "lu",
                                          "maximum_iterations": 20,
                                          "report": True,
                                          "error_on_nonconvergence": False}}

problem = NonlinearVariationalProblem(J, w, bc, H)

# handle the bounds
lower = Function(W)
upper = Function(W)

ninfty = Function(W.sub(0).collapse()); ninfty.vector()[:] = -numpy.infty
pinfty = Function(W.sub(0).collapse()); pinfty.vector()[:] =  numpy.infty

fa = FunctionAssigner(W, [V])
fa.assign(lower, [interpolate(y_lower, V)])
fa.assign(upper, [interpolate(y_upper, V)])

# set bounds and solve
problem.set_bounds(lower, upper)

solver  = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
info(solver.parameters, True)

(iter, converged) = solver.solve()


# w = Function(W)
# solve(a == L, w)
(u, c) = w.split()

plot(u, title='Final solution')
filename = 'demo_u.png'
plt.savefig(filename)
plt.close()

# Save solution in VTK format
file = File("output/demo_u.pvd")
file << u```

There are several issues with your code. I’ve addressed all of them in the following snippet.
Next time, please post the error message you are receiving and reduce the minimal code example to only contain the lines that are required to produce the error:

import matplotlib.pyplot as plt
import numpy
from dolfin import *

mesh = UnitSquareMesh(64, 64)

V = FiniteElement("CG", mesh.ufl_cell(), 1)
R = FiniteElement("R", mesh.ufl_cell(), 0)
W_elem = MixedElement([V, R])
W = FunctionSpace(mesh, W_elem)
W0 = W.sub(0).collapse()
W1 = W.sub(1).collapse()
c0 = Constant(5)

u_lower = Constant(0.0) # lower bound 
u_upper = Constant(1.0) # upper bound 

# Define variational problem
w = Function(W)
(u, c) = split(w)
(v, d) = TestFunctions(W)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("-sin(5*x[0])", degree=2)
a = inner(grad(u), grad(v)) * dx + c * v * dx + u * d * dx
L = f * v * dx + c0 * d * dx
F = a - L

# Compute solution
bc = []

J = derivative(F, w, TrialFunction(W))

snes_solver_parameters = {"nonlinear_solver": "snes",
                          "snes_solver": {"linear_solver": "lu",
                                          # "maximum_iterations": 20,
                                          # "report": True,
                                          # "error_on_nonconvergence": False
                          }}

problem = NonlinearVariationalProblem(F, w, bc, J)

# handle the bounds
lower = Function(W)
upper = Function(W)

ninfty = Function(W1);
ninfty.vector()[:] = -1e1
pinfty = Function(W1);
pinfty.vector()[:] =  1e1
u_lower = Function(W0)
u_lower.vector()[:] = 0
u_upper = Function(W0)
u_upper.vector()[:] = 1

fa = FunctionAssigner(W, [W0, W1])
fa.assign(lower, [u_lower, ninfty])
fa.assign(upper, [u_upper, pinfty])

# set bounds and solve
problem.set_bounds(lower, upper)

solver  = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
# info(solver.parameters, True)
(iter, converged) = solver.solve()


# w = Function(W)
# solve(a == L, w)
(u, c) = w.split()

plot(u, title='Final solution')
filename = 'demo_u.png'
plt.savefig(filename)
plt.close()

# Save solution in VTK format
file = File("output/demo_u.pvd")
file << u

Thank you. Next time I will be posting the error message as well.

Hi @dokken : I have a question in this reply of yours. When u is a vector-valued function, how do we set the -ve infinity and +ve infinity bounds for each component?

Thanks!

I would suggest interpolating the upper and lower bounds into each respective function.