Jan S. Hesthaven

Andreas Klöckner, Tim Warburton, Jeffrey Bridge et al.

Discontinuous Galerkin (DG) methods for the numerical solution of partial differential equations have enjoyed considerable success because they are both flexible and robust: They allow arbitrary unstructured geometries and easy control of accuracy without compromising simulation stability. Lately, another property of DG has been growing in importance: The majority of a DG operator is applied in an element-local way, with weak penalty-based element-to-element coupling. The resulting locality in memory access is one of the factors that enables DG to run on off-the-shelf, massively parallel graphics processors (GPUs). In addition, DG's high-order nature lets it require fewer data points per represented wavelength and hence fewer memory accesses, in exchange for higher arithmetic intensity. Both of these factors work significantly in favor of a GPU implementation of DG. Using a single US$400 Nvidia GTX 280 GPU, we accelerate a solver for Maxwell's equations on a general 3D unstructured grid by a factor of 40 to 60 relative to a serial computation on a current-generation CPU. In many cases, our algorithms exhibit full use of the device's available memory bandwidth. Example computations achieve and surpass 200 gigaflops/s of net application-level floating point work. In this article, we describe and derive the techniques used to reach this level of performance. In addition, we present comprehensive data on the accuracy and runtime behavior of the method.

Andreas Klöckner, Tim Warburton, Jan S. Hesthaven

We present a novel, cell-local shock detector for use with discontinuous Galerkin (DG) methods. The output of this detector is a reliably scaled, element-wise smoothness estimate which is suited as a control input to a shock capture mechanism. Using an artificial viscosity in the latter role, we obtain a DG scheme for the numerical solution of nonlinear systems of conservation laws. Building on work by Persson and Peraire, we thoroughly justify the detector's design and analyze its performance on a number of benchmark problems. We further explain the scaling and smoothing steps necessary to turn the output of the detector into a local, artificial viscosity. We close by providing an extensive array of numerical tests of the detector in use.

Babak Maboudi Afkham, Jan S. Hesthaven

While reduced-order models (ROMs) have been popular for efficiently solving large systems of differential equations, the stability of reduced models over long-time integration is of present challenges. We present a greedy approach for ROM generation of parametric Hamiltonian systems that captures the symplectic structure of Hamiltonian systems to ensure stability of the reduced model. Through the greedy selection of basis vectors, two new vectors are added at each iteration to the linear vector space to increase the accuracy of the reduced basis. We use the error in the Hamiltonian due to model reduction as an error indicator to search the parameter space and identify the next best basis vectors. Under natural assumptions on the set of all solutions of the Hamiltonian system under variation of the parameters, we show that the greedy algorithm converges with exponential rate. Moreover, we demonstrate that combining the greedy basis with the discrete empirical interpolation method also preserves the symplectic structure. This enables the reduction of the computational cost for nonlinear Hamiltonian systems. The efficiency, accuracy, and stability of this model reduction technique is illustrated through simulations of the parametric wave equation and the parametric Schrodinger equation.

Paolo Conti, Mengwu Guo, Andrea Manzoni et al.

When evaluating quantities of interest that depend on the solutions to differential equations, we inevitably face the trade-off between accuracy and efficiency. Especially for parametrized, time dependent problems in engineering computations, it is often the case that acceptable computational budgets limit the availability of high-fidelity, accurate simulation data. Multi-fidelity surrogate modeling has emerged as an effective strategy to overcome this difficulty. Its key idea is to leverage many low-fidelity simulation data, less accurate but much faster to compute, to improve the approximations with limited high-fidelity data. In this work, we introduce a novel data-driven framework of multi-fidelity surrogate modeling for parametrized, time-dependent problems using long short-term memory (LSTM) networks, to enhance output predictions both for unseen parameter values and forward in time simultaneously - a task known to be particularly challenging for data-driven models. We demonstrate the wide applicability of the proposed approaches in a variety of engineering problems with high- and low-fidelity data generated through fine versus coarse meshes, small versus large time steps, or finite element full-order versus deep learning reduced-order models. Numerical results show that the proposed multi-fidelity LSTM networks not only improve single-fidelity regression significantly, but also outperform the multi-fidelity models based on feed-forward neural networks.

