Omar Ghattas

Shape optimization under uncertainty (OUU) is computationally intensive for classical PDE-based methods due to the high cost of repeated sampling-based risk evaluation across many uncertainty realizations and varying geometries, while standard neural surrogates often fail to provide accurate and efficient sensitivities for optimization. We introduce Shape-DINO, a derivative-informed neural operator framework for learning PDE solution operators on families of varying geometries, with a particular focus on accelerating PDE-constrained shape OUU. Shape-DINOs encode geometric variability through diffeomorphic mappings to a fixed reference domain and employ a derivative-informed operator learning objective that jointly learns the PDE solution and its Fréchet derivatives with respect to design variables and uncertain parameters, enabling accurate state predictions and reliable gradients for large-scale OUU. We establish a priori error bounds linking surrogate accuracy to optimization error and prove universal approximation results for multi-input reduced basis neural operators in suitable $C^1$ norms. We demonstrate efficiency and scalability on three representative shape OUU problems, including boundary design for a Poisson equation and shape design governed by steady-state Navier-Stokes exterior flows in two and three dimensions. Across these examples, Shape-DINOs produce more reliable optimization results than operator surrogates trained without derivative information. In our examples, Shape-DINOs achieve 3-8 orders-of-magnitude speedups in state and gradient evaluations. Counting training data generation, Shape-DINOs reduce necessary PDE solves by 1-2 orders-of-magnitude compared to a strictly PDE-based approach for a single OUU problem. Moreover, Shape-DINO construction costs can be amortized across many objectives and risk measures, enabling large-scale shape OUU for complex systems.

Shape optimization under uncertainty (OUU) is computationally intensive for classical PDE-based methods due to the high cost of repeated sampling-based risk evaluation across many uncertainty realizations and varying geometries, while standard neural surrogates often fail to provide accurate and efficient sensitivities for optimization. We introduce Shape-DINO, a derivative-informed neural operator framework for learning PDE solution operators on families of varying geometries, with a particular focus on accelerating PDE-constrained shape OUU. Shape-DINOs encode geometric variability through diffeomorphic mappings to a fixed reference domain and employ a derivative-informed operator learning objective that jointly learns the PDE solution and its Fréchet derivatives with respect to design variables and uncertain parameters, enabling accurate state predictions and reliable gradients for large-scale OUU. We establish a priori error bounds linking surrogate accuracy to optimization error and prove universal approximation results for multi-input reduced basis neural operators in suitable $C^1$ norms. We demonstrate efficiency and scalability on three representative shape OUU problems, including boundary design for a Poisson equation and shape design governed by steady-state Navier-Stokes exterior flows in two and three dimensions. Across these examples, Shape-DINOs produce more reliable optimization results than operator surrogates trained without derivative information. In our examples, Shape-DINOs achieve 3-8 orders-of-magnitude speedups in state and gradient evaluations. Counting training data generation, Shape-DINOs reduce necessary PDE solves by 1-2 orders-of-magnitude compared to a strictly PDE-based approach for a single OUU problem. Moreover, Shape-DINO construction costs can be amortized across many objectives and risk measures, enabling large-scale shape OUU for complex systems.

We present a goal-oriented framework for constructing digital twins with the following properties: (1) they employ discretizations of high-fidelity PDE models governed by autonomous dynamical systems, leading to large-scale forward problems; (2) they solve a linear inverse problem to assimilate observational data to infer uncertain model components followed by a forward prediction of the evolving dynamics; and (3) the entire end-to-end, data-to-inference-to-prediction computation is carried out without approximation and in real time through a Bayesian framework that rigorously accounts for uncertainties. Several challenges must be overcome to realize this framework, including the large scale of the forward problem, the high dimensionality of the parameter space, and for a class of problems including those we target, the slow decay of the singular values of the parameter-to-observable map. Here we introduce a methodology to overcome these challenges by exploiting the autonomous structure of the forward model to decompose the solution of the inverse problem into an offline phase in which the PDE model is solved a limited number of times, and an online phase that computes the parameter inference and prediction of quantities of interest in real time, given observational data. Our goal is to apply this framework to construct digital twins for subduction zones to provide early warning for tsunamis. To this end, we show how our methodology can be used to employ seafloor pressure observations, along with the coupled acoustic-gravity wave equations, to infer the earthquake-induced seafloor motion (discretized with 10^9 parameters) and forecast the tsunami propagation. We present results of an end-to-end inference, prediction, and uncertainty quantification for a representative test problem with 10^8 inversion parameters for which goal-oriented Bayesian inference is accomplished in real time.

This work expands on a strain-based aerodynamic sensing strategy for hypersonics to account for nonlinear temperature effects in real time. The sensing strategy uses sparse strain observations to infer the aerodynamic pressure loads, which is mathematically posed as a partial differential equation (PDE)-constrained inverse problem. In previous work, this inverse problem was shown to have a closed-form solution, where offline computation of the operations requiring the PDE solution was exploited to enable real-time evaluation speeds. In this work, the temperature effects preclude the offline pre-computation acceleration because the PDE operator is nonlinearly temperature dependent. To address this challenge, the recently developed neural matrix operator (NEMO) approach is employed to account for the temperature dependence. The NEMO method explicitly incorporates the physics structure of the governing equations, and thus preserves the availability of a closed-form inverse solution that can be computed rapidly onboard the vehicle. This work further considers the tasks of estimating the temperature field from sparse temperature measurements, and compensation for the thermal strain. The overall performance is demonstrated on the Initial Concept 3.X hypersonic vehicle. The results show strong approximation performance of NEMO and corresponding inverse solution accuracy, but further work is necessary to reduce errors in the thermal strain compensation.

