How to bound the solution to a given intervall

Yoann · July 16, 2020, 2:47pm

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\
    +lamb*dot(dot(u, nabla_grad(u)), v)*dx\
    +D*(dot(grad(u[0]), grad(v[0]))+dot(grad(u[1]), grad(v[1])))*dx\
    -(alph*(p-pho_c)-a4*dot(u,u))*dot(u,v)*dx\
    +sigma*dot(grad(p),v)*dx\
    +(((p-p_1)/Dt+div(p*u))*q+D2*dot(grad(p), grad(q)))*dx\
    -D2*dot(grad(p),n)*q*ds
    

    
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

bhaveshshrimali · July 16, 2020, 5:10pm

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.

bhaveshshrimali · July 16, 2020, 5:57pm

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\
    +lamb*dot(dot(u, nabla_grad(u)), v)*dx\
    +D*(dot(grad(u[0]), grad(v[0]))+dot(grad(u[1]), grad(v[1])))*dx\
    -(alph*(p-pho_c)-a4*dot(u,u))*dot(u,v)*dx\
    +sigma*dot(grad(p),v)*dx\
    +(((p-p_1)/Dt+div(p*u))*q+D2*dot(grad(p), grad(q)))*dx\
    -D2*dot(grad(p),n)*q*ds

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)

Yoann · July 16, 2020, 8:06pm

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.

Topic		Replies	Views
Penalty formulation in Fenics variational formulation	2	888	June 30, 2021
Integral constraint variational formulation	4	455	November 16, 2020
Integral Constrains in FEniCSx General	5	947	June 16, 2023
Time out of the Solver variational formulation	1	651	March 18, 2020
Doubt about Fenics solution General	5	105	April 28, 2024

How to bound the solution to a given intervall

Related topics