High-order Extrapolation Research Articles

Energetic spectral-element time marching methods for phase-field nonlinear gradient systems

We propose two efficient energetic spectral-element methods in time for marching nonlinear gradient systems with the phase-field Allen–Cahn equation as an example: one fully implicit nonlinear method and one semi-implicit linear method. Different from other spectral methods in time using spectral Petrov-Galerkin or weighted Galerkin approximations, the presented implicit method employs an energetic variational Galerkin form that can maintain the mass conservation and energy dissipation property of the continuous dynamical system. Another advantage of this method is its superconvergence. A high-order extrapolation is adopted for the nonlinear term to get the semi-implicit method. The semi-implicit method does not have superconvergence, but can be improved by a few Picard-like iterations to recover the superconvergence of the implicit method. Numerical experiments verify that the method using Legendre elements of degree three outperforms the 4th-order implicit-explicit backward differentiation formula and the 4th-order exponential time difference Runge-Kutta method, which were known to have best performances in solving phase-field equations. In addition to the standard Allen–Cahn equation, we also apply the method to a conservative Allen–Cahn equation, in which the conservation of discrete total mass is verified. The applications of the proposed methods are not limited to phase-field Allen–Cahn equations. They are suitable for solving general, large-scale nonlinear dynamical systems.

Thermodynamic modeling of KCl-PrCl3 and KCl-LiCl-PrCl3 systems

Molten salt electrolysis can recover the actinides from spent nuclear fuels, and it involves a eutectic LiCl-KCl in molten form as an electrolyte. During reprocessing, the concentration of fission products such as La, Nd, Pr, etc., increases and affects the recovery efficiency of the electrolyte. In this work, thermodynamic modeling of KCl-PrCl3 and KCl-LiCl-PrCl3 systems was carried out using the CALPHAD (Calculation of Phase Diagrams) approach for the first time. The thermodynamic functions for the pure salts were taken from the SGTE (Scientific Group Thermodata Europe) Substances (SSUB) database. The experimental thermochemical and phase equilibria data available in the literature were used as input for the assessment of KCl-PrCl3 and KCl-LiCl-PrCl3 systems. The model parameters for the KCl-LiCl system were adjusted to include the new Gibbs energy descriptions for the pure salts. In addition, the sublattice model for the liquid phase in the LiCl-PrCl3 system was modified and reassessed to ensure the model compatibility for higher-order extrapolation. There is a good agreement between the experimental and calculated thermochemical and phase diagram data for all the optimized constituent binaries and the ternary system. This work will be beneficial for determining the solubility limit of PrCl3 in molten LiCl-KCl electrolytes and their thermodynamic properties for improving the efficiency of the pyrochemical process.

Solving SEI Model using Non-Standard Finite Difference and High Order Extrapolation with Variable Step Length

A high-level method was obtained to solve the SEI model problem involving Symmetrization measures in numerical calculations through the Implicit Midpoint Rule method (IMR). It is obtained using Non-Standard Finite Difference Schemes (NSFD) with Extrapolation techniques combined. In solving differential equation problems numerically, the Extrapolated SEI model method is able to generate more accurate results than the existing numerical method of SEI model. This study aims to investigate the accuracy and efficiency of computing between Extrapolated One-Step Active Symmetry Implicit Midpoint Rule method (1ASIMR), Extrapolated One-Step Active Symmetry Implicit Midpoint Rule method (2ASIMR), Extrapolated One-Step Passive Symmetry Midpoint Rule method (1PSIMR) and the extrapolated Two-Step Passive Symmetry Midpoint Rule method (2PSIMR). The results show that the 1ASIMR method is the most accurate method. For the determination of the efficiency of 2ASIMR and 2PSIMR methods have high efficiency. At the end of the study, the results from the numerical method obtained show that Extrapolation using Non-Standard Finite Difference has higher accuracy than the existing Implicit Midpoint Rule method.

Extremely high-order convergence in simulations of relativistic stars

