Search

The Parker conjecture, which explores whether magnetic fields in perfectly conducting plasmas can develop tangential discontinuities during magnetic relaxation, remains an open question in astrophysics. Helicity conservation provides a topological barrier during relaxation, preventing topologically nontrivial initial data relaxing to trivial solutions; preserving this mechanism discretely over long time periods is therefore crucial for numerical simulation. This work presents an energy- and helicity-preserving finite element discretization for the magneto-frictional system, for investigating the Parker conjecture. The algorithm preserves a discrete version of the topological barrier and a discrete Arnold inequality. We also discuss extensions to domains with nontrivial topology.

We derive the Helmholtz--Korteweg equation, which models acoustic waves in Korteweg fluids. We further derive a nematic variant of the Helmholtz-Korteweg equation, which incorporates an additional orientational term in the stress tensor. Its dispersion relation coincides with that arising in Virga's analysis of the Euler-Korteweg equations, which we extend to consider imaginary wave numbers and the effect of boundary conditions. In particular, our extensions allow us to analyze the effect of nematic orientation on the penetration depth of evanescent plane waves, and on the scattering of sound waves by obstacles. Furthermore, we make new, experimentally-verifiable predictions for the effect of boundary conditions for a modification of the Mullen-L\"uthi-Stephen experiment, and for the scattering of acoustic waves in nematic-Korteweg fluids by a circular obstacle.

We derive and analyze a broad class of finite element methods for numerically simulating the stationary, low Reynolds number flow of concentrated mixtures of several distinct chemical species in a common thermodynamic phase. The underlying partial differential equations that we discretize are the Stokes$\unicode{x2013}$Onsager$\unicode{x2013}$Stefan$\unicode{x2013}$Maxwell (SOSM) equations, which model bulk momentum transport and multicomponent diffusion within ideal and non-ideal mixtures. Unlike previous approaches, the methods are straightforward to implement in two and three spatial dimensions, and allow for high-order finite element spaces to be employed. The key idea in deriving the discretization is to suitably reformulate the SOSM equations in terms of the species mass fluxes and chemical potentials, and discretize these unknown fields using stable $H(\textrm{div}) \unicode{x2013} L^2$ finite element pairs. We prove that the methods are convergent and yield a symmetric linear system for a Picard linearization of the SOSM equations, which staggers the updates for concentrations and chemical potentials. We also discuss how the proposed approach can be extended to the Newton linearization of the SOSM equations, which requires the simultaneous solution of mole fractions, chemical potentials, and other variables. Our theoretical results are supported by numerical experiments and we present an example of a physical application involving the microfluidic non-ideal mixing of hydrocarbons.

We address the solution of the distributed control problem for the steady, incompressible Navier--Stokes equations. We propose an inexact Newton linearization of the optimality conditions. Upon discretization by a finite element scheme, we obtain a sequence of large symmetric linear systems of saddle-point type. We use an augmented Lagrangian-based block triangular preconditioner in combination with the flexible GMRES method at each Newton step. The preconditioner is applied inexactly via a suitable multigrid solver. Numerical experiments indicate that the resulting method appears to be fairly robust with respect to viscosity, mesh size, and the choice of regularization parameter when applied to 2D problems.

Numerical methods for the simulation of transient systems with structure-preserving properties are known to exhibit greater accuracy and physical reliability, in particular over long durations. However, there remain difficulties in devising geometric numerical integrators that preserve dissipation laws and conserve non-quadratic invariants. In this work, we propose a framework for the construction of timestepping schemes that preserve dissipation laws and conserve multiple general invariants. The framework employs finite elements in time and systematically introduces auxiliary variables; it extends to arbitrary order in time. We demonstrate the ideas by devising novel integrators that conserve (to machine precision) all known invariants of general conservative ODEs, energy-conserving finite-element discretisations of general Hamiltonian PDEs, and finite-element schemes for the compressible Navier-Stokes equations that conserve mass, momentum, and energy, and provably possess non-decreasing entropy. The approach generalises and unifies several existing ideas in the literature, including Gauss methods, the framework of Cohen & Hairer, and the energy- and helicity-conserving scheme of Rebholz., Comment: 54 pages, 13 figures

