Hello everyone,
I want to use the conjugate gradient to solve my nonlinear variational problem.
However, although I have found the scripts for linear problems using the conjugate gradient, I couldn’t find scripts for nonlinear solvers using the conjugate gradient.
Can you please help me with that?
Thank you so much!
See this for example.
Also consider an introduction to iterative methods like Saad 2003 to see why it’s a risky decision to use CG for a nonlinear problem.
Thanks for the example you shared. I think this example utilizes Newton solver, but I am looking for the conjugate gradient method for solving my nonlinear problem.
The reason that I am looking for the conjugate gradient method, and not Newton method, is that my code fails with Newton solver, and gives the following error:
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
Any help about how to use conjugate gradient and/or how to solve the current error with Newton solver is much appreciated!
It takes lots of care to get a “highly nonlinear” problem to converge. If you can you reproduce your non-converging problem with a minimal working example (\backsim 30 lines of code or fewer), we can take a look.
If not, the following tips and tricks may help: Default absolute tolerance and relative tolerance
Also, try reducing the complexity of your problem until the Newton solver converges successfully so you can find the problematic portion of your formulation.
Thanks, Nate for the tips you shared with me.
Unfortunately, I couldn’t reproduce my error with a minimal example, but I know that the problem occurs when I increase the number of meshes. The problem with the lower number of meshes is that it does not capture some small sensitive parts of the domain that I am interested in.
And when I use adaptive meshing, even for very small tolerance (1e-200), I get the message " Error estimate (0) is less than tolerance (1e-200), returning. ", and consequently, the code does not generate adaptive meshes.
Also, I have another question. I think for nonlinear problems, your suggestion is to stick with Newton solver, and not other methods (e.g., conjugate gradient, steepest descent). Is there any other method other than the Newon solver that you would recommend for nonlinear problems?
Thank you so much for your very helpful answers! 
For scalable nonlinear solvers, you could take a look at this paper as a starting point. In general it can pin down to be highly specific to the problem you are solving.
Newton combined with a proper choice of linear solver (at each Newton-iteration) and a preconditioner should usually give you out of the box good performance for many. PETScSNES is also a good place to start looking out for solvers available directly through PETSc. You could also take a look at Non linear solver to program the solver by yourself.
https://github.com/FEniCS/dolfinx/blob/master/python/test/unit/nls/test_newton.py and a corresponding dolfin-equivalent Set krylov linear solver paramters in newton solver are also good pointers on how to setup your problem/solver and preconditioner.
My personal experience: Newton with a line search and a good preconditioner (e.g., AMG for hyperelasticity) should be pretty robust and scale well for larger problems. You can use CG as the Krylov solver if your system is SPD. Direct solvers (even MUMPS, SUPERLU_DIST) will run out of memory for very large problems, and therefore unsuitable. For many problems, however, Newton is not a good choice and you have to go with a gradient-descent approach which has the disadvantage of being horribly slow.
Yet another alternative for Newton: http://www.lmm.jussieu.fr/~corrado/articles/cm-a16-FarMau.pdf, for variational phase-field, which is also a highly nonlinear and non-convex problem.
Thanks, Bhavesh for your very comprehensive response!
I’ll try to implement your very helpful suggestions on my code.
All the best!
Hi,
Thank you so much for following up on my question @nate @bhaveshshrimali
I am using SNES to solve my nonlinear variational problem, but I am getting the following error:
*** Error: Unable to solve nonlinear variational problem.
*** Reason: Unknown nonlinear solver type.
*** Where: This error was encountered inside NonlinearVariationalSolver.cpp.
Following is my minimal code:
import numpy as np
from numpy import linalg as LA
from numpy.linalg import inv
import matplotlib.pyplot as plt
from dolfin import *
from scipy import optimize
import numdifftools as nd
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
# Create mesh and define function space
mesh = RectangleMesh(Point(-1,-1),Point(1,1), 50, 50)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
################################################################################
# Define boundary condictions
################################################################################
F11=0.5
F12=0.0 + 0.1
F21=0.0 - 0.1
F22=2.0
y_BC = Expression(("F11*x[0]+F12*x[1]", "F21*x[0]+F22*x[1]"), F11=F11,F12=F12,F21=F21,F22=F22,degree=2)
tol = 1E-3
def x0_boundary(x, on_boundary):
return on_boundary and near(x[0], -1, tol) or near(x[0], 1, tol) or \
near(x[1], -1, tol) or near(x[1], 1, tol)
bcs = DirichletBC(V, y_BC, x0_boundary)
################################################################################
# Define functions
################################################################################
dy = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
y = Function(V)
B = Constant((0, 0)) # Body force per unit area
T = Constant((0, 0)) # Traction force on the boundary
y_reference = Expression(("x[0]","x[1]"),degree=3)
################################################################################
# Elasticity model: Stored strain energy density (compressible neo-Hookean model)
################################################################################
E, nu = 10.0, 0.3
mu, blk = Constant(E/(2*(1 + nu))), Constant(E/(3*(1 - 2*nu)))
def W(y): #Stored strain energy density
F = grad(y) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
Ic = tr(C)
J = det(F)
return (mu/2)*(Ic - 3) - mu*ln(J) + (blk/2 - mu/3)*(J-1)**2
################################################################################
# Total potential energy and solve
################################################################################
Pi = W(y)*dx -dot(B, (y-y_reference))*dx + dot(T, (y-y_reference))*ds
F = derivative(Pi, y, v) # first variation (derivative at u in the direction of v)
dF = derivative(F, y, dy) # Jacobian of F
problem=NonlinearVariationalProblem(F, y, bcs=bcs, J=dF)
solver=NonlinearVariationalSolver(problem)
solver.parameters['nonlinear_solver'] = 'snesncg'
solver.parameters["snes_solver"]["maximum_iterations"] = 50
solver.parameters["snes_solver"]["report"] = False
solver.solve()
To solve the problem, I tried implementing what Nate has suggested in this link: Using petsc4py.PETSc.SNES directly - #11 by niewiarowski
However, I am getting another error that says: TypeError: 'DirichletBC' object is not iterable
The code that I have used based on the mentioned link is as below:
import numpy as np
from numpy import linalg as LA
from numpy.linalg import inv
import matplotlib.pyplot as plt
from dolfin import *
from scipy import optimize
import numdifftools as nd
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
# Create mesh and define function space
mesh = RectangleMesh(Point(-1,-1),Point(1,1), 50, 50)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
################################################################################
# Define boundary condictions
################################################################################
F11=0.5
F12=0.0 + 0.1
F21=0.0 - 0.1
F22=2.0
y_BC = Expression(("F11*x[0]+F12*x[1]", "F21*x[0]+F22*x[1]"), F11=F11,F12=F12,F21=F21,F22=F22,degree=2)
tol = 1E-3
def x0_boundary(x, on_boundary):
return on_boundary and near(x[0], -1, tol) or near(x[0], 1, tol) or \
near(x[1], -1, tol) or near(x[1], 1, tol)
bcs = DirichletBC(V, y_BC, x0_boundary)
################################################################################
# Define functions
################################################################################
dy = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
y = Function(V)
B = Constant((0, 0)) # Body force per unit area
T = Constant((0, 0)) # Traction force on the boundary
y_reference = Expression(("x[0]","x[1]"),degree=3)
################################################################################
# Elasticity model: Stored strain energy density (compressible neo-Hookean model)
################################################################################
E, nu = 10.0, 0.3
mu, blk = Constant(E/(2*(1 + nu))), Constant(E/(3*(1 - 2*nu)))
def W(y): #Stored strain energy density
F = grad(y) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
Ic = tr(C)
J = det(F)
return (mu/2)*(Ic - 3) - mu*ln(J) + (blk/2 - mu/3)*(J-1)**2
################################################################################
# Total potential energy and solve
################################################################################
Pi = W(y)*dx -dot(B, (y-y_reference))*dx + dot(T, (y-y_reference))*ds
F = derivative(Pi, y, v) # first variation (derivative at u in the direction of v)
class SNESProblem():
def __init__(self, F, y, bcs):
V = y.function_space()
dy = TrialFunction(V)
self.L = F
self.a = derivative(F, y, dy)
self.bcs = bcs
self.y = y
def F(self, snes, x, F):
x = PETScVector(x)
F = PETScVector(F)
assemble(self.L, tensor=F)
for bc in self.bcs:
bc.apply(F, x)
def J(self, snes, x, J, P):
J = PETScMatrix(J)
assemble(self.a, tensor=J)
for bc in self.bcs:
bc.apply(J)
problem = SNESProblem(F, y, bcs)
import petsc4py, sys
petsc4py.init(sys.argv)
from petsc4py import PETSc
b = PETScVector() # same as b = PETSc.Vec()
J_mat = PETScMatrix()
snes = PETSc.SNES().create(MPI.comm_world)
snes.setFunction(problem.F, b.vec())
snes.setJacobian(problem.J, J_mat.mat())
snes.solve(None, problem.y.vector().vec())
I would be so thankful if you could help me to solve my nonlinear variational problem using SNES (and preferably, its conjugate gradient method (NCG)).
Thank you very much!
Mary
You can change the solver to snes by the following:
solver.parameters["nonlinear_solver"]="snes"
solver.parameters["snes_solver"]["method"] = "ncg"
As mentioned in Solving adaptive nonlinear equation in Vector Space - #11 by dokken you can inspect the parameters object to determine the keys.
For your second problem, you simply need to wrap bcs as a list, i.e. change:
to bcs= [ DirichletBC(V, y_BC, x0_boundary)]
You could also try working with PETScSNESSolver or, for even finer control, directly use petsc4py.