We present a goal-oriented framework for constructing digital twins with the following properties: (1) they employ discretizations of high-fidelity PDE models governed by autonomous dynamical systems, leading to large-scale forward problems; (2) they solve a linear inverse problem to assimilate observational data to infer uncertain model components followed by a forward prediction of the evolving dynamics; and (3) the entire end-to-end, data-to-inference-to-prediction computation is carried out without approximation and in real time through a Bayesian framework that rigorously accounts for uncertainties. Several challenges must be overcome to realize this framework, including the large scale of the forward problem, the high dimensionality of the parameter space, and for a class of problems including those we target, the slow decay of the singular values of the parameter-to-observable map. Here we introduce a methodology to overcome these challenges by exploiting the autonomous structure of the forward model to decompose the solution of the inverse problem into an offline phase in which the PDE model is solved a limited number of times, and an online phase that computes the parameter inference and prediction of quantities of interest in real time, given observational data. Our goal is to apply this framework to construct digital twins for subduction zones to provide early warning for tsunamis. To this end, we show how our methodology can be used to employ seafloor pressure observations, along with the coupled acoustic-gravity wave equations, to infer the earthquake-induced seafloor motion (discretized with 10^9 parameters) and forecast the tsunami propagation. We present results of an end-to-end inference, prediction, and uncertainty quantification for a representative test problem with 10^8 inversion parameters for which goal-oriented Bayesian inference is accomplished in real time.

This study of the effects of local geological conditions on seismic ground motion uses 1D amplification as a reference point and examines, via simple theoretical and more realistic numerical examples and observations, how 2D and 3D conditions differ from 1D estimations. Because 1D simulations cannot model basin and edge effects, 1D response tends, in general, to exhibit lower peaks and be of shorter duration than 2D and 3D results. On the other hand, due to destructive interference of different types of waves, there are sites where the response can be much smaller than predicted by 1D models.

We present approximation theories and efficient training methods for derivative-informed Fourier neural operators (DIFNOs) with applications to PDE-constrained optimization. A DIFNO is an FNO trained by minimizing its prediction error jointly on output and Fréchet derivative samples of a high-fidelity operator (e.g., a parametric PDE solution operator). As a result, a DIFNO can closely emulate not only the high-fidelity operator's response but also its sensitivities. To motivate the use of DIFNOs instead of conventional FNOs as surrogate models, we show that accurate surrogate-driven PDE-constrained optimization requires accurate surrogate Fréchet derivatives. Then, we establish (i) simultaneous universal approximation of continuously differentiable operators and their Fréchet derivatives by FNOs on compact sets, and (ii) universal approximation of continuously differentiable operators by FNOs in weighted Sobolev spaces with input measures that have unbounded supports. Our theoretical results certify the capability of FNOs for accurate derivative-informed operator learning and for the solution of PDE-constrained optimization problems. Furthermore, we develop efficient training schemes that leverage dimensionality reduction and multi-resolution techniques to significantly reduce memory and computational costs in Fréchet derivative learning. Numerical examples on nonlinear diffusion--reaction, Helmholtz, and Navier--Stokes equations demonstrate that DIFNOs are superior in sample complexity for operator learning and solving infinite-dimensional PDE-constrained inverse problems, achieving high accuracy at low training sample sizes.

We present a Bayesian inversion-based digital twin that employs acoustic pressure data from seafloor sensors, along with 3D coupled acoustic–gravity wave equations, to infer earthquake-induced spatiotemporal seafloor motion in real time and forecast tsunami propagation toward coastlines for early warning with quantified uncertainties. Our target is the Cascadia subduction zone, with one billion parameters. Computing the posterior mean alone would require 50 years on a 512 GPU machine. Instead, exploiting the shift invariance of the parameter-to-observable map and devising novel parallel algorithms, we induce a fast offline–online decomposition. The offline component requires just one adjoint wave propagation per sensor; using MFEM, we scale this part of the computation to the full El Capitan system (43,520 GPUs) with 92% weak parallel efficiency. Moreover, given real-time data, the online component exactly solves the Bayesian inverse and forecasting problems in 0.2 seconds on a modest GPU system, a ten-billion-fold speedup.

The hardware diversity in leadership-class computing facilities, alongside the immense performance boosts from today’s GPUs when computing in lower precision, incentivizes scientific HPC workflows to adopt mixed-precision algorithms and performance portability models. We present an on-the-fly framework using hipify for performance portability and apply it to FFTMatvec—an HPC application that computes matrix-vector products with block-triangular Toeplitz matrices. Our approach enables FFTMatvec, initially a CUDA-only application, to run seamlessly on AMD GPUs with excellent performance. Performance optimizations for AMD GPUs are integrated into the open-source rocBLAS library, keeping the application code unchanged. We then present a dynamic mixed-precision framework for FFTMatvec; a Pareto front analysis determines the optimal mixed-precision configuration for a desired error tolerance. Results are shown for AMD Instinct™ MI250X, MI300X, and the newly launched MI355X GPUs. The performance-portable, mixed-precision FFTMatvec is scaled to 4,096 GPUs on the OLCF Frontier supercomputer.

