Georg Stadler

Hari Sundar, Georg Stadler, George Biros

We present a comparison of different multigrid approaches for the solution of systems arising from high-order continuous finite element discretizations of elliptic partial differential equations on complex geometries. We consider the pointwise Jacobi, the Chebyshev-accelerated Jacobi and the symmetric successive over-relaxation (SSOR) smoothers, as well as elementwise block Jacobi smoothing. Three approaches for the multigrid hierarchy are compared: 1) high-order $h$-multigrid, which uses high-order interpolation and restriction between geometrically coarsened meshes; 2) $p$-multigrid, in which the polynomial order is reduced while the mesh remains unchanged, and the interpolation and restriction incorporate the different-order basis functions; and 3), a first-order approximation multigrid preconditioner constructed using the nodes of the high-order discretization. This latter approach is often combined with algebraic multigrid for the low-order operator and is attractive for high-order discretizations on unstructured meshes, where geometric coarsening is difficult. Based on a simple performance model, we compare the computational cost of the different approaches. Using scalar test problems in two and three dimensions with constant and varying coefficients, we compare the performance of the different multigrid approaches for polynomial orders up to 16. Overall, both $h$- and $p$-multigrid work well; the first-order approximation is less efficient. For constant coefficients, all smoothers work well. For variable coefficients, Chebyshev and SSOR smoothing outperforms Jacobi smoothing. While all of the tested methods converge in a mesh-independent number of iterations, none of them behaves completely independent of the polynomial order. When multigrid is used as a preconditioner in a Krylov method, the iteration number decreases significantly compared to using multigrid as a solver.

Benjamin Crestel, Georg Stadler, Omar Ghattas

Joint inversion refers to the simultaneous inference of multiple parameter fields from observations of systems governed by single or multiple forward models. In many cases these parameter fields reflect different attributes of a single medium and are thus spatially correlated or structurally similar. By imposing prior information on their spatial correlations via a joint regularization term, we seek to improve the reconstruction of the parameter fields relative to inversion for each field independently. One of the main challenges is to devise a joint regularization functional that conveys the spatial correlations or structural similarity between the fields while at the same time permitting scalable and efficient solvers for the joint inverse problem. We describe several joint regularizations that are motivated by these goals: a cross-gradient and a normalized cross-gradient structural similarity term, the vectorial total variation, and a joint regularization based on the nuclear norm of the gradients. Based on numerical results from three classes of inverse problems with piecewise-homogeneous parameter fields, we conclude that the vectorial total variation functional is preferable to the other methods considered. Besides resulting in good reconstructions in all experiments, it allows for scalable, efficient solvers for joint inverse problems governed by PDE forward models.

Chen Li, Georg Stadler

We study sparse solutions of optimal control problems governed by PDEs with uncertain coefficients. We propose two formulations, one where the solution is a deterministic control optimizing the mean objective, and a formulation aiming at stochastic controls that share the same sparsity structure. In both formulations, regions where the controls do not vanish can be interpreted as optimal locations for placing control devices. In this paper, we focus on linear PDEs with linearly entering uncertain parameters. Under these assumptions, the deterministic formulation reduces to a problem with known structure, and thus we mainly focus on the stochastic control formulation. Here, shared sparsity is achieved by incorporating the $L^1$-norm of the mean of the pointwise squared controls in the objective. We reformulate the problem using a norm reweighting function that is defined over physical space only and thus helps to avoid approximation of the random space using samples or quadrature. We show that a fixed point algorithm applied to the norm reweighting formulation leads to a variant of the well-studied iterative reweighted least squares (IRLS) algorithm, and we propose a novel preconditioned Newton-conjugate gradient method to speed up the IRLS algorithm. We combine our algorithms with low-rank operator approximations, for which we provide estimates of the truncation error. We carefully examine the computational complexity of the resulting algorithms. The sparsity structure of the optimal controls and the performance of the solution algorithms are studied numerically using control problems governed by the Laplace and Helmholtz equations. In these experiments the Newton variant clearly outperforms the IRLS method.

