Prescribed Error Tolerance Research Articles

Development and Verification of a Higher-Order Computational Fluid Dynamics Solver

Abstract Over the past two decades, higher-order methods have gained much broader use in computational science and engineering as these schemes are often more efficient per degree-of-freedom at achieving a prescribed error tolerance than lower-order methods. During this time, higher-order variants of most discretization schemes, such as finite difference methods, finite volume methods, and finite element methods, have emerged. The finite volume method is arguably the most widely used discretization technique in production-level computational fluid dynamics solvers due to its robustness and conservation properties. However, most finite volume solvers only employ a conventional second-order scheme. To leverage the benefits of higher-order methods, the higher-order finite volume method seems the most natural for those seeking to extend their legacy solvers to higher-order. Nonetheless, ensuring higher-order accuracy is maintained is quite challenging as the implementation requirements for a higher-order scheme are much greater than those of a lower-order scheme. In this work, a methodology for verifying higher-order finite volume codes is presented. The higher-order finite volume method is outlined in detail. Order verification tests are proposed for all major components, including the treatment of curved boundaries and the higher-order solution reconstruction. System-level verification tests are performed using the weak form of the method of manufactured solutions. Several canonical verification cases are also presented for the Euler and laminar Navier–Stokes equations.

Exploiting Kronecker structure in exponential integrators: Fast approximation of the action of [formula omitted]-functions of matrices via quadrature

In this article, we propose an algorithm for approximating the action of φ−functions of matrices against vectors, which is a key operation in exponential time integrators. In particular, we consider matrices with Kronecker sum structure, which arise from problems admitting a tensor product representation. The method is based on quadrature approximations of the integral form of the φ−functions combined with a scaling and modified squaring method. Owing to the Kronecker sum representation, only actions of 1D matrix exponentials are needed at each quadrature node and assembly of the full matrix can be avoided. Additionally, we derive a priori bounds for the quadrature error, which show that, as expected by classical theory, the rate of convergence of our method is supergeometric. Guided by our analysis, we construct a fast and robust method for estimating the optimal scaling factor and number of quadrature nodes that minimizes the total cost for a prescribed error tolerance. We investigate the performance of our algorithm by solving several linear and semilinear time-dependent problems in 2D and 3D. The results show that our method is accurate and orders of magnitude faster than the current state-of-the-art.

Command-Filtered Robust Adaptive NN Control With the Prescribed Performance for the 3-D Trajectory Tracking of Underactuated AUVs.

A novel robust adaptive neural network (NN) control scheme with prescribed performance is developed for the 3-D trajectory tracking of underactuated autonomous underwater vehicles (AUVs) with uncertain dynamics and unknown disturbances using new prescribed performance functions, an additional term, the radial basis function (RBF) NN, and the command-filtered backstepping approach. Different from the traditional prescribed performance functions, the new prescribed performance functions are innovatively proposed such that the time desired for the trajectory tracking errors of AUVs to reach and stay within the prescribed error tolerance band can be preset exactly and flexibly. The additional term with the Nussbaum function is designed to deal with the underactuation problem of AUVs. By means of RBF NN, the uncertain item lumped by the uncertain dynamics of AUVs and unknown disturbances is eventually transformed into a linearly parametric form with only a single unknown parameter. The developed control scheme ensures that all signals in the AUV 3-D trajectory tracking closed-loop control system are bounded. Simulation results with comparisons show the validity and the superiority of our developed control scheme.

Numerical simulation of optical propagation using sinc approximation.

This paper further develops a recently proposed method for computing the diffraction integrals of optics based on sinc series approximation by presenting a numerical implementation, parameter selection criteria based on rigorous error analysis, and example optical propagation simulations demonstrating those criteria. Unlike fast Fourier transform (FFT)-based methods that are based on Fourier series, such as the well-known angular spectrum method (ASM), the sinc method uses a basis that is naturally suited to problems on an infinite domain. As such, it has been shown that the sinc method avoids the problems of artificial periodicity inherent in the ASM. After a brief review of the method, the detailed error analysis we provide confirms its super-algebraic convergence and verifies the claim that the accuracy of the method is independent of wavelength, propagation distance, and observation plane discretization; it depends only on the accuracy of the source field approximation. Based on this analysis, we derive parameter selection criteria for achieving a prescribed error tolerance, which will be valuable to potential users. Numerical simulations of Gaussian beam and optical phased array propagation verify the high-order accuracy and computational efficiency of the proposed algorithms. To facilitate the reproduction of numerical results, we provide a Matlab code that implements our numerical approach for the Fresnel diffraction integral. For comparison, we also present numerical results obtained with the ASM as well as the band-limited angular spectrum method.

