Preconditioning for Diffuse Interface Tumor Growth Models

Jessica Bosch; Pelin Çiloğlu; Chen Greif; Martin Stoll

Source-linked AI summary

Preconditioning for Diffuse Interface Tumor Growth Models

Jessica Bosch, Pelin Çiloğlu, Chen Greif, Martin Stoll

arXiv:2608.25537v1Aug 26, 2026math.NA

TL;DR

Efficiently solving the large sparse systems arising in diffuse-interface tumor simulations remains challenging, so this paper develops block preconditioners using Schur-complement approximations for two tumor-growth models. Numerical experiments show robust behavior in two and three dimensions across model parameters, time-step sizes, and mesh refinement.

Problem
Efficient solution of the large sparse linear systems arising in diffuse-interface tumor-growth simulations lacks dedicated preconditioning strategies for Cahn–Hilliard-based models.
Method
The paper develops block preconditioners based on Schur-complement approximations for the smooth and nonsmooth GLSS models and the HKNZ model.
Results
The proposed preconditioners show robust behavior in two- and three-dimensional experiments across model parameters, time-step sizes, and mesh refinement.
Takeaways & Limitations
The methods are presented as suitable for efficiently solving large-scale diffuse-interface tumor-growth simulations.
Takeaways & Limitations
The GLSS setup assumes a quasisteady nutrient and imposes Dirichlet boundary conditions on the nutrient.

Abstract

from arXiv · show

Computational models of tumor growth form a contemporary area of study. Diffuse interface models have been proposed as an important tool for cancer modeling. In this work, we consider preconditioners for numerical simulations of two tumor growth models based on diffuse interface models. For numerical efficiency, we propose block preconditioners that rely on using effective Schur-complement approximations. We show that these methods combined with state-of-the-art adaptive finite element discretizations lead to robust simulations in two and three dimensions. In extensive numerical experiments, we show that the proposed methods show robust convergence behavior.

1. Introduction.

The paper studies efficient solution strategies for two diffuse-interface tumor-growth models, whose large sparse systems pose major computational challenges. It develops Schur-complement-based block preconditioners and reports robust two- and three-dimensional simulations.

Motivation and setting: Diffuse-interface models represent tumor and healthy tissue as phases separated by a diffuse phase-field interface, avoiding explicit tracking and accommodating topological changes.The tumor phase is characterized by u ≈ 1, the healthy phase by u ≈ −1, and interface thickness scales with ε > 0.
Models: The study considers the HKNZ model and smooth and nonsmooth GLSS models, with GLSS additionally incorporating chemotaxis and active transport.HKNZ couples a Cahn–Hilliard equation with a nutrient reaction–diffusion equation.
Computational challenge: Large sparse linear systems must be solved at every time step, including nonsymmetric saddle-point systems from Newton’s method for the smooth GLSS model.The nonsmooth GLSS model uses Moreau–Yosida regularization with a semi-smooth Newton method.
Contribution: The paper develops block preconditioners using Schur-complement approximations for the smooth and nonsmooth GLSS and HKNZ models.The authors identify these as dedicated preconditioning strategies for Cahn–Hilliard-based tumor-growth models and aim for efficient, robust performance.
Contribution: Extensive numerical experiments evaluate the proposed methods in two and three dimensions and show robust convergence behavior.The broader work combines the methods with state-of-the-art adaptive finite element discretizations.

2. Diffuse interface tumor growth models.

The section presents two diffuse-interface tumor growth models: the HKNZ model, with diffusion-governed nutrient dynamics, and the GLSS model, which additionally incorporates chemotaxis and active transport. The GLSS formulation is considered with a quasisteady nutrient and Dirichlet nutrient boundary conditions.

The nutrient phase is modeled without a surface energy or double-well potential, and its uncoupled dynamics are governed by diffusion.
HKNZ model: The HKNZ model formulates the unknown pair (u, σ) as a dissipative gradient flow and uses a double-well potential ψ(u) with minima at ±1.The model also includes a chemical potential w and mass-exchange terms motivated by linear phenomenological laws.
HKNZ model: The HKNZ growth function restricts growth to conditions with available nutrient and helps prevent growth when numerical values of u exceed magnitude 1.
GLSS model: The GLSS model adds chemotaxis and active transport through a transport mechanism absent from HKNZ, while retaining a double-well potential with minima at ±1.The parameter λ can switch off active transport as λ → 0 while preserving chemotaxis.
GLSS model: GLSS is considered with a quasisteady nutrient, neglecting its time derivative, and with Dirichlet boundary conditions imposed on the nutrient.

3. Time discretization.

The paper adopts second-order convexity-splitting time discretization for one tumor model and a smooth-potential scheme for another. Stabilization ensures energy decay under sufficiently large parameters, while the formulations enable linear solves or decoupled solution steps.

3. Time discretization.: A recently proposed second-order convexity-splitting scheme modifies Crank–Nicolson with specialized nonlinear treatment and two stabilizations.The free energy is split into convex and concave parts, treated implicitly and explicitly, respectively.
3. Time discretization.: The α1-stabilization acts as artificial diffusivity, whereas α2-stabilization represents artificial convexity.For sufficiently large α1 and α2, the time-discrete scheme is stable because its corresponding energy functional decreases over time.
3. Time discretization.: The smooth-potential time discretization treats the nonlinear term implicitly and the linear term explicitly.The resulting weak-form scheme decouples σ(n+1), which is obtained by solving a linear equation after the remaining nonlinear system.
3. Time discretization.: For the nonsmooth obstacle potential, the variational inequality is replaced by a Moreau–Yosida regularized problem with penalty parameter 0 < c ≪1.The regularized scheme retains the decoupling of σ(n+1) from the remaining system.

