Conservative Finite Difference Research Articles

Numerical simulation of a fractional Davey–Stewartson model via a conservative scheme: preservation of mass, energy and momenta

In the present manuscript, we depart from a two-dimensional form of the nonlinear Davey–Stewartson system from fluid dynamics, consisting of two partial differential equations with two unknown functions, one of them is complex-valued and the other real-valued. This model possesses four conserved quantities, namely, the mass, the energy and two momenta. In a first stage, the mathematical model is generalized to the fractional scenario considering Riesz fractional operators and three different orders. We prove mathematically that this system also possesses four quantities that are preserved in time, and which are fractional extensions of the classical mass, energy and momentum operators. Motivated by this fact, we introduce a numerical scheme for the mathematical model by making use of fractional-order centred differences. By using similar arguments as in the analysis of the continuous system, we propose some discretizations of the mass, the energy and the momenta, and we proved theoretically that these quantities are also conserved in the discrete case. In our theoretical analysis, we establish mathematically the properties of consistency for the mathematical model as well as for the discrete conserved quantities. Computationally, the finite-difference scheme is implemented using a fixed-point algorithm. Various computer simulations are performed to illustrate the conservation of the discrete quantities. Finally, we provide some interesting perspectives of research at the end of this work. This report is the first work in which a conservative finite-difference scheme for the Davey–Stewartson system is designed and rigorously analysed for the conservation properties, for both the integer-order and fractional cases.

Soliton solutions of derivative nonlinear Schrödinger equations: Conservative schemes and numerical simulation

In this paper, we numerically study soliton solutions of derivative nonlinear Schrödinger equations based on several conservative finite difference methods. All schemes own second-order accuracy with the convergence order O(τ2+h2) in the discrete L∞-norm, where h denotes the spatial step size and τ denotes the temporal step size. We show that difference schemes preserve some discrete counterparts of continuous conservation laws, and all these schemes are solvable. Extensive numerical examples with soliton solutions are carried out to verify the theoretical results. These results manifest that our schemes have potential application to soliton propagation in optical fibers.

High order conservative finite difference WENO scheme for three-temperature radiation hydrodynamics

The three-temperature (3-T) radiation hydrodynamics (RH) equations play an important role in the high-energy-density-physics fields, such as astrophysics and inertial confinement fusion (ICF). It describes the interaction between radiation and high-energy-density plasmas including electron and ion in the assumption that radiation, electron and ion are in their own equilibrium state, which means they can be characterized by their own temperatures. The 3-T RH system consists of the density, momentum and three internal energy (electron, ion and radiation) equations. In this paper, we propose a high order conservative finite difference weighted essentially non-oscillatory (WENO) scheme solving one-dimensional (1D) and two-dimensional (2D) 3-T RH equations respectively. Following our previous paper [7], we introduce the three new energy variables, and then design a finite difference scheme with both the conservative property and arbitrary high order accuracy. Based on the WENO interpolation and the strong stability preserving (SSP) high order time discretizations, taken as an example, we design a class of fifth order conservative finite difference schemes in space and third order in time. Compared with the Lagrangian method we proposed in [7], which can only reach second order accuracy for 2D 3-T RH equations if straight-line edged meshes are used, the finite difference scheme can be easily designed to arbitrary high order accuracy for multi-dimensional 3-T RH equations. The finite difference formulation is also much less expensive in multi-dimensions than finite volume schemes used in [7]. Furthermore, our method can handle fluids with large deformation easily. Numerous 1D and 2D numerical examples are presented to verify the desired properties of the high order finite difference WENO schemes such as high order accuracy, non-oscillation, conservation and adaptation to severely distorted single-material radiation hydrodynamics problems.

The thermal structure of the Colombian lithosphere: A regional and basin-scale analysis