Babak Maboudi Afkham, Jan S. Hesthaven

Reduced basis methods are popular for approximately solving large and complex systems of differential equations. However, conventional reduced basis methods do not generally preserve conservation laws and symmetries of the full order model. Here, we present an approach for reduced model construction, that preserves the symplectic symmetry of dissipative Hamiltonian systems. The method constructs a closed reduced Hamiltonian system by coupling the full model with a canonical heat bath. This allows the reduced system to be integrated with a symplectic integrator, resulting in a correct dissipation of energy, preservation of the total energy and, ultimately, in the stability of the solution. Accuracy and stability of the method are illustrated through the numerical simulation of the dissipative wave equation and a port-Hamiltonian model of an electric circuit.

Andreas Klöckner, Timothy Warburton, Jan S. Hesthaven

Discontinuous Galerkin (DG) methods for the numerical solution of partial differential equations have enjoyed considerable success because they are both flexible and robust: They allow arbitrary unstructured geometries and easy control of accuracy without compromising simulation stability. In a recent publication, we have shown that DG methods also adapt readily to execution on modern, massively parallel graphics processors (GPUs). A number of qualities of the method contribute to this suitability, reaching from locality of reference, through regularity of access patterns, to high arithmetic intensity. In this article, we illuminate a few of the more practical aspects of bringing DG onto a GPU, including the use of a Python-based metaprogramming infrastructure that was created specifically to support DG, but has found many uses across all disciplines of computational science.

Babak Maboudi Afkham, Ashish Bhatt, Bernard Haasdonk et al.

In the recent years, considerable attention has been paid to preserving structures and invariants in reduced basis methods, in order to enhance the stability and robustness of the reduced system. In the context of Hamiltonian systems, symplectic model reduction seeks to construct a reduced system that preserves the symplectic symmetry of Hamiltonian systems. However, symplectic methods are based on the standard Euclidean inner products and are not suitable for problems equipped with a more general inner product. In this paper, we generalize symplectic model reduction to allow for the norms and inner products that are most appropriate to the problem while preserving the symplectic symmetry of the Hamiltonian systems. To construct a reduced basis and accelerate the evaluation of nonlinear terms, a greedy generation of a symplectic basis is proposed. Furthermore, it is shown that the greedy approach yields a norm-bounded reduced basis. The accuracy and the stability of this model reduction technique are illustrated through the development of reduced models for a vibrating elastic beam and the sine-Gordon equation.

Fernando Henriquez, Jan S. Hesthaven

