High-order Finite Difference Method Research Articles

Hydrodynamic Coupling of Viscous and Nonviscous Numerical Wave Solutions Within the Open-Source Hydrodynamics Framework reef3d

Abstract A comprehensive understanding of the marine environment in the offshore area requires phase-resolved wave information. For far-field wave propagation, computational efficiency is crucial, as large spatial and temporal scales are involved. For the near-field extreme wave events and wave impacts, high resolution is required to resolve the flow details and turbulence. The combined use of a computationally efficient large-scale model and a high-resolution local-scale solver provides a solution that combines accuracy and efficiency. This article introduces a coupling strategy between the efficient fully nonlinear potential flow (FNPF) solver REEF3D::FNPF and the high-fidelity computational fluid dynamics (CFD) model REEF3D::CFD within the open-source hydrodynamics framework REEF3D. REEF3D::FNPF solves the Laplace equation together with the boundary conditions on a sigma-coordinate. The free surface boundary conditions are discretized using high-order finite difference methods. The Laplace equation for the velocity potential is solved with a conjugated gradient solver preconditioned with a geometric multigrid provided by the open-source library Hypre. The model is fully parallelized following the domain decomposition strategy and the message passing interface protocol. The waves calculated with the FNPF solver are used as wave generation boundary conditions for the CFD-based numerical wave tank REEF3D::CFD. The CFD model employs an interface capturing two-phase flow approach that can resolve complex wave structure interaction, including breaking wave kinematics and turbulent effects. The presented hydrodynamic coupling strategy is tested for various wave conditions and the accuracy is fully assessed.

On turbulent flow and aerodynamic noise of generic side-view mirror with cell-centred finite difference method

This paper investigates the turbulent flow and aerodynamic noise of a half-cylinder body mounted on a flat plate at using high-order cell-centred finite difference method with delayed detached-eddy simulation (DDES) and large-eddy simulation (LES). Transient flow patterns from the two simulations are found to be very different in consideration of the small-scale structures. The profiles of mean velocity, resolved turbulent kinetic energy and resolved Reynolds shear stress are found to be similar among all the simulations, indicating mean quantities are relatively insensitive to turbulence modelling and grid resolution. The power spectra density of the pressure fluctuations show that LES is more capable of resolving energies in high-frequency range than DDES. After computing the normalised wavenumber-frequency spectra of fluctuating pressure on the window, we further carried out the wavenumber-frequency decomposition to separate the acoustic and the hydrodynamic components from the pressure fluctuations. The energy distribution shows that the acoustic energy has a much slower decaying rate in the high-frequency range than the hydrodynamic energy. In addition, the space-averaged sound pressure levels of pressure fluctuations on the window indicate that the present simulation with a high-order method is able to improve the accuracy in predicting pressure spectra. Finally, we carry out proper orthogonal decomposition to extract the dominating features of the decomposed acoustic and hydrodynamic components of pressure fluctuation. Patterns of multi-scale turbulence in hydrodynamic modes and propagating wavefronts of cylinder shape in acoustic modes are identified. The present research indicates that a relatively coarse grid is still capable of resolving fluctuating quantities of energy-containing structures, and LES is suggested against DDES when near-wall aerodynamic noise is the main concern.

An application of Gaussian process modeling for high-order accurate adaptive mesh refinement prolongation

We present a new polynomial-free prolongation scheme for Adaptive Mesh Refinement (AMR) simulations of compressible and incompressible computational fluid dynamics. The new method is constructed using a multi-dimensional kernel-based Gaussian Process (GP) prolongation model. The formulation for this scheme was inspired by the GP methods introduced by A. Reyes et al. (A New Class of High-Order Methods for Fluid Dynamics Simulation using Gaussian Process Modeling, Journal of Scientific Computing, 76 (2017), 443-480; A variable high-order shock-capturing finite difference method with GP-WENO, Journal of Computational Physics, 381 (2019), 189-217). In this paper, we extend the previous GP interpolations and reconstructions to a new GP-based AMR prolongation method that delivers a high-order accurate prolongation of data from coarse to fine grids on AMR grid hierarchies. In compressible flow simulations special care is necessary to handle shocks and discontinuities in a stable manner. To meet this, we utilize the shock handling strategy using the GP-based smoothness indicators developed in the previous GP work by A. Reyes et al. We demonstrate the efficacy of the GP-AMR method in a series of testsuite problems using the AMReX library, in which the GP-AMR method has been implemented.

