Numerical Truncation Error Research Articles

Overview of Factorisation Methods in Kalman Filtering

Abstract The paper summarises and describes the most commonly used matrix factorisation methods applied in design of the Kalman filter in order to improve computational efficiency and avoid divergence issues caused by numerical round-off and truncation errors. Some forms of the Kalman filter are more prone to the growth of numerical error sand possible divergence than other implementations. In order to prevent the algorithm’s divergence additional processing is needed and this paper discusses pros and cons of different implementations and their numerical characteristics. Numerical issues still arise in finite word length implementations of algorithms, which frequently occur in embedded systems. This paper describes algorithms based on different factorisations such as Cholesky, U-D, SVD and their basic numerical properties.

A High-Order, Conservative Integrator with Local Time-Stepping

We present a family of multistep integrators based on the Adams--Bashforth methods. These schemes can be constructed for arbitrary convergence order with arbitrary step size variation. The step size can differ between different subdomains of the system. It can also change with time within a given subdomain. The methods are linearly conservative, preserving a wide class of analytically constant quantities to numerical roundoff, even when numerical truncation error is significantly higher. These methods are intended for use in solving conservative PDEs in discontinuous Galerkin formulations or in finite-difference methods with compact stencils. A numerical test demonstrates these properties and shows that significant speed improvements over the standard Adams--Bashforth schemes can be obtained.

On the automatic and a priori design of unstructured mesh resolution for coastal ocean circulation models

This study investigates the design of unstructured mesh resolution and its impact on the modeling of barotropic tides along the United States East Coast and Gulf Coast (ECGC). A discrete representation of a computational ocean domain (mesh design) is necessary due to finite computational resources and an incomplete knowledge of the physical system (e.g., shoreline and seabed topography). The selection of mesh resolution impacts both the numerical truncation error and the approximation of the system’s physical domain. To increase confidence in the design of high-resolution coastal ocean meshes and to quantify the efficacy of current mesh design practices, an automated mesh generation approach is applied to objectively control resolution placement based on a priori information such as shoreline geometry and seabed topographic features. The simulated harmonic tidal elevations for each mesh design are compared to that of a reference solution, computed on a 10.8 million vertex mesh of the ECGC region with a minimum shoreline resolution of 50 m. Our key findings indicate that existing mesh designs that use uniform resolution along the shoreline and slowly varying resolution sizes on the continental shelf inefficiently discretize the computational domain. Instead, a targeted approach that places fine resolution in narrow geometric features, along steep topographic gradients, and along pronounced submerged estuarine channels, while aggressively relaxing resolution elsewhere, leads to a mesh with an order of magnitude fewer vertices than the reference solution with comparable accuracy (within 3% harmonic elevation amplitudes in 99% of the domain).

Variable coefficient depletion equations – Reference solution based on Rational Approximation Methods

High order neutronic-depletion coupling strategies achieve better simulation performance against traditional strategies, at the cost of demanding variable coefficient depletion equations rather than constant ones to be solved. Through results comparisons of Rational Approximation solvers with various settings, different types of errors are explicitly separated and analyzed. In particular, sub-step errors of most nuclides are found to suffer from slow convergent problem; method intrinsic errors of the 14-order Chebyshev Rational Approximation Method are found to be large enough to interrupt convergence processes of sub-step errors prematurely. In addition, the fading behavior of the sub-step errors is captured and explained, and it leads to the optimal sub-step division appears to have moderate decreasing sub-step sizes; also, numerical truncation errors are shown to be dominated by method intrinsic errors. Accordingly, the extended precision 64-order Quadrature Rational Approximation Method and the Richardson Extrapolation technique are combined to obtain reference solutions.

M-PCM-OFFD: An effective output statistics estimation method for systems of high dimensional uncertainties subject to low-order parameter interactions

The evaluation of output performance statistics for systems of high-dimensional uncertain input parameters is crucial for robust real-time decision-making tasks of large-scale complex systems that operate in an uncertain environment. We develop a framework that integrates Multivariate Probabilistic Collocation Method (M-PCM) and Orthogonal Fractional Factorial Design (OFFD) to achieve an effective and scalable output statistics estimation. In this paper, we prove that when the degree of each uncertain parameter does not exceed 3 and under the widely held assumption for high-dimensional systems that the interactions among uncertain input parameters are negligible beyond certain order, the integrated M-PCM–OFFD method breaks the curse of dimensionality for correct output mean estimation by maximally reducing the number of simulations from 22m to 2log2(m+1) for a system mapping of m uncertain input parameters. In addition, the resulting reduced-size simulation set is the most robust to numerical truncation errors of simulators among all subsets of the same size in the M-PCM simulation set. The analysis also provides new insightful formal interpretations of the optimality of OFFDs.

