ABSTRACT. Model order reduction methods can reduce the runtime of forward model solves by many orders of magnitude. This has important applications in creating simulation tools with real-time response rates, accelerating the solution of stochastic PDEs and inverse problems.

A popular and straightforward way of reducing linear PDEs is the method of proper orthogonal decomposition (POD) [1]. This leads to a new set of global POD basis functions that optimally represent the solution data in a set of generated solution snapshots.

In an online phase bilinear and linear operators on the original finite element space are projected onto the pre-calculated POD basis, resulting in a small dense linear system that is fast to solve [1]. However, as a natural consequence of solving non-linear PDEs using Newton’s method, the bilinear and linear forms on the finite element space become functions of the current solution. Therefore at every Newton iteration we must re-assemble the finite element space operators and re-project onto the POD basis. This operation can become a dominant portion of the runtime of the solution of the reduced order model.

To overcome this issue we have developed a new method where we construct a reduced integration mesh that is adapted to the POD basis and the data contained in the snapshots. We generate this mesh via a greedy process in the offline phase using a simple error estimation and adaptive mesh refinement scheme on the sequence of POD basis functions. The reduced integration mesh has far fewer cells than the original mesh and therefore the linearised projected operators required at each Newton iteration can be assembled much faster.

We will show some results implemented in DOLFIN [2] using SLEPc [3] that show an order of magnitude speed up in forward runtime on a time-dependent two-dimensional non-linear reaction diffusion equation and a three-dimensional geometrically non-linear hyperelastic model.

ABSTRACT. Motivated by finite element spaces used for representation of temperature in the compatible finite element approach for numerical weather prediction, we introduce locally bounded transport schemes for (partially-)continuous finite element spaces. The underlying high-order transport scheme is constructed by injecting the partially-continuous field into an embedding discontinuous finite element space, applying a stable upwind discontinuous Galerkin (DG) scheme, and projecting back into the partially-continuous space; we call this an embedded DG scheme. We prove that this scheme is stable in L2 provided that the underlying upwind DG scheme is. We then provide a framework for applying limiters for embedded DG transport schemes. Standard DG limiters are applied during the underlying DG scheme. We introduce a new localised form of element-based flux-correction which we apply to limiting the projection back into the partially-continuous space, so that the whole transport scheme is bounded. We provide details in the specific case of tensor-product finite element spaces on wedge elements that are discontinuous P1/Q1 in the horizontal and continuous P2 in the vertical.

The framework is illustrated with numerical tests, implemented in Firedrake. The slope limiters are implemented with hand-tooled element kernels, since they can't be expressed in UFL. I will show how we make use of the Firedrake "escape hatch" to write and execute the element kernels.

ABSTRACT. Computing the solution sets of nonlinear equations is a fundamental task in applied mathematics. We consider numerical methods for computing the solutions of

f (u, λ) = 0

where f : U × R → Y is the C^1 problem residual, U and Y are isomorphic Banach spaces, u ∈ U is referred to as the solution, and λ ∈ R is referred to as the parameter. The main algorithm currently used to map out the solution sets of (1) as λ is varied is the combination of arclength continuation and branch switching. This algorithm has been enormously successful and is used in every quantitative branch of human endeavour. However, it suffers from three significant disadvantages:

• It demands the solution of expensive subproblems. In particular, it demands the computation of determinants of the Jacobian of f , to detect the presence of bifurcation points, and the nullspace of the Jacobian at those bifurcation points to switch branches. This limits  the scalability of the algorithm and restricts its applicability to crude discretisations of PDEs.

• It cannot compute disconnected bifurcation diagrams: it only attempts to compute that part of the diagram connected to the initial data. Unfortunately, bifurcation diagrams are generically disconnected.

• In fact, it can also miss connected branches associated with nonsimple bifurcation points (such as when an eigenvalue of even multiplicity crosses the origin).

In recent work we have developed an algorithm, called deflated continuation, that overcomes all three disadvantages. The algorithm only demands the solution of the Jacobian of the underlying problem, and thus scales as far as its preconditioner. Furthermore, deflated continuation can detect branches without regard to whether or how they are connected to known data.

In this talk, I will present the algorithm and its implementation in FEniCS, and give examples of its utility in discovering previously unknown solutions in varied applications such as nematic liquid crystals and Bose–Einstein condensates.

ABSTRACT. In this talk we describe a novel approximation method for the computation of solutions of high-dimensional parametric PDEs. It combines a multi-level approach, with Petrov-Galerkin discretizations approximated by means of weighted compressed sensing. We propose to decompose the solution in a hierarchy of nested finite dimensional subspaces, and approximate the details based on a discretization level dependent random sampling of the parameter space. The discretization at a certain level is then computed by means of weighted compressed sensing.

