Reliable Methods of Mathematical Modeling

1 Introduction

The development and analysis of numerical methods for partial differential equations (PDEs) is under constant progression, driven, on the one hand, by mathematical curiosity and, on the other hand, by the practical needs of the computational sciences and engineering. For the latter, the role of mathematics is to provide reliable tools to the engineers which prevent failure of numerical simulation.

For these reasons, the biennial RMMM conferences bring together scientists developing reliable methods for mathematical modeling. The topics of the conference include applied and numerical analysis, methods for the control of modeling and numerical errors, algorithmic aspects, challenging applications, and novel discretization methods for the numerical approximation of PDEs.

The 9th edition of the conference took place on September 9–13, 2019 at TU Wien, Vienna, Austria. This special issue collects selected works from participants of RMMM 2019 that are related to their presentations. The overall focus is as wide as the needs for mathematics in computational PDEs, addressing a posteriori error control [2, 5] and adaptivity [3, 8, 7, 1], reliable methods for space-time problems [4, 12, 3], non-standard numerical discretizations [2, 4, 9, 6], and iterative solvers and optimal preconditioning [8, 7, 10, 11].

2 Non-Standard Discretizations

The work [2] of Carstensen and Nataraj gives an overview on (optimal) a priori and a posteriori error analysis for the nonconforming Crouzeix–Raviart and Morley FEM, including counterexamples on the role of the data oscillations in best approximation estimates. As model problems serve the 2D Poisson problem as well as the 2D biharmonic problem, which are treated in one mathematical framework. The key argument in the analysis is a (conforming) companion operator being defined as a right-inverse of the (nonconforming) interpolation operator with additional benefits like L 2 -orthogonality. In particular, the authors advertise the use of the companion operator as a paradigm shift in the numerical analysis of nonconforming methods to circumvent, amongst others, (discrete) Helmholtz and Hodge decompositions in the a posteriori error analysis.

The work [9] of Kreuzer, Verfürth, and Zanotti considers a discontinuous Galerkin (dG) approach for the discretization of the incompressible Stokes equations. For the Stokes equation, dG discretizations have the advantage that the canonical polynomial spaces (of order ℓ for the velocity resp. ℓ - 1 for the pressure) are automatically inf-sup stable. The proposed new method is proved to be quasi-optimal and pressure-robust (i.e., the approximation quality of the velocity field is independent of a possibly large pressure error). The key tool is a new reconstruction operator that maps discontinuous polynomial test functions to conforming test functions. This new operator particularly preserves the discrete divergence and thus maps discretely divergence-free test functions to exactly divergence-free test functions. With the help of this operator, the authors are able to derive quasi-optimal and pressure robust error estimates that are also valid for non-smooth right-hand sides, which are admissible to the continuous formulation.

The work [6] of Führer, García Vera, and Heuer derives an ultra-weak variational formulation of the Timoshenko beam bending model with various boundary conditions, combining clamped, simply supported, and free ends. The unknowns of the novel formulation, which is proved to be well-posed, are the transverse deflection and the bending moment. In the limit case of vanishing thickness, the proposed formulation reduces, as expected, to the Euler–Bernoulli model. Then, a discontinuous Petrov–Galerkin (DPG) discretization is proposed and it is proved that the method is locking-free and quasi-optimal, uniformly in the thickness of the beam.

3 Space-Time Discretizations

One key concern of numerical schemes is the preservation of mathematical structure of the PDE on the discrete level. The work [4] of Egger, Habrich, and Shashkov presents and analyzes a general framework for the numerical approximation of evolution problems that allows to preserve an underlying dissipative Hamiltonian or gradient structure exactly. In particular, it is shown that using the right canonical form on the abstract level, the Galerkin discretization in space combined with a Petrov–Galerkin discretization in time preserves the relevant energy/dissipation identities. For instance, the given framework applies to the equations of magneto-quasistatics that arise in the eddy current approximation of the Maxwell equations, the Cahn–Hilliard equation that models the phase separation in binary fluids, or constrained Hamiltonian systems that arise in the modeling of multibody dynamics.

The work [12] of Zank considers a space-time discretization of parabolic evolution equations in the setting of anisotropic space-time Sobolev spaces. While classical discretizations of evolution problems rely on time marching schemes, the benefits of space-time methods are that they have the potential for space-time adaptivity as well as parallelization. One possible way to obtain a variational space-time formulation with well-posed conforming Galerkin discretization (without any CFL condition) relies on a modified Hilbert transform. Having recalled this framework, the work discusses possible numerical realizations of this modified Hilbert transform. A new series expansion based on the Legendre chi function allows to compute the entries of the space-time Galerkin matrix up to machine precision, independently of the mesh-size.