Nematic ordering describes the phenomenon where anisotropic molecules tend to locally align, like matches in a matchbox. This ordering can arise in solids (as nematic elastomers), liquids (as liquid crystals), and in gases. In the 1940s, Onsager described how nematic ordering can arise in dilute colloidal suspensions from the molecular point of view. However, the kinetic theory of nonspherical molecules has not, thus far, accounted for phenomena relating to the presence of nematic ordering. In this work we develop a kinetic theory for the behavior of rarified calamitic (rodlike) gases in the presence of nematic ordering. Building on previous work by Curtiss, we derive from kinetic theory the rate of work hypothesis that forms the starting point for Leslie--Ericksen theory. We incorporate ideas from the variational theory of nematic liquid crystals to create a moment closure that preserves the coupling between the laws of linear and angular momentum. The coupling between these laws is a key feature of our theory, in contrast to the kinetic theory proposed by {Curtiss \& Dahler}, where the couple stress tensor is assumed to be zero. This coupling allows the characterization of anisotropic phenomena arising from the nematic ordering. Furthermore, the theory leads to an energy functional that is a compressible variant of the classical Oseen--Frank energy (with a pressure-dependent Frank constant) and to a compressible analogue of the Leslie--Ericksen equations. The emergence of compressible aspects in the theory for nematic fluids enhances our understanding of these complex systems.

We present the full approximation scheme constraint decomposition (FASCD) multilevel method for solving variational inequalities (VIs). FASCD is a common extension of both the full approximation scheme (FAS) multigrid technique for nonlinear partial differential equations, due to A.~Brandt, and the constraint decomposition (CD) method introduced by X.-C.~Tai for VIs arising in optimization. We extend the CD idea by exploiting the telescoping nature of certain function space subset decompositions arising from multilevel mesh hierarchies. When a reduced-space (active set) Newton method is applied as a smoother, with work proportional to the number of unknowns on a given mesh level, FASCD V-cycles exhibit nearly mesh-independent convergence rates, and full multigrid cycles are optimal solvers. The example problems include differential operators which are symmetric linear, nonsymmetric linear, and nonlinear, in unilateral and bilateral VI problems., Comment: 25 pages, 9 figures

Liquid crystals can self-organize into a layered smectic phase. While the smectic layers are typically straight forming a lamellar pattern in bulk, external confinement may drastically distort the layers due to the boundary conditions imposed on the orientational director field. Resolving this distortion leads to complex structures with topological defects. Here, we explore the configurations adopted by two-dimensional colloidal smectics made from nearly hard rod-like particles in complex confinements, characterized by a button-like structure with two internal boundaries (inclusions): a two-holed disk and a double annulus. The topology of the confinement generates new structures which we classify in reference to previous work as generalized laminar and generalized Shubnikov states. To explore these configurations, we combine particle-resolved experiments on colloidal rods with three complementary theoretical approaches: Monte-Carlo simulation, first-principles density functional theory and phenomenological $\mathbf{Q}$-tensor modeling. This yields a consistent and comprehensive description of the structural details. In particular, we characterize a nontrivial tilt angle between the direction of the layers and symmetry axes of the confinement.

The Riesz maps of the $L^2$ de Rham complex frequently arise as subproblems in the construction of fast preconditioners for more complicated problems. In this work we present multigrid solvers for high-order finite element discretizations of these Riesz maps with the same time and space complexity as sum-factorized operator application, i.e.~with optimal complexity in polynomial degree in the context of Krylov methods. The key idea of our approach is to build new finite elements for each space in the de Rham complex with orthogonality properties in both the $L^2$- and $H(\mathrm{d})$-inner products ($\mathrm{d} \in \{\mathrm{grad}, \mathrm{curl}, \mathrm{div}\})$ on the reference hexahedron. The resulting sparsity enables the fast solution of the patch problems arising in the Pavarino, Arnold--Falk--Winther and Hiptmair space decompositions, in the separable case. In the non-separable case, the method can be applied to an auxiliary operator that is sparse by construction. With exact Cholesky factorizations of the sparse patch problems, the application complexity is optimal but the setup costs and storage are not. We overcome this with the finer Hiptmair space decomposition and the use of incomplete Cholesky factorizations imposing the sparsity pattern arising from static condensation, which applies whether static condensation is used for the solver or not. This yields multigrid relaxations with time and space complexity that are both optimal in the polynomial degree.

