Max Gunzburger

Marta D'Elia, Max Gunzburger

We analyze a nonlocal diffusion operator having as special cases the fractional Laplacian and fractional differential operators that arise in several applications. In our analysis, a nonlocal vector calculus is exploited to define a weak formulation of the nonlocal problem. We demonstrate that, when sufficient conditions on certain kernel functions hold, the solution of the nonlocal equation converges to the solution of the fractional Laplacian equation on bounded domains as the nonlocal interactions become infinite. We also introduce a continuous Galerkin finite element discretization of the nonlocal weak formulation and we derive a priori error estimates. Through several numerical examples we illustrate the theoretical results and we show that by solving the nonlocal problem it is possible to obtain accurate approximations of the solutions of fractional differential equations circumventing the problem of treating infinite-volume constraints.

Wenbin Chen, Max Gunzburger, Dong Sun et al.

We propose and study two second-order in time implicit-explicit (IMEX) methods for the coupled Stokes-Darcy system that governs flows in karst aquifers. The first is a combination of a second-order backward differentiation formula and the second-order Gear's extrapolation approach. The second is a combination of the second-order Adams-Moulton and second-order Adams-Bashforth methods. Both algorithms only require the solution of two decoupled problems at each time step, one Stokes and the other Darcy. Hence, these schemes are very efficient and can be easily implemented using legacy codes. We establish the unconditional and uniform in time stability for both schemes. The uniform in time stability leads to uniform in time control of the error which is highly desirable for modeling physical processes, e.g., contaminant sequestration and release, that occur over very long time scales. Error estimates for fully-discretized schemes using finite element spatial discretizations are derived. Numerical examples are provided that illustrate the accuracy, efficiency, and long-time stability of the two schemes.

Max Gunzburger, Nan Jiang, Michael Schneier

The definition of partial differential equation (PDE) models usually involves a set of parameters whose values may vary over a wide range. The solution of even a single set of parameter values may be quite expensive. In many cases, e.g., optimization, control, uncertainty quantification, and other settings, solutions are needed for many sets of parameter values. We consider the case of the time-dependent Navier-Stokes equations for which a recently developed ensemble-based method allows for the efficient determination of the multiple solutions corresponding to many parameter sets. The method uses the average of the multiple solutions at any time step to define a linear set of equations that determines the solutions at the next time step. To significantly further reduce the costs of determining multiple solutions of the Navier-Stokes equations, we incorporate a proper orthogonal decomposition (POD) reduced-order model into the ensemble-based method. The stability and convergence results for the ensemble-based method are extended to the ensemble-POD approach. Numerical experiments are provided that illustrate the accuracy and efficiency of computations determined using the new approach.

Max Gunzburger, Buyang Li, Jilu Wang

The stochastic time-fractional equation $\partial_t ψ-Δ\partial_t^{1-α} ψ= f + \dot W$ with space-time white noise $\dot W$ is discretized in time by a backward-Euler convolution quadrature for which the sharp-order error estimate \[ {\mathbb E}\|ψ(\cdot,t_n)-ψ_n\|_{L^2(\mathcal{O})}^2=O(τ^{1-αd/2}) \] is established for $α\in(0,2/d)$, where $d$ denotes the spatial dimension, $ψ_n$ the approximate solution at the $n^{\rm th}$ time step, and $\mathbb{E}$ the expectation operator. In particular, the result indicates optimal convergence rates of numerical solutions for both stochastic subdiffusion and diffusion-wave problems in one spatial dimension. Numerical examples are presented to illustrate the theoretical analysis.

Martin Hess, Alessandro Alla, Annalisa Quaini et al.

