Robust and scalable nonlinear solvers for finite element discretizations of biological transportation networks

TL;DR

This paper addresses robust, large-scale simulation of biological transport networks driven by gradient-flow minimization of a non-convex energy, coupled to a Poisson constraint. It introduces a fully implicit nonlinear FEM solver that uses a discontinuous conductivity space to preserve the positive semi-definiteness of the conductivity tensor during backward Euler time stepping, and it designs a Schur-complement preconditioner with AMG for scalable Newton solves. The authors validate their approach in 2D and 3D, demonstrating energy decay, network formation dynamics, and strong/weak scalability up to tens of thousands of processes; they also analyze mesh-, time-integration-, and parameter-related effects on the emergent networks. The results include the first reported three-dimensional simulations of the network formation system, highlighting the method’s ability to handle large-scale, realistic geometries and its potential to explore discretization-driven pattern formation in biological transport networks.

Abstract

We develop robust and scalable fully implicit nonlinear finite element solvers for the simulations of biological transportation networks driven by the gradient flow minimization of a non-convex energy cost functional. Our approach employs a discontinuous space for the conductivity tensor that allows us to guarantee the preservation of its positive semi-definiteness throughout the entire minimization procedure arising from the time integration of the gradient flow dynamics using a backward Euler scheme. Extensive tests in two and three dimensions demonstrate the robustness and performance of the solver, highlight the sensitivity of the emergent network structures to mesh resolution and topology, and validate the resilience of the linear preconditioner to the ill-conditioning of the model. The implementation achieves near-optimal parallel scaling on large-scale, high-performance computing platforms. To the best of our knowledge, the network formation system has never been simulated in three dimensions before. Consequently, our three-dimensional results are the first of their kind.

Figures (18)

• Figure 1: Left: leaf mesh, right: example decomposition with 12 subdomains.
• Figure 2: Parallel scalability: strong scaling (left) and weak scaling (right) results compared against the ideal case. See \ref{['sec:scalability']} for additional details.
• Figure 3: Energy (left panel), $-\frac{\partial E}{\partial t}$ (center), and time step (right), as a function of simulation time for 1024x1024 mesh. See \ref{['sec:sequence_formation']} for additional details.
• Figure 4: $\left\| \mathbb{C} \right\|$ at selected time instances (see \ref{['fig:box_sequence_logs']}) for 1024x1024 mesh in logarithmic scale. See \ref{['sec:sequence_formation']} for additional details.
• Figure 5: Energy (left panel), $-\frac{\partial E}{\partial t}$ (center), and time step (right), as a function of simulation time for 1024x1024 (blue), 2048x2048 (red), and 4096x4096 (green) meshes. See \ref{['sec:be_robust']} for additional details.