Hello,
Is there a way to set bounds on the solution of a NonLinearProblem instance in dolfinx, that is an equivalent of the set_bounds option in dolfin. The syntax was as follows.
problem = NonLinearVariationalProblem(F,u,bc)
problem.set_bounds(u_min, u_max)
Thank you
Krishna
Use the PETSc SNES interface directely, as shown in:
using
Hello,
Thank you for the reply. I tried to run test_nonlinear_pde_snes() in the test_newton.py file after adding the following few lines for bounds.
u_min = u.copy()
u_min.x.array[:] = 0.5
u_max = u.copy()
u_max.x.array[:] = 1.0
snes.setVariableBounds(u_min.vector, u_max.vector)
But I get the error that newtonls does not support bounds.
File “PETSc/SNES.pyx”, line 573, in petsc4py.PETSc.SNES.solve
petsc4py.PETSc.Error: error code 73
[0] SNESSolve() at /usr/local/petsc/src/snes/interface/snes.c:4809
[0] SNESSolve_NEWTONLS() at /usr/local/petsc/src/snes/impls/ls/ls.c:145
[0] Object is in wrong state
[0] SNES solver newtonls does not support bounds
How do I fix this?
Thank you
Read up on SNESVISetVariableBounds particularly with regard to the SNESType and Section 2.4.2 of the PETSc Manual.
Thank you for the links. I was able to figure out from the PETSc manual a list of non-linear solvers that can work with constraints.
Sorry to continue a thread already marked as solved, but I am not sure I have figured out the right way of doing this and need some more help.
Here is a MWE adapted from the unit test for nonlinear pde snes.
# Copyright (C) 2018 Garth N. Wells
#
# This file is part of DOLFINx (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
"""Unit tests for Newton solver assembly"""
import numpy as np
import ufl
from dolfinx import cpp as _cpp
from dolfinx import la
from dolfinx.fem import (Function, FunctionSpace, dirichletbc, form,
locate_dofs_geometrical)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_matrix, create_vector, set_bc)
from dolfinx.mesh import create_unit_square
from ufl import TestFunction, TrialFunction, derivative, dx, grad, inner
from mpi4py import MPI
from petsc4py import PETSc
class NonlinearPDE_SNESProblem:
def __init__(self, F, u, bc):
V = u.function_space
du = TrialFunction(V)
self.L = form(F)
self.a = form(derivative(F, u, du))
self.bc = bc
self._F, self._J = None, None
self.u = u
def F(self, snes, x, F):
"""Assemble residual vector."""
x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
x.copy(self.u.vector)
self.u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
with F.localForm() as f_local:
f_local.set(0.0)
assemble_vector(F, self.L)
apply_lifting(F, [self.a], bcs=[[self.bc]], x0=[x], scale=-1.0)
F.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(F, [self.bc], x, -1.0)
def J(self, snes, x, J, P):
"""Assemble Jacobian matrix."""
J.zeroEntries()
assemble_matrix(J, self.a, bcs=[self.bc])
J.assemble()
def test_nonlinear_pde_snes(bound):
"""Test Newton solver for a simple nonlinear PDE"""
# Create mesh and function space
mesh = create_unit_square(MPI.COMM_WORLD, 12, 15)
V = FunctionSpace(mesh, ("Lagrange", 1))
u = Function(V)
u.x.array[:] = 15
u.name = 'solution'
from dolfinx import io
with io.XDMFFile(MPI.COMM_WORLD,'u.xdmf','w') as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(u,0)
v = TestFunction(V)
F = inner(5.0, v) * dx - ufl.sqrt(u * u) * inner(
grad(u), grad(v)) * dx - inner(u, v) * dx
u_bc = Function(V)
u_bc.x.array[:] = 1.0
bc = dirichletbc(u_bc, locate_dofs_geometrical(V, lambda x: np.logical_or(np.isclose(x[0], 0.0),
np.isclose(x[0], 1.0))))
# Create nonlinear problem
problem = NonlinearPDE_SNESProblem(F, u, bc)
b = la.create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs)
J = create_matrix(problem.a)
# Create Newton solver and solve
snes = PETSc.SNES().create()
opts = PETSc.Options()
opts['snes_monitor'] = None
snes.setFromOptions()
snes.setType('vinewtonssls')
snes.setFunction(problem.F, b)
snes.setJacobian(problem.J, J)
snes.setTolerances(atol=1e-10,rtol=1.0e-9, max_it=10)
snes.getKSP().setType("preonly")
snes.getKSP().setTolerances(rtol=1.0e-9)
snes.getKSP().getPC().setType("lu")
if bound:
u_min, u_max = Function(V), Function(V)
u_min.x.array[:] = 1.0
u_max.x.array[:] = 1.5
snes.setVariableBounds(u_min.vector, u_max.vector)
snes.solve(None, u.vector)
print('Minimum in u array:',np.min(u.x.array),'\nMaximum in u array:', np.max(u.x.array))
from dolfinx import io
with io.XDMFFile(MPI.COMM_WORLD,'u.xdmf','a') as xdmf:
xdmf.write_function(u,1)
print('No bounds')
test_nonlinear_pde_snes(bound=False)
print('\nBounds set between 1.0 and 1.5')
test_nonlinear_pde_snes(bound=True)
When I run this code, I get the following results,
No bounds
0 SNES Function norm 9.083961017374e+02
1 SNES Function norm 2.603033612549e-01
2 SNES Function norm 6.870718004914e-02
3 SNES Function norm 1.678228243832e-03
4 SNES Function norm 8.126048345290e-07
Minimum in u array: 1.0
Maximum in u array: 1.3851364823572738Bounds set between 1.0 and 1.5
0 SNES Function norm 2.202572332184e+00
1 SNES Function norm 4.139002904298e-01
Minimum in u array: -0.012085010497842469
Maximum in u array: 0.03760750671222456
Considering that the final solution without bounds lies between 1.0 and 1.5, I expected the solver to converge to that solution without any issues even with those bounds set. But this doesn’t seem to be the case. Is there something wrong with the way I am setting the bounds on the problem?
Also, the two solvers that PETSc provides for problems with bounds, vinewtonssls and vinewtonrsls give different results.
Could someone please share their expertise in this matter?
Thank you
Krishna Kenja
I am no expert on this either but variational inequalities in general are tricky and not as trivial as setting the bounds to what you are getting without any constraints.
In particular, you may want to take a look at projected newton methods, e.g.
For an immediate remedy, consider replacing the solver to vinewtonrsls and the bounds to be something for reasonable:
if bound:
u_min, u_max = Function(V), Function(V)
u_min.vector.setArray(1.0)
u_max.vector.setArray(20.5) # try changing this and you will see that the solver doesn't converge
snes.setVariableBounds(u_min.vector, u_max.vector)
which does give something more reasonable:
No bounds
0 SNES Function norm 9.083961017374e+02
1 SNES Function norm 2.603033612549e-01
2 SNES Function norm 6.870718004914e-02
3 SNES Function norm 1.678228243832e-03
4 SNES Function norm 8.126048345290e-07
Minimum in u array: 1.0
Maximum in u array: 1.3851364823572738
Bounds set between 1.0 and 1.5
0 SNES Function norm 9.083961017374e+02
1 SNES Function norm 2.603033612549e-01
2 SNES Function norm 6.870718004914e-02
3 SNES Function norm 1.678228243832e-03
4 SNES Function norm 8.126048345290e-07
Minimum in u array: 1.0
Maximum in u array: 1.3851364823572738
Thank you @bhaveshshrimali for the links to the papers. I will sure to check them out. I suspected that it might be heavily problem dependent, but posted the MWE just in case I was missing something obvious/fundamental about it.
For anyone else that might land on this page, I just want to leave it here that I was able to solve the issue I was facing with my code by using a basic line search type instead of the default bt for the VI solvers.
Thank you