Reduced-order modeling (ROM) commonly refers to the construction, based on a few solutions (referred to as snapshots) of an expensive discretized partial differential equation (PDE), and the subsequent application of low-dimensional discretizations of partial differential equations (PDEs) that can be used to more efficiently treat problems in control and optimization, uncertainty quantification, and other settings that require multiple approximate PDE solutions. In this work, a ROM is developed and tested for the treatment of nonlinear PDEs whose solutions bifurcate as input parameter values change. In such cases, the parameter domain can be subdivided into subregions, each of which corresponds to a different branch of solutions. Popular ROM approaches such as proper orthogonal decomposition (POD), results in a global low-dimensional basis that does no respect not take advantage of the often large differences in the PDE solutions corresponding to different subregions. Instead, in the new method, the k-means algorithm is used to cluster snapshots so that within cluster snapshots are similar to each other and are dissimilar to those in other clusters. This is followed by the construction of local POD bases, one for each cluster. The method also can detect which cluster a new parameter point belongs to, after which the local basis corresponding to that cluster is used to determine a ROM approximation. Numerical experiments show the effectiveness of the method both for problems for which bifurcation cause continuous and discontinuous changes in the solution of the PDE.

Max Gunzburger, Buyang Li, Jilu Wang

Numerical approximation of a stochastic partial integro-differential equation driven by a space- time white noise is studied by truncating a series representation of the noise, with finite element method for spatial discretization and convolution quadrature for time discretization. Sharp-order convergence of the numerical solutions is proved up to a logarithmic factor. Numerical examples are provided to support the theoretical analysis.

Max Gunzburger, Traian Iliescu, Michael Schneier

Partial differential equations (PDEs) are often dependent on input quantities which are inherently uncertain. To quantify this uncertainty, these PDEs must be solved over a large ensemble of parameters. Even for a single realization this can a computationally intensive process. In the case of flows governed by the Navier-Stokes equations, an efficient method has been devised for computing an ensemble of solutions. To further reduce the computational cost of this method, an ensemble proper orthogonal decomposition (POD) method was recently proposed. The main contribution of this work is the introduction of POD spatial filtering for ensemble-POD methods. The POD spatial filter makes possible the construction of the Leray ensemble-POD model, which is a regularized reduced order model for the numerical simulation of convection-dominated flows. The Leray ensemble-POD model employs the POD spatial filter to smooth (regularize) the convection term in the Navier-Stokes equations and greatly diminishes the numerical inaccuracies produced by the ensemble-POD method in the numerical simulation of convection-dominated flows. Specifically, for the numerical simulation of a convection-dominated two-dimensional flow between two offset cylinders, we show that the Leray ensemble-POD method yields accurate results, whereas the ensemble-POD is highly inaccurate. The second contribution of this work is a new numerical discretization of the variable viscosity ensemble algorithm in which the average viscosity is replaced with the maximum viscosity. It is shown that this new numerical discretization is significantly more stable than those in current use. Furthermore, error estimates for the novel Leray ensemble-POD algorithm with this new numerical discretization are also proven.

James Cheung, Mauro Perego, Pavel Bochev et al.

Meshing of geometric domains having curved boundaries by affine simplices produces a polytopial approximation of those domains. The resulting error in the representation of the domain limits the accuracy of finite element methods based on such meshes. On the other hand, the simplicity of affine meshes makes them a desirable modeling tool in many applications. In this paper, we develop and analyze higher-order accurate finite element methods that remain stable and optimally accurate on polytopial approximations of domains with smooth boundaries. This is achieved by constraining a judiciously chosen extension of the finite element solution on the polytopial domain to weakly match the prescribed boundary condition on the true geometric boundary. We provide numerical examples that highlight key properties of the new method and that illustrate the optimal $H^1$ and $L^2$-norm convergence rates.

Max Gunzburger, Nan Jiang, Michael Schneier