We thus introduced an approximation scheme which is easy to implement by black-box tools such as FEniCS and any $\ell_1$ solvers. Its convergence and optimality are analyzed. Furthermore, our theoretical findings are illustrated by numerical simulations.

ABSTRACT. The FEniCS suite of tools has proven extremely useful for computing the finite element approximation to a large range of physical systems. Furthermore, the automated code generation and symbolic representation of inherent variational formulations leads to the ability to rapidly prototype exotic finite element methods. A subset of these discretisation schemes is the discontinuous Galerkin finite element method. Sadly, the DG FE formulation is often verbose and complicated, even in a linear setting.

We present a method for automatically generating the DG FE formulation of nonlinear hyperbolic and elliptic operators. The result is not only a large reduction in the complexity of code written by the user, but the underlying mathematics itself.

The automated generation of the DG FE formulation is then applied to the compressible Euler and compressible Navier-Stokes equations, where we seek steady state solutions for super-sonic flow.

ABSTRACT. The FEniCS environment offers the unique possibility to automatically generate solvers for PDEs on manifolds. For this reason, we discuss how to semi-automatically generate a solver family for geometric inverse optimization problems using FEniCS and, where possible, Dolfin-Adjoint. In particular, we are interested in reconstructing shapes of objects by measuring reflections of waves. We solve geometry reconstruction problems governed by Maxwell's equations and acoustics of up to 3.5 billion PDE state unknowns in parallel using only intrinsic FEniCS functionality. This also includes grid deformation and repair in distributed memory parallelism.

ABSTRACT. We consider a singular problem $Au = L$ and two common ways of obtaining its numerical solution: (i) Krylov subspace iterations with a near nullspace and (ii) saddle point formulation. It is known that for (i) the right hand side has to be orthogonalized with respect to the nullspace of $A$. But how should the orthogonality be measured? Should one orthogonalize $L$ or is modifying the discrete right hand side sufficient? Finally, do such distinctions matter? We shall answer these questions. In particular, using Poisson and linear elasticity equations we show that orthogonalizing the discretized continuous problem, an approach currently taken by \texttt{FEniCS}, can lead to convergence to wrong solutions. We demonstrate that correct solutions are obtained by discretizing the orthogonalized continuous problem. In case of linear elasticity the process requires an orthonormal basis of the space of rigid motions which we shall provide by an elegant construction based on physical principles. Finally we show usefulness of the constructed basis in preconditioning (ii).

ABSTRACT. The G2 General Galerkin/Direct FEM methodology and Unicorn/FEniCS realization has contributed significant new directions and results in predicting gross quantities of turbulent flow [5, 2, 4]. The question of the stability of the adjoint problem is still open and debated [8], also in the FEniCS community. In this contribution we provide an updated approach to the implementation of Unicorn, denoted Unicorn minimal, where the stabilization, continuous adjoint method and do-nothing adaptivity are all fully automated in the FEniCS-HPC branch [1]. Based on this development we carry out a detailed verification of the stability of the adjoint problem for 3D incompressible Euler flow past a cube, where we demonstrate stability in time and under adaptive mesh refinement [3]. In Unicorn-minimal we also formulate a Direct FEM variable-density variant, which can be interpreted as an implicit level-set formulation [7], where we present good benchmark results for the standard 3D MARIN dam break benchmark [6]. The implicit level-set Direct FEM approach is also extended to fixed-mesh FSI in a prototype setting.

ABSTRACT. Hybrid particle-mesh methods attempt to combine the advantages of Eulerian and Lagrangian methods: Lagrangian particles are used for the advection, whereas a Eulerian background grid is used for computing the particle interactions. Such a hybrid approach is expected to have several benefits when simulating flows involving free surfaces or material interfaces: the particles can be efficiently used to track the free surface, while the background grid can be used to solve the governing Navier-Stokes equations and impose the incompressibility constraint in a convenient manner. The prospects of such a hybrid approach were assessed by implementing the four main steps common to all hybrid particle-mesh methods: mapping the particle-properties to the background grid, solving the background grid equations, updating and advecting the particles. While particle libraries take care of the particle-to-grid mapping, the grid-to-particle mapping and the advection step, FEniCS is used for solving the governing equations at the background grid. For this, a velocity-pressure formulation of the Navier-Stokes equation in Lagrangian form is solved using the Chorin-Temam projection scheme. Although not restricted to a specific element choice, a simple admissible P0P1 element was chosen in this work . Various test cases, both for single-phase and two-phase problems, show the advantages and the disadvantages of the developed model. Model results are in good agreement with analytical results or results obtained with established numerical methods, although a few persisent, particle related issues are encountered such as the particle distribution which is slowly degrading over time. However, the hybrid approach is shown to have some attractive features: the sharp interface between materials (air-water) is for example well maintained over time by the particles. Another distinct advantage of the hybrid particle-mesh method is the implementation of kinematic boundary conditions. Since the mesh can be redefined every timestep without additional interpolation steps, the method can deal with large boundary displacements as shown for a solitary wave generated with a moving boundary. Finally, using the hybrid particle-mesh formulation, interfaces of complex topological shape are conveniently tracked by the particles. As such, the method is shown to be able to simulate a breaking wave on a submerged bar up to and including wave breaking.

