DEV: Add the possibility to specifiy a minimum number of iterations for inner (linear) solvers in Newton's method

CHANGELOG.rst

@@ -26,8 +26,14 @@ of the maintenance releases, please take a look at
 
 * :py:func:`cashocs.import_mesh` can now also directly import a Gmsh mesh file. Internally, the mesh is directly converted to xdmf and then read. At the moment, this only supports the conversion mode `physical`.
 
+* Add a new parameter :python:`min_inner_iter` to :py:func:`cashocs.newton_solve`. This parameter is used to specify the minimum amount of iterations that the inner (i.e. linear) solver should perform, regardless of termination criteria.
+
 * New configuration file parameters:
 
+  * Section StateSystem
+
+    * :ini:`min_inner_iter` is a parameter controlling the minimum number of iterations that a inner (i.e. linear) solver should perform when using Newton's method. Mainly intended to be used for inexact Newton's method.
+
   * Section LineSearch
 
     * :ini:`fail_if_not_converged` determines, whether the line search is cancelled once the state system cannot be solved or if a new iterate is tried instead.

cashocs/_pde_problems/state_problem.py

@@ -66,6 +66,7 @@ def __init__(
         self.newton_atol = self.config.getfloat("StateSystem", "newton_atol")
         self.newton_damped = self.config.getboolean("StateSystem", "newton_damped")
         self.newton_inexact = self.config.getboolean("StateSystem", "newton_inexact")
+        self.min_inner_iter = self.config.getint("StateSystem", "min_inner_iter")
         self.newton_verbose = self.config.getboolean("StateSystem", "newton_verbose")
         self.newton_iter = self.config.getint("StateSystem", "newton_iter")
 
@@ -149,6 +150,7 @@ def solve(self) -> List[fenics.Function]:
                             A_tensor=self.A_tensors[i],
                             b_tensor=self.b_tensors[i],
                             preconditioner_form=self.db.form_db.preconditioner_forms[i],
+                            min_inner_iter=self.min_inner_iter,
                         )
 
             else:

cashocs/_utils/linalg.py

@@ -24,10 +24,13 @@
 
 import fenics
 import numpy as np
+from packaging import version
+import petsc4py
 from petsc4py import PETSc
 import ufl
 
 from cashocs import _exceptions
+from cashocs import _loggers
 from cashocs._utils import forms as forms_module
 
 if TYPE_CHECKING:
@@ -323,6 +326,7 @@
     atol: Optional[float] = None,
     comm: Optional[MPI.Comm] = None,
     P: Optional[PETSc.Mat] = None,  # pylint: disable=invalid-name
+    min_iter: int = 0,
 ) -> PETSc.Vec:
     """Solves a finite dimensional linear problem.
 
@@ -344,6 +348,8 @@
             solving the linear problem. Overrides the specification in the ksp object
             and ksp_options.
         comm: The MPI communicator for the problem.
+        min_iter: The minimum number of iterations to use for the linear solver,
+            regardless of tolerances.
 
     Returns:
         The solution vector.
@@ -365,6 +371,14 @@
     options = define_ksp_options(ksp_options)
 
     setup_fieldsplit_preconditioner(fun, ksp, options)
+    if min_iter > 0:
+        if version.parse(petsc4py.__version__) < version.parse("3.20"):
+            _loggers.warning(
+                "The ability to specify a minimum number of iterations for a linear "
+                "solver is only available with PETSc >= 3.20. The specified minimum "
+                "iteration number will be ignored."
+            )
+        options["ksp_min_it"] = min_iter
     setup_petsc_options([ksp], [options])
 
     if rtol is not None:

cashocs/io/config.py

@@ -129,6 +129,9 @@ def __init__(self, config_file: Optional[str] = None) -> None:
                 "newton_inexact": {
                     "type": "bool",
                 },
+                "min_inner_iter": {
+                    "type": "int",
+                },
                 "newton_verbose": {
                     "type": "bool",
                 },
@@ -542,6 +545,7 @@ def __init__(self, config_file: Optional[str] = None) -> None:
 newton_iter = 50
 newton_damped = False
 newton_inexact = False
+min_inner_iter = 0
 newton_verbose = False
 picard_iteration = False
 picard_rtol = 1e-10

cashocs/nonlinear_solvers/newton_solver.py

@@ -57,6 +57,7 @@ def __init__(
         b_tensor: Optional[fenics.PETScVector] = None,
         is_linear: bool = False,
         preconditioner_form: ufl.Form = None,
+        min_inner_iter: int = 0,
     ) -> None:
         r"""Initializes self.
 