We extend our previous work [F. Henr'iquez and J. S. Hesthaven, arXiv:2403.02847 (2024)] to the linear, second-order wave equation in bounded domains. This technique uses two widely known mathematical tools to construct a fast and efficient method for the solution of linear, time-dependent problems: the Laplace transform (LT) and the model-order reduction (MOR) techniques, hence the name LT-MOR method. The application of the Laplace transform yields a time-independent problem parametrically depending on the Laplace variable. Following the two-phase paradigm of the reduced basis method, first in an offline stage we sample the Laplace parameter, compute the high-fidelity solution, and then resort to a Proper Orthogonal Decomposition (POD) to extract a basis of reduced dimension. Then, in an online phase, we project the time-dependent problem onto this basis and proceed to solve the evolution problem using any suitable time-stepping method. We prove exponential convergence of the reduced solution computed by the proposed method toward the high-fidelity one as the dimension of the reduced space increases. Finally, we present numerical experiments illustrating the performance of the method in terms of accuracy and, in particular, speed-up when compared to the full-order model.

Marco Sutti, Jan S. Hesthaven

We study the stability and sensitivity of an absorbing layer for the Boltzmann equation by examining the Bhatnagar-Gross-Krook (BGK) approximation and using the perfectly matched layer (PML) technique. To ensure stability, we discard some parameters in the model and calculate the total sensitivity indices of the remaining parameters using the ANOVA expansion of multivariate functions. We conduct extensive numerical experiments on two test cases to study stability and compute the total sensitivity indices, which allow us to identify the essential parameters of the model.

Max Hirsch, Federico Pichi, Jan S. Hesthaven

In this paper, we introduce the neural empirical interpolation method (NEIM), a neural network-based alternative to the discrete empirical interpolation method for reducing the time complexity of computing the nonlinear term in a reduced order model (ROM) for a parameterized nonlinear partial differential equation. NEIM is a greedy algorithm which accomplishes this reduction by approximating an affine decomposition of the nonlinear term of the ROM, where the vector terms of the expansion are given by neural networks depending on the ROM solution, and the coefficients are given by an interpolation of some "optimal" coefficients. Because NEIM is based on a greedy strategy, we are able to provide a basic error analysis to investigate its performance. NEIM has the advantages of being easy to implement in models with automatic differentiation, of being a nonlinear projection of the ROM nonlinearity, of being efficient for both nonlocal and local nonlinearities, and of relying solely on data and not the explicit form of the ROM nonlinearity. We demonstrate the effectiveness of the methodology on solution-dependent and solution-independent nonlinearities, a nonlinear elliptic problem, and a nonlinear parabolic model of liquid crystals. Code availability: https://github.com/maxhirsch/NEIM

Jan S. Hesthaven, Benjamin Peherstorfer, Benjamin Unger

This article surveys nonlinear model reduction methods that remain effective in regimes where linear reduced-space approximations are intrinsically inefficient, such as transport-dominated problems with wave-like phenomena and moving coherent structures, which are commonly associated with the Kolmogorov barrier. The article organizes nonlinear model reduction techniques around three key elements -- nonlinear parametrizations, reduced dynamics, and online solvers -- and categorizes existing approaches into transformation-based methods, online adaptive techniques, and formulations that combine generic nonlinear parametrizations with instantaneous residual minimization.

Oisín M. Morrison, Federico Pichi, Jan S. Hesthaven

This work presents a novel resolution-invariant model order reduction strategy for multifidelity applications. We base our architecture on a novel neural network layer developed in this work, the graph feedforward network, which extends the concept of feedforward networks to graph-structured data by creating a direct link between the weights of a neural network and the nodes of a mesh, enhancing the interpretability of the network. We exploit the method's capability of training and testing on different mesh sizes in an autoencoder-based reduction strategy for parametrised partial differential equations. We show that this extension comes with provable guarantees on the performance via error bounds. The capabilities of the proposed methodology are tested on three challenging benchmarks, including advection-dominated phenomena and problems with a high-dimensional parameter space. The method results in a more lightweight and highly flexible strategy when compared to state-of-the-art models, while showing excellent generalisation performance in both single fidelity and multifidelity scenarios.

Junming Duan, Qian Wang, Jan S. Hesthaven

Accurate prediction of aerodynamic forces in real-time is crucial for autonomous navigation of unmanned aerial vehicles (UAVs). This paper presents a data-driven aerodynamic force prediction model based on a small number of pressure sensors located on the surface of UAV. The model is built on a linear term that can make a reasonably accurate prediction and a nonlinear correction for accuracy improvement. The linear term is based on a reduced basis reconstruction of the surface pressure distribution, where the basis is extracted from numerical simulation data and the basis coefficients are determined by solving linear pressure reconstruction equations at a set of sensor locations. Sensor placement is optimized using the discrete empirical interpolation method (DEIM). Aerodynamic forces are computed by integrating the reconstructed surface pressure distribution. The nonlinear term is an artificial neural network (NN) that is trained to bridge the gap between the ground truth and the DEIM prediction, especially in the scenario where the DEIM model is constructed from simulation data with limited fidelity. A large network is not necessary for accurate correction as the linear model already captures the main dynamics of the surface pressure field, thus yielding an efficient DEIM+NN aerodynamic force prediction model. The model is tested on numerical and experimental dynamic stall data of a 2D NACA0015 airfoil, and numerical simulation data of dynamic stall of a 3D drone. Numerical results demonstrate that the machine learning enhanced model can make fast and accurate predictions of aerodynamic forces using only a few pressure sensors, even for the NACA0015 case in which the simulations do not agree well with the wind tunnel experiments. Furthermore, the model is robust to noise.

Federico Pichi, Beatriz Moya, Jan S. Hesthaven

The present work proposes a framework for nonlinear model order reduction based on a Graph Convolutional Autoencoder (GCA-ROM). In the reduced order modeling (ROM) context, one is interested in obtaining real-time and many-query evaluations of parametric Partial Differential Equations (PDEs). Linear techniques such as Proper Orthogonal Decomposition (POD) and Greedy algorithms have been analyzed thoroughly, but they are more suitable when dealing with linear and affine models showing a fast decay of the Kolmogorov n-width. On one hand, the autoencoder architecture represents a nonlinear generalization of the POD compression procedure, allowing one to encode the main information in a latent set of variables while extracting their main features. On the other hand, Graph Neural Networks (GNNs) constitute a natural framework for studying PDE solutions defined on unstructured meshes. Here, we develop a non-intrusive and data-driven nonlinear reduction approach, exploiting GNNs to encode the reduced manifold and enable fast evaluations of parametrized PDEs. We show the capabilities of the methodology for several models: linear/nonlinear and scalar/vector problems with fast/slow decay in the physically and geometrically parametrized setting. The main properties of our approach consist of (i) high generalizability in the low-data regime even for complex regimes, (ii) physical compliance with general unstructured grids, and (iii) exploitation of pooling and un-pooling operations to learn from scattered data.

Federico Pichi, Francesco Ballarin, Gianluigi Rozza et al.

This work deals with the investigation of bifurcating fluid phenomena using a reduced order modelling setting aided by artificial neural networks. We discuss the POD-NN approach dealing with non-smooth solutions set of nonlinear parametrized PDEs. Thus, we study the Navier-Stokes equations describing: (i) the Coanda effect in a channel, and (ii) the lid driven triangular cavity flow, in a physical/geometrical multi-parametrized setting, considering the effects of the domain's configuration on the position of the bifurcation points. Finally, we propose a reduced manifold-based bifurcation diagram for a non-intrusive recovery of the critical points evolution. Exploiting such detection tool, we are able to efficiently obtain information about the pattern flow behaviour, from symmetry breaking profiles to attaching/spreading vortices, even at high Reynolds numbers.

Przemyslaw Zielinski, Jan S. Hesthaven

Finding a reduction of complex, high-dimensional dynamics to its essential, low-dimensional "heart" remains a challenging yet necessary prerequisite for designing efficient numerical approaches. Machine learning methods have the potential to provide a general framework to automatically discover such representations. In this paper, we consider multiscale stochastic systems with local slow-fast time scale separation and propose a new method to encode in an artificial neural network a map that extracts the slow representation from the system. The architecture of the network consists of an encoder-decoder pair that we train in a supervised manner to learn the appropriate low-dimensional embedding in the bottleneck layer. We test the method on a number of examples that illustrate the ability to discover a correct slow representation. Moreover, we provide an error measure to assess the quality of the embedding and demonstrate that pruning the network can pinpoint an essential coordinates of the system to build the slow representation.

Hermes Sampedro Llopis, Allan P. Engsig-Karup, Cheol-Ho Jeong et al.

The use of model-based numerical simulation of wave propagation in rooms for engineering applications requires that acoustic conditions for multiple parameters are evaluated iteratively and this is computationally expensive. We present a reduced basis methods (RBM) to achieve a computational cost reduction relative to a traditional full order model (FOM), for wave-based room acoustic simulations with parametrized boundary conditions. In this study, the FOM solver is based on the spectral element method, however other numerical methods could be applied. The RBM reduces the computational burden by solving the problem in a low-dimensional subspace for parametrized frequency-independent and frequency-dependent boundary conditions. The problem is formulated and solved in the Laplace domain, which ensures the stability of the reduced order model based on the RBM approach. We study the potential of the proposed RBM framework in terms of computational efficiency, accuracy and storage requirements and we show that the RBM leads to 100-fold speed-ups for a 2D case with an upper frequency of 2kHz and around 1000-fold speed-ups for an analogous 3D case with an upper frequency of 1kHz. While the FOM simulations needed to construct the ROM are expensive, we demonstrate that despite this cost, the ROM has a potential of three orders of magnitude faster than the FOM when four different boundary conditions are simulated per room surface. Moreover, results show that the storage model for the ROM is relatively high but affordable for the presented 2D and 3D cases.

Mengwu Guo, Andrea Manzoni, Maurice Amendt et al.

Highly accurate numerical or physical experiments are often time-consuming or expensive to obtain. When time or budget restrictions prohibit the generation of additional data, the amount of available samples may be too limited to provide satisfactory model results. Multi-fidelity methods deal with such problems by incorporating information from other sources, which are ideally well-correlated with the high-fidelity data, but can be obtained at a lower cost. By leveraging correlations between different data sets, multi-fidelity methods often yield superior generalization when compared to models based solely on a small amount of high-fidelity data. In this work, we present the use of artificial neural networks applied to multi-fidelity regression problems. By elaborating a few existing approaches, we propose new neural network architectures for multi-fidelity regression. The introduced models are compared against a traditional multi-fidelity scheme, co-kriging. A collection of artificial benchmarks are presented to measure the performance of the analyzed models. The results show that cross-validation in combination with Bayesian optimization consistently leads to neural network models that outperform the co-kriging scheme. Additionally, we show an application of multi-fidelity regression to an engineering problem. The propagation of a pressure wave into an acoustic horn with parametrized shape and frequency is considered, and the index of reflection intensity is approximated using the multi-fidelity models. A finite element model and a reduced basis model are adopted as the high- and low-fidelity, respectively. It is shown that the multi-fidelity neural network returns outputs that achieve a comparable accuracy to those from the expensive, full-order model, using only very few full-order evaluations combined with a larger amount of inaccurate but cheap evaluations of a reduced order model.

Jim Magiera, Deep Ray, Jan S. Hesthaven et al.

Neural networks are increasingly used in complex (data-driven) simulations as surrogates or for accelerating the computation of classical surrogates. In many applications physical constraints, such as mass or energy conservation, must be satisfied to obtain reliable results. However, standard machine learning algorithms are generally not tailored to respect such constraints. We propose two different strategies to generate constraint-aware neural networks. We test their performance in the context of front-capturing schemes for strongly nonlinear wave motion in compressible fluid flow. Precisely, in this context so-called Riemann problems have to be solved as surrogates. Their solution describes the local dynamics of the captured wave front in numerical simulations. Three model problems are considered: a cubic flux model problem, an isothermal two-phase flow model, and the Euler equations. We demonstrate that a decrease in the constraint deviation correlates with low discretization errors for all model problems, in addition to the structural advantage of fulfilling the constraint.

Paolo Gatto, Jan S. Hesthaven

We introduce a preconditioner based on a hierarchical low-rank compression scheme of Schur complements. The construction is inspired by standard nested dissection, and relies on the assumption that the Schur complements can be approximated, to high precision, by Hierarchically-Semi-Separable matrices. We build the preconditioner as an approximate $LDM^t$ factorization of a given matrix $A$, and no knowledge of $A$ in assembled form is required by the construction. The $LDM^t$ factorization is amenable to fast inversion, and the action of the inverse can be determined fast as well. We investigate the behavior of the preconditioner in the context of DG finite element approximations of elliptic and hyperbolic problems.

Akil C. Narayan, Jan S. Hesthaven

We formulate and derive a generalization of an orthogonal rational-function basis for spectral expansions over the infinite or semi-infinite interval. The original functions, first presented by Wiener are a mapping and weighting of the Fourier basis to the infinite interval. By identifying the Fourier series as a biorthogonal composition of Jacobi polynomials/functions, we are able to define generalized Fourier series' which, appropriately mapped to the whole real line and weighted, form a generalization of Wiener's basis functions. It is known that the original Wiener rational functions inherit sparse Galerkin matrices for differentiation, and can utilize the fast Fourier transform (FFT) for computation of the modal coefficients. We show that the generalized basis sets also have a sparse differentiation matrix and we discuss connection problems, which are necessary theoretical developments for application of the FFT.