Accelerating CFD simulation with high order finite difference method on curvilinear coordinates for modern GPU clusters

A high fidelity flow simulation for complex geometries for high Reynolds number (Re) flow is still very challenging, requiring a more powerful HPC system. However, the development of HPC with traditional CPU architecture suffers bottlenecks due to its high power consumption and technical difficulties. Heterogeneous architecture computation is raised to be a promising solution to the challenges of HPC development. GPU accelerating technology has been utilized in low order scheme CFD solvers on the structured grid and high order scheme solvers on unstructured meshes. The high-order finite difference methods on structured grids possess many advantages, e.g., high efficiency, robustness, and low storage. However, the strong dependence among points for a high-order finite difference scheme still limits its application on the GPU platform. In the present work, we propose a set of hardware-aware technology to optimize data transfer efficiency between CPU and GPU, as well as communication efficiency among GPUs. An in-house multi-block structured CFD solver with high order finite difference methods on curvilinear coordinates is ported onto the GPU platform and obtains satisfying performance with a speedup maximum of around 2000x over a single CPU core. This work provides an efficient solution to apply GPU computing in CFD simulation with specific high order finite difference methods on current GPU heterogeneous computers. The test shows that significant accelerating effects can be achieved for different GPUs.

A high-order finite-difference scheme for time-domain modeling of time-varying seismoelectric waves

Simulation of the seismoelectric effect serves as a useful tool to capture the observed seismoelectric conversion phenomenon in porous media, thus offering promising potential in underground exploration activities to detect pore fluids such as water, oil, and gas. The static electromagnetic (EM) approximation is among the most widely used methods for numerical simulation of the seismoelectric responses. However, the static approximation ignores the accompanying electric field generated by the shear wave, resulting in considerable errors when compared to analytical results, particularly under high-salinity conditions. To mitigate this problem, we have adopted a spatial high-order finite-difference time-domain method based on Maxwell’s full equations of time-varying EM fields to simulate the seismoelectric response in 2D mode. To improve the computational efficiency influenced by the velocity differences between seismic and EM waves, different time steps are set according to the stability conditions, and the seismic feedback values of EM time nodes are obtained by linear approximation within the seismic unit time step. To improve the simulation accuracy of the seismoelectric response with the time-varying EM calculation method, finite-difference coefficients are obtained by solving the spatial high-order difference approximation based on the Taylor expansion. Our method yields consistent simulation results compared to those obtained from the analytical method under different salinity conditions, thus indicating its validity for simulating seismoelectric responses in porous media. We further apply our method to layered and anomalous body models and extend our algorithm to three dimensions. Results indicate that the time-varying EM calculation method can effectively capture the reflection and transmission phenomena of the seismic and EM wavefields at the interfaces of contrasting media. This may allow for the identification of abnormal locations, thus highlighting the capability of seismoelectric response simulation to detect subsurface properties.

A high-order numerical scheme for right Caputo fractional differential equations with uniform accuracy

<abstract><p>In this paper, we mainly study the high-order numerical scheme of right Caputo time fractional differential equations with uniform accuracy. Firstly, we construct the high-order finite difference method for the right Caputo fractional ordinary differential equations (FODEs) based on piecewise quadratic interpolation. The local truncation error of right Caputo FODEs is given, and the stability analysis of the right Caputo FODEs is proved in detail. Secondly, the time fractional partial differential equations (FPDEs) with right Caputo fractional derivative is studied by coupling the time-dependent high-order finite difference method and the spatial central second-order difference scheme. Finally, three numerical examples are used to verify that the convergence order of high-order numerical scheme is $ 3-\lambda $ in time with uniform accuracy.</p></abstract>

Solution of a One-Dimension Heat Equation Using Higher-Order Finite Difference Methods and Their Stability