We present a method to embed local electroneutrality within Onsager-Stefan-Maxwell electrolytic-transport models, circumventing their formulation as differential systems with an algebraic constraint. Flux-explicit transport laws are formulated for general multicomponent electrolytes, in which the conductivity, component diffusivities, and transference numbers relate to Stefan-Maxwell coefficients through invertible matrix calculations. A construction we call a `salt-charge basis' implements Guggenheim's transformation of species electrochemical potentials into combinations describing a minimal set of neutral components, leaving a unique combination associated with electricity. Defining conjugate component concentrations and fluxes that preserve the structures of the Gibbs function and energy dissipation retains symmetric Onsager reciprocal relations. The framework reproduces Newman's constitutive laws for binary electrolytes and the Pollard-Newman laws for molten salts; we also propose laws for salt solutions in two-solvent blends, such as lithium-ion-battery electrolytes. Finally, we simulate a potentiostatic Hull cell containing a non-ideal binary electrolyte with concentration-dependent properties.

In recent years a great deal of attention has been paid to discretizations of the incompressible Stokes equations that exactly preserve the incompressibility constraint. These are of substantial interest because these discretizations are pressure-robust, i.e. the error estimates for the velocity do not depend on the error in the pressure. Similar considerations arise in nearly incompressible linear elastic solids. Conforming discretizations with this property are now well understood in two dimensions, but remain poorly understood in three dimensions. In this work we state two conjectures on this subject. The first is that the Scott-Vogelius element pair is inf-sup stable on uniform meshes for velocity degree $k \ge 4$; the best result available in the literature is for $k \ge 6$. The second is that there exists a stable space decomposition of the kernel of the divergence for $k \ge 5$. We present numerical evidence supporting our conjectures.

We develop finite element methods for coupling the steady-state Onsager--Stefan--Maxwell equations to compressible Stokes flow. These equations describe multicomponent flow at low Reynolds number, where a mixture of different chemical species within a common thermodynamic phase is transported by convection and molecular diffusion. Developing a variational formulation for discretizing these equations is challenging: the formulation must balance physical relevance of the variables and boundary data, regularity assumptions, tractability of the analysis, enforcement of thermodynamic constraints, ease of discretization, and extensibility to the transient, anisothermal, and non-ideal settings. To resolve these competing goals, we employ two augmentations: the first enforces the mass-average constraint in the Onsager--Stefan--Maxwell equations, while its dual modifies the Stokes momentum equation to enforce symmetry. Remarkably, with these augmentations we achieve a Picard linearization of symmetric saddle point type, despite the equations not possessing a Lagrangian structure. Exploiting the structure of linear irreversible thermodynamics, we prove the inf-sup condition for this linearization, and identify finite element function spaces that automatically inherit well-posedness. We verify our error estimates with a numerical example, and illustrate the application of the method to non-ideal fluids with a simulation of the microfluidic mixing of hydrocarbons.

The fourth-order PDE that models the density variation of smectic A liquid crystals presents unique challenges in its (numerical) analysis beyond more common fourth-order operators, such as the classical biharmonic. While the operator is positive definite, the equation has a "wrong-sign" shift, making it somewhat more akin to an indefinite Helmholtz operator, with lowest-energy modes consisting of plane waves. As a result, for large shifts, the natural continuity, coercivity, and inf-sup constants degrade considerably, impacting standard error estimates. In this paper, we analyze and compare three finite-element formulations for such PDEs, based on $H^2$-conforming elements, the $C^0$ interior penalty method, and a mixed finite-element formulation that explicitly introduces approximations to the gradient of the solution and a Lagrange multiplier. The conforming method is simple but is impractical to apply in three dimensions; the interior-penalty method works well in two and three dimensions but has lower-order convergence and (in preliminary experiments) seems difficult to precondition; the mixed method uses more degrees of freedom, but works well in both two and three dimensions, and is amenable to monolithic multigrid preconditioning. Our analysis reveals different behaviours of the error bounds with the shift parameter and mesh size for the different schemes. Numerical results verify the finite-element convergence for all discretizations, and illustrate the trade-offs between the three schemes.

We develop structure-preserving finite element methods for the incompressible, resistive Hall magnetohydrodynamics (MHD) equations. These equations incorporate the Hall current term in Ohm's law and provide a more appropriate description of fully ionized plasmas than the standard MHD equations on length scales close to or smaller than the ion skin depth. We introduce a stationary discrete variational formulation of Hall MHD that enforces the magnetic Gauss's law exactly (up to solver tolerances) and prove the well-posedness and convergence of a Picard linearization. For the transient problem, we present time discretizations that preserve the energy and magnetic and hybrid helicity precisely in the ideal limit for two types of boundary conditions. Additionally, we present an augmented Lagrangian preconditioning technique for both the stationary and transient cases. We confirm our findings with several numerical experiments.

Topology optimization problems generally support multiple local minima, and real-world applications are typically three-dimensional. In previous work [I. P. A. Papadopoulos, P. E. Farrell, and T. M. Surowiec, Computing multiple solutions of topology optimization problems, SIAM Journal on Scientific Computing, (2021)], the authors developed the deflated barrier method, an algorithm that can systematically compute multiple solutions of topology optimization problems. In this work we develop preconditioners for the linear systems arising in the application of this method to Stokes flow, making it practical for use in three dimensions. In particular, we develop a nested block preconditioning approach which reduces the linear systems to solving two symmetric positive-definite matrices and an augmented momentum block. An augmented Lagrangian term is used to control the innermost Schur complement and we apply a geometric multigrid method with a kernel-capturing relaxation method for the augmented momentum block. We present multiple solutions in three-dimensional examples computed using the proposed iterative solver.

Most research on preconditioners for time-dependent PDEs has focused on implicit multi-step or diagonally-implicit multi-stage temporal discretizations. In this paper, we consider monolithic multigrid preconditioners for fully-implicit multi-stage Runge-Kutta (RK) time integration methods. These temporal discretizations have very attractive accuracy and stability properties, but they couple the spatial degrees of freedom across multiple time levels, requiring the solution of very large linear systems. We extend the classical Vanka relaxation scheme to implicit RK discretizations of saddle point problems. We present numerical results for the incompressible Stokes, Navier-Stokes, and resistive magnetohydrodynamics equations, in two and three dimensions, confirming that these relaxation schemes lead to robust and scalable monolithic multigrid methods for a challenging range of incompressible fluid-flow models., Comment: 22 pages, 9 figures. Submitted to Journal of Computational Physics on Feb 14 2022, updated to match the final published version in Journal of Computational Physics fixing any typos and adding references

We introduce a numerical technique for controlling the location and stability properties of Hopf bifurcations in dynamical systems. The algorithm consists of solving an optimization problem constrained by an extended system of nonlinear partial differential equations that characterizes Hopf bifurcation points. The flexibility and robustness of the method allows us to advance or delay a Hopf bifurcation to a target value of the bifurcation parameter, as well as controlling the oscillation frequency with respect to a parameter of the system or the shape of the domain on which solutions are defined. Numerical applications are presented in systems arising from biology and fluid dynamics, such as the FitzHugh--Nagumo model, Ginzburg--Landau equation, Rayleigh--B\'enard convection problem, and Navier--Stokes equations, where the control of the location and oscillation frequency of periodic solutions is of high interest., Comment: 22 pages, 8 figures

Viscous contact problems describe the time evolution of fluid flows in contact with a surface from which they can detach and reattach. These problems are of particular importance in glaciology, where they arise in the study of grounding lines and subglacial cavities. In this work, we propose a novel numerical method for solving viscous contact problems based on a mixed formulation with Lagrange multipliers of a variational inequality involving the Stokes equation. The advection equation for evolving the geometry of the domain occupied by the fluid is then solved via a specially-built upwinding scheme, leading to a robust and accurate algorithm for viscous contact problems. We first verify the method by comparing the numerical results to analytical results obtained by a linearised method. Then, we use this numerical scheme to reconstruct friction laws for glacial sliding with cavitation. Finally, we compute the evolution of cavities from a steady state under oscillating water pressures. The results depend strongly on the location of the initial steady state along the friction law. In particular, we find that if the steady state is located on the downsloping or rate-weakening part of the friction law, the cavity evolves towards the upsloping section, indicating that the downsloping part is unstable.

The Arnold-Winther element successfully discretizes the Hellinger-Reissner variational formulation of linear elasticity; its development was one of the key early breakthroughs of the finite element exterior calculus. Despite its great utility, it is not available in standard finite element software, because its degrees of freedom are not preserved under the standard Piola push-forward. In this work we apply the novel transformation theory recently developed by Kirby [SMAI-JCM, 4:197-224, 2018] to devise the correct map for transforming the basis on a reference cell to a generic physical triangle. This enables the use of the Arnold-Winther elements, both conforming and nonconforming, in the widely-used Firedrake finite element software, composing with its advanced symbolic code generation and geometric multigrid functionality. Similar results also enable the correct transformation of the Mardal-Tai-Winther element for incompressible fluid flow. We present numerical results for both elements, verifying the correctness of our theory., Comment: Submitted to SMAI Journal of Computational Mathematics

We analyse an energy minimisation problem recently proposed for modelling smectic-A liquid crystals. The optimality conditions give a coupled nonlinear system of partial differential equations, with a second-order equation for the tensor-valued nematic order parameter $\mathbf{Q}$ and a fourth-order equation for the scalar-valued smectic density variation $u$. Our two main results are a proof of the existence of solutions to the minimisation problem, and the derivation of a priori error estimates for its discretisation of the decoupled case (i.e., $q=0$) using the $\mathcal{C}^0$ interior penalty method. More specifically, optimal rates in the $H^1$ and $L^2$ norms are obtained for $\mathbf{Q}$, while optimal rates in a mesh-dependent norm and $L^2$ norm are obtained for $u$. Numerical experiments confirm the rates of convergence.

We present the Onsager--Stefan--Maxwell thermodiffusion equations, which account for the Soret and Dufour effects in multicomponent fluids. Unlike transport laws derived from kinetic theory, this framework preserves the structure of the isothermal Stefan--Maxwell equations, separating the thermodynamic forces that drive diffusion from the force that drives heat flow. The Onsager--Stefan--Maxwell transport-coefficient matrix is symmetric, and the second law of thermodynamics imbues it with simple spectral characteristics. This new approach allows for heat to be considered as a pseudo-species and proves equivalent to both the intuitive extension of Fick's law and the generalized Stefan--Maxwell equations popularized by Bird, Stewart, and Lightfoot. A general inversion process facilitates the unique formulation of flux-explicit transport equations relative to any choice of convective reference velocity. Stefan--Maxwell diffusivities and thermal diffusion factors are tabulated for gaseous mixtures containing helium, argon, neon, krypton, and xenon. The framework is deployed to perform numerical simulations of steady three-dimensional thermodiffusion in a ternary gas., Comment: 38 pages, 4 figures

Stokes variational inequalities arise in the formulation of glaciological problems involving contact. We consider the problem of a two-dimensional marine ice sheet with a grounding line, although the analysis presented here is extendable to other contact problems in glaciology, such as that of subglacial cavitation. The analysis of this problem and its discretisation is complicated by the nonlinear rheology commonly used for modelling ice, the enforcement of a friction boundary condition given by a power law, and the presence of rigid modes in the velocity space, which render the variational inequality semicoercive. In this work, we consider a mixed formulation of this variational inequality involving a Lagrange multiplier and provide an analysis of its finite element approximation. Error estimates in the presence of rigid modes are obtained by means of a specially-built projection operator onto the subspace of rigid modes and a Korn-type inequality. These proofs rely on the fact that the subspace of rigid modes is at most one-dimensional. Numerical results are reported to validate the error estimates.

Pavarino proved that the additive Schwarz method with vertex patches and a low-order coarse space gives a $p$-robust solver for symmetric and coercive problems. However, for very high polynomial degree it is not feasible to assemble or factorize the matrices for each patch. In this work we introduce a direct solver for separable patch problems that scales to very high polynomial degree on tensor product cells. The solver constructs a tensor product basis that diagonalizes the blocks in the stiffness matrix for the internal degrees of freedom of each individual cell. As a result, the non-zero structure of the cell matrices is that of the graph connecting internal degrees of freedom to their projection onto the facets. In the new basis, the patch problem is as sparse as a low-order finite difference discretization, while having a sparser Cholesky factorization. We can thus afford to assemble and factorize the matrices for the vertex-patch problems, even for very high polynomial degree. In the non-separable case, the method can be applied as a preconditioner by approximating the problem with a separable surrogate. We demonstrate the approach by solving the Poisson equation and a $H(\mathrm{div})$-conforming interior penalty discretization of linear elasticity in three dimensions at $p = 15$.

Many problems in engineering can be understood as controlling the bifurcation structure of a given device. For example, one may wish to delay the onset of instability, or bring forward a bifurcation to enable rapid switching between states. We propose a numerical technique for controlling the bifurcation diagram of a nonlinear partial differential equation by varying the shape of the domain. Specifically, we are able to delay or advance a given branch point to a target parameter value. The algorithm consists of solving a shape optimization problem constrained by an augmented system of equations, the Moore--Spence system, that characterize the location of the branch points. Numerical experiments on the Allen--Cahn, Navier--Stokes, and hyperelasticity equations demonstrate the effectiveness of this technique in a wide range of settings., Comment: 20 pages, 11 figures

Fourth-order differential equations play an important role in many applications in science and engineering. In this paper, we present a three-field mixed finite-element formulation for fourth-order problems, with a focus on the effective treatment of the different boundary conditions that arise naturally in a variational formulation. Our formulation is based on introducing the gradient of the solution as an explicit variable, constrained using a Lagrange multiplier. The essential boundary conditions are enforced weakly, using Nitsche's method where required. As a result, the problem is rewritten as a saddle-point system, requiring analysis of the resulting finite-element discretization and the construction of optimal linear solvers. Here, we discuss the analysis of the well-posedness and accuracy of the finite-element formulation. Moreover, we develop monolithic multigrid solvers for the resulting linear systems. Two and three-dimensional numerical results are presented to demonstrate the accuracy of the discretization and efficiency of the multigrid solvers proposed.

The magnetohydrodynamics (MHD) equations are generally known to be difficult to solve numerically, due to their highly nonlinear structure and the strong coupling between the electromagnetic and hydrodynamic variables, especially for high Reynolds and coupling numbers. In this work, we present a scalable augmented Lagrangian preconditioner for a finite element discretization of the $\mathbf{B}$-$\mathbf{E}$ formulation of the incompressible viscoresistive MHD equations. For stationary problems, our solver achieves robust performance with respect to the Reynolds and coupling numbers in two dimensions and good results in three dimensions. We extend our method to fully implicit methods for time-dependent problems which we solve robustly in both two and three dimensions. Our approach relies on specialized parameter-robust multigrid methods for the hydrodynamic and electromagnetic blocks. The scheme ensures exactly divergence-free approximations of both the velocity and the magnetic field up to solver tolerances. We confirm the robustness of our solver by numerical experiments in which we consider fluid and magnetic Reynolds numbers and coupling numbers up to 10,000 for stationary problems and up to 100,000 for transient problems in two and three dimensions.

Creating scalable, high performance PDE-based simulations requires a suitable combination of discretizations, differential operators, preconditioners and solvers. The required combination changes with the application and with the available hardware, yet software development time is a severely limited resource for most scientists and engineers. Here we demonstrate that generating simulation code from a high-level Python interface provides an effective mechanism for creating high performance simulations from very few lines of user code. We demonstrate that moving from one supercomputer to another can require significant algorithmic changes to achieve scalable performance, but that the code generation approach enables these algorithmic changes to be achieved with minimal development effort.

We perform a bifurcation analysis of the steady states of Rayleigh--B\'enard convection with no-slip boundary conditions in two dimensions using a numerical method called deflated continuation. By combining this method with an initialisation strategy based on the eigenmodes of the conducting state, we are able to discover multiple solutions to this non-linear problem, including disconnected branches of the bifurcation diagram, without the need for any prior knowledge of the solutions. One of the disconnected branches we find contains an S-shaped curve with hysteresis, which is the origin of a flow pattern that may be related to the dynamics of flow reversals in the turbulent regime. Linear stability analysis is also performed to analyse the steady and unsteady regimes of the solutions in the parameter space and to characterise the type of instabilities., Comment: 19 pages, 16 figures

We study a model system with nematic and magnetic orders, within a channel geometry modelled by an interval, $[-D, D]$. The system is characterised by a tensor-valued nematic order parameter $\mathbf{Q}$ and a vector-valued magnetisation $\mathbf{M}$, and the observable states are modelled as stable critical points of an appropriately defined free energy. In particular, the full energy includes a nemato-magnetic coupling term characterised by a parameter $c$. We (i) derive $L^\infty$ bounds for $\mathbf{Q}$ and $\mathbf{M}$; (ii) prove a uniqueness result in parameter regimes defined by $c$, $D$ and material- and temperature-dependent correlation lengths; (iii) analyse order reconstruction solutions, possessing domain walls, and their stabilities as a function of $D$ and $c$ and (iv) perform numerical studies that elucidate the interplay of $c$ and $D$ for multistability., Comment: accepted in SIAM Journal on Numerical Analysis

Mathematical modelling of ionic electrodiffusion and water movement is emerging as a powerful avenue of investigation to provide new physiological insight into brain homeostasis. However, in order to provide solid answers and resolve controversies, the accuracy of the predictions is essential. Ionic electrodiffusion models typically comprise non-trivial systems of non-linear and highly coupled partial and ordinary differential equations that govern phenomena on disparate time scales. Here, we study numerical challenges related to approximating these systems. We consider a homogenized model for electrodiffusion and osmosis in brain tissue and present and evaluate different associated finite element-based splitting schemes in terms of their numerical properties, including accuracy, convergence, and computational efficiency for both idealized scenarios and for the physiologically relevant setting of cortical spreading depression (CSD). We find that the schemes display optimal convergence rates in space for problems with smooth manufactured solutions. However, the physiological CSD setting is challenging: we find that the accurate computation of CSD wave characteristics (wave speed and wave width) requires a very fine spatial and fine temporal resolution.

A phenomenological free energy model is proposed to describe the behavior of smectic liquid crystals, an intermediate phase that exhibits orientational order and layering at the molecular scale. Advantageous properties render the functional amenable to numerical simulation. The model is applied to a number of scenarios involving geometric frustration, leading to emergent structures such as focal conic domains and oily streaks and enabling detailed elucidation of the very rich energy landscapes that arise in these problems.

We devise 3-field and 4-field finite element approximations of a system describing the steady state of an incompressible heat-conducting fluid with implicit non-Newtonian rheology. We prove that the sequence of numerical approximations converges to a weak solution of the problem. We develop a block preconditioner based on augmented Lagrangian stabilisation for a discretisation based on the Scott-Vogelius finite element pair for the velocity and pressure. The preconditioner involves a specialised multigrid algorithm that makes use of a space-decomposition that captures the kernel of the divergence and non-standard intergrid transfer operators. The preconditioner exhibits robust convergence behaviour when applied to the Navier-Stokes and power-law systems, including temperature-dependent viscosity, heat conductivity and viscous dissipation., Comment: 9 figures, 39 pages

While implicit Runge--Kutta methods possess high order accuracy and important stability properties, implementation difficulties and the high expense of solving the coupled algebraic system at each time step are frequently cited as impediments. We present IIrksome, a high-level library for manipulating UFL (Unified Form Language) expressions of semidiscrete variational forms to obtain UFL expressions for the coupled Runge--Kutta stage equations at each time step. Irksome works with the Firedrake package to enable the efficient solution of the resulting coupled algebraic systems. Numerical examples confirm the efficacy of the software and our solver techniques for various problems.

We investigate structure-preserving finite element discretizations of the steady-state Stefan--Maxwell diffusion problem which governs diffusion within a phase consisting of multiple species. An approach inspired by augmented Lagrangian methods allows us to construct a symmetric positive definite augmented Onsager transport matrix, which in turn leads to an effective numerical algorithm. We prove inf-sup conditions for the continuous and discrete linearized systems and obtain error estimates for a phase consisting of an arbitrary number of species. The discretization preserves the thermodynamically fundamental Gibbs--Duhem equation to machine precision independent of mesh size. The results are illustrated with numerical examples, including an application to modelling the diffusion of oxygen, carbon dioxide, water vapour and nitrogen in the lungs., Comment: 27 pages, 5 figures

Topology optimization problems often support multiple local minima due to a lack of convexity. Typically, gradient-based techniques combined with continuation in model parameters are used to promote convergence to more optimal solutions; however, these methods can fail even in the simplest cases. In this paper, we present an algorithm to perform a systematic exploratory search for the solutions of the optimization problem via second-order methods without a good initial guess. The algorithm combines the techniques of deflation, barrier methods and primal-dual active set solvers in a novel way. We demonstrate this approach on several numerical examples, observe mesh-independence in certain cases and show that multiple distinct local minima can be recovered.

Augmented Lagrangian preconditioners have successfully yielded Reynolds-robust preconditioners for the stationary incompressible Navier-Stokes equations, but only for specific discretizations. The discretizations for which these preconditioners have been designed possess error estimates which depend on the Reynolds number, with the discretization error deteriorating as the Reynolds number is increased. In this paper we present an augmented Lagrangian preconditioner for the Scott-Vogelius discretization on barycentrically-refined meshes. This achieves both Reynolds-robust performance and Reynolds-robust error estimates. A key consideration is the design of a suitable space decomposition that captures the kernel of the grad-div term added to control the Schur complement; the same barycentric refinement that guarantees inf-sup stability also provides a local decomposition of the kernel of the divergence. The robustness of the scheme is confirmed by numerical experiments in two and three dimensions., Comment: Fixed sign of grad-div term in (3.7)

We propose a robust and efficient augmented Lagrangian-type preconditioner for solving linearizations of the Oseen-Frank model arising in cholesteric liquid crystals. By applying the augmented Lagrangian method, the Schur complement of the director block can be better approximated by the weighted mass matrix of the Lagrange multiplier, at the cost of making the augmented director block harder to solve. In order to solve the augmented director block, we develop a robust multigrid algorithm which includes an additive Schwarz relaxation that captures a pointwise version of the kernel of the semi-definite term. Furthermore, we prove that the augmented Lagrangian term improves the discrete enforcement of the unit-length constraint. Numerical experiments verify the efficiency of the algorithm and its robustness with respect to problem-related parameters (Frank constants and cholesteric pitch) and the mesh size.

We present a mesh-independent and parameter-robust multigrid solver for the Scott-Vogelius discretisation of the nearly incompressible linear elasticity equations on meshes with a macro element structure. The discretisation achieves exact representation of the limiting divergence constraint at moderate polynomial degree. Both the relaxation and multigrid transfer operators exploit the macro structure for robustness and efficiency. For the relaxation, we use the existence of local Fortin operators on each macro cell to construct a local space decomposition with parameter-robust convergence. For the transfer, we construct a robust prolongation operator by performing small local solves over each coarse macro cell. The necessity of both components of the algorithm is confirmed by numerical experiments.

Effective relaxation methods are necessary for good multigrid convergence. For many equations, standard Jacobi and Gau{\ss}-Seidel are inadequate, and more sophisticated space decompositions are required; examples include problems with semidefinite terms or saddle point structure. In this paper we present a unifying software abstraction, PCPATCH, for the topological construction of space decompositions for multigrid relaxation methods. Space decompositions are specified by collecting topological entities in a mesh (such as all vertices or faces) and applying a construction rule (such as taking all degrees of freedom in the cells around each entity). The software is implemented in PETSc and facilitates the elegant expression of a wide range of schemes merely by varying solver options at runtime. In turn, this allows for the very rapid development of fast solvers for difficult problems., Comment: 22 pages, minor fixes in bibliography

Multigrid methods are popular solution algorithms for many discretized PDEs, either as standalone iterative solvers or as preconditioners, due to their high efficiency. However, the choice and optimization of multigrid components such as relaxation schemes and grid-transfer operators is crucial to the design of optimally efficient algorithms. It is well--known that local Fourier analysis (LFA) is a useful tool to predict and analyze the performance of these components. In this paper, we develop a local Fourier analysis of monolithic multigrid methods based on additive Vanka relaxation schemes for mixed finite-element discretizations of the Stokes equations. The analysis offers insight into the choice of "patches" for the Vanka relaxation, revealing that smaller patches offer more effective convergence per floating point operation. Parameters that minimize the two-grid convergence factor are proposed and numerical experiments are presented to validate the LFA predictions., Comment: 30 pages, 12 figures. Add new sections: multiplicative Vanka results and sensitivity of convergence factors to mesh distortion

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., Comment: 24 pages, 3 figures

In the classical theory of fluid mechanics a linear relationship between the shear stress and the symmetric velocity gradient tensor is often assumed. Even when a nonlinear relationship is assumed, it is typically formulated in terms of an explicit relation. Implicit constitutive models provide a theoretical framework that generalises this, allowing for general implicit constitutive relations. Since it is generally not possible to solve explicitly for the shear stress in the constitutive relation, a natural approach is to include the shear stress as a fundamental unknown in the formulation of the problem. In this work we present a mixed formulation with this feature, discuss its solvability and approximation using mixed finite element methods, and explore the convergence of the numerical approximations to a weak solution of the model., Comment: 4 figures. To appear in the SIAM Journal on Numerical Analysis