Schedule for: 21w5065 - Bound-Preserving Space and Time Discretizations for Convection-Dominated Problems (Online)

We consider hyperbolic systems of balance laws. Our goal is to develop well-balanced numerical methods, which respect a delicate balance between the flux and source terms and are thus capable of exactly preserving (some of the) physically relevant steady-state solutions of the studied systems. I will introduce a general approach of constructing well-balanced schemes via a flux globalization ap- proach: The source terms are incorporated into the fluxes. This results in the hyperbolic system of conser- vation laws with global fluxes. Such systems can be then integrated using Riemann-problem-solver-free numerical methods. I will show several recent non-straightforward applications of these well-balanced schemes.

Many mathematical models are equipped with an energy that is conserved or an entropy that is known to change monotonically in time. Integrators that preserve these properties discretely are usually expensive, with the best-known examples being fully-implicit Runge-Kutta methods. I will present a modification that can be applied to any integrator in order to preserve such a structural property. The resulting method can be fully explicit, or (depending on the functional) may require the solution of a scalar algebraic equation at each step. I will present examples to show the effectiveness of these “relaxation” methods, and their advantages over fully implicit methods or orthogonal projection. Examples will include applications to compressible fluid dynamics, dispersive nonlinear waves, and Hamiltonian systems.

We present a new high-resolution, non-oscillatory semi-discrete central scheme for one-dimensional two- layer shallow-water flows along channels with arbitrary cross sections and bottom topography. The scheme extends existing central semi-discrete schemes for hyperbolic conservation laws and it enjoys two properties crucial for the accurate simulation of shallow-water flows: it preserves the positivity of the water height, and it is well balanced, i.e., the source terms arising from the geometry of the channel are discretized so as to balance the non-linear hyperbolic flux gradients. The system is integrated in time using a second order Strong Stability Preserving Runge-Kutta scheme. Along with a detailed description of the scheme and proofs of these two properties, we present several numerical experiments that demonstrate the robustness of the numerical algorithm. This is joint work with Jorge Balbas.

We first propose new flux limiters that have the structure of a flux-corrected transport (FCT) algo- rithm. These limiters do not depend on the time step size and, therefore, are applicable to spatial semi- discretizations. Our limiters contain a user defined parameter that improves accuracy at a cost of a more restrictive CFL-like condition, for explicit implementations. To test the accuracy properties of these lim- iters, we consider a spatial semi-discretization of Burgers’ equation based on arbitrarily high-order WENO reconstructions and use the limiters to impose global bounds. Using the same WENO reconstructions, we obtain a full discretization based on Runge-Kutta methods. This scheme is arbitrarily high-order and ‘essentially’ non-oscillatory but does not preserve the maximum principle. To guarantee the maximum principle we combine the fluxes of the high-order scheme with those of a low-order method based on forward Euler and Local Lax-Friedrichs (LLF) fluxes. We obtain anti-diffusive fluxes that combine corrections in space-and-time to the low-order scheme. We use our proposed limiters to guarantee the solution is maximum principle preserving. Finally, we present a similar methodology using an arbitrarily high-order Singly Diagonal RK method combined with a low-order scheme based on backward Euler and LLF fluxes. The implicit scheme is maximum principle preserving for time steps of any size.

While positivity and boundedness preservation plays a key role in convection dominated problems it also plays a key in a number of situations in PDEs where discrete interpolation is required. These situations include (i) Mapping from physics to dynamics grids in weather codes (ii) discrete remappings in adaptive meshing (iii) particle to mesh remapping in particle in call an material point methods (iv) solution recreation in AMR codes to address compute node failures. Time integration error analysis shows how interpolation errors must be controlled to avoid polluting the main calculation and the connection with stage errors in Runge-Kutta methods in he case of particle methods. For interpolation itself we provide a simple constructive proof for an adaptive algorithm tensor-product grids on arbitrary spacings to preserve boundedness and positivity and show results for cases (i) and (iii) above. Reference is also made to a longer more comprehensive proof.