One-dimensional heat equation was solved for different higher-order finite difference schemes, namely, forward time and fourth-order centered space explicit method, backward time and fourth-order centered space implicit method, and fourth-order implicit Crank-Nicolson finite difference method. Higher-order schemes have complexity in computing values at the neighboring points to the boundaries. It is required there a specification of the values of field variables at some points exterior to the domain. The complexity was incorporated using Hicks approximation. The convergence and stability analysis was also computed for those higher-order finite difference explicit and implicit methods in case of solving a one dimensional heat equation. The obtained numerical results were compared with exact solutions. It is found that backward time and fourth-order centered space implicit scheme along with Hicks approximation performed well over the other mentioned higher-order approaches.

Energy stable and accurate coupling of finite element methods and finite difference methods

We introduce a hybrid method to couple continuous Galerkin finite element methods and high-order finite difference methods in a nonconforming multiblock fashion. The aim is to optimize computational efficiency when complex geometries are present. The proposed coupling technique requires minimal changes in the existing schemes while maintaining strict stability, accuracy, and energy conservation. Results are demonstrated on linear and nonlinear scalar conservation laws in two spatial dimensions.

A FFT accelerated high order finite difference method for elliptic boundary value problems over irregular domains

For elliptic boundary value problems (BVPs) involving irregular domains and Robin boundary condition, no numerical method is known to deliver a fourth order convergence and O(Nlog⁡N) efficiency, where N stands for the total degree-of-freedom of the system. Based on the matched interface and boundary (MIB) and fast Fourier transform (FFT) schemes, a new finite difference method is introduced for such problems, which involves two main components. First, a ray-casting MIB scheme is proposed to handle different types of boundary conditions, including Dirichlet, Neumann, Robin, and their mix combinations. By enclosing the concerned irregular domain by a large enough cubic domain, the ray-casting MIB scheme generates necessary fictitious values outside the irregular domain by imposing boundary conditions along the normal direction of the boundary, so that a high order central difference discretization of the Laplacian can be formed. Second, an augmented MIB formulation is built, in which Cartesian derivative jumps are reconstructed on the boundary as auxiliary variables. By treating such variables as unknowns, the discrete Laplacian can be efficiently inverted by the FFT algorithm, in the Schur complement solution of the augmented system. The accuracy and efficiency of the proposed augmented MIB method are numerically examined by considering various elliptic BVPs in two and three dimensions. Numerical results indicate that the new algorithm not only achieves a fourth order of accuracy in treating irregular domains and complex boundary conditions, but also maintains the FFT efficiency.

Axial vibration analysis of nanorods with variable density based on nonlocal elastic theory and high-order finite difference method

This paper proposes a novel high-order finite difference method (p-FDM) for predicting the axial vibration behaviour of elastic nanorods with variable density based on nonlocal elastic theory. For p-FDM, it first solves the nodal weighting coefficients for approximating the first-order derivatives of displacements at the discrete nodes by using moving least-squares approximation. Next, p-FDM recursively determines the nodal weighting coefficients for approximating the higher-order derivatives by using the obtained weighting coefficients for the first-order derivative. After verifying with exact solutions, results show that the proposed p-FDM yields accurate vibration solutions for a nanorod with either fixed–fixed ends or fixed–free ends. Only a small number of regularly or irregularly distributed grid nodes are required by a high order p-FDM to obtain a converged solution. p-FDM realises an exponential rate of convergence for frequency solution when setting p=n−2 where p is order of p-FDM and n is the number of grid nodes. New axial vibration solutions of nanorods with exponential density distribution are found by using p-FDM. It is shown that the axial frequencies are always larger for a nanorod with an exponential density distribution and no small-scale effect.

An improved CIP-based numerical model for simulating free-surface flow with adaptive mesh

