Patrick E. Farrell

Patrick E. Farrell, Tim van Beeck, Umberto Zerbinati

We analyse the nematic Helmholtz-Korteweg equation, a variant of the classical Helmholtz equation that describes time-harmonic wave propagation in calamitic fluids in the presence of nematic order. A prominent example is given by nematic liquid crystals, which can be modeled as nematic Korteweg fluids - that is, fluids whose stress tensor depends on density gradients and on a nematic director describing the orientation of the anisotropic molecules. These materials exhibit anisotropic acoustic properties that can be tuned by external electromagnetic fields, making them attractive for potential applications such as tunable acoustic resonators. We prove the existence and uniqueness of solutions to this equation in two and three dimensions for suitable (nonresonant) wave numbers and propose a convergent discretisation for its numerical solution. The discretisation of this problem is nontrivial as it demands high regularity and involves unfamiliar boundary conditions. We address these challenges by using high-order conforming finite elements and enforcing the boundary conditions with Nitsche's method. We illustrate our analysis with numerical simulations in two dimensions.

Patrick E. Farrell, Corrado Maurini

The variational approach to fracture is effective for simulating the nucleation and propagation of complex crack patterns, but is computationally demanding. The model is a strongly nonlinear non-convex variational inequality that demands the resolution of small length scales. The current standard algorithm for its solution, alternate minimization, is robust but converges slowly and demands the solution of large, ill-conditioned linear subproblems. In this paper, we propose several advances in the numerical solution of this model that improve its computational efficiency. We reformulate alternate minimization as a nonlinear Gauss-Seidel iteration and employ over-relaxation to accelerate its convergence; we compose this accelerated alternate minimization with Newton's method, to further reduce the time to solution; and we formulate efficient preconditioners for the solution of the linear subproblems arising in both alternate minimization and in Newton's method. We investigate the improvements in efficiency on several examples from the literature; the new solver is 5--6$\times$ faster on a majority of the test cases

Patrick E. Farrell, Casper H. L. Beentjes, Ásgeir Birkisson

Arclength continuation and branch switching are enormously successful algorithms for the computation of bifurcation diagrams. Nevertheless, their combination suffers from three significant disadvantages. The first is that they attempt to compute only the part of the diagram that is continuously connected to the initial data; disconnected branches are overlooked. The second is that the subproblems required (typically determinant calculation and nullspace construction) are expensive and hard to scale to very large discretizations. The third is that they can miss connected branches associated with nonsimple bifurcations, such as when an eigenvalue of even multiplicity crosses the origin. Without expert knowledge or lucky guesses, these techniques alone can paint an incomplete picture of the dynamics of a system. In this paper we propose a new algorithm for computing bifurcation diagrams, called deflated continuation, that is capable of overcoming all three of these disadvantages. The algorithm combines classical continuation with a deflation technique that elegantly eliminates known branches from consideration, allowing the discovery of disconnected branches with Newton's method. Deflated continuation does not rely on any device for detecting bifurcations and does not involve computing eigendecompositions; all subproblems required in deflated continuation can be solved efficiently if a good preconditioner is available for the underlying nonlinear problem. We prove sufficient conditions for the convergence of Newton's method to multiple solutions from the same initial guess, providing insight into which unknown branches will be discovered. We illustrate the success of the method on several examples where standard techniques fail.

Patrick E. Farrell, John W. Pearson

We propose a new preconditioner for the Ohta--Kawasaki equation, a nonlocal Cahn--Hilliard equation that describes the evolution of diblock copolymer melts. We devise a computable approximation to the inverse of the Schur complement of the coupled second-order formulation via a matching strategy. The preconditioner achieves mesh independence: as the mesh is refined, the number of Krylov iterations required for its solution remains approximately constant. In addition, the preconditioner is robust with respect to the interfacial thickness parameter if a timestep criterion is satisfied. This enables the highly resolved finite element simulation of three-dimensional diblock copolymer melts with over one billion degrees of freedom.

Patrick E. Farrell