It is well-established that the thermal state of the lithosphere strongly influences various regional and local geological processes, including crustal deformation, hydrocarbon maturation, hydrogen generation, and geothermal phenomena. Moreover, the thermal structure exhibits high sensitivity to tectonic features, a property of particular significance in Colombia, where three main tectonic plates converge, and lithospheric tearing has been documented. In this contribution, we focus on elucidating the impact of plate architecture on the thermal field in central-eastern Colombia at both shallow and deep levels. To accomplish this, we constructed a series of two-dimensional profiles and derived a numerical solution of the heat equation using the conservative finite difference method. As constraints, we incorporate an inferred distribution of rocks in the deep crust and upper mantle based on global and local lithospheric thickness models. Material parameters for the various rocks, both exposed and inferred, were obtained from the literature. Additionally, we used superficial heat flow estimates and apparent geothermal gradients compiled by the Colombian Geological Survey.Our results suggest a significant influence of the lithospheric Caldas tear on the thermal state of Colombia, with the breaking off occurring in the Nazca plate under the Eastern Cordillera range around 5°N. The modeled asthenospheric heat flow remains approximately 25 mWm−2, except in the northern Eastern Cordillera range, where the background heat flow increases rapidly to 40 mWm−2. Consequently, our model predicts partial melting in the lower crust to the north and a thermally unstable lower crust to the south of the Caldas tear. The material parameters that best fit the surface data suggest the presence of a basement moderately enriched in radioactive elements in the Eastern Llanos basin. After accounting for compaction, we also confirm a strong influence of the tectonic setting on the thermal state of sedimentary basins.

Scalable implicit solvers with dynamic mesh adaptation for a relativistic drift-kinetic Fokker–Planck–Boltzmann model

In this work we consider a relativistic drift-kinetic model for runaway electrons along with a Fokker–Planck operator for small-angle Coulomb collisions, a radiation damping operator, and a secondary knock-on (Boltzmann) collision source. We develop a new scalable fully implicit solver utilizing finite volume and conservative finite difference schemes and dynamic mesh adaptivity. A new data management framework in the PETSc library based on the p4est library is developed to enable simulations with dynamic adaptive mesh refinement (AMR), distributed memory parallelization, and dynamic load balancing of computational work. This framework and the runaway electron solver building on the framework are able to dynamically capture both bulk Maxwellian at the low-energy region and a runaway tail at the high-energy region. To effectively capture features via the AMR algorithm, a new AMR indicator prediction strategy is proposed that is performed alongside the implicit time evolution of the solution. This strategy is complemented by the introduction of computationally cheap feature-based AMR indicators that are analyzed theoretically. Numerical results quantify the advantages of the prediction strategy in better capturing features compared with nonpredictive strategies; and we demonstrate trade-offs regarding computational costs. The robustness with respect to model parameters, algorithmic scalability, and parallel scalability are demonstrated through several benchmark problems including manufactured solutions and solutions of different physics models. We focus on demonstrating the advantages of using implicit time stepping and AMR for runaway electron simulations.

Invariant conservative finite-difference schemes for the one-dimensional shallow water magnetohydrodynamics equations in Lagrangian coordinates

Invariant finite-difference schemes for the one-dimensional shallow water equations in the presence of a magnetic field for various bottom topographies are constructed. Based on the results of the group classification recently carried out by the authors, finite-difference analogues of the conservation laws of the original differential model are obtained. Some typical problems are considered numerically, for which a comparison is made between the cases of a magnetic field presence and when it is absent (the standard shallow water model). The invariance of difference schemes in Lagrangian coordinates and the energy preservation on the obtained numerical solutions are also discussed.

Uniform Error Bounds of a Conservative Compact Finite Difference Method for the Quantum Zakharov System in the Subsonic Limit Regime

Uniform Error Bounds of a Conservative Compact Finite Difference Method for the Quantum Zakharov System in the Subsonic Limit Regime

Conservative higher-order finite difference scheme for the coupled nonlinear Schrödinger equations

This paper introduces a conservative higher-order finite difference scheme for solving the coupled nonlinear Schrödinger equations. The Crank–Nicolson method is employed to discretize time derivatives and the sixth-order difference operator is used to discretize space derivatives, correspondingly, the resulting difference scheme has second-order accuracy in time and sixth-order accuracy in space. By utilizing the discrete energy method, the conservation of discrete mass and energy, the boundedness, existence and uniqueness of solution, unconditional stability and the convergence of the new scheme are proved. Then making use of Richardson extrapolation, the time accuracy is increased to the fourth order. Finally, the numerical experiments are conducted to validate the theoretical results presented in the paper.

THE FORMATION AND THE EVOLUTION OF LARGE-SCALE VORTEX STRUCTURES IN STELLAR ACCRETION DISKS

