Research

Emergency managers need to know where floodwater is, how deep it is, and how it will change over the coming hours across an entire river basin. During a flood, however, real-time measurements come from only a handful of stream gauges, and high-resolution hydrodynamic models are too costly to rerun each time new data arrive or to run as large ensembles.

We are building the components of flood digital twins that answer in real time. GeoSWE is a GPU finite-volume solver for the full 2D shallow-water equations, from a 3 m coastal county to the conterminous United States on one multi-GPU node. On the Frontier supercomputer it simulated Hurricanes Helene and Milton across Pinellas County on a 1 m grid with every building, which you can explore in 3D in Two Storms, One Metre. CLDNet, a conditional latent dynamics network surrogate of the shallow water equations, produces a 96-hour basin-wide forecast in about 29 seconds on one GPU, about 115 times faster than the simulator it is trained on. C-STRIDE, an observation-driven AI digital twin, turns short records from a few stream gauges, together with terrain and rainfall, into basin-wide maps of water depth, and extends these predictions up to a day ahead: in the Des Plaines River basin near Chicago, six gauges inform predictions over 4.2 million 30 m grid cells.

A digital twin also has a physical side: sensors on the real system, and a way for its predictions to reach the people who act on them. In Pinellas County, a low-lying peninsula on Florida's Gulf Coast, our collaborators at the University of South Florida are building this side as CRIS-HAZARD, a community resilience information system led by Prof. Barnali Dixon and funded by the National Science Foundation (award 2325631). They have installed 32 cameras on the coast, along creeks and canals, and inland, aimed at staffs marked with reflective bands so that the water level can be read from the pictures by day and by night, and residents add their own flood reports. Our part is computational: to combine these sparse observations, along with tide gauges, stream gauges, and weather stations, with real-time flood prediction, and to quantify how far the result can be trusted before it goes back to residents and emergency managers.

Bayesian inversion, optimal experimental design, and optimization under uncertainty are outer-loop tasks that solve a PDE again and again. A neural operator that replaces the solver has to be trained from a limited number of expensive simulations, and because these tasks rely on the derivatives of the solution map, it has to approximate the derivatives as well as the values.

Derivative-informed projected neural networks (DIPNet) project the high-dimensional input and output functions onto derivative-informed low-dimensional subspaces and train a parsimonious network whose size is independent of the discretization. Derivative-informed neural operators (DINO) add the derivatives to the training loss at a cost comparable to training without them, and DE-DeepONet carries the idea to DeepONet and Fourier neural operators. Both are more accurate than operators trained on function values alone, for the solution and for its derivatives, especially when training data are limited.

We also ask when such surrogates can be trusted. Our analysis of reduced basis neural operators bounds the approximation errors of the operator and of its derivatives, from the dimension reduction, the network, and the sampling. The variationally correct reduced basis neural operator is trained on a first-order system least-squares residual that is provably equivalent to the solution error, so the loss itself is a computable a posteriori error estimator. For transformer operator networks, we establish a power scaling law between the prediction error and the number of training samples, with an exponent that reflects the intrinsic dimension of the input functions.

The Tucker tensor train Taylor series (T4S) is a local alternative to global operator learning: a high-order Taylor surrogate fitted from derivative probes at a single expansion point, with each derivative tensor represented as a Tucker tensor train. Earlier, we developed a scalable randomized algorithm for converting tensors into tensor train format based on actions of the tensor as a multilinear map, which allows efficient compression of large high-order derivative tensors for nonlinear mappings that are implicitly defined through the solution of a system of equations.

Accurate prediction of complex physical systems relies on data assimilation to correct the errors of model simulations with observations. Traditional ensemble Kalman filters and the more recent ensemble score filters face significant challenges when the states are high-dimensional, the dynamics are nonlinear, and the observations are sparse, which is the usual situation in practice.

We move the assimilation into a learned low-dimensional latent space. Latent-EnSF encodes the full states and the sparse observations consistently with a coupled variational autoencoder and runs an ensemble score filter there. LD-EnSF also evolves the dynamics in the latent space, with latent dynamics networks and an LSTM encoder of the observation history, which eliminates the full-space simulations and gives speedups of orders of magnitude. LEVDA is the variational counterpart: a four-dimensional ensemble-variational smoother in the latent space of a differentiable neural dynamics surrogate, which assimilates states and unknown parameters jointly, without adjoint code, from observations at arbitrary locations and times. LD-EnFF replaces the score filter with an ensemble update based on flow matching and evaluates the observation likelihood in latent space. The benchmarks include Kolmogorov flow, tsunami propagation, and atmospheric dynamics, with noisy observations that are sparse in both space and time.

A related question is how much of a field can be recovered from a few point sensors alone. STRIDE maps a short window of sensor measurements to a latent state with a temporal encoder and reconstructs the field at arbitrary locations with an implicit neural decoder. C-STRIDE in our flood digital twin extends it with terrain and rainfall.

By exploiting the intrinsic low-dimensional geometric structure of the posterior distributions in high-dimensional Bayesian inverse problems, we developed projected variational inference methods (pSVGD and pSVN) to tackle the curse of dimensionality faced in such problems. The methods are demonstrated to be fast, accurate, and scalable with respect to the number of parameters, samples, data points, and parallel processor cores for scientific, engineering, and machine learning problems.