@@ -95,6 +96,10 @@ def __init__(
             is_linear: A boolean flag, which indicates whether the problem is actually
                 linear.
             preconditioner_form: A UFL form which defines the preconditioner matrix.
+            min_inner_iter: Minimum number of iterations for the inner (linear) solver
+                to perform. This parameter is only intended for being used with an
+                inexact Newton's method (there, in the first steps, only a small handful
+                of inner iterations are needed to satisfy the termination criterion).
 
         """
         self.nonlinear_form = nonlinear_form
@@ -116,6 +121,8 @@ def __init__(
         else:
             self.preconditioner_form = preconditioner_form
 
+        self.min_inner_iter = min_inner_iter
+
         self.verbose = verbose if not self.is_linear else False
 
         temp_derivative = derivative or fenics.derivative(self.nonlinear_form, self.u)
@@ -293,6 +300,7 @@ def solve(self) -> fenics.Function:
                 atol=self.atol / 10.0,
                 comm=self.comm,
                 P=self.P_matrix,
+                min_iter=self.min_inner_iter,
             )
 
             if self.is_linear:
@@ -440,6 +448,7 @@ def newton_solve(
     b_tensor: Optional[fenics.PETScVector] = None,
     is_linear: bool = False,
     preconditioner_form: ufl.Form = None,
+    min_inner_iter: int = 0,
 ) -> fenics.Function:
     r"""Solves a nonlinear problem with Newton\'s method.
 
@@ -475,6 +484,11 @@ def newton_solve(
             sub-problem.
         is_linear: A boolean flag, which indicates whether the problem is actually
             linear.
+        preconditioner_form: A UFL form which defines the preconditioner matrix.
+        min_inner_iter: Minimum number of iterations for the inner (linear) solver
+            to perform. This parameter is only intended for being used with an
+            inexact Newton's method (there, in the first steps, only a small handful
+            of inner iterations are needed to satisfy the termination criterion).
 
     Returns:
         The solution of the nonlinear variational problem, if converged. This overwrites
@@ -523,6 +537,7 @@ def newton_solve(
         b_tensor=b_tensor,
         is_linear=is_linear,
         preconditioner_form=preconditioner_form,
+        min_inner_iter=min_inner_iter,
     )
 
     solution = solver.solve()

cashocs/nonlinear_solvers/picard_solver.py

@@ -94,6 +94,7 @@ def picard_iteration(
     b_tensors: Optional[List[fenics.PETScVector]] = None,
     inner_is_linear: bool = False,
     preconditioner_forms: Optional[Union[List[ufl.Form], ufl.Form]] = None,
+    min_inner_iter: int = 0,
 ) -> None:
     """Solves a system of coupled PDEs via a Picard iteration.
 
@@ -122,8 +123,10 @@ def picard_iteration(
         inner_is_linear: Boolean flag, if this is ``True``, all problems are actually
             linear ones, and only a linear solver is used.
         preconditioner_forms: The list of forms for the preconditioner. The default
-                is `None`, so that the preconditioner matrix is the same as the system
-                matrix.
+            is `None`, so that the preconditioner matrix is the same as the system
+            matrix.
+        min_inner_iter: Minimum iterations the inner (linear) solver should perform
+            regardless of termination criteria.
 
     """
     is_printing = verbose and fenics.MPI.rank(fenics.MPI.comm_world) == 0
