Numerical investigation of compact grid-characteristic schemes for acoustic problems

In this paper, an acoustic approximation is considered, which allows to describe the propagation of longitudinal waves in geological media. For the two-dimensional case, a compact numerical scheme is described that provides an increased approximation order on a two-point spatial stencil. To construct it, we used the method of splitting along coordinate axes and the differential continuation of the original system of equations. A number of numerical tests have been carried out, confirming the increased order of approximation.


Introduction
Oil and natural gas are the basis of the Russian energy complex. In connection with the significant depletion of deposits, the task of searching a new one is important. One of the most popular geophysical method is the seismic exploration, based on the propagation and scattering of seismic waves in heterogeneous geological media. Different mathematical models can be used to describe its dynamic behavior, from the simplest ray approximation to complex models that take into account fragile and inelastic processes. In engineering software, due to the high computational complexity of the problem, the acoustic approximation is widely used. For the numerical solution of the govern system of equations, the following methods are used: the finite element method, finite-difference schemes, etc. [1]. Note that the original system of first-order partial differential equations in pressure-velocity variables is hyperbolic. It has a full set of real eigenvalues and eigenvectors. The mathematical property of such systems is the presence of specific characteristic curves. In the 1950s, a direct grid-characteristic method was proposed that makes it possible to carry out calculations along them [2]. Its main disadvantage was the thickening of the solution in the area of the thickening of characteristics. Subsequently, its modification was proposed -the construction of inverse characteristics with the subsequent interpolation of values at the current time layer [3]. In the scientific group at the Moscow Institute of Physics and Technology the research software was developed that makes it possible to carry out calculations by the grid-characteristic method on triangular, tetrahedral, parallelepiped, curvilinear structural and hexahedral grids [4]. By the researchers the cases of an acoustic, linear-elastic isotropic, anisotropic medium and two-continuum approximation were successfully considered [5]. To improve the accuracy of simulations, schemes on extended stencils can be used, however, they require modification of the algorithm near to the boundaries. An alternative approach is the usage of differential consequences of the original system [6,7]. In work [8] supersonic flows past simple aerodynamics forms were simulated. Lubricant transfer from media to head during heat-assisted magnetic recording writing was predicted with compact schemes [9]. For elliptic equations this approach was applied in work [10]. The compact scheme was successfully applied for incompressible Navier-Stokes equations in two-dimensional case [11]. The same approach was applied for third order boundary value problem [12]. In [13], the third order of accuracy was achieved on a compact stencil for the multidimensional transport equation. Further, in [14], a scheme was constructed for one-dimensional acoustic system. This work is an extension of the ideas formulated in that paper. The dynamic problem of a homogeneous acoustic medium is considered. A compact grid-characteristic scheme on rectangular computational grids is investigated. A numerical estimate of the approximation order is carefully carried out, a comparison of the simulation results with those obtained by the Rusanov scheme [15] is done.