In this paper, an improved CIP (constrained interpolation profile)-based numerical model, which combines with Afivo (Adaptive Finite Volume quadtree/Octree mesh, an adaptive mesh framework), has been presented to simulate free-surface flow. In the CIP-based model, the constrained interpolation profile (CIP) method, a high-order finite difference method, is employed as the base flow solver and THINC/SW (Tangent of Hyperbola for Interface Capturing with Slope Weighting) method, a volume of fluid method, is applied to capture the flow interface. In this present model, an adaptive mesh technique is introduced to reduce the computational cost and multigrid technique is used to improve the solution efficiency. In the adaptive mesh system, new spatial interpolation schemes are provided, which have better performance in mass conservation compared with schemes provided in Afivo. Lid-driven cavity flow test was carried out to validate the accuracy of the flow solver with adaptive mesh and two/three-dimensional interface benchmark tests were performed to test the reliability of the interface capturing scheme and spatial interpolation schemes. Then two/three-dimensional applications (dam-break flow case, tank sloshing case and dam-break flow over a rectangular obstacle case) were investigated to demonstrate the model accuracy. The numerical results were compared with the previous experimental and numerical results, and good agreements were observed. It is demonstrated that the present numerical model has good applicability in modeling the free-surface flow.

A High-Order Finite-Difference Method on Staggered Curvilinear Grids for Seismic Wave Propagation Applications with Topography

ABSTRACTWe developed a 3D elastic wave propagation solver that supports topography using staggered curvilinear grids. Our method achieves comparable accuracy to the classical fourth-order staggered grid velocity–stress finite-difference method on a Cartesian grid. We show that the method is provably stable using summation-by-parts operators and weakly imposed boundary conditions via penalty terms. The maximum stable timestep obeys a relationship that depends on the topography-induced grid stretching along the vertical axis. The solutions from the approach are in excellent agreement with verified results for a Gaussian-shaped hill and for a complex topographic model. Compared with a Cartesian grid, the curvilinear grid adds negligible memory requirements, but requires longer simulation times due to smaller timesteps for complex topography. The code shows 94% weak scaling efficiency up to 1014 graphic processing units.

Study on the Influence of Isolated Roughness on Hypersonic Flow over Blunt Wedge

To investigate the problem of the influence of the depth of the isolated D-type roughness element on the hypersonic flow around the blunt wedge boundary, the high-order accuracy finite difference method is used to carry out the direct numerical simulation, the influence of the depth of the isolated D-type roughness element on the interaction between the free flow and the wall is analysed, and the influence of the D-type roughness element on the wall pressure, the wall friction resistance and the wall heat flow is discussed. The results show that the strength of compression wave and expansion wave increases with the increase of the depth of D-type roughness element, but the position of compression wave and expansion wave does not change; With the increase of the depth of the D-type element, the change speed of the flow parameters will be accelerated, the peak value will increase and the trough value will decrease, and the position of the peak and trough will gradually move away from the center of the element with the increase of the depth of the element. When the depth of the D-type element is greater than or equal to 0.08, a vortex will form in the D-type element, and the length of the vortex increases with the deepening of the element. The D-type element will not change the flow state of the shear layer near it.

A hillslope-scale aquifer-model to determine past agricultural legacy and future nitrate concentrations in rivers

The long-term fate of agricultural nitrate depends on rapid subsurface transfer, denitrification and storage in aquifers. Quantifying these processes remains an issue due to time varying subsurface contribution, unknown aquifer storage and heterogeneous denitrification potential. Here, we develop a parsimonious modelling approach that uses long-term discharge and river nitrate concentration time-series combined with groundwater age data determined from chlorofluorocarbons in springs and boreholes. To leverage their informational content, we use a Boussinesq-type equivalent hillslope model to capture the dynamics of aquifer flows and evolving surface and subsurface contribution to rivers. Nitrate transport was modelled with a depth-resolved high-order finite-difference method and denitrification by a first-order law. We applied the method to three heavily nitrate loaded catchments of a crystalline temperate region of France (Brittany). We found that mean water transit time ranged 10–32 years and Damköhler ratio (transit time/denitrification time) ranged 0.12–0.55, leading to limited denitrification in the aquifer (10–20%). The long-term trajectory of nitrate concentration in rivers appears determined by flows stratification in the aquifer. The results suggest that autotrophic denitrification is controlled by the accessibility of reduced minerals which occurs at the base of the aquifer where flows decrease. One interpretation is that denitrification might be an interfacial process in zones that are weathered enough to transmit flows and not too weathered to have remaining accessible reduced minerals. Consequently, denitrification would not be controlled by the total aquifer volume and related mean transit time but by the proximity of the active weathered interface with the water table. This should be confirmed by complementary studies to which the developed methodology might be further deployed.