Gonzalo G. de Diego, Georg Stadler

Complex physical systems which exhibit fluid-like behavior are often modeled as non-Newtonian fluids. A crucial element of a non-Newtonian model is the rheology, which relates inner stresses with strain-rates. We propose a framework for inferring rheological models from data that represents the fluid's effective viscosity with a neural network. By writing the rheological law in terms of tensor invariants and tailoring the network's properties, the inferred model satisfies key physical and mathematical properties, such as isotropic frame-indifference and existence of a convex potential of dissipation. Within this framework, we propose two approaches to learning a fluid's rheology: 1) a standard regression that fits the rheological model to stress data and 2) a PDE-constrained optimization method that infers rheological models from velocity data. For the latter approach, we combine finite element and machine learning libraries. We demonstrate the accuracy and robustness of our method on land and sea ice rheologies which also depend on external parameters. For land ice, we infer the temperature-dependent Glen's law and, for sea ice, the concentration-dependent shear component of the viscous-plastic model. For these two models, we explore the effects of large data errors. Finally, we infer an unknown concentration-dependent model that reproduces Lagrangian ice floe simulation data. Our method discovers a rheology that generalizes well outside of the training dataset and exhibits both shear-thickening and thinning behaviors depending on the concentrations.

Harun Ur Rashid, Mingxin Li, Aleksandra Pachalieva et al.

Accurate characterization of subsurface heterogeneity is challenging but essential for applications such as reservoir pressure management, geothermal energy extraction and CO$_2$, H$_2$, and wastewater injection operations. This challenge becomes especially acute in extreme pressure events, which are rarely observed but can strongly affect operational risk. Traditional history matching and inversion techniques rely on expensive full-physics simulations, making it infeasible to handle uncertainty and extreme events at scale. Purely data-driven models often struggle to maintain physics consistency when dealing with sparse observations, complex geology, and extreme events. To overcome these limitations, we introduce a physics-informed machine learning method that embeds a differentiable subsurface flow simulator directly into neural network training. The network infers heterogeneous permeability fields from limited pressure observations, while training minimizes both permeability and pressure losses through the simulator, enforcing physical consistency. Because the simulator is used only during training, inference remains fast once the model is learned. In an initial test, the proposed method reduces the pressure inference error by half compared with a purely data-driven approach. We then extend the test over eight distinct data scenarios, and in every case, our method produces significantly lower pressure inference errors than the purely data-driven model. We also evaluate our method on extreme events, which represent high-consequence data in the tail of the sample distribution. Similar to the bulk distribution, the physics-informed model maintains higher pressure inference accuracy in the extreme event regimes. Overall, the proposed method enables rapid, physics-consistent subsurface inversion for real-time reservoir characterization and risk-aware decision-making.

Kewei Wang, Georg Stadler

The spherical-radial decomposition (SRD) is an efficient method for estimating probabilistic functions and their gradients defined over finite-dimensional elliptical distributions. In this work, we generalize the SRD to infinite stochastic dimensions by combining subspace SRD with standard Monte Carlo methods. The resulting method, which we call hybrid infinite-dimensional SRD (hiSRD) provides an unbiased, low-variance estimator for convex sets arising, for instance, in chance-constrained optimization. We provide a theoretical analysis of the variance of finite-dimensional SRD as the dimension increases, and show that the proposed hybrid method eliminates truncation-induced bias, reduces variance, and allows the computation of derivatives of probabilistic functions. We present comprehensive numerical studies for a risk-neutral stochastic PDE optimal control problem with joint chance state constraints, and for optimizing kernel parameters in Gaussian process regression under the constraint that the posterior process satisfies joint chance constraints.

Chen Li, Matthew Dunlop, Georg Stadler