Partial differential equations (PDE) often involve parameters, such as viscosity or density. An analysis of the PDE may involve considering a large range of parameter values, as occurs in uncertainty quantification, control and optimization, inference, and several statistical techniques. The solution for even a single case may be quite expensive; whereas parallel computing may be applied, this reduces the total elapsed time but not the total computational effort. In the case of flows governed by the Navier-Stokes equations, a method has been devised for computing an ensemble of solutions. Recently, a reduced-order model derived from a proper orthogonal decomposition (POD) approach was incorporated into a first-order accurate in time version of the ensemble algorithm. In this work, we expand on that work by incorporating the POD reduced order model into a second-order accurate ensemble algorithm. Stability and convergence results for this method are updated to account for the POD/ROM approach. Numerical experiments illustrate the accuracy and efficiency of the new approach.

Max Gunzburger, Nan Jiang, Feifei Xu

We study a turbulence closure model in which the fractional Laplacian $(-Δ)^α$ of the velocity field represents the turbulence diffusivity. We investigate the energy spectrum of the model by applying Pao's energy transfer theory. For the case $α=1/3$, the corresponding power law of the energy spectrum in the inertial range has a correction exponent on the regular Kolmogorov -5/3 scaling exponent. For this case, this model represents Richardson's particle pair-distance superdiffusion of a fully developed homogeneous turbulent flow as well as Lévy jumps that lead to the superdiffusion. For other values of $α$, the power law of the energy spectrum is consistent with the regular Kolmogorov -5/3 scaling exponent. We also propose and study a modular time-stepping algorithm in semi-discretized form. The algorithm is minimally intrusive to a given legacy code for solving Navier-Stokes equations by decoupling the local part and nonlocal part of the equations for the unknowns. We prove the algorithm is unconditionally stable and unconditionally, first-order convergent. We also derive error estimates for full discretizations of the model which, in addition to the time stepping algorithm, involves a finite element spatial discretization and a domain truncation approximation to the range of the fractional Laplacian.

Huanhuan Yang, Max Gunzburger, Lili Ju

Centroidal Voronoi tessellation (CVT)-based mesh generation is a very effective technique for creating high-quality Voronoi meshes and their dual Delaunay triangulations that often play a crucial role in applications, including ocean and atmospheric simulations using finite volume schemes. In the next generation climate models, the spacing scales change dramatically across the whole sphere and require ultra-high resolution and smooth transitions from coarse to fine grid regions. Thus fast and robust spherical CVT (SCVT) meshing algorithms become highly desirable. In this paper, we first propose a Lloyd-preconditioned limited-memory BFGS method for constructing SCVTs that is also applicable to the construction of CVTs of general domains. This method is then parallelized based on overlapping domain decomposition, enabling excellent scalability on distributed systems. Results of several computational experiments show that the new method could incur computational time costs one order of magnitude smaller compared with some existing methods for generating large-scale highly variable-resolution meshes, while also providing significantly improvements in mesh quality.

Max Gunzburger, Aretha L Teckentrup

This paper is concerned with Lagrange interpolation by total degree polynomials in moderate dimensions. In particular, we are interested in characterising the optimal choice of points for the interpolation problem, where we define the optimal interpolation points as those which minimise the Lebesgue constant. We give a novel algorithm for numerically computing the location of the optimal points, which is independent of the shape of the domain and does not require computations with Vandermonde matrices. We perform a numerical study of the growth of the minimal Lebesgue constant with respect to the degree of the polynomials and the dimension, and report the lowest values known as yet of the Lebesgue constant in the unit cube and the unit ball in up to 10 dimensions.

Olena Burkovska, Max Gunzburger

We consider linear and obstacle problems driven by a nonlocal integral operator, for which nonlocal interactions are restricted to a ball of finite radius. These type of operators are used to model anomalous diffusion and, for a special choice of the integral kernels, reduce to the fractional Laplace operator on a bounded domain. By means of a nonlocal vector calculus we recast the problems in a weak form, leading to corresponding nonlocal variational equality and inequality problems. We prove optimal regularity results for both problems, including a higher regularity of the solution and the Lagrange multiplier. Based on the regularity results, we analyze the convergence of finite element approximations for a linear problem and illustrate the theoretical findings by numerical results.