We provide a road towards obtaining gravitational waveforms from inspiraling material binaries with an accuracy viable for third-generation gravitational wave detectors, without necessarily advancing computational hardware or massively-parallel software infrastructure. We demonstrate a proof-of-principle 1 + 1-dimensional numerical implementation that exhibits up to 7th-order convergence for highly dynamic barotropic stars in curved spacetime, and numerical errors up to 6 orders of magnitude smaller than a standard method. Aside from errors associated with overshoot in high-order extrapolation from the interior of the star to the surface, there are no obvious fundamental obstacles to obtaining convergence of even higher order. The implementation uses a novel surface-tracking method, where the surface is evolved and high-order accurate boundary conditions are imposed there. Computational memory does not need to be allocated to fluid variables in the vacuum region of spacetime. We anticipate the application of this new method to full 3 + 1-dimensional simulations of the inspiral phase of compact binary systems with at least one material body. The additional challenge of a deformable surface must be addressed in multiple spatial dimensions, but it is also an opportunity to input more precise surface tension physics.

An asymptotic local solution for axi-symmetric heat conduction problems with temperature discontinuity on the boundary

An asymptotic local solution to the axisymmetric conduction problem with a discontinuous temperature boundary condition in the radial direction is presented. The Laplace's equation in cylindrical coordinates (r, z) is expanded near the discontinuity at r = 1, z = 0. The leading order term arises from the step-change in the boundary condition. The next term reflects the effect of curvature at the position of the temperature discontinuity and is a particular solution to the Poisson's equation. Subsequent higher-order terms involve general solutions to the Laplace's equation in Cartesian coordinates and depend on conditions on the far end boundaries. Three domain sizes are thus considered: (i) infinite solid cylinder (0 ≤ r < ∞, 0 ≤ z < ∞); (ii) semi-infinite solid cylinder with finite radius (0 ≤ r < a, 0 ≤ z < ∞); and (iii) finite solid cylinder (0 ≤ r ≤ a, 0 ≤ z ≤ a). Exact solutions for temperature in those three cases are used to determine the coefficients in the higher order terms of the asymptotic solution. As a result, a simple expression for local normal wall heat flux near the discontinuity is obtained and the accuracy is confirmed by comparing with the heat flux obtained through high order extrapolation of the exact temperature field near the wall.

Convergence of explicitly coupled simulation tools (co-simulations)

Abstract In engineering, it is a common desire to couple existing simulation tools together into one big system by passing information from subsystems as parameters into the subsystems under influence. As executed at fixed time points, this data exchange gives the global method a strong explicit component. Globally, such an explicit co-simulation schemes exchange time step can be seen as a step of an one-step method which is explicit in some solution components. Exploiting this structure, we give a convergence proof for such schemes. As flows of conserved quantities are passed across subsystem boundaries, it is not ensured that system-wide balances are fulfilled: the system is not solved as one single equation system. These balance errors can accumulate and make simulation results inaccurate. Use of higher-order extrapolation in exchanged data can reduce this problem but cannot solve it. The remaining balance error has been handled in past work by recontributing it to the input signal in next coupling time step, a technique labeled balance correction methods. Convergence for that method is proven. Further, the lack of stability for co-simulation schemes with and without balance correction is stated.

Limitations of Richardson Extrapolation for Kernel Density Estimation

This paper develops the process of using Richardson Extrapolation to improve the Kernel Density Estimation method, resulting in a more accurate (lower Mean Squared Error) estimate of a probability density function for a distribution of data in $R_d$ given a set of data from the distribution. The method of Richardson Extrapolation is explained, showing how to fix conditioning issues that arise with higher-order extrapolations. Then, it is shown why higher-order estimators do not always provide the best estimate, and it is discussed how to choose the optimal order of the estimate. It is shown that given n one-dimensional data points, it is possible to estimate the probability density function with a mean squared error value on the order of only $n^{-1}\sqrt{\ln(n)}$. Finally, this paper introduces a possible direction of future research that could further minimize the mean squared error.

