diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py index 805511d4a..1598cce5f 100644 --- a/src/festim/hydrogen_transport_problem.py +++ b/src/festim/hydrogen_transport_problem.py @@ -1674,37 +1674,7 @@ def create_solver(self): elif Version(dolfinx.__version__) > Version("0.9.0"): from dolfinx.fem.petsc import NonlinearProblem - if self.petsc_options is None: - # taken from https://github.com/FEniCS/dolfinx/blob/5fcb988c5b0f46b8f9183bc844d8f533a2130d6a/python/demo/demo_cahn-hilliard.py#L279C1-L286C28 - use_superlu = ( - PETSc.IntType == np.int64 - ) # or PETSc.ScalarType == np.complex64 - sys = PETSc.Sys() # type: ignore - if sys.hasExternalPackage("mumps") and not use_superlu: - linear_solver = "mumps" - elif sys.hasExternalPackage("superlu_dist"): - linear_solver = "superlu_dist" - else: - linear_solver = "petsc" - - petsc_options = { - "snes_type": "newtonls", - "snes_linesearch_type": "none", - "snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) - * 1e-2, - # TODO : make atol and rtol callable - "snes_atol": self.settings.atol, - "snes_rtol": self.settings.rtol, - "snes_max_it": self.settings.max_iterations, - "snes_divergence_tolerance": "PETSC_UNLIMITED", - "ksp_type": "preonly", - "pc_type": "lu", - "pc_factor_mat_solver_type": linear_solver, - "snes_error_if_not_converged": True, - "ksp_error_if_not_converged": True, - } - else: - petsc_options = self.petsc_options + petsc_options = self.get_petsc_options() self.solver = NonlinearProblem( self.forms, diff --git a/src/festim/problem.py b/src/festim/problem.py index 610d4a081..7f1d23b77 100644 --- a/src/festim/problem.py +++ b/src/festim/problem.py @@ -10,6 +10,7 @@ from dolfinx import fem from dolfinx.nls.petsc import NewtonSolver from packaging.version import Version +import warnings import festim as F from festim.mesh.mesh import Mesh as _Mesh @@ -118,6 +119,45 @@ def define_boundary_conditions(self): form = self.create_dirichletbc_form(bc) self.bc_forms.append(form) + def get_petsc_options(self) -> dict[str, Any]: + """ + Gets the PETSc options to pass to the NewtonProblem solver. Default + options are updated with user-provided options, if any. + + Returns: + the petsc options to pass to the NewtonProblem solver. + """ + petsc_options = get_default_petsc_options() + + # Update default PETSc options with user-provided options, if any + if self.petsc_options: + petsc_options.update(self.petsc_options) + + if self.petsc_options: + if ( + "snes_atol" in self.petsc_options + or "snes_rtol" in self.petsc_options + or "snes_max_it" in self.petsc_options + ): + warnings.warn( + "You have set one of the following PETSc options: snes_atol, " + "snes_rtol or snes_max_it. These options will be overwritten by " + "the values in festim.Settings (atol, rtol and max_iterations) to " + "ensure consistency between different versions of dolfinx. If you " + "want to set these options manually, please set them in " + "festim.Settings and not in the petsc_options dictionary." + ) + + petsc_options.update( + { + "snes_atol": self.settings.atol, + "snes_rtol": self.settings.rtol, + "snes_max_it": self.settings.max_iterations, + } + ) + + return petsc_options + def create_solver(self): """Creates the solver of the model""" @@ -159,38 +199,7 @@ def create_solver(self): elif Version(dolfinx.__version__) > Version("0.9.0"): from dolfinx.fem.petsc import NonlinearProblem - if self.petsc_options is None: - # taken from https://github.com/FEniCS/dolfinx/blob/5fcb988c5b0f46b8f9183bc844d8f533a2130d6a/python/demo/demo_cahn-hilliard.py#L279C1-L286C28 - use_superlu = ( - PETSc.IntType == np.int64 - ) # or PETSc.ScalarType == np.complex64 - sys = PETSc.Sys() # type: ignore - if sys.hasExternalPackage("mumps") and not use_superlu: - linear_solver = "mumps" - elif sys.hasExternalPackage("superlu_dist"): - linear_solver = "superlu_dist" - else: - linear_solver = "petsc" - - petsc_options = { - "snes_type": "newtonls", - "snes_linesearch_type": "none", - "snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) - * 1e-2, - # TODO : make atol and rtol callable - "snes_atol": self.settings.atol, - "snes_rtol": self.settings.rtol, - "snes_max_it": self.settings.max_iterations, - "snes_divergence_tolerance": "PETSC_UNLIMITED", - "ksp_type": "preonly", - "pc_type": "lu", - "pc_factor_mat_solver_type": linear_solver, - "snes_error_if_not_converged": True, - "ksp_error_if_not_converged": True, - } - else: - petsc_options = self.petsc_options - + petsc_options = self.get_petsc_options() self.solver = NonlinearProblem( self.formulation, self.u, @@ -285,3 +294,31 @@ def update_time_dependent_values(self): for source in self.sources: if source.value.explicit_time_dependent: source.value.update(t=t) + + +# DEFAULT PETSC OPTIONS + +# taken from https://github.com/FEniCS/dolfinx/blob/5fcb988c5b0f46b8f9183bc844d8f533a2130d6a/python/demo/demo_cahn-hilliard.py#L279C1-L286C28 +use_superlu = PETSc.IntType == np.int64 # or PETSc.ScalarType == np.complex64 +sys = PETSc.Sys() # type: ignore +if sys.hasExternalPackage("mumps") and not use_superlu: + linear_solver = "mumps" +elif sys.hasExternalPackage("superlu_dist"): + linear_solver = "superlu_dist" +else: + linear_solver = "petsc" +_DEFAULT_PETSC_OPTS = { + "snes_type": "newtonls", + "snes_linesearch_type": "none", + "snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) * 1e-2, + "snes_divergence_tolerance": "PETSC_UNLIMITED", + "ksp_type": "preonly", + "pc_type": "lu", + "pc_factor_mat_solver_type": linear_solver, + "snes_error_if_not_converged": True, + "ksp_error_if_not_converged": True, +} + + +def get_default_petsc_options(): + return _DEFAULT_PETSC_OPTS.copy()