# How to bound the solution to a given intervall

Hi,
I’m working on a PDE which I succeed to solve with Fenics but now I want to constraint a part of the solution to a given interval. In fact I am solving an équation of (u,p) where u is a 2D vector and p a density (scalar).

Before to constraint I was solving my PDE like this:

``````#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from fenics import *
from dolfin import *
from mshr import *
import numpy as np

tx=5
ty=tx
nbx=100
nby=nbx

domain=Rectangle(Point(-tx/2,-ty/2),Point(tx/2,ty/2))
mesh=RectangleMesh(Point(-tx/2,-ty/2),Point(tx/2,ty/2),nbx,nby)

V1 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
V2 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
V = VectorFunctionSpace(mesh,'P',2)
Vs= FunctionSpace(mesh,'P',2)

# Make a mixed space
TH = V1*V2
W = FunctionSpace(mesh, TH)

#Initial conditions
r=2
u_0 = Expression(('r*(rand()-0.5)', 'r*(rand()-0.5)'),degree=2,r=r)
u_alea=interpolate(u_0, W.sub(0).collapse())

p_0 = Expression('0.2', degree=0,r=r)
p_cst=interpolate(p_0, W.sub(1).collapse())

#Declare functions
w=Function(W)
w_1=Function(W)

assign(w,[u_alea,p_cst])
assign(w_1,[u_alea,p_cst])

w_t=TestFunction(W)

(u,p)=split(w)
(u_1,p_1)=split(w_1)
(v,q)=split(w_t)

#Declare constants
n = FacetNormal(mesh)

mur_g  = 'near(x[0], -2.5)'
mur_d  = 'near(x[0], 2.5)'
mur_b  = 'near(x[1], -2.5)'
mur_h  = 'near(x[1], 2.5)'
bcu_d  = DirichletBC(W.sub(0), Constant((0,1)), mur_d)
bcu_g  = DirichletBC(W.sub(0), Constant((0,-1)), mur_g)
bcu_h  = DirichletBC(W.sub(0), Constant((-1,0)), mur_h)
bcu_b  = DirichletBC(W.sub(0), Constant((1,0)), mur_b)
bcu = [bcu_d,bcu_b,bcu_g,bcu_h]

lamb=Constant(0.7) #SD
alph= Constant(1e2) #s-1
pho_c=Constant(3e-3) #SD
pho=Constant(0.1) #SD
a2=Constant(alph*(pho-pho_c)) #s-1
a4=Constant(10) #mm-2 s
D=Constant(1e-2) #mm2 s-1
D2=Constant(4e-6) #mm2 s-1
sigma=Constant(5) #mm2 s-2

dt=0.5e-2
Dt = Constant(dt)
num_step=500

#Declare the form
F = dot((u-u_1)/Dt,v)*dx\
-(alph*(p-pho_c)-a4*dot(u,u))*dot(u,v)*dx\

for i in range(0,num_step):
solve(F==0,w,bcu)
assign(w_1,w)

``````

But now I want p (the density) to remain in the interval [0,1] to be physically relevant. And I don’t want to do post solving treatment like set the negatives values to zero … Because this adds mass to my system.

The solution I am exploring is around `problem.set_bounds(wmin, wmax)`

To do so I’ve just added some code before the loop over time:

``````J = derivative(F, w, TestFunction(W))

p_lower = Constant(0.0) # lower bound
p_upper = Constant(1.0) # upper bound

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

ninfty = Function(V); ninfty.vector()[:] = -np.infty
pinfty = Function(V); pinfty.vector()[:] =  np.infty

assign(lower, [ninfty, interpolate(p_lower, Vs)])
assign(upper, [pinfty, interpolate(p_upper, Vs)])

# set bounds and solve
problem = NonlinearVariationalProblem(F, w, bcu,J=J)
problem.set_bounds(lower, upper)

solver  = NonlinearVariationalSolver(problem)
solver.solve()

``````

I get a big error, but since I never use `NonlinearVariationalProblem` I am not even sure what I am doing wrong.

If you have any hint to give me I would appreciate.

Kind regards,
Yoann

I would suggest creating your own solver. `PETScSNESSolver` allows you to provide the lower and upper bounds for your function.