Qingshan Chen, Max Gunzburger

A general framework for goal-oriented a posteriori error estimation for finite volume methods is presented. The framework does not rely on recasting finite volume methods as special cases of finite element methods, but instead directly determines error estimators from the discretized finite volume equations. Thus, the framework can be ap- plied to arbitrary finite volume methods. It also provides the proper functional settings to address well-posedness issues for the primal and adjoint problems. Numerical results are presented to illustrate the validity and effectiveness of the a posteriori error estimates and their applicability to adaptive mesh refinement.

Max Gunzburger, Michael Schneier, Clayton Webster et al.

It is shown that the computational efficiency of the discrete least-squares (DLS) approximation of solutions of stochastic elliptic PDEs is improved by incorporating a reduced-basis method into the DLS framework. The goal is to recover the entire solution map from the parameter space to the finite element space. To this end, first, a reduced-basis solution using a weak greedy algorithm is constructed, then a DLS approximation is determined by evaluating the reduced-basis approximation instead of the full finite element approximation. The main advantage of the new approach is that one only need apply the DLS operator to the coefficients of the reduced-basis expansion, resulting in huge savings in both the storage of the DLS coefficients and the online cost of evaluating the DLS approximation. In addition, the recently developed quasi-optimal polynomial space is also adopted in the new approach, resulting in superior convergence rates for a wider class of problems than previous analyzed. Numerical experiments are provided that illustrate the theoretical results.

Pavel Bochev, James Cheung, Max Gunzburger et al.

In this paper, we present a new numerical method for determining the numerical solution of interface problems to optimal accuracy with respect to the polynomial order of the Lagrangian finite element space on polytopial meshes. We introduce the notion of a virtual interface, and on this virtual interface we enforce that "extended" interface conditions are satisfied in the sense of a Dirichlet--Neumann coupling. The virtual interface framework serves to bypass geometric variational crimes incurred by the classical finite element method. Further, this approach does not require that the geometric interfaces are spatially matching. Our analysis indicates that this approach is well--posed and optimally convergent in $H^1$. Numerical experiments indicate that optimal $H^1$ and $L^2$ convergence is achieved.

Pavel Bochev, James Cheung, Max Gunzburger et al.

We present a new formulation based on the classical Dirichlet-Neumann formulation for interface coupling problems in linearized elasticity. By using Taylor series expansions, we derive a new set of interface conditions that allow our formulation to pass the linear consistency test. In addition, we propose an iterative method to determine the solution of our formulation. We demonstrate in our numerical results that we may achieve the desired piecewise linear finite element error bounds for both nonoverlapping domain decomposition problems as well as for interface coupling problems where the Lamé parameters of the structures differ.

Yuwei Geng, Olena Burkovska, Lili Ju et al.

We propose an efficient end-to-end deep learning method for solving nonlocal Allen-Cahn (AC) and Cahn-Hilliard (CH) phase-field models. One motivation for this effort emanates from the fact that discretized partial differential equation-based AC or CH phase-field models result in diffuse interfaces between phases, with the only recourse for remediation is to severely refine the spatial grids in the vicinity of the true moving sharp interface whose width is determined by a grid-independent parameter that is substantially larger than the local grid size. In this work, we introduce non-mass conserving nonlocal AC or CH phase-field models with regular, logarithmic, or obstacle double-well potentials. Because of non-locality, some of these models feature totally sharp interfaces separating phases. The discretization of such models can lead to a transition between phases whose width is only a single grid cell wide. Another motivation is to use deep learning approaches to ameliorate the otherwise high cost of solving discretized nonlocal phase-field models. To this end, loss functions of the customized neural networks are defined using the residual of the fully discrete approximations of the AC or CH models, which results from applying a Fourier collocation method and a temporal semi-implicit approximation. To address the long-range interactions in the models, we tailor the architecture of the neural network by incorporating a nonlocal kernel as an input channel to the neural network model. We then provide the results of extensive computational experiments to illustrate the accuracy, structure-preserving properties, predictive capabilities, and cost reductions of the proposed method.