Improved Variants of the Hutch++ Algorithm for Trace Estimation

This paper is concerned with two improved variants of the Hutch++ algorithm for estimating the trace of a square matrix, implicitly given through matrix-vector products. Hutch++ combines randomized low-rank approximation in a first phase with stochastic trace estimation in a second phase. In turn, Hutch++ only requires $O\left(\varepsilon^{-1}\right)$ matrix-vector products to approximate the trace within a relative error $\varepsilon$ with high probability. This compares favorably with the $O\left(\varepsilon^{-2}\right)$ matrix-vector products needed when using stochastic trace estimation alone. In Hutch++, the number of matrix-vector products is fixed a priori and distributed in a prescribed fashion among the two phases. In this work, we derive an adaptive variant of Hutch++, which outputs an estimate of the trace that is within some prescribed error tolerance with a controllable failure probability, while splitting the matrix-vector products in a near-optimal way among the two phases. For the special case of symmetric positive semi-definite matrix, we present another variant of Hutch++, called Nystr\om++, which utilizes the so called Nystr\om approximation and requires only one pass over the matrix, as compared to two passes with Hutch++. We extend the analysis of Hutch++ to Nystr\om++. Numerical experiments demonstrate the effectiveness of our two new algorithms.

A Precise Time-Domain Expanding Crack-Tip Element for 3D Planar and Axial Dynamic Cracks

Unlike two-dimensional (2D) crack problems of which the analytical symplectic eigensolutions are available, the analytical symplectic solution of three-dimensional (3D) planar and axial dynamic cracks has not been derived yet. We thus propose a trial displacement field for the 3D planar and axial crack problem using the existing solutions of the 2D crack problems and construct a new crack-tip element on that basis. The recently developed precise integration scheme for 3D planar and axial cracks under dynamic loading (Hu, X. F., Chen, W. H., Zhang, P. and Yao, W. A. [2021] “An explicit crack-tip element for stationary dynamic cracks,” Theor. Appl. Fract. Mech. 112, 102886) is further enhanced. A precise time-domain expanding algorithm is employed for the temporal discretization. All the field variables are expanded using a series expansion scheme, and the dynamic problem is then transformed into a series of quasi-static elastic problems. All the sub-elastic problems are related in a recursive manner. Taking advantage of a smart automatic truncating approach, the number of expanding terms can be determined automatically subjected to a prescribed error tolerance. A few numerical examples are investigated to demonstrate the developed method, numerical solutions indicate that the time dependent stress intensity factors (SIFs) can be calculated in a simple way and the present method has higher efficiency than the recently developed time integration scheme while having comparable accuracy.

Multilevel Monte Carlo method for topology optimization of flexoelectric composites with uncertain material properties

We present an efficient multilevel Monte Carlo (MLMC) method for the topology optimization of flexoelectric structures. A flexoelectric composite consisting of flexoelectric and purely elastic building blocks is investigated. The governing equations are solved by Non-Uniform Rational B-spline (NURBS)-based isogeometric analysis (IGA) exploiting its higher order continuity. Genetic algorithms (GA) based integer-valued optimization is used to obtain the optimal topological design. The uncertainties in the material properties and the volume fraction of the constituents are considered to quantify the uncertainty in the electromechanical coupling effect. Then, a multilevel hierarchy of computational meshes is obtained by a uniform refinement according to a geometric sequence. We estimate the growth rate of the simulation cost, in addition to the rates of decay in the expectation and the variance of the differences between the approximations over the hierarchy. Finally, we determine the minimum number of simulations required on each level to achieve the desired accuracy at different prescribed error tolerances. The results show that the proposed method reduces the computational cost in the numerical experiments without loss of the accuracy. The overall computation saving was in the range 2.0–3.5.

SHASTA