On Multistep Stabilizing Correction Splitting Methods with Applications to the Heston Model

In this note we consider splitting methods based on linear multistep methods and stabilizing corrections. To enhance the stability of the methods, we employ an idea of Bruno and Cubillos [O. P. Bruno and M. Cubillos, J. Comput. Phys., 307 (2016), pp. 476--495], who combine a high-order extrapolation formula for the explicit term with a formula of one order lower for the implicit terms. Several examples of the obtained multistep stabilizing correction methods are presented, and results on linear stability and convergence are derived. The methods are tested in the application to the well-known Heston model arising in financial mathematics and are found to be competitive with well-established one-step splitting methods from the literature.

Guided lock of a suspended optical cavity enhanced by a higher-order extrapolation.

Lock acquisition of a suspended optical cavity can be a highly stochastic process and is therefore nontrivial. Guided lock is a method to make lock acquisition less stochastic by decelerating the motion of the cavity length based on an extrapolation of the motion from an instantaneous velocity measurement. We propose an improved scheme that is less susceptible to seismic disturbances by incorporating the acceleration as a higher-order correction in the extrapolation. We implemented the new scheme in a 300-m suspended Fabry-Perot cavity and improved the success rate of lock acquisition by a factor of 30.

On the Quasi-unconditional Stability of BDF-ADI Solvers for the Compressible Navier--Stokes Equations and Related Linear Problems

The companion paper "Higher-order in time quasi-unconditionally stable ADI solvers for the compressible Navier-Stokes equations in 2D and 3D curvilinear domains", which is referred to as Part I in what follows, introduces ADI (Alternating Direction Implicit) solvers of higher orders of temporal accuracy (orders $s = 2$ to $6$) for the compressible Navier-Stokes equations in two- and three-dimensional space. The proposed methodology employs the backward differentiation formulae (BDF) together with a quasilinear-like formulation, high-order extrapolation for nonlinear components, and the Douglas-Gunn splitting. A variety of numerical results presented in Part I demonstrate in practice the theoretical convergence rates enjoyed by these algorithms, as well as their excellent accuracy and stability properties for a wide range of Reynolds numbers. In particular, the proposed schemes enjoy a certain property of "quasi-unconditional stability": for small enough (problem-dependent) fixed values of the time-step $\Delta t$, these algorithms are stable for arbitrarily fine spatial discretizations. The present contribution presents a mathematical basis for the performance of these algorithms. Short of providing stability theorems for the full BDF-ADI Navier-Stokes solvers, this paper puts forth proofs of unconditional stability and quasi-unconditional stability for BDF-ADI schemes as well as some related un-split BDF schemes, for a variety of related linear model problems in one, two and three spatial dimensions, and for schemes of orders $2\leq s\leq 6$ of temporal accuracy. Additionally, a set of numerical tests presented in this paper for the compressible Navier-Stokes equation indicate that quasi-unconditional stability carries over to the fully non-linear context.

Stability analysis of the inverse Lax–Wendroff boundary treatment for high order upwind-biased finite difference schemes

In this paper, we consider linear stability issues for one-dimensional hyperbolic conservation laws using a class of conservative high order upwind-biased finite difference schemes, which is a prototype for the weighted essentially non-oscillatory (WENO) schemes, for initial–boundary value problems (IBVP). The inflow boundary is treated by the so-called inverse Lax–Wendroff (ILW) or simplified inverse Lax–Wendroff (SILW) procedure, and the outflow boundary is treated by the classical high order extrapolation. A third order total variation diminishing (TVD) Runge–Kutta time discretization is used in the fully discrete case. Both GKS (Gustafsson, Kreiss and Sundström) and eigenvalue analyses are performed for both semi-discrete and fully discrete schemes. The two different analysis techniques yield consistent results. Numerical tests are performed to demonstrate the stability results predicted by the analysis.

On a simple way to calculate electronic resonances for polyatomic molecules.