We present an efficient and scalable algorithm for performing matrix-vector multiplications ("matvecs") for block Toeplitz matrices. Such matrices, which are shift-invariant with respect to their blocks, arise in the context of solving inverse problems governed by autonomous systems, and time-invariant systems in particular. In this article, we consider inverse problems that infer unknown parameters from observational data of a linear time-invariant dynamical system given in the form of partial differential equations (PDEs). Matrix-free Newton-conjugate-gradient methods are often the gold standard for solving these inverse problems, but they require numerous actions of the Hessian on a vector. Matrix-free adjoint-based Hessian matvecs require solution of a pair of linearized forward/adjoint PDE solves per Hessian action, which may be prohibitive for large-scale inverse problems. Time invariance of the forward PDE problem leads to a block Toeplitz structure of the discretized parameter-to-observable (p2o) map defining the mapping from inputs (parameters) to outputs (observables) of the PDEs. This block Toeplitz structure enables us to exploit two key properties: (1) compact storage of the p2o map and its adjoint; and (2) efficient fast Fourier transform (FFT)-based Hessian matvecs. The proposed algorithm is mapped onto large multi-GPU clusters and achieves more than 80 percent of peak bandwidth on NVIDIA A100 GPUs. Excellent weak scaling is shown for up to 48 A100 GPUs. For the targeted problems, the implementation executes Hessian matvecs within fractions of a second, which is orders of magnitude faster than can be achieved by conventional matrix-free Hessian matvecs via forward/adjoint PDE solves.

Design problems in aerospace engineering often require numerous evaluations of expensive-to-evaluate high-fidelity models, resulting in prohibitive computational costs. One way to address the computational cost is through building surrogates, such as deep neural networks (DNNs). However, DNNs may only be an effective surrogate when sufficient evaluations of the high-fidelity model are required such that the up-front training cost is amortized, or in situations that require real-time responses (such as interactive visualizations). Typically, the data requirements for adequately accurate training of DNNs are often impractical for engineering applications. To alleviate this issue, the proposed work utilizes output dimensionality reduction along with information from multiple models of varying fidelities and cost to develop accurate projection-enabled multifidelity neural networks (MF-NNs) with limited training samples. The dimensionality reduction leads to a more parsimonious network and the multifidelity aspect adds more training data from lower-cost, lower-fidelity models. Three approaches for MF-NNs that leverage proper orthogonal decomposition based projections are introduced: (i) pre-training method, (ii) additive method, and (iii) multi-step method. The MF-NN is applied to approximate the optimal design of 2D aerodynamic airfoils given the performance and design requirements. The MF-NN leads to ~27% computational cost reduction compared to single-fidelity neural networks at the same accuracy (90%), with the multi-step approach performing the best for this application.

Sustained advances in the mathematics of modeling and simulation have resulted in the capability today for routine simulation of a number of large scale complex DOE-relevant systems. As remarkable as this capability for solving the so-called forward problem is, it is typically only the first step-an inner loop within an outer loop that explores the simulation model's parameter space and decision space to characterize uncertainty in the model's predictions, learn unknown model parameters from data, design the most informative experiments, determine optimal control strategies, and create optimal designs. Broadly, what unifies all of these outer loop problems is that they are, in one form or another, optimization problems over parameter/control/design space that are constrained by complex uncertain models. To fully realize the power of scientific simulation as a basis for scientific discovery, technological innovation, and rational decision-making, it is imperative to move beyond simulation to tackle the outer loop of optimization for learning from data, experimental design, and control with complex uncertain models. When the models under consideration are large-scale and complex, and when the optimization variable and uncertain parameter spaces are high (or infinite) dimensional, this constitutes a grand challenge of the highest order, and is intractable with conventional methods. To overcome these challenges, the AEOLUS Center was established to develop a unified mathematical, computational, and statistical framework for (1) Learning predictive models from complex data via Bayesian inference and optimization, and (2) Optimizing experiments, processes, and designs using the resulting uncertain models. These problems are intractable with conventional methods, for several reasons: (1) The simulation problems that govern the inner loops of the optimization problems are expensive to execute (due to severe nonlinearity, heterogeneity, multiphysics/multiscale coupling); (2) The optimization variable and uncertain parameter spaces are high dimensional, often stemming from discretizations of infinite dimensional fields such as initial conditions, sources, or material properties. We argue that the key to overcoming these challenges is to develop new mathematical, computational, and statistical methods that exploit the structure of the Bayesian inference and optimization problems mediated by their underlying complex uncertain models. This structure includes the regularity, sparsity, geometry, low intrinsic dimensionality, and multifidelity nature of the maps from uncertain parameter/optimization variable spaces to the specific objectives targeted: Bayesian inference, optimal experimental design, and optimal control design. Black box methods developed as generic tools are incapable of exploiting this structure. To be successful, we must create, integrate, and cross-fertilize ideas across multiple areas of applied math--including approximation theory, Bayesian inference, data science, experimental design, information theory, machine learning, model reduction, optimal control theory, parallel algorithms, PDE-constrained optimization, randomized algorithms, stochastic optimization, and uncertainty quantification--all while exploiting the structure of the problems at hand. With this goal in mind, we have marshaled a team of leading authorities in these areas. While the methods we develop will be broadly applicable across a wide spectrum of DOE problems in which experiments inform models and the systems those models describe must be optimized under uncertainty, we have chosen a specific area, advanced manufacturing and materials, to drive our work. AMM is characterized by complex models across multiple scales, and is a rich source of challenging problems in inference, experimental design, and optimal control, requiring multifaceted and integrated advances in applied mathematics. As such, AMM serves as an excellent vehicle to motivate and demonstrate the advances in applied mathematics developed by our center.