A key requirement for efficient general purpose approximate computing is an amalgamation of flexible hardware design and intelligent application tuning, which together can leverage the appropriate amount of approximation that the applications engender and reap the best efficiency gains from them. To achieve this, we have identified three important features to build better general-purpose cross-layer approximation systems: ① individual per-operation (“spatio-temporally fine-grained”) approximation, ② hardware-cognizant application tuning for approximation, ③ systemwide approximation-synergy. We build an efficient general purpose approximation system called SHASTA: Synergic HW-SW Architecture for Spatio-Temporal Approximation, to achieve these goals. 1 First, in terms of hardware, SHASTA approximates both compute and memory—SHASTA proposes (a) a form of timing approximation called Slack-control Approximation, which controls the computation timing of each approximation operation and (b) a Dynamic Pre-L1 Load Approximation mechanism to approximate loads prior to cache access. These hardware mechanisms are designed to achieve fine-grained spatio-temporally diverse approximation. Next, SHASTA proposes a Hardware-cognizant Approximation Tuning mechanism to tune an application’s approximation to achieve the optimum execution efficiency under the prescribed error tolerance. The tuning mechanism is implemented atop a gradient descent algorithm and, thus, the application’s approximation is tuned along the steepest error vs. execution efficiency gradient. Finally, SHASTA is designed with a full-system perspective, which achieves Synergic benefits across its optimizations, building a closer-to-ideal general purpose approximation system. SHASTA is implemented on top of an OOO core and achieves mean speedups/energy savings of 20%–40% over a non-approximate baseline for greater than 90% accuracy—these benefits are substantial for applications executing on a traditional general purpose processing system. SHASTA can be tuned to specific accuracy constraints and execution metrics and is quantitatively shown to achieve 2–15× higher benefits, in terms of performance and energy, compared to prior work.

Randomized matrix approximation to enhance regularized projection schemes in inverse problems

The authors consider a randomized solution to ill-posed operator equations in Hilbert spaces. In contrast to statistical inverse problems, where randomness appears in the noise, here randomness arises in the low-rank matrix approximation of the forward operator, which results in using a Monte Carlo method to solve the inverse problems. In particular, this approach follows the paradigm of the study N. Halko et al 2011 SIAM Rev. 53 217–288, and hence regularization is performed based on the low-rank matrix approximation. Error bounds for the mean error are obtained which take into account solution smoothness and the inherent noise level. Based on the structure of the error decomposition the authors propose a novel algorithm which guarantees (on the mean) a prescribed error tolerance. Numerical simulations confirm the theoretical findings.

Online adaptive basis refinement and compression for reduced-order models via vector-space sieving

In many applications, projection-based reduced-order models (ROMs) have demonstrated the ability to provide rapid approximate solutions to high-fidelity full-order models (FOMs). However, there is no a priori assurance that these approximate solutions are accurate; their accuracy depends on the ability of the low-dimensional trial basis to represent the FOM solution. As a result, ROMs can generate inaccurate approximate solutions, e.g., when the FOM solution at the online prediction point is not well represented by training data used to construct the trial basis. To address this fundamental deficiency of standard model-reduction approaches, this work proposes a novel online-adaptive mechanism for efficiently enriching the trial basis in a manner that ensures convergence of the ROM to the FOM, yet does not incur any FOM solves. The mechanism is based on the previously proposed adaptive h-refinement method for ROMs (carlberg, 2015), but improves upon this work in two crucial ways. First, the proposed method enables basis refinement with respect to any orthogonal basis (not just the Kronecker basis), thereby generalizing the refinement mechanism and enabling it to be tailored to the physics characterizing the problem at hand. Second, the proposed method provides a fast online algorithm for periodically compressing the enriched basis via an efficient proper orthogonal decomposition (POD) method, which does not incur any operations that scale with the FOM dimension. These two features allow the proposed method to serve as (1) a failsafe mechanism for ROMs, as the method enables the ROM to satisfy any prescribed error tolerance online (even in the case of inadequate training), and (2) an efficient online basis-adaptation mechanism, as the combination of basis enrichment and compression enables the basis to adapt online while controlling its dimension.

A Gauss–Jacobi Kernel Compression Scheme for Fractional Differential Equations