We interpret the enriched Galerkin (EG) method as generalization of standard finite elements (contin- uous Galerkin, CG) and of the discontinuous Galerkin (DG) method by combining the continuous and the discontinuous trial spaces of CG and DG, and by using the DG bilinear and linear forms. Then, we introduce algebraic flux correction schemes for the standard enriched (P1⊕P0 and Q1⊕P0) Galerkin discretizations of the linear advection equation. Here, the piecewise-constant component stabilizes the continuous Galerkin approximation without introducing free parameters. However, violations of discrete maximum principles are possible in the neighborhood of discontinuities and steep fronts. To keep the cell averages and the degrees of freedom of the continuous P1/Q1 component in the admissible range, we limit the fluxes and element contributions, the complete removal of which would correspond to first-order upwinding. Finally, we discuss a further generalization of the enriched Galerkin method. The key feature of this step is an adaptive two-mesh approach that, in addition to the standard enrichment of a conforming finite element discretization via discontinuous degrees of freedom, allows to subdivide selected (e.g. troubled) mesh cells in a non-conforming fashion and to use further discontinuous enrichment on this finer submesh. Here, we prove stability and sharp a priori error estimates for a linear advection equation under appropriate assumptions.

Since the celebrated Lax Wendroff convergence theorem, published in 1960 in CPAM, every one knows what should be the structure of a finite volume/finite difference scheme so that one can have a reasonable hope of convergence towards the ‘true’ entropy solution. The proof can easily be adapted to schemes like the discontinuous Galerkin ones, thought it becomes a bit less clear. There are however many schemes that does not fit clearly in that framework: continuous finite element methods, for example. Though it is relatively easy to prove a variant of the Lax Wendroff for them, this does not answer the question of the engineer: show me explicitely the flux. There are other questions related to conservation. We all know that it is forbiden by the Law to discretize the non conservative version of a conservative system (for example, the Euler equations in primitive variables), and there are many counter examples. However, to which extend is that statement true? If one has an additional conservation law satisfied by the system (for example entropy conservation for smooth solutions, or kinetic momentum preservation), how can we modify a given ‘good’ scheme so that the modified one will satisfy all constraints? In this talk, which sumarizes [1,2,3,4,5,6], I will try to show the boundaries of these statements, and provide example of schemes, some already known, some more recent, that contradict, in some sense, the standard beliefs. But not too much. References 1. R. Abgrall. Some remarks about conservation for residual distribution schemes. Computational Meth- ods in Applied Mathematics, 18(3):327–350, 2018. 2. R. Abgrall and S. Tokareva. Staggered grid residual distribution scheme for Lagrangian hydrodynamics. SIAM J. Scientific Computing, 39(5):A2317–A2344, 2017. 3. R. Abgrall. A general framework to construct schemes satisfying additional conservation relations, application to entropy conservative and entropy dissipative schemes. J. Comput. Phys, 372(1), 2020. 4. Nathaniel Morgan Rémi Abgrall, Konstantin Lipnikov and Svetlana Tokareva. Multidimensional Stag- gered Grid Residual Distribution Scheme for Lagrangian Hydrodynamics. SIAM J. Sci. Comput., 42(1):A343–A370, 2020. 5. R. Abgrall. A combination of Residual Distribution and the Active Flux formulations or a new class of schemes that can combine several writings of the same hyperbolic problem: application to the 1D Euler equations, 2021. https://arxiv.org/abs/2011.12572 6. R. Abgrall, P. Öffner, and H. Ranocha. Reinterpretation and Extension of Entropy Correction Terms for Resi

Discontinuous Galerkin (DG) methods are among the most widely used numerical discretization tech- niques for solving partial differential equations (PDEs). Their local conservation property, inherent stability, and, favorable scalability in parallel make these schemes attractive for many applications. There are how- ever, many shock-dominated examples, for which even DG methods fail to produce stable approximations. To overcome this shortcoming, we developed an algebraic flux limiter, which blends a provably property- preserving low order method with a corresponding high order DG target scheme. This monolithic convex limiter is primarily ]used to impose local (and global) bounds on numerical approximations, but extensions for incorporating entropy inequalities, as well as relaxation of the constraints in smooth regions are also possible. In my talk, I will discuss the details of the approach, which include the sparsification of the low order method, stabilization of the numerical flux, as well as the design of the monolithic limiter. Sequential limiting for products of unknowns and the preservation of global constraints, such as nonnegativity of pres- sure will also be touched upon and similarities to comparable schemes will be put into context. All presented numerical results were obtained with a code that is based on the open source C++ library MFEM. The performance of the method will be evaluated by considering a variety of classical benchmarks for scalar conservation laws, as well as the systems of shallow water and Euler equations.