We present LazyDINO, a transport map variational inference method for fast, scalable, and efficiently amortized solutions of high-dimensional nonlinear Bayesian inverse problems with expensive parameter-to-observable (PtO) maps. Our method consists of an offline phase in which we construct a derivative-informed neural surrogate of the PtO map using joint samples of the PtO map and its Jacobian. During the online phase, when given observational data, we seek rapid posterior approximation using surrogate-driven training of a lazy map [Brennan et al., NeurIPS, (2020)], i.e., a structure-exploiting transport map with low-dimensional nonlinearity. The trained lazy map then produces approximate posterior samples or density evaluations. Our surrogate construction is optimized for amortized Bayesian inversion using lazy map variational inference. We show that (i) the derivative-based reduced basis architecture [O'Leary-Roseberry et al., Comput. Methods Appl. Mech. Eng., 388 (2022)] minimizes the upper bound on the expected error in surrogate posterior approximation, and (ii) the derivative-informed training formulation [O'Leary-Roseberry et al., J. Comput. Phys., 496 (2024)] minimizes the expected error due to surrogate-driven transport map optimization. Our numerical results demonstrate that LazyDINO is highly efficient in cost amortization for Bayesian inversion. We observe one to two orders of magnitude reduction of offline cost for accurate posterior approximation, compared to simulation-based amortized inference via conditional transport and conventional surrogate-driven transport. In particular, LazyDINO outperforms Laplace approximation consistently using fewer than 1000 offline samples, while other amortized inference methods struggle and sometimes fail at 16,000 offline samples.

Estimating the parameters of compact binaries which coalesce and produce gravitational waves is a challenging Bayesian inverse problem. Gravitational-wave parameter estimation lies within the class of multifidelity problems, where a variety of models with differing assumptions, levels of fidelity, and computational cost are available for use in inference. In an effort to accelerate the solution of a Bayesian inverse problem, cheaper surrogates for the best models may be used to reduce the cost of likelihood evaluations when sampling the posterior. Importance sampling can then be used to reweight these samples to represent the true target posterior, incurring a reduction in the effective sample size. In cases when the problem is high dimensional, or when the surrogate model produces a poor approximation of the true posterior, this reduction in effective samples can be dramatic and render multifidelity importance sampling ineffective. We propose a novel method of tempered multifidelity importance sampling in order to remedy this issue. With this method the biasing distribution produced by the low-fidelity model is tempered, changing the chi-squared divergence between the two distributions and thereby affecting the efficiency of importance sampling. There is an optimal temperature that maximizes the efficiency in this setting, and we propose a low-cost strategy for approximating this optimal temperature using samples from the untempered distribution. In this paper, we motivate this method by applying it to Gaussian target and biasing distributions. Finally, we apply it to a series of problems in gravitational wave parameter estimation and demonstrate improved efficiencies when applying the method to real gravitational wave detections.

This work develops an efficient real-time inverse formulation for inferring the aerodynamic surface pressures on a hypersonic vehicle from sparse measurements of the structural strain. The approach aims to provide real-time estimates of the aerodynamic loads acting on the vehicle for ground and flight testing, as well as guidance, navigation, and control applications. Specifically, the approach targets hypersonic flight conditions where direct measurement of the surface pressures is challenging due to the harsh aerothermal environment. For problems employing a linear elastic structural model, the inference problem can be posed as a least-squares problem with a linear constraint arising from a finite element discretization of the governing elasticity partial differential equation. Due to the linearity of the problem, an explicit solution is given by the normal equations. Precomputation of the resulting inverse map enables rapid evaluation of the surface pressure and corresponding integrated quantities, such as the force and moment coefficients. The inverse approach additionally allows for uncertainty quantification, providing insights for theoretical recoverability and robustness to sensor noise. Numerical studies demonstrate the estimator performance for reconstructing the surface pressure field, as well as the force and moment coefficients, for the Initial Concept 3.X (IC3X) conceptual hypersonic vehicle.

This work develops an efficient real-time inverse formulation for inferring the aerodynamic surface pressures on a hypersonic vehicle from sparse measurements of the structural strain. The approach aims to provide real-time estimates of the aerodynamic loads acting on the vehicle for ground and flight testing, as well as guidance, navigation, and control applications. Specifically, the approach targets hypersonic flight conditions where direct measurement of the surface pressures is challenging due to the harsh aerothermal environment. For problems employing a linear elastic structural model, we show that the inference problem can be posed as a least-squares problem with a linear constraint arising from a finite element discretization of the governing elasticity partial differential equation. Due to the linearity of the problem, an explicit solution is given by the normal equations. Pre-computation of the resulting inverse map enables rapid evaluation of the surface pressure and corresponding integrated quantities, such as the force and moment coefficients. The inverse approach additionally allows for uncertainty quantification, providing insights for theoretical recoverability and robustness to sensor noise. Numerical studies demonstrate the estimator performance for reconstructing the surface pressure field, as well as the force and moment coefficients, for the Initial Concept 3.X (IC3X) conceptual hypersonic vehicle.