We prove a new theorem relating the number of distinct eigenvalues of a matrix after perturbation to the prior number of distinct eigenvalues, the rank of the update, and the degree of nondiagonalizability of the matrix. In particular, a rank one update applied to a diagonalizable matrix can at most double the number of distinct eigenvalues. The theorem applies to both symmetric and nonsymmetric matrices and perturbations, of arbitrary magnitudes. An an application, we prove that in exact arithmetic the number of Krylov iterations required to exactly solve a linear system involving a diagonalizable matrix can at most double after a rank one update.

Patrick E. Farrell, Michael Neilan, Charles Parker et al.

We present new convergence estimates for the iterated penalty method applied to structure-preserving discretizations of linear generalized saddle point systems. The method may be viewed as an Uzawa iteration on an augmented Lagrangian formulation of the system. As a by-product, we obtain sharper stability estimates for penalized/perturbed saddle point problems. Three model finite element applications show agreement with the theory.

Ganghui Zhang, Boris D. Andrews, Patrick E. Farrell

Geometric flows, where an immersed manifold evolves in time according to its own geometry, exhibit important structural properties. For example, surface diffusion dissipates surface area while conserving volume; it is desirable to preserve these properties on discretization. This has motivated a substantial body of research on structure-preserving discretizations for these flows, albeit at low order in time. In this work, we present the first discretization of geometric curvature flows (curve shortening/mean curvature flow and curve/surface diffusion) that preserves the evolution of area and volume at arbitrary order in space and time. The key idea is to introduce auxiliary variables in a particular way so that the derivation of the area dissipation law can be replicated after discretization with continuous Petrov--Galerkin in time. These auxiliary variables are indicated by a general strategy for structure-preservation in time that applies to many other problems. The proposed scheme also preserves mesh quality in the same manner as the minimal deformation rate strategy. We demonstrate its structure-preserving properties and high-order convergence on several benchmark examples.

Aaron Baier-Reinio, Patrick E. Farrell, Charles W. Monroe

We present a broad family of high-order finite element algorithms for simulating the flow of electroneutral electrolytes. The governing partial differential equations that we solve are the electroneutral Navier--Stokes--Onsager--Stefan--Maxwell (NSOSM) equations, which model momentum transport, multicomponent diffusion and electrical effects within the electrolyte. Our algorithms can be applied in the steady and transient settings, in two and three spatial dimensions, and under a variety of boundary conditions. Moreover, we allow for the material parameters (e.g. viscosity, diffusivities, thermodynamic factors and density) to be dependent on the local state of the mixture and thermodynamically non-ideal. The flexibility of our approach requires us to address subtleties that arise in the governing equations due to the interplay between boundary conditions and the equation of state. We demonstrate the algorithms in various physical configurations, including (i) electrolyte flow around a microfluidic rotating disk electrode and (ii) the flow in a Hull cell of a cosolvent electrolyte mixture used in lithium-ion batteries.

Pablo D. Brubeck, Patrick E. Farrell, Robert C. Kirby et al.

We present new finite elements for solving the Riesz maps of the de Rham complex on triangular and tetrahedral meshes at high order. The finite elements discretize the same spaces as usual, but with different basis functions, so that the resulting matrices have desirable properties. These properties mean that we can solve the Riesz maps to a given accuracy in a $p$-robust number of iterations with $\mathcal{O}(p^6)$ flops in three dimensions, rather than the naïve $\mathcal{O}(p^9)$ flops. The degrees of freedom build upon an idea of Demkowicz et al., and consist of integral moments on an equilateral reference simplex with respect to a numerically computed polynomial basis that is orthogonal in two different inner products. As a result, the interior-interface and interior-interior couplings are provably weak, and we devise a preconditioning strategy by neglecting them. The combination of this approach with a space decomposition method on vertex and edge star patches allows us to efficiently solve the canonical Riesz maps at high order. We apply this to solving the Hodge Laplacians of the de Rham complex with novel augmented Lagrangian preconditioners.

Patrick E. Farrell, Mingdong He, Kaibo Hu et al.