So the way forward could be to create a subclass of `NonlinearProblem`, calculate the jacobian as you have above, and then use PETScSNESSolver to solve the problem like here, along with providing the `bound` vectors. Basically the array of dofs for your `lower` and `upper` functions (namely `lower.vector()` and `upper.vector()`). See also here

Edit: It does seem that you can do the same with `NonlinearVariationalProblem`, I can take a look at it sometime later today.

1 Like

This should instead be

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

You can stack the bounds for your vectors directly, without having to use interpolate

``````import numpy as np

lowerVec = np.hstack((-1.e9*np.ones(W.sub(0).dim()),np.zeros(W.sub(1).dim())))
upperVec = np.hstack((1.e9*np.ones(W.sub(0).dim()),np.ones(W.sub(1).dim())))

lower.vector()[:] = lowerVec
upper.vector()[:] = upperVec
``````

I set the lower and upper bounds of `u` to some absurdly low and high values (didn’t want to use inf), but I assume you must be having an idea apriori as to what those must be.

A cleaned up version (assuming the variational formulation is correct):

``````
from dolfin import *
import numpy as np
tx=5
ty=tx
nbx=100
nby=nbx

mesh=RectangleMesh(Point(-tx/2,-ty/2),Point(tx/2,ty/2),nbx,nby)
V1 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
V2 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
V = VectorFunctionSpace(mesh,'P',2)
Vs= FunctionSpace(mesh,'P',2)
TH = V1*V2
W = FunctionSpace(mesh, TH)

#Initial conditions
r=2
u_0 = Expression(('r*(rand()-0.5)', 'r*(rand()-0.5)'),degree=2,r=r)
u_alea=interpolate(u_0, W.sub(0).collapse())

p_0 = Expression('0.2', degree=0,r=r)
p_cst=interpolate(p_0, W.sub(1).collapse())

#Declare functions
w=Function(W)
w_1=Function(W)

assign(w,[u_alea,p_cst])
assign(w_1,[u_alea,p_cst])

w_t=TestFunction(W)

(u,p)=split(w)
(u_1,p_1)=split(w_1)
(v,q)=split(w_t)

#Declare constants
n = FacetNormal(mesh)

mur_g  = 'near(x[0], -2.5)'
mur_d  = 'near(x[0], 2.5)'
mur_b  = 'near(x[1], -2.5)'
mur_h  = 'near(x[1], 2.5)'
bcu_d  = DirichletBC(W.sub(0), Constant((0,1)), mur_d)
bcu_g  = DirichletBC(W.sub(0), Constant((0,-1)), mur_g)
bcu_h  = DirichletBC(W.sub(0), Constant((-1,0)), mur_h)
bcu_b  = DirichletBC(W.sub(0), Constant((1,0)), mur_b)
bcu = [bcu_d,bcu_b,bcu_g,bcu_h]

lamb=Constant(0.7) #SD
alph= Constant(1e2) #s-1
pho_c=Constant(3e-3) #SD
pho=Constant(0.1) #SD
a2=Constant(alph*(pho-pho_c)) #s-1
a4=Constant(10) #mm-2 s
D=Constant(1e-2) #mm2 s-1
D2=Constant(4e-6) #mm2 s-1
sigma=Constant(5) #mm2 s-2

dt=0.5e-2
Dt = Constant(dt)
num_step=500

#Declare the form
F = dot((u-u_1)/Dt,v)*dx\
-(alph*(p-pho_c)-a4*dot(u,u))*dot(u,v)*dx\

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

p_lower = Constant(0.0) # lower bound
p_upper = Constant(1.0) # upper bound

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

lowerVec = np.hstack((-1.e9*np.ones(W.sub(0).dim()),np.zeros(W.sub(1).dim())))
upperVec = np.hstack((1.e9*np.ones(W.sub(0).dim()),np.ones(W.sub(1).dim())))
lower.vector()[:] = lowerVec
upper.vector()[:] = upperVec

problem = NonlinearVariationalProblem(F, w, bcu,J=J)
problem.set_bounds(lower, upper)

solver  = NonlinearVariationalSolver(problem)
solver.parameters.update({"nonlinear_solver":"snes"})

for i in range(num_step):
solver.solve()
w_1.assign(w)
``````
2 Likes

Thank you so much this works perfectly fine!
And indeed no need to use inf. I will try to read about the PETScSNESSolver and the parameters to use.