Mathematical model and numerical method
Acoustic equations in 2D can be written in generalized form as follows: where In order to simplify the system and use one-dimensional method of characteristics, we will apply coordinate splitting: on each iterative time step: (i) Find q n+1/2 by solving q t + A 1 q x = 0; (ii) Find˜ q n+1 by solving q t + A 2 q y = 0 with initial values after step 1; (iii) Find q n+1 by solving q t = f with initial values after step 2.
Let's consider the following one-dimensional equation: Using the grid-characteristic method, we can apply spectral decomposition to the constant matrix A 1 : where the matrices used are defined as follows: We can introduce the Riemann invariants If we multiply equation (5) by the matrix Ω x and use the substitution (10), we obtain the following diagonalized system of independent equations: This equation can be solved by utilizing the following characteristic property: the right side of which is determined via interpolation of the invariant function ω at the known time step t in the coordinate point x − λdt. Now let us design a scheme to improve the interpolation order.
Assuming the unknown vector-function q to be sufficiently smooth, we can differentiate (5) with respect to the spatial coordinate x and change the order of differentiation, obtaining We can consider q x to be another unknown and start solving the system (13) the same way as we did with the initial system (5), by introducing the Riemann invariants w defined by The problems of interpolating ω(t, x − λdt) and w(t, x − λdt) can be tackled simultaneously. We can, for example, use two stencil points for each function (corresponding to mesh points x m and x m−1 for positive λ and x m and x m+1 for negative λ), thus determining a third order interpolation polynomial.
After the interpolation procedure we know ω and w values for the next time step t + dt, and we can use the inverse transformations defined below to calculate the values of the original unknowns: In order to be able to accurately perform the similar operations at the second step of the splitting procedure, we need to know the values of q y and q xy as well. So, assuming smoothness, we can differentiate both (5) and (13) with respect to y, thus obtaining the system We can consider q y and q xy to be just some unknown functions satisfying the given equations, and, since the second function is the x-derivative of the first one, we can use the same procedure to find q n+1/2 y and q n+1/2 xy by applying transformations (10) and (14), then interpolating, and then applying inverse transformation (15) Similarly, using the grid-characteristic method, we utilize spectral decomposition of the constant matrix A 2 : where the matrices are defined as follows: and matrix Λ is the same as in the definition (9). We introduce the Riemann invariants Here the same symbols ω and w are used in order to demonstrate the commonality of the concept for a one-dimensional equation (both for the X and Y splitting steps), although the functions itself are different, meaning that Ω x q = Ω y q.
Multiplying (19) by Ω y and using the substitution (23), we obtain the following diagonalized system of independent equations: which is solved using ω i (t + dt, y) = ω i (t, y − λ i dt).
Assuming the smoothness of q, we can differentiate (19) with respect to y, which yields ( q y ) t + A 2 ( q y ) y = 0.
This equation can also be solved with respect to unknown function r = q y using the gridcharacteristic method: Riemann invariants w will be defined as the inverse transformation will be given by and the interpolation of values w i (t, y − λ i dt) will be performed together with the interpolation of ω i (t, y − λ i dt). In order for the described procedure to work correctly, the values ω and w, calculated from q and q y , must be taken from the (n + 1/2)-th step (i.e. after the X-step). That is the reason we had to solve equation (17) during the X-step. And to solve it accurately (with higher order of approximation) we used its differential continuation (18).
Now the corresponding procedure needs to be done during the Y -step. Differentiating (19) with respect to x and xy, we obtain ( q xy ) t + A 2 ( q xy ) y = 0.
These two equations are solved using the grid-characteristic method described above with the interpolation procedure performed simultaneously for both of them.
To sum up, here is the general scheme of each time step:

Results of simulations
One of the benefits of the grid-characteristic method for one-dimensional hyperbolic equations is the preserve of the interpolation order as an approximation order of the final scheme. For example, the described procedure (differential continuation) leads to the third accuracy order. However, in two-dimensional case the coordinate procedure may influence on it. In this work we numerically estimated the scheme accuracy. A square area with dimensions of 17 by 17 m was considered. The P-wave velocity was 10 m/s, and the density was 20 kg/m 3 . As an initial perturbation, an analytically specified plane wave of unit amplitude was used. It had a spatial distribution equals to the sin 4 (ξ), which ensures the necessary smoothness of the solution and its derivatives. It is known that the analytical solution of this problem is the propagation of a P-wave with a constant velocity perpendicular to a given wavefront.
Firstly, the P-wave was initiated along X axis. Secondly, the P-wave was initiated along Y axis. In these cases, the split method reduces to the pure 1D problem. Due to this fact, the third order of approximation was obtained for all used schemes with L 1 and L inf norms.
As a next test, we set the two P-waves tilted on the 45 degrees and minus 45 degrees. This symmetry with the usage of periodic boundary condition on all boundaries eliminates all numerical artifacts (see figure 1). The full period was simulated, and results were summarized at Table 1 and Table 2. To exclude the influence of splitting directions on the results, the P-wave was initiated tilted on 30 degrees to the vertical axis (see figure 2). Unfortunately, in this case it is not possible to set periodic boundary conditions with zero boundary artifacts (see figure 3). To filter them the computational domain was increased so that numerical artifacts arising near the boundaries could become to the estimation region. Table 3 and Table 4 present the results of estimating the order of convergence of algorithms using the Rusanov scheme (extended stencil) and a compact two-point scheme. In general, the other order of coordinate steps can be used. To estimate its influence on the scheme approximation, we tested both A y A x A y A x and A x A y A y A x (see Table 5 and Table 6) algorithms. It can be seen that both of them demonstrate the achievement of the second order of approximation, which is lower than in one-dimensional case. The explanation of this fact is the usage of the splitting method.

Conclusions
In this work, an acoustic approximation is considered, which allows one to describe the propagation of P-waves in geological media. For the two-dimensional case, a compact scheme is described that provides an increased order of approximation on an narrow spatial stencil. To construct it, the method of splitting along coordinate directions and differential continuation of the original system of equations were used. A number of numerical experiments have been carried out to confirm the high order of approximation. This work was carried out with the financial support of the Russian Science Foundation, project no. 19-71-10060.