diff --git a/.cspell_dict.txt b/.cspell_dict.txt index 795a5ff..d7c3dd3 100644 --- a/.cspell_dict.txt +++ b/.cspell_dict.txt @@ -10,6 +10,7 @@ Barnafi Bartolucci basix bcast +bddc behaviour bestel Bestel @@ -78,16 +79,20 @@ Jilberto Karabelas Kasra Kluwer +kwonlyargs Laszlo LDRB levelname libgl libxrender linalg +linesearch linspace lmbda Luca +MATIS Matteo +maxstep mechanicsproblem ment meshtags @@ -101,8 +106,11 @@ nabla ndarray neohookean Neumann +newtonls +newtontr Niederer Nolte +nullspace numpy Oreste Orovio @@ -130,6 +138,7 @@ SELLIER Severi Shrier Silvano +snes Sorine Stefano stica @@ -139,11 +148,13 @@ tomek Tomek TRAME transmural +unassembled ureg Usyk varepsilon Venant Verlag +vinewtonrsls Virag viscoelastic viscoelasticity 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" ] 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. 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 c6acf7a..e863f9f 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() @@ -113,8 +120,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: @@ -138,8 +145,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, name="u") - self.u_old = dolfinx.fem.Function(self.u_space, name="u_old") + self.u = self.Function(self.u_space, name="u") + self.u_old = self.Function(self.u_space, name="u_old") self.u_test = ufl.TestFunction(self.u_space) self.du = ufl.TrialFunction(self.u_space) self.u_full = dolfinx.fem.Function(self.u_space, name="u_full") @@ -161,8 +168,24 @@ 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_atol": 1e-6, + "snes_rtol": 1e-10, + "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": "l2", # "mat_mumps_icntl_24": 1, # Zero pivot detection - # "mat_mumps_icntl_25": 0, # Which null space to extract + # "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 @@ -184,8 +207,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) @@ -204,8 +227,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: @@ -381,7 +404,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)) @@ -456,8 +479,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() @@ -471,8 +498,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() @@ -525,14 +556,73 @@ 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, + ) + + 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 + 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( + R, + K, + u, + bcs=bcs, + max_iterations=25, + petsc_options=petsc_options, + ) def update_fields(self): pass @@ -542,23 +632,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 @@ -570,10 +648,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(rtol=1e-10, atol=1e-6) + self.update_fields() - return ret + return converged class DynamicProblem(StaticProblem): @@ -586,14 +671,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 @@ -601,7 +686,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): 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}")