We consider Bayesian inverse problems wherein the unknown state is assumed to be a function with discontinuous structure a priori. A class of prior distributions based on the output of neural networks with heavy-tailed weights is introduced, motivated by existing results concerning the infinite-width limit of such networks. We show theoretically that samples from such priors have desirable discontinuous-like properties even when the network width is finite, making them appropriate for edge-preserving inversion. Numerically we consider deconvolution problems defined on one- and two-dimensional spatial domains to illustrate the effectiveness of these priors; MAP estimation, dimension-robust MCMC sampling and ensemble-based approximations are utilized to probe the posterior distribution. The accuracy of point estimates is shown to exceed those obtained from non-heavy tailed priors, and uncertainty estimates are shown to provide more useful qualitative information.

Tobin Isaac, Noemi Petra, Georg Stadler et al.

The majority of research on efficient and scalable algorithms in computational science and engineering has focused on the forward problem: given parameter inputs, solve the governing equations to determine output quantities of interest. In contrast, here we consider the broader question: given a (large-scale) model containing uncertain parameters, (possibly) noisy observational data, and a prediction quantity of interest, how do we construct efficient and scalable algorithms to (1) infer the model parameters from the data (the deterministic inverse problem), (2) quantify the uncertainty in the inferred parameters (the Bayesian inference problem), and (3) propagate the resulting uncertain parameters through the model to issue predictions with quantified uncertainties (the forward uncertainty propagation problem)? We present efficient and scalable algorithms for this end-to-end, data-to-prediction process under the Gaussian approximation and in the context of modeling the flow of the Antarctic ice sheet and its effect on sea level. The ice is modeled as a viscous, incompressible, creeping, shear-thinning fluid. The observational data come from InSAR satellite measurements of surface ice flow velocity, and the uncertain parameter field to be inferred is the basal sliding parameter. The prediction quantity of interest is the present-day ice mass flux from the Antarctic continent to the ocean. We show that the work required for executing this data-to-prediction process is independent of the state dimension, parameter dimension, data dimension, and number of processor cores. The key to achieving this dimension independence is to exploit the fact that the observational data typically provide only sparse information on model parameters. This property can be exploited to construct a low rank approximation of the linearized parameter-to-observable map.

Tobin Isaac, Georg Stadler, Omar Ghattas

Motivated by the need for efficient and accurate simulation of the dynamics of the polar ice sheets, we design high-order finite element discretizations and scalable solvers for the solution of nonlinear incompressible Stokes equations. We focus on power-law, shear thinning rheologies used in modeling ice dynamics and other geophysical flows. We use nonconforming hexahedral meshes and the conforming inf-sup stable finite element velocity-pressure pairings $\mathbb{Q}_k\times \mathbb{Q}^\text{disc}_{k-2}$ or $\mathbb{Q}_k \times \mathbb{P}^\text{disc}_{k-1}$. To solve the nonlinear equations, we propose a Newton-Krylov method with a block upper triangular preconditioner for the linearized Stokes systems. The diagonal blocks of this preconditioner are sparse approximations of the (1,1)-block and of its Schur complement. The (1,1)-block is approximated using linear finite elements based on the nodes of the high-order discretization, and the application of its inverse is approximated using algebraic multigrid with an incomplete factorization smoother. This preconditioner is designed to be efficient on anisotropic meshes, which are necessary to match the high aspect ratio domains typical for ice sheets. We develop and make available extensions to two libraries---a hybrid meshing scheme for the p4est parallel AMR library, and a modified smoothed aggregation scheme for PETSc---to improve their support for solving PDEs in high aspect ratio domains. In a numerical study, we find that our solver yields fast convergence that is independent of the element aspect ratio, the occurrence of nonconforming interfaces, and of mesh refinement, and that depends only weakly on the polynomial finite element order. We simulate the ice flow in a realistic description of the Antarctic ice sheet derived from field data, and study the parallel scalability of our solver for problems with up to 383M unknowns.