Radial basis function generated finite difference (RBF-FD) method for solving PDEs is a generalization of the classical finite difference method to scattered nodes. In this work we use an RBF-FD method to solve a nonlinear conservation law in the conservative form:
is a solution vector offunctions, defined over an open and bounded domain , where is a spatial position and is a position in time, where is the final simulation time. We have that:
where is a flux term, and is denoted as the Jacobian matrix when or the velocity field when . In this paper we focus on both, scalar conservation laws () and systems of conservation laws ().
Since is nonlinear in , the analytic solutions of (1) with appropriate boundary conditions are discontinuous as , even if the initial condition is smooth. The points of discontinuities are referred to as shocks. The occurence of shocks is problematic when they are numerically approximated, as it is known that the approximations in this case contain oscillations, see Figure 1. In order to avoid the oscillations which gradually evolve into instabilities while advancing the numerical solution in time, it is important to filter these oscillations out.
1.1 Radial basis function generated finite difference methods
In this paper we define and use an oversampled RBF-FD method for time-dependent PDEs, where an explicit time stepping method is used to advance the solution in time. The discretization is constructed using a sequence of local interpolation problems on stencils, such that the differential operators become rectangular matrices. The matrices then form a rectangular system of equations, which is reformulated to a square system of equations by using an appropriate projection technique. A similar formulation for time-dependent problems was in the context of a global RBF method made inRodrigoPlatte16. The more standard, collocation RBF-FD methods Tolstykh02; FBW16; Lehto17; Talat18; Berljavac21; Shankar17 are different to the oversampled RBF-FD methods in the sense that the differential operators are discretized such that the resulting system of equations is square by construction. The motivation to formulate an oversampled RBF-FD method for time-dependent problems stems from a study of an oversampled RBF-FD method with a least-squares projection for an elliptic model problem ToLaHe21, where we found benefits compared with the collocation RBF-FD method in terms of stability and approximation error.
A standard way of stabilizing hyperbolic time-dependent problems within the RBF community is to augment the collocation scheme using a hyperviscosity term that damps high-frequency oscillations Shankar_hypervi1; Shankar_timedep; FornbergLehto; Flyer12, and makes the scheme stable in time. However, the hyperviscosity stabilization is ineffective in the presence of shocks, as is shown later in this paper, for example in Section 6, and then an additional type of stabilization has to be included to treat the shocks.
1.2 Stabilization of RBF-FD methods using a residual-based artificial viscosity term
A famous approach to control the Gibbs phenomenon and also other types of spurious oscillations, was introduced in the 1950s by von Neumann and Richtmyer NeumannRichtmyer. The authors expanded a numerical discretization using a vanishing upwind viscosity term (also called a first-order viscosity term), which diffuses the oscillations, and give a so-called viscous solution. The solution is diffused throughout the whole computational domain, and not only in the regions with oscillations (shock formations). This can over diffuse the solution in general, including many interesting waves such as a contact discontinuity and a rarefaction, which would remain unresolved.
To avoid an excessive smearing of the solution, we localize the action of the viscosity term to regions with shocks. We formulate a residual based viscosity stabilization (RV) for the collocation and the oversampled RBF-FD methods. RV is mainly used within the finite element community. Its parent is the streamline diffusion (SD) method supplemented with an artificial viscosity term to form the SD+RV method SDandRV, where the viscosity coefficient is determined from the residual of the PDE and has a shock capturing feature. The SD+RV method is difficult to implement and requires a coupled space-time discretization. To make a simplification, the SD terms were removed by Nazarov in Nazarov13, who only kept the RV term, and observed sufficient localized stabilization effects. In the same paper Nazarov13, Nazarov proved convergence of the RV stabilized finite element discretization that employs an implicit time stepping algorithm to advance the solution in time. The RV stabilization is related to the well known entropy viscosity method of Guermond et al. Guermond11, since the solution of a scalar conservation law can be one of the entropy functionals, however, this is not the case for the system of conservation laws. The RV method has also been used for stabilizing systems of conservation laws and discretizations other than the classical finite element methods in NazarovHoffman13; stiernstrom2021; MarrasNazarov15; Lu_spectral_rv,
The objective of the RV stabilization is to detect the position where the shock starts to form, by means of the magnitude of the numerical residual in (1). In the regions of where the numerical solution starts to produce an overshoot over a shock, the residual is expected to be large, whereas in the regions where the solution is free of oscillations, the residual is expected to be small. The residual is then used as a weight for first-order viscosity term, which is included in the discretization. As a consequence, the weighted viscosity term is large only in the regions with shocks and potentially leads to a better approximation error.
1.3 Other stabilization techniques
A family of essentially non-oscillatory (ENO) methods introduced in ENO are based on grid discretization of a PDE and stabilize the Gibbs phenomenon by using a reconstruction step based on interpolation with the least oscillatory stencil. Weighted ENO (WENO) methods introduced in WENO use a convex combination of the candidate interpolation stencils to reconstruct a non-oscillatory function. There exist many subsequent papers on (W)ENO methods, see e.g., Monkeberg3; Monkeberg1; Monkeberg2, where the authors develop (W)ENO methods using the RBF interpolants and by that allow approximations on unstructured grids.
Another class of stabilization methods are flux/slope limiters: these methods limit the derivatives of the solution to physically relevant magnitudes, see e.g., limiter_minmod; limiter_vanLeer; limiter_Kurganov. Employing limiters on unstructured grids is a non-trivial task.
1.4 Outline of the paper
In Section 2 we formulate an oversampled RBF-FD method for solving a scalar conservation law, and the Euler system of equations that models compressed gas dynamics. In Section 3 we derive a residual viscosity stabilization in the context of the RBF-FD method, for scalar conservation laws and Euler’s system of PDEs. In Section 4 we make a comparison between the oversampled and the collocation RBF-FD methods when solving a linear conservation law (linear advection problem), and also numerically verify the consistency of the RV stabilization. In Section 5 and Section 6 we present numerical results for solving Burger’s and the Kurganov-Petrova-Popov equations respectively. In Section 7 we solve several benchmark problems for the Euler equations. Final remarks are given in Section 8.
2 The oversampled RBF-FD discretization
In this section we describe the oversampled RBF-FD discretization that we use to solve (1). In particular, we describe our choice of point sets spread over , the formation of evaluation and differentiation matrices using the RBF-FD method, and the discretization of a time-dependent PDE where an explicit method is used to advance the solution in time.
2.1 Point sets
The domain is discretized using two point sets:
the interpolation point set for generating the cardinal functions (left plot in Figure 2),
In the present work, we generate the point set such that the mean distance between the points is set to a given value . Depending on the application, the algorithms which we use are the 2D node generator from FBF15_nodes, Gmsh Gmsh or DistMesh Distmesh. The point set is generated by placing points in each Voronoi region centered around every , , so that the relation between the number of and points is , where we define as the oversampling parameter.
2.2 Evaluation and differentiation matrices
We use the RBF-FD method in order to generate the RBF-FD trial space, which is spanned by a set of compactly supported cardinal basis functions . These functions are constructed through a sequence of interpolation problems over stencils. An example of a stencil over is given in Figure 2.
The cardinal basis functions are used to compute a PDE solution and a derivative , in the form:
where is a linear differential operator and are the solution nodal values which are the unknowns.
Using , we reformulate the relations in (3) using a semi-discrete matrix notation:
Here is the vector of nodal values and , are semi-discrete matrices constructed upon the RBF-FD evaluation/differentiation weights, where each set of weights is generated specifically for a point . An algorithmic discussion on how to generate the evaluation and differentiation weights is given in tominecDiaphragm2D. A mathematical formulation for computing these weights is given in tominec2020unfitted; ToLaHe21. In this paper we use these approaches to generate a vector of weights for every , such that the weights are exact for a cubic polyharmonic spline basis and a monomial basis of degree , which is a concept introduced in Barnett15 and studied in FFBB16; BFFB17; Bayona19; Jancic21. In our computations we use the stencil size for , and for . To form evaluation and differentiation matrices, we sample (4) in each , . This gives the following two relations:
where is a rectangular evaluation matrix of size and is a rectangular differentiation matrix of size . The two matrices can easily be computed in MATLAB by using the code provided in Tominec_rbffdcode2021.
2.3 Discretization of a conservation law
where are the -th unknown nodal values of the functions that we are solving the conservation law for, and where we used an interpolant of the flux , , to compute the divergence term in (1). Now we sample the equation above in every in order to obtain a system of equations with unknowns:
The matrix-vector formulations for the two problems that we consider in this paper are collected in the subsections below.
2.4 Scalar conservation law in matrix-vector format
Discussion in this section is kept general in the sense that it applies to for scalar conservation laws. In this case we have that and . Below, we rewrite the system (6) using the matrices defined in (5) and (4):
Here , and are the rectangular evaluation and differentiation matrices of size defined in Section 2.2. The vectors and of size correspond to the components of the (non)linear flux, thus their elements are and respectively. Now we multiply (7) with , the approximate mean distance between the points to arrive at:
where and . The multiplication by is a norm scaling introduced in ToLaHe21, which makes it possible to couple our discrete problem to an equivalent continuous problem ToLaHe21.
The rectangular system of equations (8) can not be solved using an explicit time-stepping technique in the form it has, as the system is rectangular and the time derivative on the left hand side is not explicitly defined. For this reason we project the residual of (8) onto the column space of by multiplication with , where for simplicity, we also drop the notation and obtain:
Looking at (9), the components of the left hand side matrix are:
In other words, every component of is a discrete inner product , implying that is a mass matrix in this inner product. The same type of inner product also applies to the components of the right-hand-side matrix products. Thus, our projection using has a tight connection with a classical Galerkin projection. A key difference is that our type of discretization uses discrete (summation) instead of a continuous (exact integration) inner product.
Now we solve (9) for the time derivative by inverting the mass matrix , in order to get a system of ODEs ready to be advanced in time:
Instead of directly inverting , we solve the system for iteratively, using the conjugate gradient method (function pcg() in MATLAB).
To (11) we also add two stabilization terms; (i) the term to stabilize the system of ODEs in time using a hyperviscosity operator, (ii) the term to treat the discontinuities which appear when using a nonlinear flux . The scheme is then:
This is a discretization of the scalar conservation law in conservative form. In this work we formulate the hyperviscosity term as:
where is the norm scaling and is the hyperviscosity scaling. Our choice of is justified by making the term consistent with the PDE up to the -th order, and then increasing that number to in order to allow a larger time step when using an explicit time-stepping algorithm. Term is constructed using the residual viscosity concept described in Section 3.
Note that when the flux is linear, then the conservative and the non-conservative forms of a conservation laws exhibit the same numerical properties, and the ODE system (12) transforms to:
Here and correspond to the components of the velocity field . More specifically, they represent rectangular identity matrices of size with components , , .
At this point we can use an explicit time-stepping algorithm to advance the solution of (12) from to . Throughout this paper we use an explicit classical Runge-Kutta 4 time-stepping algorithm. In the numerical experiment section of the paper we solve problems with Dirichlet-type boundary conditions on the inflow boundary. To impose these conditions we use a so-called injection method, that is, we step the discretization (12) and then in every time step overwrite the correct subset of boundary elements of with the corresponding values of the Dirichlet boundary condition. After the nodal solution at the final time is obtained, we evaluate the solution at points as:
2.5 System of conservation laws in matrix-vector format: compressible Euler equations
Consider compressible Euler equations in two spatial dimensions. The unknowns and the flux are in this case given by:
is an identity matrix,is the density, is the momentum, is the total energy and is the pressure of the fluid. The relation between the momentum and the velocity field is . When forming a system of equations we have unknown functions, but only equations. We close the system of equations using an ideal gas condition, by defining the pressure variable as: where is the temperature of the fluid and is the adiabatic gas constant. In this case the expanded form of the nonlinear conservation law (1) is:
We discretize the above system by introducing unknown nodal values for all unknown functions . The oversampled RBF-FD discretization analogous to (12) is then:
with differentiation matrices defined in Section 2.2, and stabilizers and described in Section 2.4. The Dirichlet-type boundary conditions for a given unknown function are imposed using the injection method, just as described for the scalar conservation law in Section 2.4. After the system of ODEs is advanced to the final time , each unknown set of nodal values is evaluated at the points, analogously to (15).
3 Residual based artificial viscosity stabilization of shocks for nonlinear conservation laws
Solutions to nonlinear conservation laws (1) are discontinuous, when the flux does not include any physical viscous forces. Numerical approximation of discontinuous functions suffers from the Gibbs phenomenon. Based on Nazarov13; NazarovHoffman13; NazarovLarcher17 we formulate a residual based viscosity stabilization of the solution in the context of the RBF-FD method.
The numerical scheme is stabilized using a viscosity term , where is a spatially variable coefficient, computed such that the viscous term is prevalently active in the regions of discontinuities. The viscosity term is below discretized by using the differentiation matrices introduced in Section 2.2:
The definition of for scalar conservation laws and the Euler system, is given in the following two subsections.
3.1 Definition of the viscosity coefficient for scalar conservation laws
The discrete problem which we aim to stabilize is (12). The discussion is centered around the stabilization term defined in (20). For each node , , we define the viscosity coefficient vector included inside the definition of as:
Here and are the residual and the upwind coefficients respectively. Let be a Voronoi cell with a corresponding Voronoi center , and an evaluation point (see the right plot in Figure 2). Then the definitions of and are Nazarov13:
Here we define as the minimum pairwise distance between points in a patch centered around , where the patch consists of points closest to stiernstrom2021. The temporal solution at a point comes from the relation . The user-defined constant is of size and is independent of . As , then , and . On the other hand, as , then , and . The term is the residual normalization defined by:
where is the mean of the temporal solution . The role of the normalization is to unify the physical units of and . The residual is defined as:
where is the flux divergence of the solution at time computed using the RBF-FD differentiation matrices introduced in Section 2.2. Furthermore, is an evaluation matrix also introduced in Section 2.2, and the term is an approximation of the time derivative at using the already computed solution from previous time points . We construct in two different ways: (i) when the solution is advanced in time using a constant we use the backward-differentiation formulae (BDF), (ii) when the solution is advanced in time by a variable , , then we construct an approximation to a time derivative using the polynomial basis, in each time step. In any case, the approximation order in (i) or (ii) has to match the approximation order of the spatial discretization. At the beginning of the simulation we can only use a few points from the past temporal solutions. When a sufficient amount of temporal solutions have been computed we are allowed to use a high-order formula. We then write that:
where is a differentiation operator of order , of order two, and so on. In the case (i) we use the BDF4 formula to define , for consecutively as:
In the case (ii) we construct , by computing a set of differentiation weights through a 1D polynomial interpolant , where is the number of the time points in which we already know the solution , and is the desired order of the approximation. Then we sample in to obtain a system of equations , where and , , are the interpolation matrix and the unknown coefficients respectively. A time-derivative of is , an equivalent vector formulation is , where . Plugging the computed into we have . Here is a final vector of weights for computing a derivative in time when the time grid is non-uniform. The MATLAB code to compute is provided in A.
Note that the time derivative approximated using the approaches (i) and (ii), is used only with the purpose to compute the residual, once the solution at has already been computed. The time stepping scheme that is actually advancing the solution in time is chosen differently.
3.2 Definition of the viscosity coefficient for the Euler system
Here we formulate the viscosity coefficient for the Euler system of equations (18). The coefficient is computed using the same relation as in (21), however, the definitions of and in (21) are different. The upwind viscosity coefficient at time and a point is defined as:
where the term in parenthesis is the local wave speed computed from the eigenvalues of, and where , are the velocities in horizontal and vertical directions respectively, and is the temperature of the fluid, and the internodal distance is defined in the scope of (22). The coefficient is defined by:
Here , , and are defined in (19), and the time derivative approximation is the same as defined in (25) and (26). Each normalization in (28) follows the same definition as in (23). For example, is defined as:
where is the mean density at and is a Voronoi cell centered around .
4 Numerical study I: linear advection problem
Throughout this section we compare collocation and oversampled RBF-FD methods, and investigate how do the hyperviscosity term (13) and the residual viscosity term (20) influence the numerical solution. Consider is a linear flux with time-dependent velocity field given by:
The boundary of is defined using polar coordinates , where: for . We use the explicit classical Runge-Kutta 4 method to advance the solution in time. The constant time step is chosen as:
where we also used that , since is in (31) chosen as the rotational velocity field. The term is defined in the scope of (22). The constant CFL is the Courant-Friedrichs-Lax number, specified in each subsection separately. The relative approximation errors in - and -norm at final time are computed as:
where is the exact solution at , sampled at points.
4.1 Smooth initial condition: approximation error and stability properties
Oversampled methods are normally associated with improved stability properties. In this section we investigate whether this is also true when solving hyperbolic problems, i.e., if there is a positive effect of oversampling on the eigenvalue spectrum of the discretized advection operator.