Réunion d'été SMC 2023

Ottawa, 2 - 5 juin 2023


Méthodes numériques pour les équations différentielles partielles
Org: Yves Bourgault et Diane Guignard (University of Ottawa)

FRANCIS AZNARAN, University of Oxford
Transformations for Piola-mapped elements  [PDF]

Many finite element spaces are not preserved by the standard pullback to the reference cell. Robust implementation therefore requires studying the relation between degrees of freedom under pushforward, in order to obtain the correct bases on a generic physical triangle [Kirby, 2018]. In this work, we extend this transformation theory to vector- and tensor-valued elements mapped by the contravariant Piola transform. We apply this theory, and describe its efficient implementation in Firedrake, for the the Mardal–Tai–Winther elements discretising $H(\mathrm{div})$ for Stokes–Darcy flow, and the conforming and nonconforming Arnold--Winther elements discretising $H(\mathrm{div};\mathbb{S})$ for the stress-displacement formulation of linear elasticity.

In particular, the Arnold–Winther elements were the first to stably enforce exact symmetry of the Cauchy stress tensor; we demonstrate how they may be efficiently mapped, while the few prior implementations are either custom-made for specific numerical experiments, or require the explicit element-by-element construction of the basis. Our novel implementation of these exotic elements composes inexpensively and automatically with the rest of the Firedrake code stack. We also demonstrate the effectiveness of appropriate multigrid smoothers for this system, prove convergence of Nitsche's method for the weak enforcement of traction conditions, and provide a uniform construction of all standard reference-to-physical Piola pullback maps using the finite element exterior calculus.

BAPTISTE BERLIOUX, Polytechnique Montréal
A-stable and high order nonlinear time integration methods based on deferred correction schemes  [PDF]

The starting point of our methods is the A-stable BDF (Backward Differentiation Formulae) schemes of order 1 and 2. Our methods consist in gradually increasing the order of these schemes from deferred corrections, while keeping the unconditional stability. The dominant terms of the local truncation error are substituted from the initial expression of the BDF1 and BDF2 methods. The substituted term is then the correction, interpolated by a Newton polynomial from the solution of the lower order method, supposed to be the most accurate. For example, the DC2/BDF1 method substitutes the term $-\ddot{u}(t^n)\frac{k_n}{2!}$ (first dominant term of the truncation error of the BDF1 scheme) in the expression of the BDF1 scheme. Concerning $\ddot{u}(t^n)$, it is interpolated from the numerical solution $u_1^n$, resulting from the BDF1 solution at time $t^n$, in order to make it discrete. By construction of the correction, the DC2/BDF1 method is then of order 2.

By applying the same strategy, we gradually increase the order of the methods and can expect to reach very high orders (DC6/BDF, DC7/BDF, etc.).

But, because of the interdependence of the methods, the restriction of the order allows limiting the complexity of implementation as well as the computational cost. Thus, in order to limit this one, we aim to build an algorithm with adaptive time steps in order to minimize the number of time steps. Before that, it is necessary to verify the numerical behaviour of the methods in different time step configurations.

JEAN DETEIX, Université Laval
A projection scheme for the Navier–Stokes/Allan–Cahn model  [PDF]

Projection methods are known for their efficicency for Newtonian fluids. However their use in the context of variable density or viscosity is problematic (but not impossible). We propose in this work a time-discrete formulation of the coupled Navier– Stokes/Allen–Cahn equations based on a projection method. The scheme is based on two ingredients: a projection method for heterogeneous fluid and the concept of coupled projection scheme.

We first establish the well posedness and stability of the time-discrete formulation. Next we propose an iterative schemes for the actual approximation of solutions.In the last part of our presentation, using finite element for spatial discretization, we estimate the order of accuracy in time and illustrate the validity of the scheme through numerical results.

QIWEI FENG, University of Alberta
High Order Finite Difference Methods for Interface Problems  [PDF]