The work [3] of Dubuis, Picasso, and Wittwer considers a PDE system, which couples the Navier–Stokes equations with the Newton laws describing the motion of a rigid body within an incompressible viscous fluid under the action of gravity. Numerical simulations of such systems face several difficulties because of possible collisions between the rigid body and the container boundary. The paper proposes an error indicator driven adaptive space-time algorithm, which makes it possible to compute quasi-collisions with high accuracy. The algorithm is validated numerically by simulating an analytical result, namely the non-collision in finite time for a rigid smooth disk falling toward a bottom wall in an incompressible Newtonian fluid.

4 Analysis of Adaptive Methods

An important aspect of adaptive algorithms is the overall computational cost due to the inherent cumulative nature (since any solution depends on the full adaptive history and computation). The work [7] of Giani, Grubišić, Heltai, and Mulita considers the interplay of adaptive mesh-refinement with an algebraic iterative eigenvalue solver. More precisely, the inexact algebraic iterative eigenvalue solvers on previous mesh levels are used to design an efficient solver on the last adaptively refined mesh. The astonishing observation of the proposed S-AFEM algorithm is that it is enough to smooth out only the high frequency error of the numerical approximation in order for the computed a posteriori error estimator to generate a sequence of quasi-optimal meshes. Overall, the computational costs for adaptive mesh refinement are thus greatly reduced, because no intermediate discrete solutions appear to be necessary to drive the mesh-refinement.

Many practically relevant PDEs are nonlinear so that the numerical approximation usually relies on an appropriate linearization of the discrete equations. The work [8] of Heid, Praetorius, and Wihler considers strongly-monotone elliptic PDEs. Based on an a posteriori error estimator, the authors formulate an adaptive iterative linearized finite element method (AILFEM), which steers the adaptive mesh-refinement (to resolve possible singularities of the unknown solution) as well as the linearization so that in each step only one linear system is solved. Under natural assumptions, it is proved that the proposed AILFEM algorithm leads to linear convergence in each step with optimal rates with respect to the computational cost. Numerical experiments underline the theory and show that AILFEM leads to optimal decay of the error with respect to the measured computational time.

In many situations, the relevant output of a numerical simulation is not the (accurate approximation of the) solution u of a PDE, but rather a functional value Q ⁢ ( u ) of some goal functional (also called quantity of interest). The work [5] of Endtmayer, Langer, and Wick proposes a new algorithm for the dual weighted residual method for goal-oriented FEM. The weights of the a posteriori error estimator are computed by comparing the primal and dual solution from the ansatz space with the corresponding ones from a higher order ansatz space. To lower the computational cost, the algorithm substitutes the solution of primal and dual problem in the higher-order space by an interpolation of the corresponding solutions from the lower-order space in some steps, based on an a posteriori criterion. For this procedure, the authors derive lower and upper bounds for the goal error. The new methodology is validated by considering a variety of model problems like the Poisson problem, a (regularized) p-Laplacian, and the Navier–Stokes problem.

A related question is addressed in the work [1] of Becker, Innerberger, and Praetorius, which is concerned with the thorough convergence analysis of a goal-oriented adaptive algorithm. Considering a slightly simpler model problem than [5], namely a general linear and elliptic PDE and a quadratic goal-functional, the work formulates an adaptive strategy which steers the linearization of the goal-functional as well as the local mesh-refinement for the discretization of the primal and the (linearized) dual problem. Under canonical assumptions on the a posteriori error estimators (which are met, for instance, for standard residual error estimators), it is proved that the adaptive strategy leads to linear convergence with respect to the adaptive level and, eventually, also to optimal convergence rates with respect to the number of the degrees of freedom.

5 Solvers and Preconditioning

The work [10] of Miraci, Papež, and Vohralik proposes an adaptive geometric multigrid method for the iterative solution of systems arising from the discretization of linear symmetric elliptic PDEs. First, using one V-cycle (“full-smoothing” substep) the proposed adaptive multigrid solver chooses adaptively the optimal step-size as well as the type of smoothing per level (weighted restricted additive or additive Schwarz), generalizing previous work of the authors. Second, another V-cycle (“adaptive-smoothing” substep) concentrates smoothing only to marked regions with (estimated) high error identified through a bulk-chasing criterion. The authors prove that the complete multigrid algorithm yields a uniform and p-robust contraction. Moreover, the proposed solver comes with a built-in a posteriori error estimator which provides a lower and upper bound for the unknown algebraic solver error (with known constant 1 for the lower bound).

Finally, the work [11] of van Venetië and Stevenson deals with the definition and analysis of preconditioners for problems involving operators of negative orders as met, e.g., in the boundary element method for weakly-singular integral operators. The work combines the ideas of an abstract framework developed by the authors and a stable multilevel splitting to obtain a preconditioner whose evaluation is of linear complexity. The result appears to be the first preconditioner based on the concept of operator preconditioning that can be applied in linear complexity and yields uniformly bounded condition numbers on locally refined meshes for the class of problems under consideration. For symmetric operators like the Laplace single-layer integral operator, this ensures that, e.g., the preconditioned conjugate gradient solver is uniformly contractive independently of the mesh.