Magnetic relaxation drives plasma toward lower-energy equilibria under helicity constraints. In ideal magnetohydrodynamics (MHD), helicity is locally conserved, while resistive theories such as Taylor relaxation preserve only global helicity. This distinction has important implications for structure-preserving numerical methods. We compare three finite element formulations: an unconstrained scheme that does not conserve helicity, a mixed method based on finite element exterior calculus that preserves all local helicities, and a Lagrange multiplier approach that enforces only global helicity conservation. Numerical experiments on braided and knotted magnetic fields show that local helicity preservation prevents spurious reconnection and maintains nontrivial topology in ideal MHD or magneto-friction, whereas enforcing only global helicity allows further relaxation through local reconnection. Numerical results on magnetic knots and braids are provided. These results clarify how different levels of discrete helicity constraints influence magnetic relaxation and equilibrium structure in numerical computation.

Kars Knook, Aaron Baier-Reinio, Patrick E. Farrell

The Onsager-Stefan-Maxwell (OSM) equations are an important model of mass transport in multicomponent flows with multiple chemical species. They describe the coupling of diffusive fluxes between species, accounting for their interactions through frictional and thermodynamic driving forces. In this work we propose an augmented Lagrangian preconditioner and prove its discretization-robustness for a Picard linearization of the stationary OSM equations in the isobaric, isothermal, ideal gaseous setting. For the Newton linearization we employ the augmented Lagrangian preconditioner as a block diagonal smoother inside a monolithic geometric multigrid iteration and combine with vertex star Schwarz methods. This strategy is shown to be applicable in a wide variety of settings which incorporate cross-diffusion, nonideal mixing, thermal, pressure, convective, and electrochemical effects. We demonstrate robustness or mild dependence with respect to mesh refinement and polynomial degree of the proposed monolithic preconditioning strategy for different types of multicomponent flows in several applications: cross-diffusion in the human airways, separation of gases under a temperature gradient, nonideal mixing of benzene and cyclohexane, and electrolytic transport in a Hull cell undergoing electroplating.

Patrick E. Farrell, Matteo Croci, Thomas M. Surowiec

Variational inequalities can in general support distinct solutions. In this paper we study an algorithm for computing distinct solutions of a variational inequality, without varying the initial guess supplied to the solver. The central idea is the combination of a semismooth Newton method with a deflation operator that eliminates known solutions from consideration. Given one root of a semismooth residual, deflation constructs a new problem for which a semismooth Newton method will not converge to the known root, even from the same initial guess. This enables the discovery of other roots. We prove the effectiveness of the deflation technique under the same assumptions that guarantee locally superlinear convergence of a semismooth Newton method. We demonstrate its utility on various finite- and infinite-dimensional examples drawn from constrained optimization, game theory, economics and solid mechanics.

Matteo Croci, Michael B. Giles, Marie E. Rognes et al.

When solving stochastic partial differential equations (SPDEs) driven by additive spatial white noise, the efficient sampling of white noise realizations can be challenging. Here, we present a new sampling technique that can be used to efficiently compute white noise samples in a finite element method and multilevel Monte Carlo (MLMC) setting. The key idea is to exploit the finite element matrix assembly procedure and factorize each local mass matrix independently, hence avoiding the factorization of a large matrix. Moreover, in a MLMC framework, the white noise samples must be coupled between subsequent levels. We show how our technique can be used to enforce this coupling even in the case of non-nested mesh hierarchies. We demonstrate the efficacy of our method with numerical experiments. We observe optimal convergence rates for the finite element solution of the elliptic SPDEs of interest in 2D and 3D and we show convergence of the sampled field covariances. In a MLMC setting, a good coupling is enforced and the telescoping sum is respected.

Patrick E. Farrell, Johan E. Hake, Simon W. Funke et al.

Mathematical models that couple partial differential equations (PDEs) and spatially distributed ordinary differential equations (ODEs) arise in biology, medicine, chemistry and many other fields. In this paper we discuss an extension to the FEniCS finite element software for expressing and efficiently solving such coupled systems. Given an ODE described using an augmentation of the Unified Form Language (UFL) and a discretisation described by an arbitrary Butcher tableau, efficient code is automatically generated for the parallel solution of the ODE. The high-level description of the solution algorithm also facilitates the automatic derivation of the adjoint and tangent linearization of coupled PDE-ODE solvers. We demonstrate the capabilities of the approach on examples from cardiac electrophysiology and mitochondrial swelling.