We present the Onsager--Stefan--Maxwell thermodiffusion equations, which account for the Soret and Dufour effects in multicomponent fluids by treating heat as a pseudo-component. Unlike transport laws derived from kinetic theory, this framework prese
rves 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 proves equivalent to both the intuitive extension of Ficks law and the generalized Stefan--Maxwell equations popularized by Bird, Stewart, and Lightfoot. A general inversion process allows 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.
Stokes variational inequalities arise in the formulation of glaciological problems involving contact. Two important examples of such problems are that of the grounding line of a marine ice sheet and the evolution of a subglacial cavity. In general, r
igid modes are present in the velocity space, rendering 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 novel technique involving metric projections onto closed convex cones. Numerical results are reported to validate the error estimates and demonstrate the advantages of using a mixed formulation in a glaciological application.
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.
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 propos
e 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 bifurcation point to a given parameter value, often to within machine precision. 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 bifurcation points. Numerical experiments on the Allen--Cahn, Navier--Stokes, and hyperelasticity equations demonstrate the effectiveness of this technique in a wide range of settings.
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 Nitsches 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 an
d 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 state solutions of Rayleigh--Benard convection with no-slip boundary conditions in two dimensions using a numerical method called deflated continuation. By combining this method with an initialisation s
trategy 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 of any prior knowledge of the dynamics. One of the disconnected branches we find contains a s-shape bifurcation with hysteresis, which is the origin of the 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.
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.
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, t
he 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.