Anthony Gruber, Max Gunzburger, Lili Ju et al.

The popularity of deep convolutional autoencoders (CAEs) has engendered new and effective reduced-order models (ROMs) for the simulation of large-scale dynamical systems. Despite this, it is still unknown whether deep CAEs provide superior performance over established linear techniques or other network-based methods in all modeling scenarios. To elucidate this, the effect of autoencoder architecture on its associated ROM is studied through the comparison of deep CAEs against two alternatives: a simple fully connected autoencoder, and a novel graph convolutional autoencoder. Through benchmark experiments, it is shown that the superior autoencoder architecture for a given ROM application is highly dependent on the size of the latent space and the structure of the snapshot data, with the proposed architecture demonstrating benefits on data with irregular connectivity when the latent space is sufficiently large.

Anthony Gruber, Max Gunzburger, Lili Ju et al.

A dimension reduction method based on the "Nonlinear Level set Learning" (NLL) approach is presented for the pointwise prediction of functions which have been sparsely sampled. Leveraging geometric information provided by the Implicit Function Theorem, the proposed algorithm effectively reduces the input dimension to the theoretical lower bound with minor accuracy loss, providing a one-dimensional representation of the function which can be used for regression and sensitivity analysis. Experiments and applications are presented which compare this modified NLL with the original NLL and the Active Subspaces (AS) method. While accommodating sparse input data, the proposed algorithm is shown to train quickly and provide a much more accurate and informative reduction than either AS or the original NLL on two example functions with high-dimensional domains, as well as two state-dependent quantities depending on the solutions to parametric differential equations.

Max Gunzburger, Nan Jiang, Zhu Wang

We consider settings for which one needs to perform multiple flow simulations based on the Navier-Stokes equations, each having different values for the physical parameters and/or different initial condition data, boundary conditions data, and/or forcing functions. For such settings, we propose a second-order time accurate ensemble-based method that to simulate the whole set of solutions, requires, at each time step, the solution of only a single linear system with multiple right-hand-side vectors. Rigorous analyses are given proving the conditional stability and error estimates for the proposed algorithm. Numerical experiments are provided that illustrate the analyses.

Max Gunzburger, Nan Jiang, Zhu Wang

Many applications of computational fluid dynamics require multiple simulations of a flow under different input conditions. In this paper, a numerical algorithm is developed to efficiently determine a set of such simulations in which the individually independent members of the set are subject to different viscosity coefficients, initial conditions, and/or body forces. The proposed scheme applied to the flow ensemble leads to need to solve a single linear system with multiple right-hand sides, and thus is computationally more efficient than solving for all the simulations separately. We show that the scheme is nonlinearly and long-term stable under certain conditions on the time-step size and a parameter deviation ratio. Rigorous numerical error estimate shows the scheme is of first-order accuracy in time and optimally accurate in space. Several numerical experiments are presented to illustrate the theoretical results.

Guannan Zhang, Weidong Zhao, Clayton Webster et al.

We propose a novel numerical approach for nonlocal diffusion equations [8] with integrable kernels, based on the relationship between the backward Kolmogorov equation and backward stochastic differential equations (BSDEs) driven by Lèvy processes with jumps. The nonlocal diffusion problem under consideration is converted to a BSDE,for which numerical schemes are developed and applied directly. As a stochastic approach, the proposed method does not require the solution of linear systems, which allows for embarrassingly parallel implementations and also enables adaptive approximation techniques to be incorporated in a straightforward fashion. Moreover, our method is more accurate than classic stochastic approaches due to the use of high-order temporal and spatial discretization schemes. In addition, our approach can handle a broad class of problems with general nonlinear forcing terms as long as they are globally Lipchitz continuous. Rigorous error analysis of the new method is provided as several numerical examples that illustrate the effectiveness and efficiency of the proposed approach.