In this talk, we will discuss design of structure-preserving central-upwind finite volume methods for shallow water models in domains with irregular geometry and for shallow water models with uncertainty. Shallow water models are widely used in many scientific and engineering applications related to modeling of water flows in rivers, lakes and coastal areas. Shallow water equations are examples of hyperbolic systems of balance laws and such mathematical models can present a significant challenge for the construction of accurate and efficient numerical algorithms. We will show that the developed structure-preserving central-upwind schemes for shallow water equations deliver high-resolution, can handle complicated geometry, and satisfy necessary stability conditions. We will illustrate the performance of the designed methods on a number of challenging numerical tests. Current and future research will be discussed as well.

In this talk we present our work in updating the high-order finite element-based ALE remap method in the MARBL multi-physics code from LLNL [1] for high performance on GPU platforms. MARBL is a multi- material hydrodynamics code based on a three phase Arbitrary-Lagrangian-Eulerian (ALE) framework: evolution of physical conservation laws within a moving material (Lagrangian) frame; mesh optimization, and field remap. The remap step corresponds to the transfer of state variables from the initial mesh to the optimized mesh and can often dominate the run time for typical calculations. For high fidelity simulations it is imperative that the remap procedure be accurate, preserve physical quantities, and be monotonic (not introduce new extrema). Furthermore, for scalable performance in large scale calculations it is imperative that the characteristics of the algorithm are well suited for modern computing platforms (e.g. GPU based architectures). This work is based on adopting a matrix-free algorithmic approach. The current remap algorithm in MARBL is based on the work of Anderson and co-authors [1, 2] which introduce a high-order approach based on concepts from flux corrected transport (FCT) and a discontinuous Galerkin (dG) discretization for the advection equation but requires full matrix assembly due to its algebraic nature. Methods based on full matrix assembly are known to have poor performance as the order of the method is increased. In addition, there are involved memory motion operations which do not work well on GPU architectures. The new matrix-free framework we have been developing combines the residual distribution schemes of Hajduk and co-authors [3,4] with a high-order DG scheme using the clip scale strategy of Anderson et al. in [5]. Lastly, we describe the algorithmic tailoring for the GPU and present performance comparisons between the different frameworks. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LLNL-ABS-824642.

In the first part of the talk, the SUPG method for continuous piecewise linear finite elements is considered. Numerical solutions computed with these methods are known to possess spurious oscillations in a vicinity of layers. In [JKS11], a general approach for optimizing the stabilization parameter was presented. This talk will address an open problem stated in this paper: the restriction of the optimization to subregions were the choice of the stabilization parameter is essential. In this way, a reduction of the dimension of the space for optimization is achieved. Suitable algorithms are discussed and numerical studies are presented. The second part of the talk deals with discontinuous Galerkin (DG) finite element methods. These methods are known to be stable and to compute sharp layers in the convection-dominated regime, but also to show large spurious oscillations. Post-processing methods for reducing spurious oscillations are discussed, which re- place the DG solution in a vicinity of layers by a constant or linear approximation. A survey of methods that are available in the literature is presented and several generalizations and modifications are proposed. Numerical studies illustrate the behavior of these methods. Details can be found in [FJ21]. This talk presents joint work with Ulrich Wilbrandt (WIAS) and Derk Frerichs (WIAS). References FJ21 Derk Frerichs and Volker John. On reducing spurious oscillations in dis- continuous Galerkin (DG) methods for steady-state convection-diffusion equations. J. Comput. Appl. Math., 393:113487, 20, 2021. JKS11 Volker John, Petr Knobloch, and Simona B. Savescu. A posteriori op- timization of parameters in stabilized methods for convection-diffusion problems—Part I. Comput. Methods Appl. Mech. Engrg., 200(41-44):2916– 2929, 2011.

