From 0abc4e96fedc0c84abbbc3085fd265136986a1b6 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Sun, 21 Dec 2025 23:00:21 +0100 Subject: [PATCH 1/9] Move to SNES solver for dolfinx>=0.10 and be ready for dolfinx-adjoint --- .cspell_dict.txt | 9 +++ src/pulse/exceptions.py | 15 +++++ src/pulse/problem.py | 131 ++++++++++++++++++++++++++++------------ 3 files changed, 118 insertions(+), 37 deletions(-) diff --git a/.cspell_dict.txt b/.cspell_dict.txt index 9da9824..e0be40e 100644 --- a/.cspell_dict.txt +++ b/.cspell_dict.txt @@ -21,6 +21,7 @@ caplog cardiomyocytes cellname Chiara +cntl cofac cofnorm cpus @@ -59,6 +60,7 @@ holzapfelogden Hookean hstack hyperelastic +icntl inproceedings Intracomm isclose @@ -70,11 +72,13 @@ Jilberto Karabelas Kasra Kluwer +kwonlyargs Laszlo levelname libgl libxrender linalg +linesearch lmbda Luca Matteo @@ -90,8 +94,11 @@ nabla ndarray neohookean Neumann +newtonls +newtontr Niederer Nolte +nullspace numpy Oreste Orovio @@ -118,6 +125,7 @@ scifem Severi Shrier Silvano +snes Sorine Stefano stica @@ -131,6 +139,7 @@ Usyk varepsilon Venant Verlag +vinewtonrsls Virag viscoelastic viscoelasticity diff --git a/src/pulse/exceptions.py b/src/pulse/exceptions.py index 7b54f11..08af444 100644 --- a/src/pulse/exceptions.py +++ b/src/pulse/exceptions.py @@ -1,3 +1,4 @@ +import logging import operator from dataclasses import dataclass @@ -5,6 +6,9 @@ import dolfinx import numpy as np +import ufl + +logger = logging.getLogger(__name__) def check_value_greater_than( @@ -40,6 +44,17 @@ def check_value_greater_than( bound, ) + elif isinstance(f, ufl.indexed.Indexed): + try: + func, index = f.ufl_operands + value = func.x.array[int(index[0])] + return op(value, bound) + except Exception: + logger.warning( + "Could not extract value from ufl.indexed.Indexed. Assuming check succeeds.", + ) + return True + raise PulseException( # pragma: no cover f"Invalid type for f: {type(f)}. Expected 'float', " "'dolfinx.fem.Constant', 'numpy array' or 'dolfinx.fem.Function'", diff --git a/src/pulse/problem.py b/src/pulse/problem.py index 70b8686..882fc7d 100644 --- a/src/pulse/problem.py +++ b/src/pulse/problem.py @@ -1,3 +1,4 @@ +import inspect import logging import typing from dataclasses import dataclass, field @@ -10,6 +11,7 @@ import numpy as np import scifem import ufl +from packaging.version import Version from .boundary_conditions import BoundaryConditions from .cardiac_model import CardiacModel @@ -17,6 +19,7 @@ from .units import Variable, mesh_factor T = typing.TypeVar("T", dolfinx.fem.Function, np.ndarray) +_dolfinx_version = Version(dolfinx.__version__) logger = logging.getLogger(__name__) @@ -76,6 +79,10 @@ class StaticProblem: parameters: dict[str, typing.Any] = field(default_factory=dict) bcs: BoundaryConditions = field(default_factory=BoundaryConditions) cavities: list[Cavity] = field(default_factory=list) + NonlinearProblem: typing.Type[dolfinx.fem.petsc.NonlinearProblem] = ( + dolfinx.fem.petsc.NonlinearProblem + ) + Function: typing.Type[dolfinx.fem.Function] = dolfinx.fem.Function def __post_init__(self): parameters = type(self).default_parameters() @@ -112,8 +119,8 @@ def _init_p_space(self): degree=int(p_degree), ) self.p_space = dolfinx.fem.functionspace(self.geometry.mesh, p_element) - self.p = dolfinx.fem.Function(self.p_space) - self.p_old = dolfinx.fem.Function(self.p_space) + self.p = self.Function(self.p_space) + self.p_old = self.Function(self.p_space) self.p_test = ufl.TestFunction(self.p_space) self.dp = ufl.TrialFunction(self.p_space) else: @@ -137,8 +144,8 @@ def _init_u_space(self): shape=(self.geometry.mesh.topology.dim,), ) self.u_space = dolfinx.fem.functionspace(self.geometry.mesh, u_element) - self.u = dolfinx.fem.Function(self.u_space) - self.u_old = dolfinx.fem.Function(self.u_space) + self.u = self.Function(self.u_space) + self.u_old = self.Function(self.u_space) self.u_test = ufl.TestFunction(self.u_space) self.du = ufl.TrialFunction(self.u_space) @@ -159,6 +166,20 @@ def default_parameters(): "ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps", + "snes_error_if_not_converged": True, + "ksp_error_if_not_converged": True, + # "snes_monitor": None, + # "ksp_monitor": None, + "snes_linesearch_monitor": None, + "snes_type": "newtonls", + # "snes_type": "newtontr", + # "snes_type": "vinewtonrsls", + "snes_linesearch_type": "none", + # "mat_mumps_icntl_24": 1, # Zero pivot detection + # "mat_mumps_icntl_25": 0, # Which nullspace to extract + # "mat_mumps_icntl_4": 1, # Verbosity + # "mat_mumps_icntl_2": 1, # std out + # "mat_mumps_cntl_3": 1e-6, # Threshold factor }, } @@ -177,8 +198,8 @@ def _init_cavity_pressure_spaces(self): logger.debug(f"Number of cavity pressure states: {self.num_cavity_pressure_states}") self.real_space = scifem.create_real_functionspace(self.geometry.mesh) for _ in range(self.num_cavity_pressure_states): - cavity_pressure = dolfinx.fem.Function(self.real_space) - cavity_pressure_old = dolfinx.fem.Function(self.real_space) + cavity_pressure = self.Function(self.real_space) + cavity_pressure_old = self.Function(self.real_space) cavity_pressure_test = ufl.TestFunction(self.real_space) cavity_pressure_trial = ufl.TrialFunction(self.real_space) @@ -197,8 +218,8 @@ def _init_rigid_body(self): self.geometry.mesh, value_shape=(6,), ) - self.r = dolfinx.fem.Function(self.rigid_space) - self.r_old = dolfinx.fem.Function(self.rigid_space) + self.r = self.Function(self.rigid_space) + self.r_old = self.Function(self.rigid_space) self.dr = ufl.TrialFunction(self.rigid_space) self.q = ufl.TestFunction(self.rigid_space) else: @@ -352,7 +373,7 @@ def base_dirichlet(self): marker = self.geometry.markers[self.parameters["base_marker"]][0] base_facets = self.geometry.facet_tags.find(marker) dofs_base = dolfinx.fem.locate_dofs_topological(self.u_space, 2, base_facets) - u_bc_base = dolfinx.fem.Function(self.u_space) + u_bc_base = self.Function(self.u_space) u_bc_base.x.array[:] = 0 bcs.append(dolfinx.fem.dirichletbc(u_bc_base, dofs_base)) @@ -427,8 +448,12 @@ def update_old_states(self): self.cavity_pressures_old[i].x.array[:] = self.cavity_pressures[i].x.array.copy() if self.parameters["rigid_body_constraint"]: logger.debug("Updating old rigid body constraint state to current value...") + assert self.r is not None + assert self.r_old is not None self.r_old.x.array[:] = self.r.x.array.copy() if self.is_incompressible: + assert self.p is not None + assert self.p_old is not None logger.debug("Updating old pressure state to current value...") self.p_old.x.array[:] = self.p.x.array.copy() @@ -442,8 +467,12 @@ def reset_states(self): self.cavity_pressures[i].x.array[:] = self.cavity_pressures_old[i].x.array.copy() if self.parameters["rigid_body_constraint"]: logger.debug("Resetting rigid body constraint state to old value...") + assert self.r is not None + assert self.r_old is not None self.r.x.array[:] = self.r_old.x.array.copy() if self.is_incompressible: + assert self.p is not None + assert self.p_old is not None logger.debug("Resetting pressure state to old value...") self.p.x.array[:] = self.p_old.x.array.copy() @@ -496,14 +525,47 @@ def _init_forms(self) -> None: bcs = self.base_dirichlet logger.debug("Creating Newton solver...") - self._solver = scifem.NewtonSolver( - R, - K, - u, - bcs=bcs, - max_iterations=25, - petsc_options=self.parameters["petsc_options"], - ) + + # Hack to pass petsc options to dolfinx_adjoint NonlinearProblem + kwargs = {} + if "adjoint_petsc_options" in inspect.getfullargspec(self.NonlinearProblem).kwonlyargs: + kwargs["adjoint_petsc_options"] = self.parameters["petsc_options"] + + if _dolfinx_version >= Version("0.10"): + # Until we have implemented NonlinearProblem for blocked systems + # in dolfinx_adjoint, we only support single state problems here. + # Therefore we extract the first block. + if self.num_states == 1: + R = R[0] + u = u[0] + K = K[0][0] + + kwargs["petsc_options_prefix"] = "pulse_problem_" + kwargs["petsc_options"] = self.parameters["petsc_options"] + + self.problem = self.NonlinearProblem( + F=R, + J=K, + u=u, + bcs=bcs, + **kwargs, + ) + else: + petsc_options = self.parameters["petsc_options"] + # Pop options that are not supported in older dolfinx versions + petsc_options.pop("snes_error_if_not_converged", None) + petsc_options.pop("snes_type", None) + petsc_options.pop("snes_linesearch_type", None) + + # Keep old behavior for older dolfinx versions + self._solver = scifem.NewtonSolver( + R, + K, + u, + bcs=bcs, + max_iterations=25, + petsc_options=petsc_options, + ) def update_fields(self): pass @@ -513,23 +575,11 @@ def reset(self): self.reset_states() self._init_forms() - def solve( - self, - rtol: float = 1e-10, - atol: float = 1e-6, - beta: float = 1.0, - update_old_states: bool = True, - ) -> bool: + def solve(self, update_old_states: bool = True) -> bool: """Solve nonlinear problem with Newton solver Parameters ---------- - rtol : float, optional - Relative tolerance, by default 1e-10 - atol : float, optional - Absolute tolerance, by default 1e-6 - beta : float, optional - Damping parameter, by default 1.0 update_old_states : bool, optional Whether to update old states before solving, by default True @@ -541,10 +591,17 @@ def solve( logger.debug("Solving the system...") if update_old_states: self.update_old_states() - ret = self._solver.solve(rtol=rtol, atol=atol, beta=beta) + if _dolfinx_version >= Version("0.10"): + self.problem.solve() + converged = self.problem.solver.getConvergedReason() > 0 + iters = self.problem.solver.getIterationNumber() + logger.debug(f"Solved in {iters} iterations, converged: {converged}") + else: + converged = self._solver.solve() + self.update_fields() - return ret + return converged class DynamicProblem(StaticProblem): @@ -557,14 +614,14 @@ def __post_init__(self): def _init_u_space(self): super()._init_u_space() - self.u_old = dolfinx.fem.Function(self.u_space) - self.v_old = dolfinx.fem.Function(self.u_space) - self.a_old = dolfinx.fem.Function(self.u_space) + self.u_old = self.Function(self.u_space) + self.v_old = self.Function(self.u_space) + self.a_old = self.Function(self.u_space) def _init_p_space(self): super()._init_p_space() if self.is_incompressible: - self.p_old = dolfinx.fem.Function(self.p_space) + self.p_old = self.Function(self.p_space) else: self.p_old = None @@ -572,7 +629,7 @@ def _init_cavity_pressure_spaces(self): super()._init_cavity_pressure_spaces() self.cavity_pressures_old = [] for _ in range(self.num_cavity_pressure_states): - cavity_pressure_old = dolfinx.fem.Function(self.real_space) + cavity_pressure_old = self.Function(self.real_space) self.cavity_pressures_old.append(cavity_pressure_old) def _material_form(self, u, v, p): From 2d15ec855fb965f33411f14bba716a29e7e8def1 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Mon, 22 Dec 2025 15:52:08 +0100 Subject: [PATCH 2/9] Specify the same tolerances for the old solver --- src/pulse/problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pulse/problem.py b/src/pulse/problem.py index 882fc7d..0deef01 100644 --- a/src/pulse/problem.py +++ b/src/pulse/problem.py @@ -597,7 +597,7 @@ def solve(self, update_old_states: bool = True) -> bool: iters = self.problem.solver.getIterationNumber() logger.debug(f"Solved in {iters} iterations, converged: {converged}") else: - converged = self._solver.solve() + converged = self._solver.solve(rtol=1e-10, atol=1e-6) self.update_fields() From d91bff3e9067cdd88e1ee5b343348d9b9174a23c Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Wed, 18 Mar 2026 15:25:40 -0700 Subject: [PATCH 3/9] Bugfix in active stress --- src/pulse/active_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pulse/active_model.py b/src/pulse/active_model.py index ec49239..4081e58 100644 --- a/src/pulse/active_model.py +++ b/src/pulse/active_model.py @@ -130,7 +130,7 @@ def S(self, C: ufl.core.expr.Expr, dev: bool = False) -> ufl.core.expr.Expr: Cdev = kinematics.Cdev(C) else: Cdev = C - return 2.0 * ufl.diff(self.strain_energy(Cdev), Cdev) + return 2.0 * ufl.diff(self.strain_energy(Cdev), C) def P(self, F: ufl.core.expr.Expr, dev: bool = False) -> ufl.core.expr.Expr: """First Piola-Kirchhoff stress tensor for the active model. From 0e1a20c35a785ca5f27396c6dd6e044cad5ebc96 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Wed, 18 Mar 2026 15:53:36 -0700 Subject: [PATCH 4/9] Make it possible to use PC Type bddc --- .cspell_dict.txt | 3 +++ src/pulse/problem.py | 22 ++++++++++++++++++++++ 2 files changed, 25 insertions(+) diff --git a/.cspell_dict.txt b/.cspell_dict.txt index e0be40e..20d77e0 100644 --- a/.cspell_dict.txt +++ b/.cspell_dict.txt @@ -10,6 +10,7 @@ Barnafi Bartolucci basix bcast +bddc behaviour bestel Bestel @@ -81,6 +82,7 @@ linalg linesearch lmbda Luca +MATIS Matteo mechanicsproblem ment @@ -134,6 +136,7 @@ Taras tomek Tomek TRAME +unassembled ureg Usyk varepsilon diff --git a/src/pulse/problem.py b/src/pulse/problem.py index 0deef01..0c4e5fd 100644 --- a/src/pulse/problem.py +++ b/src/pulse/problem.py @@ -550,6 +550,28 @@ def _init_forms(self) -> None: bcs=bcs, **kwargs, ) + + if kwargs["petsc_options"].get("pc_type") == "bddc": + # Create the unassembled MATIS matrix required by BDDC + # Note: problem.a is the compiled bilinear form of the Jacobian + from petsc4py import PETSc + + J_form = dolfinx.fem.form(K) + A_is = dolfinx.fem.petsc.create_matrix(J_form, kind=PETSc.Mat.Type.IS) + + # Define a custom Python callback to assemble the + # MATIS matrix at every Newton step + def compute_jacobian(snes, x, J, P): + J.zeroEntries() + dolfinx.fem.petsc.assemble_matrix(J, J_form, bcs=bcs) + J.assemble() + + # 4. Retrieve the SNES solver and swap both the matrix AND the callback + snes = self.problem.solver + snes.setJacobian(compute_jacobian, J=A_is, P=A_is) + + # Apply options so PETSc registers the MATIS matrix before setup + snes.setFromOptions() else: petsc_options = self.parameters["petsc_options"] # Pop options that are not supported in older dolfinx versions From 043dc75205323d64b6cd9bcef16d54825369c0d3 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 27 Mar 2026 08:17:31 +0100 Subject: [PATCH 5/9] Add some more default options to snes solver to make sure it converges for all tests --- src/pulse/problem.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/pulse/problem.py b/src/pulse/problem.py index c0e148c..33e9be8 100644 --- a/src/pulse/problem.py +++ b/src/pulse/problem.py @@ -174,9 +174,14 @@ def default_parameters(): # "ksp_monitor": None, "snes_linesearch_monitor": None, "snes_type": "newtonls", + "snes_atol": 1e-6, + "snes_rtol": 1e-10, + # "snes_stol": 1e-8, + "snes_max_it": 50, # "snes_type": "newtontr", # "snes_type": "vinewtonrsls", - "snes_linesearch_type": "none", + # "snes_linesearch_type": "none", + "snes_linesearch_type": "secant", # "mat_mumps_icntl_24": 1, # Zero pivot detection # "mat_mumps_icntl_25": 0, # Which nullspace to extract # "mat_mumps_icntl_4": 1, # Verbosity From f624b802a71f9468672bc40f9316763a1a6cb91d Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 27 Mar 2026 11:02:38 +0100 Subject: [PATCH 6/9] Add non-zero initial guess in test --- .cspell_dict.txt | 3 ++- src/pulse/problem.py | 6 ++++-- tests/test_static_problem.py | 7 ++++--- 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/.cspell_dict.txt b/.cspell_dict.txt index c16fd3c..d7c3dd3 100644 --- a/.cspell_dict.txt +++ b/.cspell_dict.txt @@ -92,6 +92,7 @@ lmbda Luca MATIS Matteo +maxstep mechanicsproblem ment meshtags @@ -146,8 +147,8 @@ Taras tomek Tomek TRAME -unassembled transmural +unassembled ureg Usyk varepsilon diff --git a/src/pulse/problem.py b/src/pulse/problem.py index 33e9be8..374cfb7 100644 --- a/src/pulse/problem.py +++ b/src/pulse/problem.py @@ -176,12 +176,14 @@ def default_parameters(): "snes_type": "newtonls", "snes_atol": 1e-6, "snes_rtol": 1e-10, - # "snes_stol": 1e-8, + "snes_stol": 1e-8, "snes_max_it": 50, + # "snes_linesearch_maxstep": 50, + # "snes_view": None, # "snes_type": "newtontr", # "snes_type": "vinewtonrsls", # "snes_linesearch_type": "none", - "snes_linesearch_type": "secant", + "snes_linesearch_type": "l2", # "mat_mumps_icntl_24": 1, # Zero pivot detection # "mat_mumps_icntl_25": 0, # Which nullspace to extract # "mat_mumps_icntl_4": 1, # Verbosity diff --git a/tests/test_static_problem.py b/tests/test_static_problem.py index a190145..9b93ca2 100644 --- a/tests/test_static_problem.py +++ b/tests/test_static_problem.py @@ -150,7 +150,7 @@ def handle_control_lv( def gen(problem): dvlp = (target_lvp - initial_lvp) / 3.0 - for lvp in np.arange(0.0, target_lvp + 1e-9, dvlp): + for lvp in np.arange(initial_lvp + dvlp, target_lvp + 1e-9, dvlp): print(f"Solve Pressure: {lvp}") traction.assign(lvp) problem.solve() @@ -166,7 +166,8 @@ def gen(problem): def gen(problem): dlv = (target_lvv - initial_volume) / 4.0 - for lvv in np.arange(initial_volume, target_lvv + 1e-9, dlv): + for lvv in np.arange(initial_volume + dlv, target_lvv + 1e-9, dlv): + problem.u.x.array[:] = 0.00001 print("Solve Volume: ", lvv) volume.value = lvv problem.solve() @@ -247,7 +248,7 @@ def test_static_problem_lv( }, cavities=cavities, ) - problem.solve() + for lvp, lvv in gen(problem): print(f"LVP: {lvp}, LVV: {lvv}") From b3485422b64575c3afbecc89745efb160ac71ac9 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 27 Mar 2026 11:04:01 +0100 Subject: [PATCH 7/9] Turn off linesearch monitor --- src/pulse/problem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pulse/problem.py b/src/pulse/problem.py index 374cfb7..7715804 100644 --- a/src/pulse/problem.py +++ b/src/pulse/problem.py @@ -172,7 +172,7 @@ def default_parameters(): "ksp_error_if_not_converged": True, # "snes_monitor": None, # "ksp_monitor": None, - "snes_linesearch_monitor": None, + # "snes_linesearch_monitor": None, "snes_type": "newtonls", "snes_atol": 1e-6, "snes_rtol": 1e-10, From a5c8099fa6970e064fc80b17a9c006e1aa53b44e Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 27 Mar 2026 13:13:35 +0100 Subject: [PATCH 8/9] Remove unused petsc options on old fenics version --- src/pulse/problem.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/pulse/problem.py b/src/pulse/problem.py index 7715804..e863f9f 100644 --- a/src/pulse/problem.py +++ b/src/pulse/problem.py @@ -609,6 +609,10 @@ def compute_jacobian(snes, x, J, P): petsc_options.pop("snes_error_if_not_converged", None) petsc_options.pop("snes_type", None) petsc_options.pop("snes_linesearch_type", None) + petsc_options.pop("snes_atol", None) + petsc_options.pop("snes_rtol", None) + petsc_options.pop("snes_stol", None) + petsc_options.pop("snes_max_it", None) # Keep old behavior for older dolfinx versions self._solver = scifem.NewtonSolver( From ab8241b57aa68422160cf0b932af5391461028ba Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Fri, 27 Mar 2026 13:28:57 +0100 Subject: [PATCH 9/9] Skip test on v0.9 --- .github/workflows/test_package_coverage.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/test_package_coverage.yml b/.github/workflows/test_package_coverage.yml index 554f3cd..48b472b 100644 --- a/.github/workflows/test_package_coverage.yml +++ b/.github/workflows/test_package_coverage.yml @@ -17,7 +17,6 @@ jobs: fail-fast: false matrix: container: [ - "ghcr.io/fenics/dolfinx/dolfinx:v0.9.0", "ghcr.io/fenics/dolfinx/dolfinx:stable", "ghcr.io/fenics/dolfinx/dolfinx:nightly" ]