A scheme for approximating the kernel w of the fractional $$\alpha $$ -integral by a linear combination of exponentials is proposed and studied. The scheme is based on the application of a composite Gauss–Jacobi quadrature rule to an integral representation of w. This results in an approximation of w in an interval $$[\delta ,T]$$ , with $$0<\delta $$ , which converges rapidly in the number J of quadrature nodes associated with each interval of the composite rule. Using error analysis for Gauss–Jacobi quadratures for analytic functions, an estimate of the relative pointwise error is obtained. The estimate shows that the number of terms required for the approximation to satisfy a prescribed error tolerance is bounded for all $$\alpha \in (0,1)$$ , and that J is bounded for $$\alpha \in (0,1)$$ , $$T>0$$ , and $$\delta \in (0,T)$$ .

Finite-error metrological bounds on multiparameter Hamiltonian estimation

Estimation of multiple parameters in an unknown Hamiltonian is investigated. We present upper and lower bounds on the time required to complete the estimation within a prescribed error tolerance $\ensuremath{\delta}$. The lower bound is given on the basis of the Cram\'er-Rao inequality, where the quantum Fisher information is bounded by the squared evolution time. The upper bound is obtained by an explicit construction of estimation procedures. By comparing the cases with different numbers of Hamiltonian channels, we also find that the few-channel procedure with adaptive feedback and the many-channel procedure with entanglement are equivalent in the sense that they require the same amount of time resource up to a constant factor.

An FC-based spectral solver for elastodynamic problems in general three-dimensional domains

This paper presents a spectral numerical algorithm for the solution of elastodynamics problems in general three-dimensional domains. Based on a recently introduced “Fourier continuation” (FC) methodology for accurate Fourier expansion of non-periodic functions, the proposed approach possesses a number of appealing properties: it yields results that are essentially free of dispersion errors, it entails mild CFL constraints, it runs at a cost that scales linearly with the discretization sizes, and it lends itself easily to efficient parallelization in distributed-memory computing clusters. The proposed algorithm is demonstrated in this paper by means of a number of applications to problems of isotropic elastodynamics that arise in the fields of materials science and seismology. These examples suggest that the new approach can yield solutions within a prescribed error tolerance by means of significantly smaller discretizations and shorter computing times than those required by other methods.

An error estimator for real-time simulators based on model order reduction

Model order reduction is one of the most appealing choices for real-time simulation of non-linear solids. In this work a method is presented in which real time performance is achieved by means of the off-line solution of a (high dimensional) parametric problem that provides a sort of response surface or computational vademecum. This solution is then evaluated in real-time at feedback rates compatible with haptic devices, for instance (i.e., more than 1 kHz). This high dimensional problem can be solved without the limitations imposed by the curse of dimensionality by employing proper generalized decomposition (PGD) methods. Essentially, PGD assumes a separated representation for the essential field of the problem. Here, an error estimator is proposed for this type of solutions that takes into account the non-linear character of the studied problems. This error estimator allows to compute the necessary number of modes employed to obtain an approximation to the solution within a prescribed error tolerance in a given quantity of interest.

An efficient algorithm based on splitting for the time integration of the Schrödinger equation

We present a practical algorithm based on symplectic splitting methods intended for the numerical integration in time of the Schrödinger equation when the Hamiltonian operator is either time-independent or changes slowly with time. In the later case, the evolution operator can be effectively approximated in a step-by-step manner: first divide the time integration interval in sufficiently short subintervals, and then successively solve a Schrödinger equation with a different time-independent Hamiltonian operator in each of these subintervals. When discretized in space, the Schrödinger equation with the time-independent Hamiltonian operator obtained for each time subinterval can be recast as a classical linear autonomous Hamiltonian system corresponding to a system of coupled harmonic oscillators. The particular structure of this linear system allows us to construct a set of highly efficient schemes optimized for different precision requirements and time intervals. Sharp local error bounds are obtained for the solution of the linear autonomous Hamiltonian system considered in each time subinterval. Our schemes can be considered, in this setting, as polynomial approximations to the matrix exponential in a similar way as methods based on Chebyshev and Taylor polynomials. The theoretical analysis, supported by numerical experiments performed for different time-independent Hamiltonians, indicates that the new methods are more efficient than schemes based on Chebyshev polynomials for all tolerances and time interval lengths. The algorithm we present automatically selects, for each time subinterval, the most efficient splitting scheme (among several new optimized splitting methods) for a prescribed error tolerance and given estimates of the upper and lower bounds of the eigenvalues of the discretized version of the Hamiltonian operator.