We propose a simple method for calculation of low-lying shape electronic resonances of polyatomic molecules. The method introduces a perturbation potential and requires only routine bound-state type calculations in the real domain of energies. Such a calculation is accessible by most of the free or commercial quantum chemistry software. The presented method is based on the analytical continuation in a coupling constant model, but unlike its previous variants, we experience a very stable and robust behavior for higher-order extrapolation functions. Moreover, the present approach is independent of the correlation treatment used in quantum many-electron computations and therefore we are able to apply Coupled Clusters (CCSD-T) level of the correlation model. We demonstrate these properties on determination of the resonance position and width of the (2)Πu temporary negative ion state of diacetylene using CCSD-T level of theory.

Numerical solution of Burgers’ equation with high order splitting methods

In this work, high order splitting methods have been used for calculating the numerical solutions of Burgers’ equation in one space dimension with periodic, Dirichlet, Neumann and Robin boundary conditions. However, splitting methods with real coefficients of order higher than two necessarily have negative coefficients and cannot be used for time-irreversible systems, such as Burgers’ equations, due to the time-irreversibility of the Laplacian operator. Therefore, the splitting methods with complex coefficients and extrapolation methods with real and positive coefficients have been employed. If we consider the system as the perturbation of an exactly solvable problem (or one that can be easily approximated numerically), it is possible to employ highly efficient methods to approximate Burgers’ equation. The numerical results show that both the methods with complex time steps having one set of coefficients real and positive, say ai∈R+ and bi∈C+, and high order extrapolation methods derived from a lower order splitting method produce very accurate solutions of Burgers’ equation.

Theoretical and numerical analysis of nonconforming grid interface for unsteady flows

A framework is proposed for the spectral analysis of the numerical scheme applied on a nonconforming grid interface between two structured blocks for the two-dimensional advection equation, in a cell-centered finite-volume formalism. The conservative flux computation on the nonconforming grid interface is based on the sum of partial fluxes computed with the same numerical scheme used for a standard cell interface. This framework is used to analyze the effect of grid refinement or coarsening on the stability of the second-order centered scheme. New theoretical results are given and compared to numerical results. Considering the convection of a two-dimensional isentropic vortex, the refinement/coarsening are shown to be the cause of instabilities, poor accuracy and reflection of high-frequency waves. A new approach to compute partial fluxes, which is based both on a high-order extrapolation that accounts for the local metric and on a Riemann solver, is then proposed to reduce spurious modes.

Benchmarks and numerical methods for the simulation of boiling flows

Comparisons of different numerical methods suited to the simulations of phase changes are presented in the framework of interface capturing computations on structured fixed computational grids. Due to analytical solutions, we define some reference test-cases that every numerical technique devoted to phase change should succeed. Realistic physical properties imply some drastic interface jump conditions on the normal velocity or on the thermal flux. The efficiencies of Ghost Fluid and Delta Function Methods are compared to compute the normal velocity jump condition. Next, we demonstrate that high order extrapolation methods on the thermal field allow performing accurate and robust simulations for a thermally controlled bubble growth. Finally, some simulations of the growth of a rising bubble are presented, both for a spherical bubble and a deformed bubble.

Quasi-a priori truncation error estimation and higher order extrapolation for non-linear partial differential equations

In this paper, we show how to accurately estimate the local truncation error of partial differential equations in a quasi-a priori way. We approximate the spatial truncation error using the -estima...

Comprehensive Interpretation of a Three-Point Gauss Quadrature with Variable Sampling Points and Its Application to Integration for Discrete Data

This study examined the characteristics of a variable three-point Gauss quadrature using a variable set of weighting factors and corresponding optimal sampling points. The major findings were as follows. The one-point, two-point, and three-point Gauss quadratures that adopt the Legendre sampling points and the well-known Simpson’s 1/3 rule were found to be special cases of the variable three-point Gauss quadrature. In addition, the three-point Gauss quadrature may have out-of-domain sampling points beyond the domain end points. By applying the quadratically extrapolated integrals and nonlinearity index, the accuracy of the integration could be increased significantly for evenly acquired data, which is popular with modern sophisticated digital data acquisition systems, without using higher-order extrapolation polynomials.