Computational models governed by partial differential equations (PDEs) are frequently used by engineers to optimize the performance of various physical systems through decisions relating to their configuration (optimal design) and operation (optimal control).However, the ability to make optimal choices is often hindered by uncertainty, such as uncertainty in model parameters (e.g., material properties) and operating conditions (e.g., forces on a structure).The need to account for these uncertainties in order to arrive at robust and risk-informed decisions thus gives rise to problems of optimization under uncertainty (OUU) (D.

Rock strength has long been linked to lithospheric deformation and seismicity. However, independent constraints on the related elastic heterogeneity are missing, yet could provide key information for solid Earth dynamics. Using coseismic Global Navigation Satellite Systems (GNSS) data for the 2011 M9 Tohoku-oki earthquake in Japan, we apply an inverse method to infer elastic structure and fault slip simultaneously. We find compliant material beneath the volcanic arc and in the mantle wedge within the partial melt generation zone inferred to lie above ~100 km slab depth. We also identify low-rigidity material closer to the trench matching seismicity patterns, likely associated with accretionary wedge structure. Along with traditional seismic and electromagnetic methods, our approach opens up avenues for multiphysics inversions. Those have the potential to advance earthquake and volcano science, and in particular once expanded to InSAR type constraints, may lead to a better understanding of transient lithospheric deformation across scales.

Estimating the parameters of compact binaries which coalesce and produce gravitational waves is a challenging Bayesian inverse problem. Gravitational-wave parameter estimation lies within the class of multifidelity problems, where a variety of models with differing assumptions, levels of fidelity, and computational cost are available for use in inference. In an effort to accelerate the solution of a Bayesian inverse problem, cheaper surrogates for the best models may be used to reduce the cost of likelihood evaluations when sampling the posterior. Importance sampling can then be used to reweight these samples to represent the true target posterior, incurring a reduction in the effective sample size. In cases when the problem is high dimensional, or when the surrogate model produces a poor approximation of the true posterior, this reduction in effective samples can be dramatic and render multifidelity importance sampling ineffective. We propose a novel method of tempered multifidelity importance sampling in order to remedy this issue. With this method the biasing distribution produced by the low-fidelity model is tempered, allowing for potentially better overlap with the target distribution. There is an optimal temperature which maximizes the efficiency in this setting, and we propose a low-cost strategy for approximating this optimal temperature using samples from the untempered distribution. In this paper, we motivate this method by applying it to Gaussian target and biasing distributions. Finally, we apply it to a series of problems in gravitational wave parameter estimation and demonstrate improved efficiencies when applying the method to real gravitational wave detections.

We propose an operator learning approach to accelerate geometric Markov chain Monte Carlo (MCMC) for solving infinite-dimensional Bayesian inverse problems (BIPs). While geometric MCMC employs high-quality proposals that adapt to posterior local geometry, it requires repeated computations of gradients and Hessians of the log-likelihood, which becomes prohibitive when the parameter-to-observable (PtO) map is defined through expensive-to-solve parametric partial differential equations (PDEs). We consider a delayed-acceptance geometric MCMC method driven by a neural operator surrogate of the PtO map, where the proposal exploits fast surrogate predictions of the log-likelihood and, simultaneously, its gradient and Hessian. To achieve a substantial speedup, the surrogate must accurately approximate the PtO map and its Jacobian, which often demands a prohibitively large number of PtO map samples via conventional operator learning methods. In this work, we present an extension of derivative-informed operator learning [O'Leary-Roseberry et al., J. Comput. Phys., 496 (2024)] that uses joint samples of the PtO map and its Jacobian. This leads to derivative-informed neural operator (DINO) surrogates that accurately predict the observables and posterior local geometry at a significantly lower training cost than conventional methods. Cost and error analysis for reduced basis DINO surrogates are provided. Numerical studies demonstrate that DINO-driven MCMC generates effective posterior samples 3--9 times faster than geometric MCMC and 60--97 times faster than prior geometry-based MCMC. Furthermore, the training cost of DINO surrogates breaks even compared to geometric MCMC after just 10--25 effective posterior samples.

We consider the Bayesian calibration of models describing the phenomenon of block copolymer (BCP) self-assembly using image data produced by microscopy or X-ray scattering techniques. To account for the random long-range disorder in BCP equilibrium structures, we introduce auxiliary variables to represent this aleatory uncertainty. These variables, however, result in an integrated likelihood for high-dimensional image data that is generally intractable to evaluate. We tackle this challenging Bayesian inference problem using a likelihood-free approach based on measure transport together with the construction of summary statistics for the image data. We also show that expected information gains (EIGs) from the observed data about the model parameters can be computed with no significant additional cost. Lastly, we present a numerical case study based on the Ohta--Kawasaki model for diblock copolymer thin film self-assembly and top-down microscopy characterization. For calibration, we introduce several domain-specific energy- and Fourier-based summary statistics, and quantify their informativeness using EIG. We demonstrate the power of the proposed approach to study the effect of data corruptions and experimental designs on the calibration results.

This work develops an efficient inverse formulation for inferring the aerodynamic surface pressures on a hypersonic vehicle from sparse measurements of the structural strain. The approach aims to provide accurate estimates of the aerodynamic loads acting on the vehicle for real-time guidance, navigation, and control with quantifiable uncertainty. Specifically, the approach targets hypersonic flight conditions where direct measurement of the surface pressures is challenging due to the harsh aerothermal environment. For problems employing a linear elastic structural model, we show that the inference problem can be posed as a least-squares problem with a linear constraint arising from a finite element discretization of the governing elasticity partial differential equation. Due to the linearity of the problem, an explicit solution is given by the normal equations. Pre-computation of resulting inverse map enables rapid evaluation of the surface pressure and corresponding integrated quantities, such as the force and moment coefficients. The inverse approach additionally allows for uncertainty quantification, providing insights for optimal sensor placement and robustness to sensor noise. Numerical studies demonstrate the estimator performance for reconstructing the surface pressure field, as well as the force and moment coefficients, for the Initial Concept 3.X (IC3X) conceptual hypersonic vehicle.

SUMMARY Analysis of tectonic and earthquake-cycle associated deformation of the crust can provide valuable insights into the underlying deformation processes including fault slip. How those processes are expressed at the surface depends on the lateral and depth variations of rock properties. The effect of such variations is often tested by forward models based on a priori geological or geophysical information. Here, we first develop a novel technique based on an open-source finite-element computational framework to invert geodetic constraints directly for heterogeneous media properties. We focus on the elastic, coseismic problem and seek to constrain variations in shear modulus and Poisson’s ratio, proxies for the effects of lithology and/or temperature and porous flow, respectively. The corresponding nonlinear inversion is implemented using adjoint-based optimization that efficiently reduces the cost function that includes the misfit between the calculated and observed displacements and a penalty term. We then extend our theoretical and numerical framework to simultaneously infer both heterogeneous Earth’s structure and fault slip from surface deformation. Based on a range of 2-D synthetic cases, we find that both model parameters can be satisfactorily estimated for the megathrust setting-inspired test problems considered. Within limits, this is the case even in the presence of noise and if the fault geometry is not perfectly known. Our method lays the foundation for a future reassessment of the information contained in increasingly data-rich settings, for example, geodetic GNSS constraints for large earthquakes such as the 2011 Tohoku-oki M9 event, or distributed deformation along plate boundaries as constrained from InSAR.

In this paper we propose a variant of enriched Galerkin methods for second order elliptic equations with over-penalization of interior jump terms. The bilinear form with interior over-penalization gives a non-standard norm which is different from the discrete energy norm in the classical discontinuous Galerkin methods. Nonetheless we prove that optimal a priori error estimates with the standard discrete energy norm can be obtained by combining a priori and a posteriori error analysis techniques. We also show that the interior over-penalization is advantageous for constructing preconditioners robust to mesh refinement by analyzing spectral equivalence of bilinear forms. Numerical results are included to illustrate the convergence and preconditioning results.

Analysis of tectonic and earthquake-cycle associated deformation of the crust can provide valuable insights into the underlying deformation processes including fault slip. How those processes are expressed at the surface depends on the lateral and depth variations of rock properties. The effect of such variations is often tested by forward models based on a priori geological or geophysical information. Here, we first develop a novel technique based on an open-source finite-element computational framework to invert geodetic constraints directly for heterogeneous media properties. We focus on the elastic, coseismic problem and seek to constrain variations in shear modulus and Poisson's ratio, proxies for the effects of lithology and/or temperature and porous flow, respectively. The corresponding non-linear inversion is implemented using adjoint-based optimization that efficiently reduces the cost functional that includes the misfit between the calculated and observed displacements and a penalty term. We then extend our theoretical and numerical framework to simultaneously infer both heterogeneous Earth's structure and fault slip from surface deformation. Based on a range of 2-D synthetic cases, we find that both model parameters can be satisfactorily estimated for the megathrust setting-inspired test problems considered. Within limits, this is the case even in the presence of noise and if the fault geometry is not perfectly known. Our method lays the foundation for a future reassessment of the information contained in increasingly data-rich settings, e.g. geodetic GNSS constraints for large earthquakes such as the 2011 Tohoku-oki M9 event, or distributed deformation along plate boundaries as constrained from InSAR.

We present an efficient matrix-free point spread function (PSF) method for approximating operators that have locally supported non-negative integral kernels. The method computes impulse responses at scattered points, and interpolates these impulse responses to approximate integral kernel entries. Impulse responses are computed by applying the operator to Dirac comb batches of point sources, which are chosen via an ellipsoid packing procedure. Evaluation of kernel entries allows us to construct a hierarchical matrix approximation of the operator, which is used for further matrix computations. We illustrate the end-to-end method on a blur problem, then use the method to build preconditioners for the Hessian in two inverse problems governed by partial differential equations (PDEs): inversion for the basal friction coefficient in an ice sheet flow problem and for the initial condition in an advective-diffusive transport problem. While for many ill-posed inverse problems the Hessian of the data misfit term exhibits a low rank structure, and hence a low rank approximation is suitable, for many problems of practical interest the numerical rank of the Hessian is still large. But Hessian impulse responses typically become more local as the numerical rank increases, which benefits the PSF method. Numerical results reveal that the PSF preconditioner clusters the spectrum of the preconditioned Hessian near one, yielding roughly 5x-10x reductions in the required number of PDE solves, as compared to regularization preconditioning and no preconditioning. We also present a numerical study for the influence of various parameters (that control the shape of the impulse responses) on the effectiveness of the advection-diffusion Hessian approximation. The results show that the PSF-based preconditioners are able to form good approximations of high-rank Hessians using a small number of operator applications.

Identifying parameters of computational models from experimental data, or model calibration, is fundamental for assessing and improving the predictability and reliability of computer simulations. In this work, we propose a method for Bayesian calibration of models that predict morphological patterns of diblock copolymer (Di-BCP) thin film self-assembly while accounting for various sources of uncertainties in pattern formation and data acquisition. This method extracts the azimuthally-averaged power spectrum (AAPS) of the top-down microscopy characterization of Di-BCP thin film patterns as summary statistics for Bayesian inference of model parameters via the pseudo-marginal method. We derive the analytical and approximate form of a conditional likelihood for the AAPS of image data. We demonstrate that AAPS-based image data reduction retains the mutual information, particularly on important length scales, between image data and model parameters while being relatively agnostic to the aleatoric uncertainties associated with the random long-range disorder of Di-BCP patterns. Additionally, we propose a phase-informed prior distribution for Bayesian model calibration. Furthermore, reducing image data to AAPS enables us to efficiently build surrogate models to accelerate the proposed Bayesian model calibration procedure. We present the formulation and training of two multi-layer perceptrons for approximating the parameter-to-spectrum map, which enables fast integrated likelihood evaluations. We validate the proposed Bayesian model calibration method through numerical examples, for which the neural network surrogate delivers a fivefold reduction of the number of model simulations performed for a single calibration task.

We propose a novel machine learning framework for solving optimization problems governed by large-scale partial differential equations (PDEs) with high-dimensional random parameters. Such optimization under uncertainty (OUU) problems may be computational prohibitive using classical methods, particularly when a large number of samples is needed to evaluate risk measures at every iteration of an optimization algorithm, where each sample requires the solution of an expensive-to-solve PDE. To address this challenge, we propose a new neural operator approximation of the PDE solution operator that has the combined merits of (1) accurate approximation of not only the map from the joint inputs of random parameters and optimization variables to the PDE state, but also its derivative with respect to the optimization variables, (2) efficient construction of the neural network using reduced basis architectures that are scalable to high-dimensional OUU problems, and (3) requiring only a limited number of training data to achieve high accuracy for both the PDE solution and the OUU solution. We refer to such neural operators as multi-input reduced basis derivative informed neural operators (MR-DINOs). We demonstrate the accuracy and efficiency our approach through several numerical experiments, i.e. the risk-averse control of a semilinear elliptic PDE and the steady state Navier--Stokes equations in two and three spatial dimensions, each involving random field inputs. Across the examples, MR-DINOs offer $10^{3}$--$10^{7} \times$ reductions in execution time, and are able to produce OUU solutions of comparable accuracies to those from standard PDE based solutions while being over $10 \times$ more cost-efficient after factoring in the cost of construction.

We explore using neural operators, or neural network representations of nonlinear maps between function spaces, to accelerate infinite-dimensional Bayesian inverse problems (BIPs) with models governed by nonlinear parametric partial differential equations (PDEs). Neural operators have gained significant attention in recent years for their ability to approximate the parameter-to-solution maps defined by PDEs using as training data solutions of PDEs at a limited number of parameter samples. The computational cost of BIPs can be drastically reduced if the large number of PDE solves required for posterior characterization are replaced with evaluations of trained neural operators. However, reducing error in the resulting BIP solutions via reducing the approximation error of the neural operators in training can be challenging and unreliable. We provide an a priori error bound result that implies certain BIPs can be ill-conditioned to the approximation error of neural operators, thus leading to inaccessible accuracy requirements in training. To reliably deploy neural operators in BIPs, we consider a strategy for enhancing the performance of neural operators, which is to correct the prediction of a trained neural operator by solving a linear variational problem based on the PDE residual. We show that a trained neural operator with error correction can achieve a quadratic reduction of its approximation error, all while retaining substantial computational speedups of posterior sampling when models are governed by highly nonlinear PDEs. The strategy is applied to two numerical examples of BIPs based on a nonlinear reaction--diffusion problem and deformation of hyperelastic materials. We demonstrate that posterior representations of the two BIPs produced using trained neural operators are greatly and consistently enhanced by error correction.

Directed self-assembly (DSA) of block-copolymers (BCPs) is one of the most promising developments in the cost-effective production of nanoscale devices. The process makes use of the natural tendency for BCP mixtures to form nanoscale structures upon phase separation. The phase separation can be directed through the use of chemically patterned substrates to promote the formation of morphologies that are essential to the production of semiconductor devices. Moreover, the design of substrate pattern can formulated as an optimization problem for which we seek optimal substrate designs that effectively produce given target morphologies. In this paper, we adopt a phase field model given by a nonlocal Cahn--Hilliard partial differential equation (PDE) based on the minimization of the Ohta--Kawasaki free energy, and present an efficient PDE-constrained optimization framework for the optimal design problem. The design variables are the locations of circular- or strip-shaped guiding posts that are used to model the substrate chemical pattern. To solve the ensuing optimization problem, we propose a variant of an inexact Newton conjugate gradient algorithm tailored to this problem. We demonstrate the effectiveness of our computational strategy on numerical examples that span a range of target morphologies. Owing to our second-order optimizer and fast state solver, the numerical results demonstrate five orders of magnitude reduction in computational cost over previous work. The efficiency of our framework and the fast convergence of our optimization algorithm enable us to rapidly solve the optimal design problem in not only two, but also three spatial dimensions.

We address the solution of large-scale Bayesian optimal experimental design (OED) problems governed by partial differential equations (PDEs) with infinite-dimensional parameter fields. The OED problem seeks to find sensor locations that maximize the expected information gain (EIG) in the solution of the underlying Bayesian inverse problem. Computation of the EIG is usually prohibitive for PDE-based OED problems. To make the evaluation of the EIG tractable, we approximate the (PDE-based) parameter-to-observable map with a derivative-informed projected neural network (DIPNet) surrogate, which exploits the geometry, smoothness, and intrinsic low-dimensionality of the map using a small and dimension-independent number of PDE solves. The surrogate is then deployed within a greedy algorithm-based solution of the OED problem such that no further PDE solves are required. We analyze the EIG approximation error in terms of the generalization error of the DIPNet and show they are of the same order. Finally, the efficiency and accuracy of the method are demonstrated via numerical experiments on OED problems governed by inverse scattering and inverse reactive transport with up to 16,641 uncertain parameters and 100 experimental design variables, where we observe up to three orders of magnitude speedup relative to a reference double loop Monte Carlo method.

Analysis of geodetic and seismological data helps constrain earthquake dynamics and the physics of lithospheric deformation. Here, we discuss a new modeling approach based on an open-source finite-element framework to invert surface deformation data for constitutive laws and their parameters, such as the Poisson’s ratio or shear modulus in the crust and mantle wedge.These inversions can be realized by using adjoint-based optimization methods which efficiently reduce the misfit between the calculated and observed displacements. To quantify the associated model uncertainties, we extend the inverse approach to a Bayesian inference problem. Since the data are usually informative only in a few directions in parameter space, we use a low-rank Laplace approximation of the posterior distribution to make the inverse problem computationally tractable. The mean and the posterior covariance are approximated by the solution of the inverse problem (MAP point) and the inverse of the Hessian of the negative log posterior evaluated at the MAP point, respectively. We show how smoothly varying parameter fields can be reconstructed satisfactorily from noisy data.To improve the spatial resolution of the inverse solution we solve a Bayesian optimal experimental design problem to find the best station configuration by maximizing the expected information gain, defined as the Kullback-Leibler divergence between posterior and prior distributions. We show how and why the optimal network improves the material property inference more than evenly spaced stations. Based on our previous work on inverting for fault slip without Green’s function computations, we combine the two inversion schemes to jointly infer both model parameters, the coseismic slip, and material properties distribution. Lastly, we test this numerical forward/inverse framework with an application, the 2011 Tohoku-oki M9 earthquake. Both continuous land-based and six offshore acoustic GNSS stations located around the earthquake epicenter are inverted to jointly estimate the shear modulus and the fault slip during the megathrust event.Our results demonstrates the potential of our computational framework and the general approach for inferring constitutive laws to evaluate sensitivity to parameters, and define strategies to improve our understanding of relevant parameters for earthquake dynamics.  

Bayesian inference provides a systematic framework for integration of data with mathematical models to quantify the uncertainty in the solution of the inverse problem. However, the solution of Bayesian inverse problems governed by complex forward models described by partial differential equations (PDEs) remains prohibitive with black-box Markov chain Monte Carlo (MCMC) methods. We present hIPPYlib-MUQ, an extensible and scalable software framework that contains implementations of state-of-the art algorithms aimed to overcome the challenges of high-dimensional, PDE-constrained Bayesian inverse problems. These algorithms accelerate MCMC sampling by exploiting the geometry and intrinsic low-dimensionality of parameter space via derivative information and low rank approximation. The software integrates two complementary open-source software packages, hIPPYlib and MUQ. hIPPYlib solves PDE-constrained inverse problems using automatically-generated adjoint-based derivatives, but it lacks full Bayesian capabilities. MUQ provides a spectrum of powerful Bayesian inversion models and algorithms, but expects forward models to come equipped with gradients and Hessians to permit large-scale solution. By combining these two complementary libraries, we created a robust, scalable, and efficient software framework that realizes the benefits of each and allows us to tackle complex large-scale Bayesian inverse problems across a broad spectrum of scientific and engineering disciplines. To illustrate the capabilities of hIPPYlib-MUQ, we present a comparison of a number of MCMC methods available in the integrated software on several high-dimensional Bayesian inverse problems. These include problems characterized by both linear and nonlinear PDEs, various noise models, and different parameter dimensions. The results demonstrate that large (∼ 50×) speedups over conventional black box and gradient-based MCMC algorithms can be obtained by exploiting Hessian information (from the log-posterior), underscoring the power of the integrated hIPPYlib-MUQ framework.