We also work with other transports of probability distributions. Projected Wasserstein gradient descent (pWGD) projects the parameters into a low-dimensional subspace to overcome the curse of dimensionality of kernel density estimation, and a convex semidefinite program gives the optimal approximation of the Wasserstein gradient direction in a class of two-layer neural networks. PDPO computes action-minimizing paths between probability densities by representing the path as the pushforward of a reference density through a parametric map; three to five control points of a spline in parameter space suffice for multimodal and high-dimensional problems.

For Bayesian inverse problems governed by PDEs at scale, hIPPyMFEM provides the full workflow on GPUs, from the MAP point and a low-rank Laplace approximation of the posterior to MCMC. In a geothermal application it solves the forward problem with 1.08 billion state unknowns on 30 GPUs, and infers the conductivity of the rock on a 128³ mesh from temperatures logged in 60 boreholes.

Experiments and sensors are expensive, so we would like to choose the most informative ones before any data are collected. Bayesian optimal experimental design does this by maximizing the expected information gain. This involves a double integration over the unknown parameters and the possible data, a Bayesian inverse problem for every candidate design, and, when sensors are selected from a set of candidates, a combinatorial optimization.

We make this tractable for models governed by PDEs with infinite-dimensional parameters. The low-rank structure of the parameter-to-observable map and an offline–online decomposition remove the PDE solves from the design optimization. Derivative-informed neural operators then supply the MAP points and the eigenvalues that the optimality criteria need: for selecting sensors in a three-dimensional nonlinear convection–diffusion–reaction problem with tens of thousands of parameters, the speedup is more than 1000×, counting both offline construction and online evaluation.

When experiments are run in sequence, each design should adapt to the data already seen. A derivative-informed latent attention neural operator (LANO) chooses the times at which to take MRI images to monitor tumor growth, with an amortized speedup of 180×. A design policy trained by policy-gradient reinforcement learning places sensors online to track a contaminant source, about 100× faster than high-fidelity finite element methods, and discovers a strategy of moving the sensors upstream. For MRI data acquisition with binary sampling masks, a conditional normalizing flow is trained jointly with a probabilistic formulation of the design.

A design that is optimal under nominal conditions can fail once manufacturing errors, uncertain material properties, or a random environment enter. Optimizing a risk measure instead requires a large number of samples at every iteration of the optimization, each one an expensive PDE solve, while the random parameters and the design variables are both high-dimensional.

We develop Taylor approximations and neural operators that reduce this cost by orders of magnitude. A Taylor approximation-based approximate Newton method designs an acoustic metamaterial cloak that hides an obstacle under manufacturing errors and incident waves from multiple angles and frequencies; it is scalable, in terms of the number of PDE solves, with respect to both the uncertain parameter dimension (up to half a million) and the optimization variable dimension (over one million). Gaussian mixture Taylor approximations keep the accuracy when the input variances are large, and estimate the conditional value-at-risk to within 1% with on the order of 10 PDE solves, comparable to Monte Carlo with on the order of 10,000 samples. Multi-input reduced basis derivative-informed neural operators (MR-DINO) solve risk-averse control problems for semilinear elliptic and Navier–Stokes equations with random field inputs 10³ to 10⁷ times faster in execution, and more than 10 times cheaper once the cost of training is included.

On the theory side, we establish performance bounds for approximation-based algorithms: the optimality gap of the computed controls becomes arbitrarily small as the inaccuracies of discretization, sampling, smoothing, and penalization vanish. Scalable algorithms for PDE-constrained optimization under uncertainty are available in our open-source library SOUPy.

Digital twins, Bayesian inversion, and the training of surrogates all rest on forward models that are fast at the scale of the real system: a continent at 30 m resolution, or a PDE with a billion unknowns. We write these codes for GPUs and release them as open source.

GeoSWE solves the full 2D shallow-water equations on NVIDIA and AMD GPUs. Its compressed active-cell mesh packs the cells that matter into flat arrays, so memory and work scale with the flooded landscape rather than its bounding rectangle. A 72-hour Hurricane Helene scenario over the conterminous United States at 30 m, with 8.88 billion active cells, runs on the eight GPUs of one node in 10.8 hours. On a 214-million-cell county benchmark against TRITON, SynxFlow, and SERGHEI, it has the lowest wall time and GPU memory at every GPU count. On OLCF Frontier, 1.024 trillion cells on 128 nodes advance at 162 ms per step, a weak-scaling efficiency of 96.7% against a single GPU.

hIPPyMFEM runs Bayesian inverse problems governed by PDEs on GPUs by differentiable programming: the residual, the Jacobian, the adjoint, and every Hessian block come from one JAX function by automatic differentiation, with no derivative coded by hand. Two Newton-CG steps at 287 million unknowns take 38 seconds on 16 GPUs, 161 times faster than hIPPYlibx on 32 CPU cores, and at 1.09 billion unknowns they take 134 seconds on 24 GPUs. The whole Bayesian workflow, from the MAP point and a low-rank Laplace posterior to samples and pointwise variance, takes 10 minutes on four GPUs for a problem with 36 million unknowns.