Collaborators: Cody J. Balos, David J. Gardner, Daniel R. Reynolds, and Steven Roberts
The SUNDIALS library of time integrators and nonlinear solvers has recently increased its support for large-scale GPU-based systems through new data structures and solver package interfaces. This talk will overview the SUNDIALS packages and recently added capabilities then show results on newly deployed systems. Several applications have taken advantage of these capabilities to improve their time integration performance, and results from a selection of these, including combustion, phase field modeling, and cosmology, will be shown.
Prepared by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LLNL-ABS-844590.
In this talk we briefly review some basic PDE models that are used to model phase separation in materials science. They have since become important tools in image processing and over the last years semi-supervised learning strategies could be implemented with these PDEs at the core. The main ingredient is the graph Laplacian that stems from a graph representation of the data. This matrix is large and typically dense. We illustrate some of its crucial features and show how to efficiently work with the graph Laplacian. In particular, we need some of its eigenvectors and for this the Lanczos process needs to be implemented efficiently. Here, we suggest the use of the NFFT method for evaluating the matrix vector products without even fully constructing the matrix. We illustrate the performance on several examples.
Data assimilation is used in a variety of disciplines and applications to combine mathematical models and numerical simulations with measurement data to provide a more accurate representation of the actual processes. In this way, potential model errors can be corrected, taking into account possible measurement errors, while still respecting the basic laws of physics.
In this talk, we investigate two different aspects of data assimilation. First, we propose a strategy for constructing adaptive time grids in variational data assimilation using a posteriori error estimators, with the goal of taking measurements at time points where particularly important dynamical processes occur. This part is joint work with Jannis Marquardt. Second, we present an application problem where the goal is to identify potential damages in fiber metal laminates using numerical models and measured data by means of data assimilation techniques. For this, we use Ensemble Kalman Filter methods and apply model reduction techniques to accelerate the underlying simulations. This part is joint work with Saddam Hijazi and is supported by the DFG (FOR 3022, project number 418311604).
In many experiments in the natural sciences, the quantity of interest is not directly accessible, but needs to be reconstructed from measurements of its impact. Mathematically, the identification of a cause from its impact means solving an inverse problem, where cause and impact are connected via a mathematical model. Often, inverse problems that arise in applications rely on a physical model that is either somewhat idealized or partly unknown. This is due to the high complexity of the actual physical processes that are involved in an experiment. By consequence, this modelling inexactness has to be taken into account during regularization of the inverse problem to reduce artifacts in the reconstruction.
In this talk, we discuss the role of modelling errors and present strategies to compensate them. We focus on tomographic X-ray imaging on the nano-scale, where inexactness in the model is introduced due to vibrations of the investigated object as well as fluctuations in the X-ray beam during data acquisition.
The goal of climate neutrality is driving the transformation of energy systems worldwide. Energy markets, in particular, undergo fundamental changes due to the integration of new technologies, including renewable energies and fuels, and necessary adjustments in regulatory frameworks, such as CO2 pricing. Equilibrium problems offer a powerful framework to capture these dynamic developments in energy systems, enabling valuable insights into key research questions related to energy market design and climate policy.
This talk provides an introduction to equilibrium concepts, relevant theory, and their applications in energy markets. We begin by exploring foundational concepts, including market equilibrium problems that model competitive market behavior and Stackelberg games that capture hierarchical decision-making structures frequently observed in energy systems. In the theory part of the talk, we delve into key mathematical properties of equilibria, such as conditions for the existence and uniqueness. Finally, the practical impact of equilibrium modeling is demonstrated using selected examples from energy markets.
Gaussian Processes (GPs) and kernel-based methods, such as Kernel Ridge Regression (KRR), are widely used for prediction making in Machine Learning. Both methods rely on positive definite kernels and are closely related. However, this relationship results in shared limitations, particularly their cubic computational complexity with respect to the number of training samples. This makes them impractical for large data sets, as the primary computational bottleneck arises from solving a linear system of equations during model fitting.
To address this challenge while preserving predictive performance, various approaches have been proposed, including the Nyström method and Sparse Variational Gaussian Processes. Both offer approximate solutions to the linear system of equations arising during training. In this work, we introduce an alternative Randomized Low-Rank Approximation (RLoRA) technique which leverages randomized techniques to identify a suitable low-rank subspace of the kernel matrix. In contrast, to the Nystr¨om method, which draws random samples from the columns and rows of the kernel matrix, RLoRA uses random samples of the column space to build a subspace. This subspace is assumed to capture most of the kernel matrix’ range.
We demonstrate the efficacy of our approach using quantum chemistry data sets, showing that RLoRA consistently outperforms the Nyström approximation in Kernel Ridge Regression in terms of predictive accuracy. With RLoRA, we aim to make kernel-based methods suitable for large data sets and easily accessible in practice.
The problem of modelling processes with partial differential equations posed on random domains arises in various applications like biology or engineering. We study uncertainty quantification for partial differential equations subject to domain uncertainty, where we parameterize the random domain using the model recently considered by Chernov and Lê (Comput. Math. Appl., 2024, and SIAM J. Numer. Anal., 2024) as well as Harbrecht, Schmidlin, and Schwab (Math. Models Methods Appl. Sci., 2024) in which the input random field is assumed to belong to a Gevrey smoothness class. This approach has the advantage of being substantially more general than models which assume a particular parametric representation of the input random field such as a Karhunen–Loève series expansion. As model problems we consider both the Poisson equation as well as the heat equation and design randomly shifted lattice quasi-Monte Carlo (QMC) cubature rules for the computation of response statistics subject to domain uncertainty. We show that these QMC rules exhibit dimension-independent, faster-than-Monte Carlo cubature convergence rates in this framework. Our theoretical results will be illustrated by numerical examples. This is a joint work with Ana Djurdjevac, Vesa Kaarnioja and Claudia Schillings.
The Fast Newton Transform (FNT) addresses the computational bottleneck that arises in solving high-dimensional problems such as 6d Boltzmann, Fokker-Planck, or Vlaslov equations, multi-body Hamiltonian systems, and the inference of governing equations in complex self-organizing systems. Specifically, the challenge lies in numerically computing function expansions and their derivatives fast, while achieving high approximation power. The FNT is a Newton interpolation algorithm with runtime complexity O(Nnm), where N is the dimension of the downward closed polynomial space, n its degree and m the spatial dimension. We select subgrids from tensorial Leja-ordered Chebyshev-Lobatto grids based on downward-closed sets. This significantly reduces the number of coefficients, N≪(n + 1)m, while achieving optimal geometric approximation rates for a class of analytic functions known as Bos–Levenberg–Trefethen functions. Specifically, we investigate ℓp-multi-index sets, where the Euclidean degree (p=2) turns out to be the pivotal choice for mitigating the curse of dimensionality. Furthermore, the differentiation matrices in Newton basis are sparse, enabling the implementation of fast pseudo-spectral methods on flat spaces, polygonal domains, and regular manifolds. Numerical experiments validate the algorithm’s superior runtime performance over state-of-the-art approaches.
Joint work with Damar Wicaksono & Michael Hecht.
In this talk, we introduce a mesh-free two-level hybrid Tucker tensor format for the approximation of multivariate functions. This new format combines the product Chebyshev interpolation with the ALS-based Tucker decomposition of the coefficients tensor. The benefits of this tensor approximation are two-fold. On the one hand, it allows to avoid the rank-structured approximation of functional tensors defined on large spatial grids, while on the other side, this leads to the Tucker decomposition of the core coefficient tensor with nearly optimal ε-rank parameters, which are shown to be much smaller than both the polynomial degree of the Chebyshev interpolant and the potential number of spatial grid points in the commonly used grid-based discretizations. We discuss the error and complexity estimates of the presented method and demonstrate its efficiency on the demanding example of multi-particle interaction potentials generated by the 3D Newton kernel.
Joint work with Peter Benner, Venera Khoromskaia and Boris N. Khoromskij.
Accurately detecting the coarse-grained coordinates is an essential task in the model discovery of complex systems, such as molecular dynamics. In this work, we show how kernel-based approximation to the Koopman generator – the kgEDMD algorithm – can be used to identify implied timescales and meta stable sets in stochastic dynamical systems, and to learn a coarsegrained dynamics on reduced variables, which retains the essential kinetic properties of the full model. The centerpiece of this study is a learning method to identify a state-dependent effective diffusion in coarse-grained space by leveraging the kgEMD model for the Koopman generator. By combining this method with force matching, a stochastic differential equation governing the effective dynamics can be inferred. Using a two-dimensional model system and molecular dynamics simulation data of alanine dipeptide, we demonstrate that the proposed method successfully and robustly recovers the essential thermodynamic and kinetic properties of the full model.
Machine Learning (ML) has seen great use in the recent years in order to mitigate high cost of solving mathematical models that study reality. Depending on the data they are trained on, ML models can predict either high fidelity (HF) or low fidelity (LF) outputs. The former requires a high amount of computational resources to generate training data while the latter is less accurate. Multifidelity (MF) methods involve the use of statistical techniques to combine HF and LF data resulting in high-accuracy lowcost ML models.
This talk will first introduce the recently developed Multifidelity Machine Learning (MFML) method built with kernel ridge regression aimed at prediction Quantum Chemistry (QC) properties. A cross validation scheme retaining the homogeneous structure of MF training data will be discussed as a model assessment metric. A novel dataadapted optimal combination for MF data, Optimized MFML (o-MFML), will be demonstrated as a robust and accurate successor to MFML. The effect of homogeneous and non-homogeneous MF training data on model accuracy shall be studied for both MFML and o-MFML for the prediction of QC properties including excitation energies which are essential to core life processes such as photosynthesis. Analyses covering data hierarchy and training set structures will be presented as the foundation to the novel high-efficiency Γ-curve technique. Finally, application of these diverse methods for QC such as excitation energy simulations of photosynthetic systems and ground state energy prediction atmospherically relevant molecules will be briefly presented as evidence of real world effectiveness of the herein presented MF methods.
Parameter-dependent partial differential equations (PPDE) are used to model various real-world phenomena. Techniques such as finite volumes, finite elements and finite differences are used to solve the PPDE numerically. Solving the resulting high-dimensional systems of equations often leads to high computational efforts. Here model order reduction methods might help. The basic procedure consists of replacing a high-dimensional model with a reduced model.
A wave equation modelling wave propagation in damaged material is considered. A damage is modelled as a local reduction in the wave propagation velocity; causing the reflection of a propagating wave at the damage. The wave equation is discretized in space and time using finite elements and the Crank-Nicoloson method, respectively.
For the model order reduction of the resulting high-dimensional model, a reduced basis method with a POD-Greedy approach is implemented. The POD-Greedy approach is used creating a reduced basis. Based on the reduced basis and the high-dimensional model, the reduced model is created using the RB-Galerkin method. Model order reduction is difficult for a hyperbolic equation such as the considered wave equation. Therefore, strategies are investigated to simplify the model order reduction process with suitable parameter choices.
We present a continuous and a discontinuous linear Finite Element method based on a predictor-corrector scheme for the numerical approximation of the Ericksen-Leslie equations, a model for nematic liquid crystal flow including a non-convex unit-sphere constraint. As predictor step we propose a linear semi-implicit Finite Element discretization which naturally offers a local orthogonality relation between the approximate director field and its time derivative. Afterwards an explicit discrete projection onto the unit-sphere constraint is applied without increasing the modeled energy. We compare the established method of using a discrete inner product usually referred to as mass-lumping for a globally continuous Finite Element discretization against a piecewise constant discontinuous Galerkin approach. Discrete well-posedness results and energy laws are established. Conditional convergence of subsequences of the approximate solutions to energy-variational solutions of the Ericksen-Leslie equations is shown for a time-step restriction. Computational studies indicate the efficiency of the proposed linearization and the improved accuracy by including a projection step in the algorithm.
Image-based, patient-specific modelling of hemodynamics can improve diagnostic capabilities and provide complementary insights to better understand the hemodynamic treatment outcomes. However, computational fluid dynamics simulations remain relatively costly in a clinical context. Moreover, projection-based reducedorder models and purely data-driven surrogate models struggle due to the high variability of anatomical shapes in a population. A possible solution is shape registration: a reference template geometry is designed from a cohort of available geometries, which can then be diffeomorphically mapped onto it. This provides a natural encoding that can be exploited by machine learning architectures and, at the same time, a reference computational domain in which efficient dimension-reduction strategies can be performed. We compare state-of-the-art graph neural network models with recent data assimilation strategies for the prediction of physical quantities and clinically relevant biomarkers in the context of aortic coarctation.
Joint work with Felipe Galarce, Jia-Jie Zhu, Jan Brüning, Leonid Goubergrits, and Alfonso Caiazzo-
Simulations of fracture propagation are commonly performed using phase field methods, where the discrete fracture is replaced by a continuous variable, the so-called phase field. The evolution of the phase field is then interpreted as a propagation of the fracture. However, the introduction of a continuous variable, that replaces the discrete fracture, causes the appearance of additional parameters. These parameters affect the propagation of fracture and need to be tuned in order to yield satisfactory results.
In this talk, we present an alternative approach to determine brittle fracture propagation which is based on shape optimization algorithms. Here, the fracture is interpreted as part of the boundary of the domain, avoiding the above-mentioned challenges. First, we present the theoretical framework and the optimization problem to be solved. Then, we show numerical results for different benchmark problems in a quasi-static setting.
Joint work with Winnifried Wollner, and Kathrin Welker.
Semiconductor lasers are pivotal components in modern technologies, spanning medical procedures, manufacturing, and autonomous systems like LiDARs. Understanding their operation and developing simulation tools are paramount for advancing such technologies. In this talk, we present a mathematical PDE model for an edge-emitting laser, combining charge transport and light propagation. The charge transport will be described by a drift-diffusion model and the light propagation by the Helmholtz equation. We discretize the drift-diffusion model using the finite volume method and solve it with a Newton solver. Moreover, we discretize the Helmholtz model using the finite element method and solve the resulting matrix eigenvalue problem with ARPACK. We discuss a coupling strategy for both models and showcase numerical simulations for a benchmark example. Our implementation is done in Julia, a language that combines the ease of use of e.g. MATLAB with the performance and efficiency of C++. Finally, we validate our results against a well-established software. We then briefly outline our ongoing efforts to prove the existence of solutions for the coupled PDE system, specifically involving the Helmholtz equation.
Radial basis function finite difference (RBF-FD) discretization has recently emerged as an alternative to classical finite difference or finite element discretization of (systems) of partial differential equations. After an introduction to the RBF-FD I describe how to discretize the steady state Oseen equations in an RBF-FD setting. Particularly, I show different approaches to deal with the pressure constraint. In our numerical results, we focus on RBF-FD discretizations based on polyharmonic splines (PHS) with polynomial augmentation. We illustrate the convergence of the method for different degrees of polynomial augmentation, viscosities and domains. In particular I show why the error in the velocity increases when the viscosity parameter is decresed.
In Bayesian inverse problems, the computation of the posterior distribution can be computationally demanding, especially in many-query settings such as filtering, where a new posterior distribution must be computed many times. When doing inference on a function, the parameter is moreover infinite-dimensional. In this work we consider some computationally efficient approximations of the posterior distribution for linear Gaussian inverse problems defined on separable Hilbert spaces. We measure the quality of these approximations using, among others, the Kullback-Leibler and Rényi divergence of the approximate posterior with respect to the true posterior and investigate their optimality properties. We construct optimal approximations to the posterior mean and covariance separately, and, for the Kullback-Leibler divergence, we also consider optimal joint approximation of mean and covariance. The approximation method exploits low dimensional behaviour of the update from prior to posterior, originating from a combination of prior smoothing, forward smoothing, measurement error and limited number of observations, analogous to the results of [Spantini et al. (SIAM J. Sci. Comput. 2015] for finite dimensional parameter spaces. Since the data is only informative on a low dimensional subspace of the parameter space, the approximation class we consider for the posterior covariance consists of suitable low rank updates of the prior. In the Hilbert space setting, care must be taken, such as when inverting covariance operators. We address such challenges by using the Feldman-Hajek theorem for Gaussian measures.
Reduced rank extrapolation (RRE) is an acceleration method typically used to accelerate the iterative solution of nonlinear systems of equations using a fixed-point process. In this context, the iterates are vectors generated from a fixed-point mapping function. However, when considering the iterative solution of large-scale matrix equations, the iterates are low-rank matrices generated from a fixed-point process for which, generally, the mapping function changes in each iteration. To enable acceleration of the iterative solution for these problems, we propose two novel generalizations of RRE. First, we show how to effectively compute RRE for sequences of low-rank matrices. Second, we derive a formulation of RRE that is suitable for fixed-point processes for which the mapping function changes each iteration. We demonstrate the potential of the methods on several numerical examples involving the iterative solution of large-scale Lyapunov and Riccati matrix equations.
The calculation of exact derivatives is a cornerstone of scientific computing. It enables deep learning, the solution of non-linear partial differential equations, the simulation of models in quantum field theory or even the quantification of the sensitivity of a portfolio value to a changing volatility, to name just a few applications. The calculation of derivatives quickly becomes very challenging when derivatives of complex algorithms or higher derived tensors are needed. Therefore, there is a need for practical algorithmic differentiation (AD) software packages to automate the calculation for the scientific computing community. One of the most comprehensive AD packages is ADOL-C, which is written in C++. Recently, ADOL-C has been extended with an interface to the Julia programming language. This newly developed package ADOLC.jl combines the diverse functionalities of ADOL-C with the ease of use
of Julia. It allows the use of drivers available in ADOL-C for the appropriate calculation of derivation information, even of very high order.
Beside a short introduction to AD, the calculation of derivatives with the help of ADOLC.jl will be shown in this talk.