4. Semi-smooth Newton method.

For each regularization parameter c in a sequence tending to zero, the system is solved using a semi-smooth Newton algorithm. The method uses Newton differentiability despite max/min nonsmoothness and is presented with finite-element discretization for the regularized setting.

4. Semi-smooth Newton method.: For each c in a specified sequence c → 0, the system (3.3a)–(3.3c) is solved using a semi-smooth Newton algorithm.The systems are written compactly before applying the iteration.
4. Semi-smooth Newton method.: Because max- and min-operators prevent Fréchet differentiability, Fc instead satisfies the weaker notion of Newton differentiability.This property supports a superlinear convergence result for the semi-smooth Newton iteration.
4. Semi-smooth Newton method.: The operator Gc(u, w) serves as a Newton-derivative for Fc, using χM(u), the characteristic function of M(u) := {x ∈ Ω: |u(x)| > 1}.The construction follows analogous operators from the cited references.
4. Semi-smooth Newton method.: The finite element method is applied to the regularized system, while smooth-potential cases use analogous procedures and receive their fully discrete linear system later.The methodology is presented only for the regularized setting because the procedures are similar.

5. Finite element approximation.

The section develops finite element discretizations for the tumor-growth models, using continuous piecewise multilinear spaces on polyhedral domains. It derives the resulting algebraic systems and records matrix properties relevant to their solution.

Finite element spaces: The spatial discretization assumes a polyhedral domain and triangulation with maximal element size h, approximating H1(Ω) by continuous piecewise multilinear functions.In two dimensions, the functions are bilinear Q1 functions with standard nodal basis functions.
Discrete systems: The discretized tumor-growth equations are formulated in finite element spaces with nodal basis representations, yielding matrix systems for each semismooth Newton step.The coefficient vectors represent the finite element unknowns and previous-step solutions.
Discrete systems: After solving for (uh, wh), the discretized third equation determines σh through a separate linear system with prescribed Dirichlet boundary data.Boundary values are incorporated by extending the basis with fixed coefficients interpolating σB on ∂Ω.
Matrix properties: The mass matrices M and M0 are diagonal symmetric positive definite, ˆK0 is symmetric positive definite, and K and ˆK are symmetric positive semidefinite.The section further notes that K has a one-dimensional kernel spanned by the constant vector 1.
Linear solvers: Because the σh system is a decoupled second-order problem, it is solved using a smoothed aggregation algebraic multigrid preconditioner with MINRES.The AMG method approximates the coefficient matrix ˆK0 + C/2 (Im + u)M0.

6. Preconditioning.

Section 6 develops block preconditioners for diffuse-interface tumor-growth systems by simplifying and approximating their Schur complements. The proposed methods use lower-triangular block structures and approximate operator inversions tailored to the smooth and nonsmooth formulations.

Schur-complement approximations: The preconditioners are constructed by neglecting a block and simplifying the exact Schur complement into a practical approximation.The exact Schur complement is omitted, one term is dropped, and a practical approximation is introduced.
Implementation: Applying the preconditioners requires approximate inversion of the principal blocks and Schur-complement factors.For one formulation, this involves the blocks A1, A2, and a factor of ˆS; another requires the block A and two factors of ˆS.
Block preconditioners: Lower triangular block preconditioners are proposed for the resulting formulations.The section presents lower triangular block structures after deriving the Schur-complement approximations.
Implementation: In two dimensions, selected operators use ILU with fill factor 10 and drop tolerance 10−4, while three-dimensional cases use AMG.The operators 2 ˆ K+B and √τA+M use these ILU settings in two dimensions and AMG in three dimensions.
Formulations and solvers: The preconditioners support smooth HKNZ and (non)smooth GLSS formulations and can be used with a range of Krylov subspace solvers.Both formulations employ left preconditioning; the HKNZ formulation is solved directly at each time step without an outer nonlinear iteration.

7. Numerical results.

Numerical experiments evaluate the proposed preconditioners for HKNZ and GLSS tumor-growth models using finite-element discretizations, adaptive refinement, and iterative nonlinear solvers. Across two- and three-dimensional simulations, iteration counts remain robust under parameter, mesh, and long-time changes, including challenging nonsmooth cases.

Solver setup: HKNZ simulations use MINRES with rtol = 10^-10, while GLSS Newton systems use Bi-CGSTAB with rtol = 10^-15.The GLSS preconditioner is also compatible with GMRES and QMR.
HKNZ results: HKNZ MINRES iteration counts remain bounded across mesh resolution, proliferation rate, and interface thickness in two-dimensional experiments.Iterations increase slightly at higher proliferation rates while retaining robust convergence.
HKNZ results: In three dimensions, the block preconditioner remains robust under mesh refinement, proliferation, and interface-thickness changes, with iteration counts nearly constant over 20 time steps.The simulations include a two-tumor configuration whose tumors eventually merge while nutrient levels are depleted inside the tumor.
Smooth GLSS results: For smooth GLSS simulations, the Newton method converges in only two iterations throughout almost the entire simulation, while Bi-CGSTAB performance remains robust across time-step sizes and initial conditions.The reported experiments compare solver statistics over the first 20 time steps.
Smooth GLSS results: Bi-CGSTAB iteration counts remain within a relatively small range as chemotaxis, proliferation, and nutrient-consumption parameters vary in the smooth GLSS model.The result indicates robustness of the proposed preconditioner with respect to physical model parameters.
Nonsmooth GLSS results: For nonsmooth GLSS, Newton converges within at most four iterations, while the maximum Bi-CGSTAB count reaches 251 over 10,000 time steps and unpreconditioned Bi-CGSTAB fails within 500 iterations.The adaptive mesh averages approximately 20,173 degrees of freedom, and the phase field stays close to obstacle constraints.