ABSTRACT. The float glass process is the standard industrial scale process for manufacturing flat glass, see Pilkington [1]. It is a practical example of a multicomponent system composed of molten glass, molten tin and nitrogen. As a mathematical model describing the flow of such immiscible chemically noninteracting substances we use a Cahn-Hilliard-Navier-Stokes type model proposed by Boyer and Lapuerta [2], which conceptually belongs to the class of so-called diffuse interface models. The main challenge in the computations is the very fine spatial resolution required for capturing the dynamics of the interface between the components. This places high demands on hardware as well as on numerical methods.

The most challenging part of the solution process is the efficient resolution of the incompressible Navier-Stokes system with variable density and viscosity. Several numerical approaches to the problem have been suggested, but only few of them have been tested beyond academic examples. Building on the FEniCS Project [3] we have redesigned and implemented some of the most promising methods using Krylov subspace solvers with a field split based "pressure convection-diffusion" (PCD) preconditioner, see Elman et al. [4]. The implementation uses the interface provided by DOLFIN's PETScKrylovSolver to set up PCD operators. The action of the preconditioner is implemented using subroutines provided by PETSc and petsc4py. We will discuss the efficiency of the proposed methods in the large scale computing setting.

References:

[1] L.A.B. Pilkington: Review lecture. The float glass process, Proc. R. Soc. A-Math. Phys. Eng. Sci., 314(1516) (1969), 1–25.

[2] F. Boyer, C. Lapuerta: Study of a three component Cahn-Hilliard flow model, ESAIM: Mathematical Modelling and Numerical Analysis, 40 (2006, 7), 653–687.

[3] M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, G. N. Wells: The FEniCS Project Version 1.5, Archive of Numerical Software, 3(100) (2015).

[4] H.C. Elman, D.J. Silvester, A.J. Wathen, Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics. Numerical Mathematics and Scientific Computation, Oxford University Press, 2014.

ABSTRACT. In a resonator cavity, the electromagnetic field is governed by the vector Helmholtz equation. For example, the electric foeld E can be taken into account. The problem resembles an extended eigenvalue problem Ax−λBx = 0. Normally, one would choose to solve for λ = ω^2_res. In this work, ω_res is known from a measurement (k^2_0). The eigenvalue corresponds to the permittivity of the material under test (MUT) of the resonator, as previously demonstrated using an FDFD approach. Using SLEPcEigenSolver() the real-valued eps_r of the MUT is calculated. Exploiting the cylindrical symmetry of the problem greatly enhances precision and allows for real-time use. The code shown solely using the E_phi component, returns results for TM0xx modes. TE/TMmxx modes (m > 0) are calculated using an extended code, based on [2]. Losses are calculated via a perturbation approach.