Confidence-based reasoning in stochastic constraint programming

In this work we introduce a novel approach, based on sampling, for finding assignments that are likely to be solutions to stochastic constraint satisfaction problems and constraint optimisation problems. Our approach reduces the size of the original problem being analysed; by solving this reduced problem, with a given confidence probability, we obtain assignments that satisfy the chance constraints in the original model within prescribed error tolerance thresholds. To achieve this, we blend concepts from stochastic constraint programming and statistics. We discuss both exact and approximate variants of our method. The framework we introduce can be immediately employed in concert with existing approaches for solving stochastic constraint programs. A thorough computational study on a number of stochastic combinatorial optimisation problems demonstrates the effectiveness of our approach.

Novel adaptive time stepping method and its application to soil seismic liquefaction analysis

In solid–fluid coupled and dynamic analysis, the temporal discretization error related to the time step size is unavoidable. To improve the calculation efficiency, a novel adaptive time stepping procedure based on a finite element and finite difference (FEM–FDM) coupled method is proposed for soil seismic liquefaction analysis. The core concept of this adaptive stepping method is the mixed displacement and pore water pressure error estimation, which is obtained by an embedded error estimator, and a time stepping strategy, which is operated according to the relation between the current mixed error and the prescribed error tolerance. Two numerical examples were performed to validate the proposed method. It is shown that under the same condition of mesh size and other numerical parameters, the time step size obviously affects the calculation results; using adaptive stepping method is economical, robust and has the same degree of accuracy as compared with the fixed stepping method in the soil earthquake liquefaction analysis.

Adaptive h‐refinement for reduced‐order models

SummaryThis work presents a method to adaptively refine reduced‐order models a posteriori without requiring additional full‐order‐model solves. The technique is analogous to mesh‐adaptive h‐refinement: it enriches the reduced‐basis space online by ‘splitting’ a given basis vector into several vectors with disjoint support. The splitting scheme is defined by a tree structure constructed offline via recursive k‐means clustering of the state variables using snapshot data. The method identifies the vectors to split online using a dual‐weighted‐residual approach that aims to reduce error in an output quantity of interest. The resulting method generates a hierarchy of subspaces online without requiring large‐scale operations or full‐order‐model solves. Further, it enables the reduced‐order model to satisfy any prescribed error tolerance regardless of its original fidelity, as a completely refined reduced‐order model is mathematically equivalent to the original full‐order model. Experiments on a parameterized inviscid Burgers equation highlight the ability of the method to capture phenomena (e.g., moving shocks) not contained in the span of the original reduced basis. Copyright © 2014 John Wiley &amp; Sons, Ltd.

Poisson Preconditioning for Self-Adjoint Elliptic Problems

In this paper, we formulate interface problem and Neumann elliptic boundary value problem into a form of linear operator equations with self-adjoint positive definite operators. We prove that in the discrete level the condition number of these operators is independent of the mesh size. Therefore, given a prescribed error tolerance, the classical conjugate gradient algorithm converges within a fixed number of iterations. The main computation task at each iteration is to solve a Dirichlet Poisson boundary value problem in a rectangular domain, which can be furnished with fast Poisson solver. The overall computational complexity is essentially of linear scaling.

Higher-Order Linear-Time Unconditionally Stable Alternating Direction Implicit Methods for Nonlinear Convection-Diffusion Partial Differential Equation Systems

We introduce a class of alternating direction implicit (ADI) methods, based on approximate factorizations of backward differentiation formulas (BDFs) of order p≥2, for the numerical solution of two-dimensional, time-dependent, nonlinear, convection-diffusion partial differential equation (PDE) systems in Cartesian domains. The proposed algorithms, which do not require the solution of nonlinear systems, additionally produce solutions of spectral accuracy in space through the use of Chebyshev approximations. In particular, these methods give rise to minimal artificial dispersion and diffusion and they therefore enable use of relatively coarse discretizations to meet a prescribed error tolerance for a given problem. A variety of numerical results presented in this text demonstrate high-order accuracy and, for the particular cases of p=2,3, unconditional stability.