Explaining the causes of angular momentum transfer in accretion stellar disks is an important astrophysical problem, since it is the process that determines the rate of accretion of matter onto the central gravitating body. Previously, within the framework of a two-dimensional approach, it was shown that the introduction of small perturbations into the flow of disk matter leads to the appearance of shear instability. This process is accompanied by the development of large-scale vortex structures. Their movement and evolution lead to a redistribution of angular momentum in the accretion disk. The action of the described mechanism was previously studied numerically only within a two-dimensional approximation, so the goal of the current work is to carry out full-scale three-dimensional modeling. The processes under study are described within the framework of the system of equations of ideal gas dynamics. The article briefly describes the method for their numerical integration, which is based on a conservative finite-difference scheme and the solution of the Riemann problem. The initial data is a stationary toroidal gas state surrounded by a matter with low density and pressure. At the next step, small perturbations of one of the gas-dynamic variables are introduced. The modeling and analysis of the results of numerical calculations show the emergence of vortex structures in the shear flow of a three-dimensional accretion disk. Their movement is accompanied by a redistribution of matter and angular momentum in the volume of the disk, leading to accretion of matter onto the central body.

Hermes-3: Multi-component plasma simulations with BOUT++

A new open source tool for fluid simulation of multi-component plasmas is presented, based on a flexible software design that is applicable to scientific simulations in a wide range of fields. Hermes-3 is built on plasma simulation framework BOUT++, consolidating earlier SD1D and Hermes models into a single code that can be configured at run-time to solve plasma models in 1D, 2D or 3D, either for transport (steady-state) or turbulent (time-evolving) problems, with an arbitrary number of ion and neutral species. We describe the improved numerical algorithms and software design that have been implemented in Hermes-3.To demonstrate the capabilities of this tool, applications relevant to the boundary of tokamak plasmas are presented: 1D simulations of diveror plasmas evolving equations for all charge states of neon and deuterium; 2D transport simulations of tokamak equilibria in single-null X-point geometry with plasma ion and neutral atom species; and simulations of the time-dependent propagation of plasma filaments (blobs).Hermes-3 is publicly available on Github under the GPL-3 open source license. The repository includes documentation and a suite of unit, integrated and convergence tests. Program summaryProgram Title: Hermes-3CPC Library link to program files:https://doi.org/10.17632/rtkvc3r743.1Developer's repository link:https://github.com/bendudson/hermes-3Licensing provisions: GPLv3Programming language: C++14Nature of problem: Simulation of dynamics and steady-state solutions of multiple ion and neutral species in magnetised plasmas using drift-reduced fluid plasma models in 1 to 3 spatial dimensions. Focused on but not restricted to the simulation of turbulence in the boundary of tokamaks.Solution method: Hermes-3 implements a system of flexible software components, that are configured at run-time and coupled by passing a state object between them (the Encapsulate State and Command patterns). Hermes-3 builds on BOUT++ [1], employing the Method of Lines with implicit and explicit time-integration methods in curvilinear coordinates on block-structured meshes, making use of libraries including SUNDIALS, PETSc and Hypre. Hermes-3 implements drift-reduced plasma fluid equations using conservative finite difference methods, and atomic processes that couple plasma and neutral fluids.

Hybrid conservative central/WENO finite difference scheme for two-dimensional detonation problems

In the present work, a hybrid finite difference scheme is proposed and tested for the study of two-dimensional detonation problems. The hybrid scheme consists of a fifth-order weighted essentially nonoscillatory (WENO) method for discontinuous regions and a fourth or sixth-order robust conservative central finite difference scheme for the remaining domain. An arbitrary shock sensor is used to identify the discontinuities. The proposed shock sensor is based on the absolute value of the density gradient with an arbitrary threshold computed at each time step using a fast image segmentation technique. The sensor is complemented with a Ducros sensor to avoid the selection of vortices as discontinuities. The hybrid scheme is tested for several benchmark examples, including a two-dimensional transverse detonation wave, detonation wave diffraction, and a multiple obstacles test case. The obtained results show that the proposed sensor detects discontinuities adequately and the hybrid schemes give the same results as the fifth-order WENO in faster computational time.

Conservative second-order finite difference method for Camassa–Holm equation with periodic boundary condition