Non-linear discretizations are necessary for convection-diffusion-reaction equations for obtaining accurate solutions that satisfy the discrete maximum principle (DMP). Algebraic stabilizations, also known as Algebraic Flux Correction (AFC) schemes, belong to the very few finite element discretizations that satisfy this property. Results regarding the convergence of the scheme [1] and efficient solution of the nonlinear system of equations [2] have been obtained recently.. The talk is devoted to the proposal of a new residual based a posteriori error estimator for AFC schemes. We derive a global upper bound in the energy norm of the system which is independent of the choice of the limiter in the AFC scheme. We also derive a global upper bound by combining the estimators from [3] and the AFC schemes. Numerical simulations in 2d are presented which support the analytical findings. References 1. Gabriel R. Barrenechea, Volker John, and Petr Knobloch. Analysis of algebraic flux correction schemes. SIAM J. Numer. Anal., 54(4):2427–2451, 2016. 2. Abhinav Jha and Volker John. A study of solvers for nonlinear AFC discretizations of convection- diffusion equations. Comput. Math. Appl., 78(9):3117–3138, 2019. 3. Volker John and Julia Novo. A robust SUPG norm a posteriori error estimator for stationary convection- diffusion equations. Comput. Methods Appl. Mech. Engrg., 255:289–305, 2013.

In this talk we will describe a second-order continuous finite element technique for solving hyperbolic systems in the arbitrary Lagrangian Eulerian framework (ALE). The main property of the method presented is that, provided the user-defined ALE velocity is reasonable, the approximate solution produced by the algorithm is formally second-order accurate in space, is conservative and preserves as many convex invariant sets of the hyperbolic system as desired by the user, by using a convex limiting technique. The time stepping is explicit, the approximation in space is done with continuous finite elements. One of the key issues in the definition of the ALE motion is the mesh optimization process in order to locally optimize the mesh quality and avoid the generation of invalid elements with negative determinant. We divide the optimization strategy in two stages: (i) a high order reconstruction of the Lagrangian velocity; (ii) a smoothing of the mesh. The second-order accuracy is numerically shown to hold in the maximum norm for smooth solutions of the compressible Euler equations and compared to the same method in Eulerian coordinates.

The objective of this talk is to present a fully-discrete approximation technique for the compressible Navier-Stokes equations. The method is implicit-explicit, second-order accurate in time and space, and guaranteed to be invariant domain preserving. The restriction on the time-step size is the standard hyperbolic CFL condition. One key originality of the method is that it is guaranteed to be conservative and invariant domain preserving under the standard hyperbolic CFL condition. The method is numerically illustrated on the OAT15a airfoil in the critical transonic regime at Re=3 millions. This is a joint work with M. Kronbichler, M. Maier, B. Popov and I. Tomas.

We introduce a relaxation technique for solving the Serre Equations for dispersive water waves. The novelty of this technique is the reformulation of the Serre Equations into a hyperbolic system which allows for explicit time stepping in the numerical method. We then propose a second-order approximation of the model using continuous finite elements that is well-balanced and positivity preserving. The method is then numerically validated and illustrated by comparison with laboratory experiments.

Patankar-type scheme are linearly implicit time integration methods constructed to satisfy positivity properties of certain ordinary differential equations. Since they are outside of the class of general linear methods, they can have superior positivity preserving properties. However, standard notions of stability do not apply. For example, classical linear stability analysis cannot be used since the schemes do not commute with diagonalization. Hence, new concepts of stability need to be introduced. We provide preliminary investigations of stability of Patankar-type schemes. In particular, we demon- strate problematic behavior of these methods that can lead to undesired oscillations or order reduction. Extreme cases of the latter manifest as spurious steady states. We investigate stability properties of vari- ous classes of Patankar-type schemes based on classical Runge-Kutta methods, strong stability preserving Runge-Kutta methods, and deferred correction schemes. This project is joint work with Davide Torlo and Philipp O ̈ffner.

Since the Cauchy problem for the complete Euler system is in general ill- posed in the class of admissible entropy weak solutions, one searches alternatives and here the concept of dissipative weak solutions seems quite promising to analyze this system analytically and numerically. In [1], the authors have studied the convergence properties of a class of entropy dissipative finite volume schemes for the barotropic and complete compressible Euler equations in the multidimensional case and could prove suitable stability and consistency properties to ensure convergence of the FV schemes via dissipative measure- valued solutions. In a series of paper, the theory has been further developed for several (classical) FV schemes (of maximum order two) and have been tested numerically, cf. [2,3]. In this talk, we consider as well convergence via dissipative weak solutions for the Euler equation, but focus on high-order finite element based methods, in particular on a specific discontinuous Galerkin schemes. For the convergence proof, we need certain properties like the preservation of several phyiscal quantities and some entropy estimates. We demonstrate how we ensure these and prove convergence of our DG scheme via dissipative weak solutions. In numerical simulations, we verify our theoretical findings. References 1. E. Feireisl, M. Luk ́aˇcov ́a-Medvid’ova ́ and H. Mizerov ́a. Convergence of finite volume schemes for the Euler equations via dissipative measure-valued solutions, Foundations of Computational Mathematics, 20 (4) pp. 1–44, 2019. 2. E. Feireisl, M. Luk ́aˇcov ́a-Medvid’ova ́ and H. Mizerov ́a. A finite volume scheme for the Euler system inspired by the two velocities approach, Numerische Mathematik, 144(1), pp. 89–132, 2020. 3. E. Feireisl, M. Luk ́aˇcov ́a-Medvid’ova ́, B. She and Y. Wang. Computing oscillatory solutions of the Euler system via K-convergence, Mathematical Models and Methods in Applied Sciences, pp. 1–40, 2021.