Interface problem arises in many applications such as modeling of underground waste disposal, oil reservoir, composite material, and many others. The coefficient $a$, the source term $f$, the wave number $\textsf{k}$, the solution $u$, and the flux $a\nabla u\cdot \vec{n}$ are possibly discontinuous across the interface curve $\Gamma$ in such problems. To obtain the reasonable numerical solution, the higher order numerical scheme is desirable. Firstly, we propose a sixth order compact 9-point finite difference method (FDM) for the Poisson interface problem with the singular source. For the elliptic interface problem with the discontinuous and piecewise smooth coefficient, we propose a high order compact 9-point FDM and a high order local calculation for the approximation of the solution $u$ and the gradient $\nabla u$ respectively. Furthermore, we derive the compact 9-point FDM with high accuracy order and/or M-matrix property for the elliptic cross-interface problem. Finally, we propose a sixth order compact FDM for the Helmholtz Equation with the singular source. To reduce the pollution effect, we propose a new pollution minimization strategy that is based on the average truncation error of plane waves. All above schemes are developed on the uniform Cartesian grid in a rectangular domain. Our numerical experiments confirm the flexibility and the expected order accuracy in $l_2$ and $l_{\infty}$ norms of the proposed schemes. Except the Helmholtz Equation, we prove the corresponding convergence rate for the proposed schemes using the discrete maximum principle. This is joint work with Bin Han, Peter Minev and Michelle Michelle.

ANDRÉ FORTIN, Université Laval
A discontinuous Galerkin method for stiff ODEs and DDEs  [PDF]

We present a very high order implicit discontinuous Galerkin method for the solution of stiff ordinary differential equations (ODEs) and delay differential equations (DDEs). The proposed method is based on Legendre orthogonal polynomials of degree $k$ and is shown to converge at order $k+1$ in $L^2$ and sup norms. We show how the error can be estimated allowing a very efficient control of the time step. We also propose a breaking point detection algorithm for DDEs. We will present numerous examples (with $k=10$ and thus converging at order 11) to illustrate the efficiency of the method: stiff ODES with discontinuous right-hand sides, stiff DDEs with discontinuous solutions (breaking points), neutral DDEs, DDEs where classical solutions cease to exist, etc.

Error Assessment for a Finite Element - Neural Network Approach Applied to Parametric PDEs  [PDF]

A parametric PDE \begin{align}\label{pde} \mathcal{F}(u(x;\mu);\mu)= 0 \ \ &\ x \in \Omega,\ \ \mu \in \mathcal{P}, \end{align} where $\Omega$ denotes the physical domain and $\mathcal{P}$ the parameters domain, is considered. Depending on the differential operator $\mathcal{F}$, numerical methods used to approximate the solution of \eqref{pde} may be time consuming, which is an issue in the many query context (e.g. if an inverse problem has to be solved) or if the solution is needed in real time. To overcome this issue, we aim to build a neural network that approximates the map $(x;\mu) \mapsto u(x;\mu)$. For this purpose, we compute numerical approximations $u_h(x;\mu_i)$ of $u(x;\mu_i)$, $i=1,\ldots,M$, and use them as a train set for the neural network. We are then interested in the error between $u$ and the approximation given by the network, that we denote by $u_\mathcal{N}$. More precisely, we want to estimate $||u-u_{\mathcal{N}}||_{L^2(\Omega \times \mathcal{P})}$, which can be split as \begin{align*} ||u-u_{\mathcal{N}}||_{L^2(\Omega \times \mathcal{P})}\leq ||u-u_{h}||_{L^2(\Omega \times \mathcal{P})}+||u_h-u_{\mathcal{N}}||_{L^2(\Omega \times \mathcal{P})}. \end{align*} In the presentation, we discuss how the two error terms can be estimated, to what extend we can ensure that they are balanced and we present numerical results for a model problem. If time permits, we will also discuss how the neural network can be used to solve an inverse problem.

Neural network-based discontinuity tracking for hyperbolic conservation laws  [PDF]

We develop neural network-based algorithms for accurately solving weak solutions to hyperbolic conservation laws. The principle is to compute the solution in space-time subdomains defined by the curves of discontinuity, constructed from the Rankine-Hugoniot jump conditions. The proposed approach allows to efficiently consider an arbitrary number of entropic shock waves, shock wave generation, as well as wave interactions. Some numerical experiments are presented to illustrate the strengths of the algorithms.

LUIS MORA, University of Waterloo
On the Strictly Uniform Exponential Decay of a Mixed-FEM Discretization for the Wave Equation with Boundary Dissipation  [PDF]