High-order finite difference method for the Schrödinger equation on deforming domains

A high-order finite difference discretisation of the Schrödinger equation on domains that deform in time is presented. For the domain deformation, a time-dependent coordinate transformation is used. The utilisation of summation-by-parts finite difference operators, combined with boundary conditions imposed weakly by a penalty technique, leads to a provably stable scheme. The convergence properties of the scheme are verified by numerical computations. The capabilities of the numerical scheme are demonstrated by simulations of Berry phases in deforming quantum billiards.

Non-explicit large eddy simulations of turbulent channel flows from [formula omitted] up to [formula omitted]

This numerical study based on high-order finite-difference schemes presents LES-NWR (Large Eddy Simulation with near-wall resolution) of turbulent channel flows up to Reτ=5200 using non-explicit approaches for which numerical dissipation is introduced via the discretisation of viscous terms of the Navier-Stokes equations. These models are cheaper than explicit LES models as no extra terms are needed in the equations to account for the contribution of the unresolved scales. The Approximate Deconvolution Method with Relaxation Term (ADM-RT) approach is also assessed and the present LES data are compared with reference Direct Numerical Simulations (DNS) data for first and second order moments as well as for the turbulent kinetic energy budget. Even if the viscous sublayer is not resolved, the proposed non-explicit LES approaches are in excellent agreement with the reference DNS data and to a certain extent with the ADM-RT model, for a fraction of the cost of the DNS. The proposed non-explicit models, for which the possibility to locally adjust the added numerical dissipation is investigated, are straightforward to implement and come with a negligible additional computational cost while the ADM-RT model is 30% more expensive than the non-explicit models.The parameters of the models are defined before the simulations and no modifications of the parameters are needed when the Reynolds number and the mesh resolution are changed. It is shown that a modulation of the magnitude of the numerical dissipation in time and in space is not necessarily needed, at least for the mesh resolutions and Reynolds numbers considered in the present study. The main conclusion is that non-explicit models can replace advantageously explicit models when high-order finite-difference methods are used. They can generate accurate LES-NWR of turbulent channel flows over a wide range of Reynolds numbers at a fraction of a cost of DNS.

A higher-order hybrid finite difference method based on grid equidistribution for fourth-order singularly perturbed differential equations

A higher-order hybrid finite difference method based on grid equidistribution for fourth-order singularly perturbed differential equations

A GPU-Accelerated Mixed-Precision WENO Method for Extremal Black Hole and Gravitational Wave Physics Computations

We develop and use a novel mixed-precision weighted essentially non-oscillatory (WENO) method for solving the Teukolsky equation, which arises when modeling perturbations of Kerr black holes. We show that WENO methods outperform higher-order finite-difference methods, standard in the discretization of the Teukolsky equation, due to the need to add dissipation for stability purposes in the latter. In particular, as the WENO scheme uses no additional dissipation, it is well suited for scenarios requiring long-time evolution such as the study of price tails and gravitational wave emission from extreme mass ratio binaries. In the mixed-precision approach, the expensive computation of the WENO weights is performed in reduced floating-point precision that results in a significant speedup factor of \(\approx 3.3\). In addition, we use state-of-the-art Nvidia general-purpose graphics processing units and cluster parallelism to further accelerate the WENO computations. Our optimized WENO solver can be used to quickly generate accurate results of significance in the field of black hole and gravitational wave physics. We apply our solver to study the behavior of the Aretakis charge—a conserved quantity, that if detected by a gravitational wave observatory like LIGO/Virgo would prove the existence of extremal black holes.

High-Order Bound-Preserving Finite Difference Methods for Multispecies and Multireaction Detonations

High-Order Bound-Preserving Finite Difference Methods for Multispecies and Multireaction Detonations

The design of new high-order group iterative method in the solution of two-dimensional fractional cable equation

In this article, a new high-order explicit group iterative scheme is developed for the solution of the two-dimensional time fractional cable equation. The proposed scheme is derived from the Crank-Nicolson (C-N) high-order compact finite difference method, where the Caputo discretization and C-N high-order approximations are used for the time fractional and space derivative respectively. The stability and convergence of the proposed scheme are also established. Finally, two numerical examples are presented to show the accuracy and efficiency of the scheme.