Massive identification of asteroids in three-body resonances

An essential role in the asteroidal dynamics is played by the mean motion resonances. Two-body planet–asteroid resonances are widely known, due to the Kirkwood gaps. Besides, so-called three-body mean motion resonances exist, in which an asteroid and two planets participate. Identification of asteroids in three-body (namely, Jupiter–Saturn–asteroid) resonances was initially accomplished by Nesvorný and Morbidelli (Nesvorný D., Morbidelli, A. [1998]. Astron. J. 116, 3029–3037), who, by means of visual analysis of the time behaviour of resonant arguments, found 255 asteroids to reside in such resonances. We develop specialized algorithms and software for massive automatic identification of asteroids in the three-body, as well as two-body, resonances of arbitrary order, by means of automatic analysis of the time behaviour of resonant arguments. In the computation of orbits, all essential perturbations are taken into account. We integrate the asteroidal orbits on the time interval of 100,000yr and identify main-belt asteroids in the three-body Jupiter–Saturn–asteroid resonances up to the 6th order inclusive, and in the two-body Jupiter–asteroid resonances up to the 9th order inclusive, in the set of ∼250,000 objects from the “Asteroids – Dynamic Site” (AstDyS) database. The percentages of resonant objects, including extrapolations for higher-order resonances, are determined. In particular, the observed fraction of pure-resonant asteroids (those exhibiting resonant libration on the whole interval of integration) in the three-body resonances up to the 6th order inclusive is ≈0.9% of the whole set; and, using a higher-order extrapolation, the actual total fraction of pure-resonant asteroids in the three-body resonances of all orders is estimated as ≈1.1% of the whole set.

Efficient implementation of high order inverse Lax–Wendroff boundary treatment for conservation laws

In [20], two of the authors developed a high order accurate numerical boundary condition procedure for hyperbolic conservation laws, which allows the computation using high order finite difference schemes on Cartesian meshes to solve problems in arbitrary physical domains whose boundaries do not coincide with grid lines. This procedure is based on the so-called inverse Lax–Wendroff (ILW) procedure for inflow boundary conditions and high order extrapolation for outflow boundary conditions. However, the algebra of the ILW procedure is quite heavy for two dimensional (2D) hyperbolic systems, which makes it difficult to implement the procedure for order of accuracy higher than three. In this paper, we first discuss a simplified and improved implementation for this procedure, which uses the relatively complicated ILW procedure only for the evaluation of the first order normal derivatives. Fifth order WENO type extrapolation is used for all other derivatives, regardless of the direction of the local characteristics and the smoothness of the solution. This makes the implementation of a fifth order boundary treatment practical for 2D systems with source terms. For no-penetration boundary condition of compressible inviscid flows, a further simplification is discussed, in which the evaluation of the tangential derivatives involved in the ILW procedure is avoided. We test our simplified and improved boundary treatment for Euler equations with or without source terms representing chemical reactions in detonations. The results demonstrate the designed fifth order accuracy, stability, and good performance for problems involving complicated interactions between detonation/shock waves and solid boundaries.

Analysis of a new Kolgan-type scheme motivated by the shallow water equations

This paper presents the analysis of two different finite volume schemes for hyperbolic conservation laws: the Kolgan high-resolution scheme, and a new Kolgan-type scheme in which the high-order extrapolation of the conservative variables is used just in the upwind contribution of the numerical flux and source terms. Both schemes are compared in terms of the local truncation error, the stability conditions and the C-property. The schemes are applied to different hyperbolic conservation equations, including the one-dimensional scalar transport equation, the Burgers equation and the 2D shallow water equations, in order to compute the observed order of accuracy and to verify the C-property. When applied to the 2D shallow water equations, the new approach avoids spurious oscillations in the solution without the need of using high-order corrections in the definition of the bed slope source term.