ABSTRACT. \title{FEM-BEM Coupling, Maxwell's Equations, and BEM++}

\author{\textbf{Matthew Scroggs}, University College London, matthew.scroggs.14@ucl.ac.uk, \\ Timo Betcke, University College London, t.betcke@ucl.ac.uk, \\ Erik Burman, University College London, e.burman@ucl.ac.uk}

\maketitle\thispagestyle{empty}

Keywords: \emph{Boundary Element Method, BEM++, Coupling, Maxwell}\\

For homogeneous exterior problems, the boundary element method is a useful tool for solving PDEs, as it only requires discretizing over the finite boundary of the domain. For problems containing small, localized inhomogeneities, the boundary and finite element methods can be coupled to efficiently find the solutions.

BEM++ (http://www.bempp.org) is an open source boundary element method library whose development is based at UCL. It provides a simple Python interface based on a fast C++ computational kernel. At FEniCS '15, I spoke about using BEM++ with FEniCS for Helmholtz FEM-BEM coupled problems. Recently, the development of BEM++ has centred around Maxwell's equations.

In this talk, I will present some of the newest features of BEM++. These will include a stabilized version of Helmholtz FEM-BEM coupling \cite{hiptmeury}, multiplicative Calder\'on preconditioners for Maxwell electric field problems \cite{precon}, and progress towards implementing FEM-BEM coupling for Maxwell problems, again using FEniCS. I hope to present our first working Maxwell FEM-BEM examples.

\begin{thebibliography}{00} \addcontentsline{toc}{chapter}{References} \bibitem{precon} F.P. Andriulli, K. Cools, H. Bagci, F. Olyslager, A. Buffa, S. Christiansen, E. Michielssen: A multiplicative Calderon preconditioner for the electric field integral equation, \textit{IEEE Transactions on antennas and propagation}, 56(8,1), 2008, 2398-2412. \bibitem{buffa} A. Buffa, R. Hiptmair: Galerkin Boundary Element Methods for Electromagnetic Scattering, in \textit{Computational Methods in Wave Propagation} (M. Ainsworth), Springer, 2003. \bibitem{hiptmeury} R. Hiptmair, P. Meury: Stabilized FEM-BEM coupling for Helmholtz transmission problems, SIAM journal on numerical analysis, 44(5), 2006, 2107-2130. \bibitem{johned} C. Johnson, J.C. N\'ed\'elec: On the coupling of boudary integral and finite-element methods, \textit{Mathematics of computation}, 35(152), 1980, 1063-1079. \bibitem{bempp} W. \'Smigaj, T. Betcke, S. Arridge, J. Phillips, M. Schweiger: Solving boundary integral problems with BEM++, \textit{ACM Transactions on mathematical software}, 41(2), 2015, 6:1-6:40. \end{thebibliography}

ABSTRACT. Archlength continuation is a classical numerical approach [3] to study parameter dependent nonlinear PDEs. It establishes the equilibrium path of steady state solutions as well it detects bifurcation points. Combination of the flexibility of FEniCS with the rigorously justified prin- ciples of predictor-corrector techniques provides a very powerful method. As all linear algebra operations are performed by PETSc, the developed code is completely parallel. Furthermore, the discontinuous Galerkin elements in FEniCS enable one to efficiently handle problems with derivatives higher order than two (instead of the classical Hermite elements). These features together imply that a promising target of the new method is simulating nonlinear problems with the biharmonic operator included. As an application we present a widely discussed problem from the solid mechanics literature: wrinkling of ultrathin films [2], [1], [4]. We demonstrate, that on the one hand the results obtained by the new method is in perfect agreement with the findings of previous studies [4], but on the other hand in efficiency and flexibility it outperforms them much. The research was supported by the Janos Bolyai Research Scholarship of the Hungarian Academy of Sciences [SA]

ABSTRACT. We present multimesh methods for finite element simulations of flow and computation of view for geometries described as collections of overlapping, non-matching meshes, with each mesh describing a separate object embedded into a common fixed background mesh. The computation of flow is based on a cut finite element method for the Stokes equations, and the computation of view is based on a new proposed measure of the view from a given location within the geometry described by the overlapping meshes. We consider the application to the design of settlement layouts, where the challenge is to position a number of objects (houses) in a background mesh (a landscape), so as to achieve good wind conditions and views for all houses in the settlement. The implementation is based on the latest multimesh functionality in FEniCS. A number of improvements and extensions to FEniCS have been developed for the current project. A code snippet for using multimesh in FEniCS can be seen below.

ABSTRACT. The Einstein-Vlasov system describes a large collection of collissionless particles interacting via the mean gravitational field, where gravity is modeled by general relativity. Solutions to this system provide models in astrophysics and are useful in understanding gravitational collapse. We present numerical solutions of these equations which are far-from spherically symmetric in the sense that the particle distributions take flattened and toroidal shapes (for example Figure~\ref{fig1}), and the solutions have non-zero net angular momentum. In addition, certain families of solutions are found to contain ergoregions \cite{Ames:2016vu}.

The solutions are obtained by making use of an ansatz, after which the Einstein-Vlasov system can be written as a system of semi-linear equations of the form \begin{equation} \Delta u - \frac C\rho \nabla \rho \cdot \nabla u = \Phi(u), \end{equation} where $\rho$ is the radial coordinate, $\Delta, \nabla$ refer to the Euclidean operators in $(\rho, z)$, $C$ is a real number, and $\Phi(u)$ is a nonlinear integral function of the unknown fields. This system of equations is solved numerically using an iteration, in which, at each step the total mass of the configuration is fixed and a linear system of PDEs is solved. FEniCS is used in solving the linear PDEs. The desired solutions have certain decay at spatial infinity, and so it is necessary to take a large computational domain in order to obtain accurate solutions. Mesh refinement thus becomes important in order to reduce computational costs. This talk will include a discussion of the properties of the solutions obtained as well as the numerical methods.