@@ -188,6 +191,7 @@ def picard_iteration(
                 b_tensor=b_tensor,
                 is_linear=inner_is_linear,
                 preconditioner_form=preconditioner_form_list[j],
+                min_inner_iter=min_inner_iter,
             )
 
     if is_printing:

docs/source/user/demos/optimal_control/doc_config.rst

@@ -110,6 +110,14 @@ Additionally, we have the boolean parameter :ini:`newton_inexact`, defined via
 
 which sets up an inexact Newton method for solving nonlinear problems in case this is :ini:`newton_inexact = True`. The default is :ini:`newton_inexact = False`.
 
+For the inexact Newton's method, it can be beneficial to set a minimum number of iterations that have to be performed regardless of termination criteria. This is sensible, as the tolerance for the early iterations is usually very low, meaning that the linear solver converges in a handful of iterations. This leads to unwanted behavior, as the Jacobian is re-assembled very often initially. The parameter
+
+.. code-block:: ini
+
+    min_inner_iter = 0
+
+can be used to set the minimum amount of iterations that a linear solver should perform (when used in Newton's method). If this is :ini:`0` (which is the default), there is no minimum number of iterations prescribed.
+
 The parameter 
 
 .. code-block:: ini
@@ -633,6 +641,8 @@ in the following.
       - if :ini:`newton_damped = True`, damping is enabled
     * - :ini:`newton_inexact = False`
       - if :ini:`newton_inexact = True`, an inexact Newton's method is used
+    * - :ini:`min_inner_iter = 0`
+      - Specifies the minimum amount of iterations the inner (linear) solver has to perform in a Newton method
     * - :ini:`newton_verbose = False`
       - :ini:`newton_verbose = True` enables verbose output of Newton's method
     * - :ini:`picard_iteration = False`

docs/source/user/demos/shape_optimization/doc_config.rst

@@ -161,6 +161,14 @@ Additionally, we have the boolean parameter :ini:`newton_inexact`, defined via
 
 which sets up an inexact Newton method for solving nonlinear problems in case this is :ini:`newton_inexact = True`. The default is :ini:`newton_inexact = False`.
 
+For the inexact Newton's method, it can be beneficial to set a minimum number of iterations that have to be performed regardless of termination criteria. This is sensible, as the tolerance for the early iterations is usually very low, meaning that the linear solver converges in a handful of iterations. This leads to unwanted behavior, as the Jacobian is re-assembled very often initially. The parameter
+
+.. code-block:: ini
+
+    min_inner_iter = 0
+
+can be used to set the minimum amount of iterations that a linear solver should perform (when used in Newton's method). If this is :ini:`0` (which is the default), there is no minimum number of iterations prescribed.
+
 Next, we have the parameter
 
 .. code-block:: ini
@@ -1129,6 +1137,8 @@ in the following.
       - if :ini:`newton_damped = True`, damping is enabled
     * - :ini:`newton_inexact = False`
       - if :ini:`newton_inexact = True`, an inexact Newton's method is used
+    * - :ini:`min_inner_iter = 0`
+      - Specifies the minimum amount of iterations the inner (linear) solver has to perform in a Newton method
     * - :ini:`newton_verbose = False`
       - :ini:`newton_verbose = True` enables verbose output of Newton's method
     * - :ini:`picard_iteration = False`

pyproject.toml

@@ -47,6 +47,7 @@ dependencies = [
 	"numpy",
 	"typing_extensions",
     "matplotlib",
+    "packaging",
 ]
 
 [project.urls]
@@ -136,5 +137,6 @@ module = [
 	"meshio",
 	"h5py",
     "matplotlib.*",
+    "packaging.*",
 ]
 ignore_missing_imports = true