We consider continuous finite element approximations of hyperbolic problems and modify them to satisfy relevant inequality constraints. The proposed approaches apply limiters to fluxes that represent the difference between a high-order target scheme and a low-order property-preserving approximation of Lax-Friedrichs type. In the first step of the limiting procedure, the given target fluxes are adjusted in a way that guarantees preservation of local and/or global bounds. In the second step, additional limiting is performed, if necessary, to ensure the validity of fully discrete and/or semi-discrete entropy inequalities. The limiter-based entropy fixes impose entropy-conservative or entropy-dissipative bounds on entropy production by antidiffusive fluxes and Runge-Kutta time discretizations. We present three algorithms developed for this purpose. The semi- discrete (SD) fix is based on Tadmor’s entropy stability theory and constrains the spatial semi-discretization. The fully discrete explicit (FDE) fix incorporates temporal entropy production into the flux constraints, which makes them more restrictive. The fully discrete implicit (FDI) fix performs iterative flux correction under SD-type constraints in the final Runge-Kutta stage. The effectiveness of these fixes is verified in numerical experiments for scalar equations and systems. References: 1. D. Kuzmin and M. Quezada de Luna, Algebraic entropy fixes and convex limiting for continuous finite element discretizations of scalar hyperbolic conservation laws. Computer Methods Appl. Mech. Engrg. 372 (2020) 113370. 2. D. Kuzmin, H. Hajduk and A. Rupp, Limiter-based entropy stabilization of semi-discrete and fully discrete schemes for nonlinear hyperbolic problems. Preprint arXiv:2107.11283 [math.NA], July 2021.

Chemotaxis models describe the evolution of biological migration processes. In a migration process, or- ganisms (or a group of cells) migrate in response to a chemical stimulus, which either attracts or repels them. In this work we focus on Keller-Segel equations. A model that despite its biologically inaccurate results, it is interesting and challenging from the mathematical point of view. Keller-Segel model solutions satisfy lower bounds, and enjoy an energy law. The mathematical interest lies in developing numerical discretiza- tions that yield solutions that preserve these properties, with the aim to contribute and give insights into discretizations for more complex and realistic chemotaxis models. Positivity-preservation is very important in these models to ensure physically meaningful results and proper evolution of the solution. Otherwise, it might lead to negative concentrations of cells or chemo-attractant. This is especially challenging in this model because, depending on the particular initial conditions specified, it might lead to a blow-up of the solution. Recently, Guti ́errez-Santacreu and Rodriguez-Galv ́an [1] have published a numerical scheme for strictly acute meshes, that also yields solutions satisfying lower bounds and an energy law in the discrete sense. In the present work, we aim to extend these results to general meshes using an artificial diffusion stabilization method based on [3]. Numerical results also show that the artificially added diffusion does not significantly smear the solution, while it enforces lower bounds even for solutions that blow-up in a finite time. References 1. J.V. Gutierr ́ez-Santacreu, J.R. Rodriguez-Galv ́an, Analysis of a fully discrete approximation for the classical Keller-Segel model: Lower and a priori bounds, Comput. Math. Appl. 85 (2021) 69–81. 2. S. Badia, J. Bonilla, J.V. Gutierr ́ez-Santacreu. Solving the Keller-Segel equations with finite element approximations over general meshes, In preparation. 3. Badia, S., Bonilla, J. Monotonicity-preserving finite element schemes based on differentiable nonlinear stabilization, Comput. Methods Appl. Mech. Engrg. 313 (2017) 133–15