Enhancing the Performance of the Distributed Gauss-Newton Optimization Method by Reducing the Effect of Numerical Noise and Truncation Error With Support-Vector Regression

Summary Numerical optimization is an integral part of many history-matching (HM) workflows. However, the optimization performance can be affected negatively by the numerical noise existent in the forward models when the gradients are estimated numerically. As an unavoidable part of reservoir simulation, numerical noise refers to the error caused by the incomplete convergence of linear or nonlinear solvers or truncation errors caused by different timestep cuts. More precisely, the allowed solver tolerances and allowed changes of pressure and saturation imply that simulation results no longer smoothly change with changing model parameters. For HM with linear-distributed Gaussian-Newton (L-DGN), caused by the discontinuity of simulation results, the sensitivity matrix computed by linear interpolation might be less accurate, which might result in slow convergence or, even worse, failure of convergence. Recently, we have developed an HM workflow by integrating the support-vector regression (SVR) with the distributed-Gaussian-Newton (DGN) method optimization method referred to as SVR-DGN. Unlike L-DGN that computes the sensitivity matrix with a simple linear proxy, SVR-DGN computes the sensitivity matrix by taking the gradient of the SVR proxies. In this paper, we provide theoretical analysis and case studies to show that SVR-DGN can compute a more-accurate sensitivity matrix than L-DGN, and SVR-DGN is insensitive to the negative influence of numerical noise. We also propose a cost-saving training procedure by replacing bad-training points, which correspond to relatively large values of the objective function, with those training-data points (simulation data) that have smaller values of the objective function and are generated at most-recent iterations for training the SVR proxies. Both the L-DGN approach and the newly proposed SVR-DGN approach are tested first with a 2D toy problem to show the effect of numerical noise on their convergence performance. We find that their performance is comparable when the toy problem is free of numerical noise. As the numerical-noise level increases, the performance of the L-DGN degrades sharply. By contrast, the SVR-DGN performance is quite stable. Then, both methods are tested using a real-field HM example. The convergence performance of the SVR-DGN is quite robust for both the tight and loose numerical settings, whereas the performance of the L-DGN degrades significantly when loose numerical settings are applied.

Meshfree simulations of ultrasound vector flow imaging using smoothed particle hydrodynamics

Before embarking on a series of in vivo tests, design of ultrasound-flow-imaging modalities are generally more efficient through computational models as multiple configurations can be tested methodically. To that end, simulation models must generate realistic blood flow dynamics and Doppler signals. The current in silico ultrasound simulation techniques suffer mainly from uncertainty in providing accurate trajectories of moving ultrasound scatterers. In mesh-based Eulerian methods, numerical truncation errors from the interpolated velocities, both in the time and space dimensions, can accumulate significantly and make the pathlines unreliable. These errors can distort beam-to-beam inter-correlation present in ultrasound flow imaging. It is thus a technical issue to model a correct motion of the scatterers by considering their interaction with boundaries and neighboring scatterers. We hypothesized that in silico analysis of emerging ultrasonic imaging modalities can be implemented more accurately with meshfree approaches. We developed an original fluid-ultrasound simulation environment based on a meshfree Lagrangian CFD (computational fluid dynamics) formulation, which allows analysis of ultrasound flow imaging. This simulator combines smoothed particle hydrodynamics (SPH) and Fourier-domain linear acoustics (SIMUS = simulator for ultrasound imaging). With such a particle-based computation, the fluid particles also acted as individual ultrasound scatterers, resulting in a direct and physically sound fluid-ultrasonic coupling. We used the in-house algorithms for fluid and ultrasound simulations to simulate high-frame-rate vector flow imaging. The potential of the particle-based method was tested in 2D simulations of vector Doppler for the intracarotid flow. The Doppler-based velocity fields were compared with those issued from SPH. The numerical evaluations showed that the vector flow fields obtained by vector Doppler components were in good agreement with the original SPH velocities, with relative errors less than 10% and 2% in the cross-beam and axial directions, respectively. Our results showed that SPH-SIMUS coupling enables direct and realistic simulations of ultrasound flow imaging. The proposed coupled algorithm has also the advantage to be 3D compatible and parallelizable.

Numerical Modeling and Validation of a Novel 2D Compositional Flooding Simulator Using a Second-Order TVD Scheme