Many approximation methods, such as standard finite elements, fail to preserve the decay rate of dissipative wave problems. Strictly preservation of the exponential stability by a first-order mixed finite element approximation method is shown for the one-dimensional wave equation with a partially reflective boundary. We use the multiplier method to analyze the continuous system's exponential stability, expressing the exponential decay rate and amplitude as functions of the physical parameters and boundary dissipation gain. An equivalent analysis is applied to prove that the energy of the approximated model is exponentially stable, and also to provide a similar bound in terms of the physical parameters.

PARIDE PASSELI, Ecole Polytechnique Fédérale de Lausanne
Anisotropic Adaptive Finite Elements for a p-Laplace Like Problem. An Application to Aluminium Electrolysis  [PDF]

Simulating Aluminium Electrolysis is a complex multi-scale and multi-physics task. Here we are interested in the fluid-flow problem describing the movements of liquid aluminium and electrolytic bath. In order to find a trade off between computational time and accuracy, Adaptive Finite Elements with large aspect ratio are considered. As a test case the p-Laplacian like problem $-\nabla\cdot((\mu+|\nabla u|^{p-2})\nabla u)=f$ , when $\mu\ge 0$ and $p\ge 2$ is considered. Using the anisotropic setting of [1, 2] and the quasi-norm techniques in [3], an anisotropic a posteriori error estimate is proved. A mesh adaptive strategy is presented. Numerical experiments show the sharpness of the estimator on both fixed and adapted meshes. Finally the developed strategy is applied to Aluminium Electrolysis.


[1] L. Formaggia, and S. Perotto. New anisotropic a priori error estimates. Numer. Math., Vol. 89(4), pp. 641-667, 2001.

[2] L. Formaggia, and S. Perotto. Anisotropic error estimates for elliptic problems. Numer. Math., Vol. 94(1), pp. 67-92, 2003.

[3] W. Liu, and N. Yan. Quasi-norm local error estimators for p-laplacian. SIAM J. NUMER. ANAL., Vol. 30(1), pp. 100-127, 2001.

SANDER RHEBERGEN, University of Waterloo
Space-time HDG for the advection-diffusion equation on time-dependent domains in the limit of small diffusion  [PDF]

This work is in collaboration with Yuan Wang. The time-dependent advection-diffusion equation on a time-dependent domain $D(t) \subset R^d$ is given by:

\begin{equation*} \partial_t \theta + {\bf u} \cdot \nabla \theta - \nu \nabla^2 \theta = f, \qquad {\bf x} \in D(t),\ t \in (0,T]. \qquad (1) \end{equation*}

Here $\theta$ is the quantity of interest, ${\bf u}$ is a flow field, and $\nu>0$ a diffusion parameter.

The space-time framework facilitates the discretization of PDEs on moving domains: a time-dependent PDE is re-formulated as a ``stationary'' PDE in $(d+1)$ space-time which is then discretized by a finite element method.

In [1] we introduced the space-time hybridizable discontinuous Galerkin method for (1). In [2] we analyzed this discretization assuming a sufficiently large diffusion parameter $\nu$. In this talk I will present a new analysis of the space-time HDG method focusing on the small diffusion limit ($\nu \ll 1$).

[1] S. Rhebergen and B. Cockburn, Space-time hybridizable discontinuous Galerkin method for the advection-diffusion equation on moving and deforming meshes, in The Courant-Friedrichs-Lewy (CFL) condition, 80 years after its discovery, ed. C.A. de Moura and C.S. Kubrusly, pp. 45-63 (2013).

[2] K.L.A. Kirk, T. Horvath, A. Cesmelioglu and S. Rhebergen, Analysis of a space-time hybridizable discontinuous Galerkin method for the advection-diffusion problem on time-dependent domains, SIAM J. Numer. Anal., 57, 4, pp. 1677-1696 (2019).

SETH TAYLOR, McGill University
A characteristic mapping method for incompressible hydrodynamics on a rotating sphere  [PDF]

We present a semi-Lagrangian method for diffeomorphism approximation and its application to incompressible hydrodynamics on the sphere. The method approximates the flow of a velocity field using a spatio-temporal discretization formed by a composition of submaps. This technique substitutes the effects of spatial refinement with the operation of composition by adaptively growing the temporal discretization. In turn, the method has the capacity of accurately and sparsely representing the generation of fine scales globally using only a linear increase in the degrees of freedom. The design and analysis of the method is presented and supported through numerical experimentation on some canonical geophysical flows.

© Société mathématique du Canada : http://www.smc.math.ca/