Search

Convection-diffusion-reaction equations are a class of second-order partial differential equations widely used to model phenomena involving the change of concentration/population of one or more substances/species distributed in space. Understanding and preserving their stability properties in numerical simulation is crucial for accurate predictions, system analysis, and decision-making. This work presents a comprehensive framework for constructing fully discrete Lyapunov-consistent discretizations of any order for convection-diffusion-reaction models. We introduce a systematic methodology for constructing discretizations that mimic the stability analysis of the continuous model using Lyapunov's direct method. The spatial algorithms are based on collocated discontinuous Galerkin methods with the summation-by-parts property and the simultaneous approximation terms approach for imposing interface coupling and boundary conditions. Relaxation Runge-Kutta schemes are used to integrate in time and achieve fully discrete Lyapunov consistency. To verify the properties of the new schemes, we numerically solve a system of convection-diffusion-reaction partial differential equations governing the dynamic evolution of monomer and dimer concentrations during the dimerization process. Numerical results demonstrated the accuracy and consistency of the proposed discretizations. The new framework can enable further advancements in the analysis, control, and understanding of general convection-diffusion-reaction systems.

We propose a parallel (distributed) version of the spectral proper orthogonal decomposition (SPOD) technique. The parallel SPOD algorithm distributes the spatial dimension of the dataset preserving time. This approach is adopted to preserve the non-distributed fast Fourier transform of the data in time, thereby avoiding the associated bottlenecks. The parallel SPOD algorithm is implemented in the PySPOD (https://github.com/MathEXLab/PySPOD) library and makes use of the standard message passing interface (MPI) library, implemented in Python via mpi4py (https://mpi4py.readthedocs.io/en/stable/). An extensive performance evaluation of the parallel package is provided, including strong and weak scalability analyses. The open-source library allows the analysis of large datasets of interest across the scientific community. Here, we present applications in fluid dynamics and geophysics, that are extremely difficult (if not impossible) to achieve without a parallel algorithm. This work opens the path toward modal analyses of big quasi-stationary data, helping to uncover new unexplored spatio-temporal patterns.

MPI is the most widely used interface for high-performance computing (HPC) workloads. Its success lies in its embrace of libraries and ability to evolve while maintaining backward compatibility for older codes, enabling them to run on new architectures for many years. In this paper, we propose a new level of MPI compatibility: a standard Application Binary Interface (ABI). We review the history of MPI implementation ABIs, identify the constraints from the MPI standard and ISO C, and summarize recent efforts to develop a standard ABI for MPI. We provide the current proposal from the MPI Forum's ABI working group, which has been prototyped both within MPICH and as an independent abstraction layer called Mukautuva. We also list several use cases that would benefit from the definition of an ABI while outlining the remaining constraints.

We study temporal step size control of explicit Runge-Kutta methods for compressible computational fluid dynamics (CFD), including the Navier-Stokes equations and hyperbolic systems of conservation laws such as the Euler equations. We demonstrate that error-based approaches are convenient in a wide range of applications and compare them to more classical step size control based on a Courant-Friedrichs-Lewy (CFL) number. Our numerical examples show that error-based step size control is easy to use, robust, and efficient, e.g., for (initial) transient periods, complex geometries, nonlinear shock capturing approaches, and schemes that use nonlinear entropy projections. We demonstrate these properties for problems ranging from well-understood academic test cases to industrially relevant large-scale computations with two disjoint code bases, the open source Julia packages Trixi.jl with OrdinaryDiffEq.jl and the C/Fortran code SSDC based on PETSc.

Nonlinear entropy stability analysis is used to derive entropy stable no-slip wall boundary conditions for the Eulerian model proposed by Sv\"{a}rd (Physica A: Statistical Mechanics and its Applications, 2018). and its spatial discretization based on entropy stable collocated discontinuous Galerkin operators with the summation-by-parts property for unstructured grids. A set of viscous test cases of increasing complexity are simulated using both the Eulerian and the classic compressible Navier-Stokes models. The numerical results obtained with the two models are compared, and differences and similarities are then highlighted., Comment: The datasets generated during and/or analysed during the current study are available in the Zenodo repository, http://doi.org/10.5281/zenodo.5041436

We develop error-control based time integration algorithms for compressible fluid dynamics (CFD) applications and show that they are efficient and robust in both the accuracy-limited and stability-limited regime. Focusing on discontinuous spectral element semidiscretizations, we design new controllers for existing methods and for some new embedded Runge-Kutta pairs. We demonstrate the importance of choosing adequate controller parameters and provide a means to obtain these in practice. We compare a wide range of error-control-based methods, along with the common approach in which step size control is based on the Courant-Friedrichs-Lewy (CFL) number. The optimized methods give improved performance and naturally adopt a step size close to the maximum stable CFL number at loose tolerances, while additionally providing control of the temporal error at tighter tolerances. The numerical examples include challenging industrial CFD applications.

Recently, relaxation methods have been developed to guarantee the preservation of a single global functional of the solution of an ordinary differential equation. Here, we generalize this approach to guarantee local entropy inequalities for finitely many convex functionals (entropies) and apply the resulting methods to the compressible Euler and Navier-Stokes equations. Based on the unstructured $hp$-adaptive SSDC framework of entropy conservative or dissipative semidiscretizations using summation-by-parts and simultaneous-approximation-term operators, we develop the first discretizations for compressible computational fluid dynamics that are primary conservative, locally entropy stable in the fully discrete sense under a usual CFL condition, explicit except for the parallelizable solution of a single scalar equation per element, and arbitrarily high-order accurate in space and time. We demonstrate the accuracy and the robustness of the fully-discrete explicit locally entropy-stable solver for a set of test cases of increasing complexity.

In computational fluid dynamics, the demand for increasingly multidisciplinary reliable simulations, for both analysis and design optimization purposes, requires transformational advances in individual components of future solvers. At the algorithmic level, hardware compatibility and efficiency are of paramount importance in determining viability at exascale and beyond. However, equally important (if not more so) is algorithmic robustness with minimal user intervention, which becomes progressively more challenging to achieve as problem size and physics complexity increase. We numerically show that low and high order entropy stable discontinuous spatial discretizations based on summation-by-part operators and simultaneous-approximation-terms technique provides an essential step toward a truly enabling technology in terms of reliability and robustness for both under-resolved turbulent flow simulations and flows with discontinuities., Comment: 39 pages. arXiv admin note: substantial text overlap with arXiv:1911.03682 and text overlap with arXiv:1905.03007 by other authors

Computational fluid dynamics and aerodynamics, which complement more expensive empirical approaches, are critical for developing aerospace vehicles. During the past three decades, computational aerodynamics capability has improved remarkably, following advances in computer hardware and algorithm development. However, for complex applications, the demands on computational fluid dynamics continue to increase in a quest to gain a few percent improvements in accuracy. Herein, we numerically demonstrate that optimizing the metric terms which arise from smoothly mapping each cell to a reference element, lead to a solution whose accuracy is practically never worse and often noticeably better than the one obtained using the widely adopted Thomas and Lombard metric terms computation (Geometric conservation law and its application to flow computations on moving grids, AIAA Journal, 1979). Low and high-order accurate entropy stable schemes on distorted, high-order tensor product elements are used to simulate three-dimensional inviscid and viscous compressible test cases for which an analytical solution is known., Comment: 22 pages and one appendix section

In this paper, the entropy conservative/stable algorithms presented by Del Rey Fernandez and coauthors [18,16,17] for the compressible Euler and Navier-Stokes equations on nonconforming p-refined/coarsened curvilinear grids is extended to h/p refinement/coarsening. The main difficulty in developing nonconforming algorithms is the construction of appropriate coupling procedures across nonconforming interfaces. Here, a computationally simple and efficient approach based upon using decoupled interpolation operators is utilized. The resulting scheme is entropy conservative/stable and element-wise conservative. Numerical simulations of the isentropic vortex and viscous shock propagation confirm the entropy conservation/stability and accuracy properties of the method (achieving ~ p + 1 convergence) which are comparable to those of the original conforming scheme [4,35]. Simulations of the Taylor-Green vortex at Re = 1,600 and turbulent flow past a sphere at Re = 2,000 show the robustness and stability properties of the overall spatial discretization for unstructured grids. Finally, to demonstrate the entropy conservation property of a fully-discrete explicit entropy stable algorithm with h/p refinement/coarsening, we present the time evolution of the entropy function obtained by simulating the propagation of the isentropic vortex using a relaxation Runge-Kutta scheme., Comment: 37 pages. arXiv admin note: text overlap with arXiv:1909.12536, arXiv:1909.12546

The entropy conservative, curvilinear, nonconforming, p-refinement algorithm for hyperbolic conservation laws of Del Rey Fernandez et al. (2019), is extended from the compressible Euler equations to the compressible Navier-Stokes equations. A simple and flexible coupling procedure with planar interpolation operators between adjoining nonconforming elements is used. Curvilinear volume metric terms are numerically approximated via a minimization procedure and satisfy the discrete geometric conservation law conditions. Distinct curvilinear surface metrics are used on the adjoining interfaces to construct the interface coupling terms, thereby localizing the discrete geometric conservation law constraints to each individual element. The resulting scheme is entropy conservative/stable, element-wise conservative, and freestream preserving. Viscous interface dissipation operators are developed that retain the entropy stability of the base scheme. The accuracy and stability properties of the resulting numerical scheme are shown to be comparable to those of the original conforming scheme (achieving ~p+1 convergence) in the context of the viscous shock problem, the Taylor-Green vortex problem at a Reynolds number of Re=1,600, and a subsonic turbulent flow past a sphere at Re = 2,000., Comment: 16 pages

The framework of inner product norm preserving relaxation Runge-Kutta methods (David I. Ketcheson, \emph{Relaxation Runge-Kutta Methods: Conservation and Stability for Inner-Product Norms}, SIAM Journal on Numerical Analysis, 2019) is extended to general convex quantities. Conservation, dissipation, or other solution properties with respect to any convex functional are enforced by the addition of a {\em relaxation parameter} that multiplies the Runge-Kutta update at each step. Moreover, other desirable stability (such as strong stability preservation) and efficiency (such as low storage requirements) properties are preserved. The technique can be applied to both explicit and implicit Runge-Kutta methods and requires only a small modification to existing implementations. The computational cost at each step is the solution of one additional scalar algebraic equation for which a good initial guess is available. The effectiveness of this approach is proved analytically and demonstrated in several numerical examples, including applications to high-order entropy-conservative and entropy-stable semi-discretizations on unstructured grids for the compressible Euler and Navier-Stokes equations.

We present a novel technique for the imposition of non-linear entropy conservative and entropy stable solid wall boundary conditions for the compressible Navier-Stokes equations in the presence of an adiabatic wall, or a wall with a prescribed heat entropy flow. The procedure relies on the formalism and mimetic properties of diagonal-norm, summation-by-parts, and simultaneous-approximation-term operators, and is a generalization of previous works on discontinuous interface coupling [1] and solid wall boundary conditions [2]. Using the method of lines, a semi-discrete entropy estimate for the entire domain is obtained when the proposed numerical imposition of boundary conditions are coupled with an entropy-conservative or entropy-stable discrete interior operator. The resulting estimate mimics the global entropy estimate obtained at the continuous level. The boundary data at the wall are weakly imposed using a penalty flux approach and a simultaneous-approximation-term technique for both the conservative variables and the gradient of the entropy variables. Discontinuous spectral collocation operators (mass lumped nodal discontinuous Galerkin operators), on high-order unstructured grids, are used for the purpose of demonstrating the robustness and efficacy of the new procedure for weakly enforcing boundary conditions. Numerical simulations confirm the non-linear stability of the proposed technique, with applications to three-dimensional subsonic and supersonic flows. The procedure described is compatible with any diagonal-norm summation-by-parts spatial operator, including finite element, finite difference, finite volume, discontinuous Galerkin, and flux reconstruction schemes., Comment: 34 pages

We present a new method for performing global redistributions of multidimensional arrays essential to parallel fast Fourier (or similar) transforms. Traditional methods use standard all-to-all collective communication of contiguous memory buffers, thus necessary requiring local data realignment steps intermixed in-between redistribution and transform steps. Instead, our method takes advantage of subarray datatypes and generalized all-to-all scatter/gather from the MPI-2 standard to communicate discontiguous memory buffers, effectively eliminating the need for local data realignments. Despite generalized all-to-all communication of discontiguous data being generally slower, our proposal economizes in local work. For a range of strong and weak scaling tests, we found the overall performance of our method to be on par and often better than well-established libraries like MPI-FFTW, P3DFFT, and 2DECOMP&FFT. We provide compact routines implemented at the highest possible level using the MPI bindings for the C programming language. These routines apply to any global redistribution, over any two directions of a multidimensional array, decomposed on arbitrary Cartesian processor grids (1D slabs, 2D pencils, or even higher-dimensional decompositions). The high level implementation makes the code easy to read, maintain, and eventually extend. Our approach enables for future speedups from optimizations in the internal datatype handling engines within MPI implementations.

The phase-field crystal equation, a parabolic, sixth-order and nonlinear partial differential equation, has generated considerable interest as a possible solution to problems arising in molecular dynamics. This is because the phase-field crystal model can capture atomic-scale effects at time-scales that are orders of magnitude larger than what molecular dynamics simulations can afford presently. Nonetheless, solving this equation is not a trivial task, as a non-increasing free energy and mass conservation need to be verified for the numerical solution to be valid. This work focuses on these issues, and proposes a novel algorithm that guarantees mass conservation, unconditional energy stability and is second-order accurate in time. This is achieved through a convex-concave splitting of the nonlinearity present in the equation, along with the use of a stabilization term that bounds possible increases in free energy. We present numerical results that validate our mathematical proofs, and show two and three dimensional simulations involving crystal growth that showcase the robustness of the method., Comment: 43 pages, 13 figures

We present PetIGA, a code framework to approximate the solution of partial differential equations using isogeometric analysis. PetIGA can be used to assemble matrices and vectors which come from a Galerkin weak form, discretized with Non-Uniform Rational B-spline basis functions. We base our framework on PETSc, a high-performance library for the scalable solution of partial differential equations, which simplifies the development of large-scale scientific codes, provides a rich environment for prototyping, and separates parallelism from algorithm choice. We describe the implementation of PetIGA, and exemplify its use by solving a model nonlinear problem. To illustrate the robustness and flexibility of PetIGA, we solve some challenging nonlinear partial differential equations that include problems in both solid and fluid mechanics. We show strong scaling results on up to 4096 cores, which confirm the suitability of PetIGA for large scale simulations.

In this paper we study how the use of a more continuous set of basis functions affects the cost of solving systems of linear equations resulting from a discretized Galerkin weak form. Specifically, we compare performance of linear solvers when discretizing using $C^0$ B-splines, which span traditional finite element spaces, and $C^{p-1}$ B-splines, which represent maximum continuity. We provide theoretical estimates for the increase in cost of the matrix-vector product as well as for the construction and application of black-box preconditioners. We accompany these estimates with numerical results and study their sensitivity to various grid parameters such as element size $h$ and polynomial order of approximation $p$. Finally, we present timing results for a range of preconditioning options for the Laplace problem. We conclude that the matrix-vector product operation is at most $\slfrac{33p^2}{8}$ times more expensive for the more continuous space, although for moderately low $p$, this number is significantly reduced. Moreover, if static condensation is not employed, this number further reduces to at most a value of 8, even for high $p$. Preconditioning options can be up to $p^3$ times more expensive to setup, although this difference significantly decreases for some popular preconditioners such as Incomplete LU factorization.