The aim of this paper is to present the latter and develop a numerical simulator aimed at solving a 2D domain porous medium, using the compositional approach to simulate chemical flooding processes. The simulator consists in a two-phase, multicomponent system solved by the IMplicit in Pressure, Explicit in Concentration (IMPEC) approach, which can be operated under an iterative/non-iterative condition on each time-step. The discretization of the differential equations is done using a fully second order of accuracy, along with a Total Variation Diminishing (TVD) scheme with a flux limiter function. This allowed reducing the artificial diffusion and dispersion on the transport equation, improving the chemical species front tracking, decreasing the numerical influence on the recovery results. The new model was validated against both commercial and academic simulators and moreover, the robustness and stability were also tested, showing that the iterative IMPEC is fully stable, behaving as an implicit numerical scheme. The non-iterative IMPEC is conditionally stable, with a critical time-step above which numerical spurious oscillations begin to appear until the system numerically crashes. The results showed a good correspondence in different grid sizes, being largely affected by the time-step, with caused a decrease in the recovery efficiency in the iterative scheme, and the occurrence of numerical oscillations in the non-iterative one. Numerically speaking, the second-order scheme using a flux splitting TVD discretization proved to be a good approach for compositional reservoir simulation, decreasing the influence of numerical truncation errors on the results when compared to traditional, first-order linear schemes. Along with these studies, secondary recoveries in constant and random permeability fields are simulated before employing them in tertiary recovery processes.

Uncertainty propagation in the N-body problem using Dromo elements

Reliable and efficient uncertainty propagation is crucial for the task of monitoring possible impacts of Near Earth Asteroids with our planet. It is well known that a switch of the primary body can greatly reduce the numerical truncation error in the case of planetary flybys. In the present work, the advantages of performing a primary body switch in the uncertainty propagation problem are explored. For this purpose, first we present the linear uncertainty propagation using Dromo formulation, which has been shown in previous works to have a satisfactory performance when propagating the orbit uncertainty of Near Earth Asteroids. Next, we introduce the concept of primary body switch for the uncertainty propagation problem. The algorithm is based on the following procedure. We sample the initial orbit uncertainty distribution and linearly propagate the samples, considering the N-body gravitational influence. A primary body switch, which is a nonlinear mapping, is performed for all the samples at a threshold distance from the approaching planet, the Earth. The orbit of the samples is then linearly propagated with respect to the Earth using the Dromo formulation until the threshold distance is reached again. Finally, the propagation center is changed back to the heliocentric frame and the linear propagation continues. We apply the proposed method to an extensive set of asteroids that approach the Earth. Results suggest that the average error of the linear propagation can be reduced up to a factor of 30 when compared to a purely heliocentric linear propagation using Dromo elements.

Implementation and Comparison of High-Resolution Spatial Discretization Schemes for Solving Two-Fluid Seven-Equation Two-Pressure Model

As compared to the two-fluid single-pressure model, the two-fluid seven-equation two-pressure model has been proved to be unconditionally well-posed in all situations, thus existing with a wide range of industrial applications. The classical 1st-order upwind scheme is widely used in existing nuclear system analysis codes such as RELAP5, CATHARE, and TRACE. However, the 1st-order upwind scheme possesses issues of serious numerical diffusion and high truncation error, thus giving rise to the challenge of accurately modeling many nuclear thermal-hydraulics problems such as long term transients. In this paper, a semi-implicit algorithm based on the finite volume method with staggered grids is developed to solve such advanced well-posed two-pressure model. To overcome the challenge from 1st-order upwind scheme, eight high-resolution total variation diminishing (TVD) schemes are implemented in such algorithm to improve spatial accuracy. Then the semi-implicit algorithm with high-resolution TVD schemes is validated on the water faucet test. The numerical results show that the high-resolution semi-implicit algorithm is robust in solving the two-pressure two-fluid two-phase flow model; Superbee scheme and Koren scheme give two highest levels of accuracy while Minmod scheme is the worst one among the eight TVD schemes.

Quantifying spatial distribution of spurious mixing in ocean models

Numerical mixing is inevitable for ocean models due to tracer advection schemes. Until now, there is no robust way to identify the regions of spurious mixing in ocean models. We propose a new method to compute the spatial distribution of the spurious diapycnic mixing in an ocean model. This new method is an extension of available potential energy density method proposed by Winters and Barkan (2013). We test the new method in lock-exchange and baroclinic eddies test cases. We can quantify the amount and the location of numerical mixing. We find high-shear areas are the main regions which are susceptible to numerical truncation errors. We also test the new method to quantify the numerical mixing in different horizontal momentum closures. We conclude that Smagorinsky viscosity has less numerical mixing than the Leith viscosity using the same non-dimensional constant.

