Set krylov linear solver paramters in newton solver

variational formulation
1

Hi,

I wonder if I wanted to set the linear solver to a krylov solver such as gmres, cg, how can I set its parameters such as absolute tolerance or maximum iteration, note these parameters are not the same as those in newton solver in that they set parameters on when to stop the krylov iteration rather than when to stop the newton solver.

I mean something like

solver = NewtonSolver()
solver.parameters[“linear_solver”] = “gmres”
#solver.parameters[“linear_solver”][‘maximum_iterations’] = 1000 #error
solver.parameters[“preconditioner”] = “ilu”
solver.parameters[“maximum_iterations”] = 20

Many thanks

2

For example

solver.parameters["newton_solver"]["krylov_solver"]["absolute_tolerance"] =

See also:
https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1006.html#working-with-linear-solvers

info(parameters, verbose=True)

3

I recommend designing your own NewtonSolver and setting the options using PETSc. Here’s a quick example based on the nonlinear poisson demo.

import matplotlib.pyplot as plt
from dolfin import *


class Problem(NonlinearProblem):
    def __init__(self, J, F, bcs):
        self.bilinear_form = J
        self.linear_form = F
        self.bcs = bcs
        NonlinearProblem.__init__(self)

    def F(self, b, x):
        assemble(self.linear_form, tensor=b)
        for bc in self.bcs:
            bc.apply(b, x)

    def J(self, A, x):
        assemble(self.bilinear_form, tensor=A)
        for bc in self.bcs:
            bc.apply(A)


class CustomSolver(NewtonSolver):
    def __init__(self):
        NewtonSolver.__init__(self, mesh.mpi_comm(),
                              PETScKrylovSolver(), PETScFactory.instance())

    def solver_setup(self, A, P, problem, iteration):
        self.linear_solver().set_operator(A)

        PETScOptions.set("ksp_type", "gmres")
        PETScOptions.set("ksp_monitor")
        PETScOptions.set("pc_type", "ilu")

        self.linear_solver().set_from_options()


mesh = UnitSquareMesh(32, 32)

V = FunctionSpace(mesh, "CG", 1)
g = Constant(1.0)
bcs = [DirichletBC(V, g, "near(x[0], 1.0) and on_boundary")]
u = Function(V)
v = TestFunction(V)
f = Expression("x[0]*sin(x[1])", degree=2)
F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx
J = derivative(F, u)

problem = Problem(J, F, bcs)
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())

plt.figure()
plot(u, title="Solution")

plt.figure()
plot(grad(u), title="Solution gradient")

plt.show()
8 Likes
4

Hi @nate,

is there a way to use the Krylov Solver and its precondition feature for solving nonlinear problem?
I did the common solver setup (in the bellow), but it is not working.

solve(F == 0, sol, boundaries, J=Jac
         ,solver_parameters={"newton_solver":{"krylov_solver":'default', "relative_tolerance":tol, "absolute_tolerance":tol, "maximum_iterations":max_iter}}
         ,form_compiler_parameters={"cpp_optimize": True})

Thank you in advance.

Thomas

5

If you want to use a separate preconditioner, you need to call PETScKrylovSolver.set_operators(A, P) where P is your preconditioner matrix. I don’t think there’s an interface to that through the generic solve() syntactic sugar.

Like I’ve said in the past, solve() is a gross oversimplification which detracts from the underlying operations that need to be performed.

I strongly recommend designing your own Newton solver as above, or use a PETScSNESSolver, or use PETSc directly through petsc4py.

3 Likes
6

Hi Nate,
Is it possibel to set the PETScKrylovSolver to not throw error when the desired tolerance has not been met ? I only managed to find KSPConvergedSkip here but setting

class CustomSolver(NewtonSolver):
    def __init__(self):
        NewtonSolver.__init__(self, mesh.mpi_comm(),
                              PETScKrylovSolver(), PETScFactory.instance())

    def solver_setup(self, A, P, problem, iteration):
        self.linear_solver().set_operator(A)

        PETScOptions.set("KSPConvergedSkip") #this has no effect