In this paper, we propose two momentum-preserving finite difference schemes for solving one-dimensional Camassa–Holm equation with periodic boundary conditions. A two-level nonlinear difference scheme and a three-level linearized difference scheme are constructed by using the method of order reduction. For nonlinear scheme, we combine mid-point rule and a specific difference operator, which ensures that our obtained scheme is of second-order convergence in both temporal and spatial directions. For linearized scheme, we apply a linear implicit Crank–Nicolson scheme in the temporal direction, then unique solvability and momentum conservation are analysed in detail. Numerical experiments are provided for Camassa–Holm equation admitting different types of solutions, which demonstrate the convergence order and accuracy of the proposed methods coincide with theoretical analysis. Moreover, numerical results show that the nonlinear scheme exhibits better accuracy for mass conservation, while the linearized scheme is more time-saving in computation.

A conservative splitting high-order finite difference method for coupled Gross–Pitaevskii equations in 2D

A conservative splitting high-order finite difference method for coupled Gross–Pitaevskii equations in 2D

Analysis of two conservative fourth-order compact finite difference schemes for the Klein-Gordon-Zakharov system in the subsonic limit regime

We propose two conservative fourth-order compact finite difference (CFD4C) schemes and give the rigorous error analysis for the Klein-Gordon-Zakharov system (KGZS) with ε∈(0,1] being a small parameter. In the case 0<ε≪1, i.e., the subsonic limit regime, the solution of KGZS propagates waves with wavelength O(ε) in time and O(1) in space, respectively, with amplitude at O(εα†) with α†=min⁡{α,β+1,2}, where α and β describe the incompatibility between the initial data of the KGZS and the limiting Klein-Gordon equation (KGE) as ε→0+ and satisfy α≥0,β+1≥0. Different from the standard KGZS, the oscillation in time becomes the main difficulty in constructing numerical schemes and giving the error analysis for KGZS in this regime. We prove that the schemes conserve the discrete energies. By defining a discrete ‘energy’, applying cut-off technique and considering the convergence rate of the KGZS to the limiting KGE, we obtain two independent error estimates with the bounds at O(h4/ε1−α⁎+τ2/ε3−α†) and O(h4+τ2+εα⁎) where α⁎=min⁡{1,α,1+β}, τ is time step and h is mesh size. Hence, we get uniformly accurate error estimates with the bounds at O(h4+τ2/(4−α†)) when α≥1 and β≥0, and O(h4α⁎+τ2α⁎/3) when 0<α<1 or −1<β<0. While for α=0 or β=−1, the meshing requirements of CFD4Cs are h=O(ε1/4) and τ=O(ε3/2) for 0<ε≪1. The compact schemes provide much better spatial resolution compared with the usual second order finite difference schemes. Our numerical results support energy-conserving properties and the results of the convergence analysis.

Reformulated Dissipation for the Free-Stream Preserving of the Conservative Finite Difference Schemes on Curvilinear Grids

Reformulated Dissipation for the Free-Stream Preserving of the Conservative Finite Difference Schemes on Curvilinear Grids

Influences of conservative and non-conservative Lorentz forces on energy conservation properties for incompressible magnetohydrodynamic flows

In the analysis of magnetohydrodynamic (MHD) flow, the Lorentz force significantly affects energy properties because the work generated by the Lorentz force changes the kinetic and magnetic energies. Therefore, the Lorentz force and energy conversion should be predicted accurately. Some energy conservation schemes have been proposed and validated. However, the influences of the Lorentz force discretization on conservation and conversion of energy have not yet been clarified. In this study, a conservative finite difference method is constructed for incompressible MHD flows considering the induced magnetic field. We compare the difference in energy conservation properties among three methods of calculating the Lorentz force. The Lorentz forces are calculated in conservative and non-conservative forms, and both compact and wide-range interpolations of magnetic flux density are used to calculate the non-conservative Lorentz force. The compact interpolation method proposed in this study can perform conversions between conservative and non-conservative forms of the Lorentz force even when using the finite difference method. The present numerical method improves the conservation of transport quantity. Five models were analyzed, and the accuracy and convergence of the present numerical method were verified. From the viewpoint of the conservation of the total energy in an ideal inviscid periodic MHD flow, we consider that the calculation using compact interpolation for the Lorentz force is appropriate. This method preserves the total energy even on non-uniform grids. Moreover, the divergence-free condition of the magnetic flux density is discretely satisfied even without the correction of the magnetic flux density. The present numerical method can capture the Hartmann layer in the propagation of an Alfvén wave and accurately predict the tendency of energy attenuation in the analysis of a Taylor decaying vortex under magnetic fields. Analysis of the Orszag–Tang vortex reveals energy dissipation processes and the generation of high current densities. The present numerical method has excellent energy conservation properties and can accurately predict energy conversion. Therefore, this method can contribute to understanding complex unsteady MHD flows.