Effect of Numerical Truncation Error on Implicit Finite Difference Methods in Groundwater Transport Models

A correction for truncation errors associated with implicit finite difference method (FDM) of the groundwater transport equation with reaction term, which is generally used in groundwater transport models, is developed here from a Taylor series analysis. An application example is formulated to illustrate the effect of truncation errors on the numerical solution of implicit FDM. The study compares the effect of truncation error on numerical model accuracy of implicit FDM and explicit FDM for groundwater transport equation with reaction term. The explicit FDM constitutes large deviation from analytical solution without truncation error correction. The relative error analysis reveals that error reduces from 75% to 30% after truncation error correction. Therefore, we can conclude that numerical truncation error correction has significant impact on accuracy of groundwater transport models based on explicit FDM than those based on implicit FDM.

On the accuracy and precision of numerical waveforms: effect of waveform extraction methodology

We present a new set of 95 numerical relativity simulations of non-precessing binary black holes (BBHs). The simulations sample comprehensively both black-hole spins up to spin magnitude of 0.9, and cover mass ratios 1–3. The simulations cover on average 24 inspiral orbits, plus merger and ringdown, with low initial orbital eccentricities . A subset of the simulations extends the coverage of non-spinning BBHs up to mass ratio q = 10. Gravitational waveforms at asymptotic infinity are computed with two independent techniques: extrapolation and Cauchy characteristic extraction. An error analysis based on noise-weighted inner products is performed. We find that numerical truncation error, error due to gravitational wave extraction, and errors due to the Fourier transformation of signals with finite length of the numerical waveforms are of similar magnitude, with gravitational wave extraction errors dominating at noise-weighted mismatches of . This set of waveforms will serve to validate and improve aligned-spin waveform models for gravitational wave science.

Constructing reference metrics on multicube representations of arbitrary manifolds

Reference metrics are used to define the differential structure on multicube representations of manifolds, i.e., they provide a simple and practical way to define what it means globally for tensor fields and their derivatives to be continuous. This paper introduces a general procedure for constructing reference metrics automatically on multicube representations of manifolds with arbitrary topologies. The method is tested here by constructing reference metrics for compact, orientable two-dimensional manifolds with genera between zero and five. These metrics are shown to satisfy the Gauss–Bonnet identity numerically to the level of truncation error (which converges toward zero as the numerical resolution is increased). These reference metrics can be made smoother and more uniform by evolving them with Ricci flow. This smoothing procedure is tested on the two-dimensional reference metrics constructed here. These smoothing evolutions (using volume-normalized Ricci flow with DeTurck gauge fixing) are all shown to produce reference metrics with constant scalar curvatures (at the level of numerical truncation error).

Time-Space Decoupled Model with a Variable-Coefficient Dispersive Condition to Simulate Tsunamis over Slowly Varying Topography

AbstractIn this study, a new numerical scheme with variable-coefficient dispersive conditions was developed to simulate tsunamis over slowly varying topography. The new model was developed using a novel time-space decoupled numerical scheme that used the classic wave equation. By adapting this model to simulate tsunamis, a variable-coefficient dispersive condition was derived to consider the dispersive effects using a numerical truncation error from the numerical scheme. The stability condition and range of the variable coefficient were analyzed using a stability analysis. Using the specific dispersive condition and the modified time-space decoupled scheme, a uniform grid system can be used to simulate a dispersive tsunami over slowly varying topography. Numerical experiments for wave propagation with Gaussian humps and dipolar sources in uniform and varying topography demonstrate that the scheme is efficient for simulating tsunami propagation over slowly varying water depths with dispersive effects.

Dispersion phenomena in symmetric pre-stressed layered elastic structures

The dispersion relation associated with a symmetric three layer structure, composed of compressible, pre-stressed elastic layers, is derived. This mathematically elaborate transcendental equation gives phase speed as an implicit function of wave number. Numerical solutions are established to show a wide range of dispersion behaviour which is delicately dependent on the material parameters and pre-stress in each layer. Particularly interesting behaviour is observed within the short wave (high wave number) regime, with six possible cases of short wave liming behaviour shown possible. Within each of these, a short wave asymptotic analysis is carried out, resulting in a set of approximations which provide explicit relationships between phase speed and wave number. It is envisaged that these approximations may prove helpful to approximate numerical truncation errors associated with impact response, as well as providing excellent first approximations for particularly (numerically) challenging sets of material parameters.

Generalized Bondi-Sachs equations for characteristic formalism of numerical relativity