I am aware that globally setting

parameters['krylov_solver']['error_on_nonconvergence'] = False

does that, but is it possible to set this while creating a custom solver ?

7

I think

CustomSolver().parameters.update({"krylov_solver" : 
                                 {"error_on_nonconvergence": False}
                                  })

should do the trick. But apparently, it still fails

8

Hi, Nate. Thanks for the explanation. I would like to add a separate preconditioner to my newton solver, however, I am not sure how to pass the P into the Newton Solver as shown in the poisson demo. Do have know how this can be done? Thank you in advance!

9

Here is a silly example which should only be considered for its syntax. The preconditioner I’ve designed here is deliberately not a good approximation of the Jacobian.

import matplotlib.pyplot as plt
from dolfin import *


class Problem(NonlinearProblem):
    def __init__(self, J, J_pc, F, bcs):
        self.bilinear_form = J
        self.preconditioner_form = J_pc
        self.linear_form = F
        self.bcs = bcs
        NonlinearProblem.__init__(self)

    def F(self, b, x):
        pass

    def J(self, A, x):
        pass

    def form(self, A, P, b, x):
        assemble(self.linear_form, tensor=b)
        assemble(self.bilinear_form, tensor=A)
        assemble(self.preconditioner_form, tensor=P)
        for bc in self.bcs:
            bc.apply(b, x)
            bc.apply(A)
            bc.apply(P)


class CustomSolver(NewtonSolver):
    def __init__(self):
        NewtonSolver.__init__(self, mesh.mpi_comm(),
                              PETScKrylovSolver(), PETScFactory.instance())

    def solver_setup(self, A, P, problem, iteration):
        self.linear_solver().set_operators(A, P)

        PETScOptions.set("ksp_type", "minres")
        PETScOptions.set("pc_type", "hypre")
        PETScOptions.set("ksp_view")

        self.linear_solver().set_from_options()


mesh = UnitSquareMesh(32, 32)

V = FunctionSpace(mesh, "CG", 1)
g = Constant(1.0)
bcs = [DirichletBC(V, g, "near(x[0], 1.0) and on_boundary")]
u = Function(V)
v = TestFunction(V)
f = Expression("x[0]*sin(x[1])", degree=2)
F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx
J = derivative(F, u)

du = TrialFunction(V)
J_pc = dot(grad(du), grad(v))*dx

problem = Problem(J, J_pc, F, bcs)
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())

plt.figure()
plot(u, title="Solution")

plt.figure()
plot(grad(u), title="Solution gradient")

plt.show()
2 Likes
10

Thank you so much. This is exactly what I need for the current problem.

11

Hi, Nate, thanks again, the customized newton solver is working well. However, I didn’t manage to find a way to make the linear solver not throw an error and stop the program in the customized newton solver, the closest option I found is the following which is not working:
image
Do you know how this can be done?
Since in the time-dependent problem that I am working on, it uses an adaptive time stepping, and the nonconvergence (due to linear solver or newton solver) is sometimes caused by a big time step, so I have to reduce the time step and continue solving the problem when nonconvergence occurs.
Thank you very much in advance!

12

Looks like old DOLFIN intercepts a PETSc error and then throws its own should the KSP iterations exceed their limit. Quick solution is to handle it with python:

try:
    custom_solver.solve(problem, u.vector())
except Exception as e:
    # Handle the error within your time stepping algorithm gracefully...

Consider using a PETScSNESSolver and its underlying SNES and KSP objects to get more information about the lack of convergence. Also consider foregoing DOLFIN’s implementations of Newton and Krylov solvers altogether and just use a SNES and KSP directly.

13

Hi @nate ,

for each newton iteration, we need to solve the linear system until it converges. I want to record this iteration number for each newton iteration. I can find it on terminal. But is there any way to get these numbers and output to a txt file?

Thank you in advance.