Control of nonlinear bulk deformation and large shear strain on first-order phase transformation kinetics

Phase transformations play a key role in numerous coupled natural processes, and they are important for many industrial applications. However, the kinetics of phase transformations in coupled chemo-mechanical systems undergoing large mechanical deformations still needs to be better quantified. Here, we study the phase transformation kinetics of a two-phase binary mixture using the diffuse interface approach. We couple a Cahn–Hilliard type model with a mechanical model for a compressible viscous flow. The bulk compressibility is a nonlinear function of the pressure, and the shear viscosity is a nonlinear function of the concentration. The mechanical coupling is achieved by employing a pressure-dependent mechanical mixing term in the equation for the Gibbs energy. We derive a dimensionless system of equations which we solve numerically with a pseudo-transient method using conservative finite differences for discretization. We perform numerical simulations in 1D and 2D model setups considering far-field simple shear and pure shear. For a chemo-mechanically coupled system, we show that the velocity of the phase boundary is a linear function of the degree of metastability and, hence, confirm the hypothesis of “normal growth.” A stronger mechanical coupling and a larger volumetric effect of the chemical reaction result in lower phase boundary velocities. The 2D results show a significant impact of the mechanical coupling and the far-field deformation on the orientation and kinetics of the phase transformations. Under far-field simple shear and pure shear in 2D, the phase transformations generate string-like patterns. The orientation of these patterns is controlled by the applied far-field deformation and orientations differ by 45 degrees between simple shear and pure shear.

A nonlinear heat equation arising from automated-vehicle traffic flow models

In this paper, a new nonlinear heat equation is studied that arises as a model of the collective behaviour of automated vehicles. The properties of the solutions of this equation are studied by introducing the appropriate notion of a weak solution that requires certain entropy-like conditions. To obtain an approximation of the solution of the nonlinear heat equation, a new conservative first-order finite difference scheme is proposed that respects the corresponding entropy conditions, and certain links between the weak solution and the numerical scheme are provided. Finally, a traffic simulation scenario and a comparison with the Lighthill–Whitham–Richards (LWR) model are provided, illustrating the benefits of the use of automated vehicles.

Conservative second-order accurate finite-difference scheme for the coupled Maxwell-Dirac equations

The recent development of nanoelectronic devices that incorporate Dirac materials has highly increased the need for adequate simulation and modelling tools. This paper introduces an accurate, multiphysics finite-difference time-domain method to solve the pertinent Maxwell-Dirac equations. The stability criterion for the Dirac equation with electromagnetic fields is derived, which reduces to the Courant-Friedrichs-Lewy condition in the absence of electromagnetic fields. Validation examples show the second-order accuracy of the novel fully coupled Maxwell-Dirac scheme and illustrate that total norm and energy are caonserved within a relative error of order 10−4. The method is applied to a ZrTe5 waveguide and it is found that even at low field strengths, the charge carriers can be accelerated to 80% of the Fermi velocity. Furthermore, the flexibility of the advocated method allows for the seamless integration into existing computational electromagnetics frameworks and the possible extension to higher-order schemes.

An efficient fourth-order accurate conservative scheme for Riesz space fractional Schrödinger equation with wave operator

In this study, we investigate a high-order accurate conservative finite difference scheme by utilizing a fourth-order fractional central finite difference method for the two-dimensional Riesz space-fractional nonlinear Schrödinger wave equation. The conservation laws of the discrete difference scheme are shown. Meanwhile, the exactness, uniqueness, and prior estimate of the numerical solution are rigorously established. Then, it is proved that the proposed scheme is unconditionally convergent in the discrete L2 and Hγ/2 norm, where γ is a fractional order. Furthermore, we demonstrate that when the fractional order γ and the spatial grid number J increase, the block-Toeplitz coefficient matrix generated by the spatial discretization becomes ill-conditioned. As a result, we adopt an effective linearized iteration method for the nonlinear system, allowing it to be solved efficiently by the Krylov subspace solver with an appropriate circulant preconditioner, in which the fast Fourier transform is applied to speed up the computational cost at each iterative step. Finally, numerical experiments are presented to validate the theoretical findings and the efficiency of the fast algorithm.