The Cauchy formalism of numerical relativity has been successfully applied to simulate various dynamical spacetimes without any symmetry assumption. But discovering how to set a mathematically consistent and physically realistic boundary condition is still an open problem for Cauchy formalism. In addition, the numerical truncation error and finite region ambiguity affect the accuracy of gravitational wave form calculation. As to the finite region ambiguity issue, the characteristic extraction method helps much. But it does not solve all of the above issues. Besides the above problems for Cauchy formalism, the computational efficiency is another problem. Although characteristic formalism of numerical relativity suffers the difficulty from caustics in the inner near zone, it has advantages in relation to all of the issues listed above. Cauchy-characteristic matching (CCM) is a possible way to take advantage of characteristic formalism regarding these issues and treat the inner caustics at the same time. CCM has difficulty treating the gauge difference between the Cauchy part and the characteristic part. We propose generalized Bondi-Sachs equations for characteristic formalism for the Cauchy-characteristic matching end. Our proposal gives out a possible same numerical evolution scheme for both the Cauchy part and the characteristic part. And our generalized Bondi-Sachs equations have one adjustable gauge freedom which can be used to relate the gauge used in the Cauchy part. Then these equations can make the Cauchy part and the characteristic part share a consistent gauge condition. So our proposal gives a possible new starting point for Cauchy-characteristic matching.

Cascades, backscatter and conservation in numerical models of two‐dimensional turbulence

AbstractThe equations governing atmospheric flow imply transfers of energy and potential enstrophy between scales. Accurate simulation of turbulent flow requires that numerical models, which have finite resolution and truncation errors, adequately capture these interscale transfers, particularly between resolved and unresolved scales. It is therefore important to understand how accurately these transfers are modelled in the presence of scale‐selective dissipation or other forms of subgrid model. Here, the energy and enstrophy cascades in numerical models of two‐dimensional turbulence are investigated using the barotropic vorticity equation.Energy and enstrophy transfers in spectral space due to truncated scales are calculated for a high‐resolution reference solution and for several explicit and implicit subgrid models at coarser resolution. The reference solution shows that enstrophy and energy are removed from scales very close to the truncation scale and energy is transferred (backscattered) into the large scales. Some subgrid models are able to capture the removal of enstrophy from small scales, though none are scale‐selective enough; however, none are able to capture accurately the energy backscatter.We propose a scheme that perturbs the vorticity field at each time step by the addition of a particular vorticity pattern derived by filtering the predicted vorticity field. Although originally conceived as a parametrization of energy backscatter, this scheme is best interpreted as an energy ‘fixer’ that attempts to repair the damage to the energy spectrum caused by numerical truncation error and an imperfect subgrid model. The proposed scheme improves the energy and enstrophy behaviour of the solution and, in most cases, slightly reduces the root mean square vorticity errors.

NUMERICAL STUDY OF THE INTERACTION BETWEEN POSITIVE AND NEGATIVE MASS OBJECTS IN GENERAL RELATIVITY

Based on Baumgarte–Shapiro–Shibata–Nakamura formalism and moving puncture method, we demonstrate the first numerical evolutions of the interaction between positive and negative mass objects. Using the causal property of general relativity, we set our computational domain around the positive mass black hole while excluding the region around the naked singularity introduced by the negative mass object. Besides the usual Sommerfeld numerical boundary condition, an approximate boundary condition is proposed for this nonasymptotically-flat computational domain. Careful checks show that either boundary condition introduces smaller error than the numerical truncation errors. This is consistent with the causal property of general relativity. Except for the numerical truncation error and round-off error, our method gives an exact solution to the full Einstein's equation for a portion of spacetime with two objects whose masses have opposite signs. So our method opens the door for numerical explorations with negative mass objects. Based on this method, we investigate the Newtonian limit of spacetime with two objects whose masses have opposite sign. Our result implies that this spacetime does have a Newtonian limit which corresponds to a negative mass particle chasing a positive mass particle. This result sheds some light on an interesting debate about the Newtonian limit of a spacetime with positive and negative point masses.

Evaluation of a general CFD-solver for a micro-scale urban flow

In this work, we carry out CFD-simulations of the Mock Urban Setting Test (MUST) case in the framework of the COST action 732. The main objective is to assess the quality assurance protocol and guidelines constructed in the COST action 732, and to study the role of the numerical truncation error as part of the total error. The computations are performed with a general flow solver FINFLO. The steady-state RANS approach is used with k-ε and k-ω SST turbulence models. In the near-wall region, standard Wall-Function Approach (WFA) is employed and the tracer is modelled as a passive scalar. The CFD results are compared to wind-tunnel data.