Journal of Computational Physics 508 (2024) 112981

Contents lists available at ScienceDirect

Journal of Computational Physics
journal homepage: www.elsevier.com/locate/jcp

A non-oscillatory ﬁnite volume scheme using a weighted
smoothed reconstruction
Davoud Mirzaei a,∗ , Navid Soodbakhsh b,ca Division of Scientiﬁc Computing, Department of Information Technology, Uppsala University, Swedenb Department of Applied Mathematics and Computer Science, Faculty of Mathematics and Statistics, University of Isfahan, 81746-73441, Isfahan,

Iranc School of Mathematics, Institute for Research in Fundamental Sciences (IPM), Tehran, 19395-5746, IranA R T I C L E

I N F O

Keywords:
Conservation laws
Finite volume method
Radial basis functions
Polyharmonic splines
Weighted smoothed reconstruction (WSR)A B S T R A C T
In this research article, we introduce a high-order and non-oscillatory ﬁnite volume method
in combination with radial basis function approximations and use it for the solution of scalarconservation laws on unstructured meshes. This novel approach departs from conventional
non-oscillatory techniques, which often require the use of multiple stencils to achieve smoothreconstructions. Instead, the new method uses a single central stencil and hinges on an
approximate interpolation methodology called the weighted smoothed reconstruction (WSR),
with a foundation on polyharmonic spline interpolation. Through some numerical experiments,
we demonstrate the eﬃciency and accuracy of the new approach. It reduces the computational
cost and performs well in capturing shocks and sharp solution fronts.

1. Introduction
Many physics and engineering problems give rise to time-dependent partial diﬀerential equations (PDEs), and one of the mostsigniﬁcant and extensively used among them is the hyperbolic conservation law. A characteristic feature of conservation laws is thatthe solution to the problem may develop discontinuities (shocks) or sharp fronts, even when the initial condition is smooth. These
discontinuities or sharp fronts result in non-physical oscillations in the numerical solution. The purpose of numerical techniques isto precisely capture these discontinuities (shocks) while avoiding unphysical oscillations and achieving a solution with a high orderof convergence.
Over recent decades, the ﬁnite volume method (FVM) has been established as a well-known and commonly used numerical
technique for solving conservation laws. FVM holds a rich history, and has been evolved into various forms and methodologies[31,34,38,51]. The main structure behind almost all ﬁnite volume methods is the spatial discretization of the domain into a set ofcells (control volumes), temporal discretization of time interval, reconstructing the solution in each control volume using the given
cell average values (average of the solution in a control volume) at a certain time, and then updating the cell average values in
the next time step by calculating the numerical ﬂux at the cell interfaces. These methods are either linear or non-linear and sharesimilarities with ﬁnite diﬀerence schemes for hyperbolic PDEs, but are also adaptable to unstructured grids.
According to the Godunov theorem a linear scheme must be either only ﬁrst-order accurate or oscillatory. Some nonlinear numerical schemes have been developed towards high-order non-oscillatory methods, including artiﬁcial viscosity schemes, total variation* Corresponding author.
E-mail address: davoud.mirzaei@it.uu.se (D. Mirzaei).

https://doi.org/10.1016/j.jcp.2024.112981
Received 6 October 2023; Received in revised form 29 February 2024; Accepted 28 March 2024
Available online 3 April 2024
0021-9991/© 2024 The Author(s).
Published by Elsevier Inc.
(http://creativecommons.org/licenses/by/4.0/).

This is an open access article under the CC BY license

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhshdiminishing (TVD) schemes, total variation bounded (TVB) schemes, essentially non-oscillatory (ENO) and weighted essentiallynon-oscillatory (WENO) schemes, monotonicity-preserving schemes, and boundary variation diminishing (BVD) schemes. There are
typically three approaches employed in these techniques: Godunov’s approach, the use of artiﬁcial viscosity, or the use of a high-orderreconstruction. Standard numerical methods of Godunov-type are well introduced in LeVeque [34] and Toro [51], Hesthaven [17],and references therein in the context of FVM. A ﬁnite diﬀerence artiﬁcial viscosity method was proposed by Von Neumann and Richtmyer [52] in the 1950s. Then it becomes the predominant approach in the ﬁnite element method (FEM) community for solving ﬂuid
dynamic problems [14,21,28,39]. As an application of the artiﬁcial viscosity and hyperviscosity methods for stabilizing RBF-basedalgorithms refer to [41,50].
The earliest attempts at achieving higher than ﬁrst-order reconstructions date back to the application of ﬂux and slope limitermethods to obtain second-order accurate schemes in FVM [34]. Subsequently, higher-order ENO reconstructions were developed
and widely used for approximating hyperbolic PDEs in computational ﬂuid dynamic problems [1,15,16,46–48]. The primary idea
behind the ENO reconstruction includes determining sets of stencils surrounding a control volume, then computing a reconstruction
on each of these stencils, and then picking the smoothest (least oscillatory) reconstruction as the desired reconstruction on thecontrol volume. The ENO reconstruction on a control volume can be improved by using a weighted sum of diﬀerent reconstructions
on diﬀerent stencils. The weights are chosen based on the smoothness of the reconstructions, such that smoother (least oscillatory)
reconstructions have higher weights. This technique is known as WENO reconstruction (see for example [10,20,27,35,42,43,55,57]).
A comprehensive review of this approach in FDM and FVM frameworks, using polynomial reconstructions, can be found in [44].
Eﬀorts to combine WENO schemes with RBF interpolations were made earlier in [2,3,23].
Another class of reconstruction techniques are central essentially non-oscillatory (CENO) methods that aim to overcome some
restriction of ENO/WENO methods by using a central reconstruction through combining a high-order 𝑘-exact and a monotonicity
preserving piecewise linear reconstructions [24].
In this article, we introduce a novel non-oscillatory ﬁnite volume method that adeptly captures shocks and sharp fronts in thesolution using just one central stencil. This places it within the CENO category, although it fundamentally diﬀers from the existing
CENO methods. Our approach employs a single central stencil alongside a generalized spline smoothing technique, called weighted
smoothed reconstruction (WSR), to suppress non-physical oscillations in the vicinity of shocks or sharp fronts. This stands in contrast
to the WENO reconstruction method, which involves the selection of several stencils (central, forward, and backward), computationof corresponding reconstructions, and forming a weighted reconstruction by assigning speciﬁc weights to diﬀerent reconstructions.
To achieve a non-oscillatory reconstruction, we extend the approximate interpolation method presented in [54] to a weighted
smoothed approximation with varying smoothing parameters. This WSR approach is based on radial basis function (RBF) approximation and is integrated into the reconstruction step of the ﬁnite volume method (FVM). When comparing it to other techniques,
we denote this method as FVM-WSR. It is applicable to unstructured meshes and can be adapted for various kernel-based interpolations. To ensure stability and accuracy, we employ polyharmonic spline basis functions and use their scalability properties
to construct stable local interpolation matrices. In the temporal domain, we use strong stability-preserving Runge-Kutta (SSPRK)
methods [6,11,12,46].
The structure of this paper is as follows: In Section 2, we provide a brief overview of kernel-based interpolation and polyharmonic
spline RBF approximation. In Section 3, we outline the fundamental principles of FVM for conservation laws and brieﬂy discuss theWENO reconstruction. In Section 4, we explore the approximate interpolation and the smoothing spline technique, and then develop
the generalized version WSR that will be employed in subsequent sections. In Section 5 we focus on the adaptation of the WSR
technique for use in the FVM framework, and ﬁnally in Section 6, we demonstrate the performance and accuracy of FVM-WSR
in solving various types of scalar conservation laws, including linear advection, Burgers’, and the Kurganov-Petrova-Popov (KKP)
equations, both with smooth and discontinuous initial conditions.
2. RBF interpolation
In this section, we provide a short overview of standard interpolation using conditionally positive deﬁnite kernels. For moredetails see, for example, [9,53].
Assume that ℝ𝑑 represents the set of all real vectors in dimension 𝑑 , and ‖ ⋅ ‖ stands for the Euclidean norm onℝ𝑑). Also, assume(𝑚+𝑑that ℙ𝑑𝑚 is the space of polynomials of degree at most 𝑚 (order 𝑚 + 1) on ℝ𝑑 . The dimension of this space is 𝑄 =.
𝑑Deﬁnition 2.1. A univariate, continuous and even function 𝜙 is conditionally positive deﬁnite on ℝ𝑑 with respect to polynomial
space ℙ𝑑𝑚 (or of order 𝑚 + 1) if
𝑁 ∑𝑁∑
𝑖=1 𝑗=1𝛼𝑖 𝛼𝑗 𝜙(‖𝑥𝑖 − 𝑥𝑗 ‖) > 0

for any 𝑁 and any set of pairwise distinct points 𝑥1 , 𝑥2 , … , 𝑥𝑁 ∈ ℝ𝑑 and any real values 𝛼1 , 𝛼2 , … , 𝛼𝑁 that satisfy the side condition𝑁∑𝑗=1

𝛼𝑗 𝑝 (𝑥𝑗 ) = 0,

for all 𝑝 ∈ ℙ𝑑𝑚 .
2

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh
Table 1
Typical RBFs for 𝑟 = ‖𝑥‖2 , 𝑥 ∈ ℝ𝑑 .
RBF

𝜙(𝑟)

Parameters

Order

Gaussian
Sobolev spline
IMQ
MQ
Polyharmonic splines(PHS)
Wendland’s functionsexp(𝜀2 𝑟2 )𝑟𝛽−𝑑∕2 𝐾𝛽−𝑑∕2 (𝜀𝑟)(1 + 𝜀2 𝑟2 )−𝛽(1 + 𝜀2 𝑟2 )𝛽{
𝑟2𝑘−𝑑 log 𝑟𝑟2𝑘−𝑑
𝜙𝑑,𝑘 (𝜀𝑟)𝜀>0
𝛽 > 𝑑∕2, 𝜀 > 0
𝛽 > 0, 𝜀 > 0
𝛽 > 0, 𝛽 ∉ ℕ, 𝜀 > 0
𝑘 ∈ ℕ, 𝑑 even𝑘 ∈ ℕ, 𝑑 odd
𝑘 ∈ ℕ0 , 𝑑 ∈ ℕ, 𝜀 > 0

0
0
0
⌈𝛽⌉
𝑘 − ⌈𝑑∕2⌉ + 1
0

Function 𝜙 is called positive deﬁnite if it is conditionally positive deﬁnite of order zero. Such a (conditionally) positive deﬁnite
function is called a radial basis function (RBF).
Some typical examples of RBF 𝜙(𝑟) for 𝑟 = ‖𝑥‖2 , 𝑥 ∈ ℝ𝑑 , are listed in Table 1. The positive parameter 𝜀 is called the shapeparameter. Polyharmonic spline (PHS) kernels are examples of 𝜀-free RBFs. In the second row of the table, 𝐾𝛽 is the modiﬁed Bessel
function of the second kind of order 𝛽 . For explicit formula of Wendland’s functions 𝜙𝑑,𝑘 (𝜀𝑟) see [53, Chapter 9].
The RBF interpolation of a continuous function 𝑓 ∶ Ω ⊂ ℝ𝑑 → ℝ on a discrete set 𝑋 = {𝑥1 , … , 𝑥𝑁 } ⊂ Ω is expressed as follows

𝑠𝑓 ,𝑋 (𝑥) =

𝑁
∑
𝑗=1

𝛼𝑗 𝜙(‖𝑥 − 𝑥𝑗 ‖) +

𝑄
∑
𝑗=1

𝛽𝑗 𝑝𝑗 (𝑥)(2.1)

where {𝑝1 , … , 𝑝𝑄 } forms a basis for ℙ𝑑𝑚 . Imposing the interpolation conditions 𝑠𝑓 ,𝑋 (𝑥𝑖 ) = 𝑓 (𝑥𝑖 ) for 𝑖 = 1, … , 𝑁 , and applying theabove side condition imply that the vectors 𝛼 = (𝛼1 , … , 𝛼𝑁 )𝑇 and 𝛽 = (𝛽1 , … , 𝛽𝑄 )𝑇 satisfy the system of equations

[

𝐴
𝑃𝑇

𝑃
0

][ ] [ ]
𝛼
𝑓
= 𝑋
𝛽
0

(2.2)

where

𝐴 = (𝜙(‖𝑥𝑗 − 𝑥𝑖 ‖)) ∈ ℝ𝑁×𝑁 ,𝑃 = (𝑝𝑗 (𝑥𝑖 )) ∈ ℝ𝑁×𝑄 ,

(2.3)
𝑇

𝑓𝑋 = (𝑓 (𝑥1 ), … , 𝑓 (𝑥𝑁 )) .
From deﬁnition, the matrix 𝐴 is symmetric and positive deﬁnite on the null space of 𝑃 𝑇 . We also assume that 𝑁 ⩾ 𝑄 and 𝑋 is ℙ𝑑𝑚 unisolvent set to guarantee the solvability of system (2.2). For a (strictly) positive deﬁnite RBF, the symmetric matrix 𝐴 is positivedeﬁnite and the interpolation problem is solvable in the absence of polynomial terms.
An important class of RBFs that will be used in this paper is polyharmonic splines (PHS)𝜙𝑑,𝑘 (𝑟) ∶={ 2𝑘−𝑑log 𝑟, 𝑘 ∈ ℕ, 𝑑 even𝑟𝑟2𝑘−𝑑 ,𝑘 ∈ ℕ, 𝑑 odd(2.4)

for 2𝑘 − 𝑑 > 0. The PHS function 𝜙𝑑,𝑘 is (up to a sign) conditionally positive deﬁnite of order 𝑚 + 1 = 𝑘 − ⌈𝑑∕2⌉ + 1. The native spaceassociated to this RBF is the Beppo-Levi space

BL𝑘 (ℝ𝑑 ) ∶= {𝑓 ∈ 𝐶(ℝ𝑑 ) ∶ 𝐷𝛼 𝑓 ∈ 𝐿2 (ℝ𝑑 ), ∀𝛼 ∈ ℕ𝑑0 with |𝛼| = 𝑘}provided that it is considered as a conditionally positive deﬁnite kernel of order 𝑘. Note that a conditionally positive deﬁnite function
with respect to a polynomial space of a certain order is also conditionally positive deﬁnite with respect to higher order polynomial
spaces. We refer the reader to [9,53] for more details about (conditionally) positive deﬁnite functions and their approximation
properties.
3. Finite volume method
In this section, we provide a brief description about ﬁnite volume method (FVM) for solving the scalar conservation law. For a
more detailed explanation, interested readers are referred to [17,34,51].
3.1. The scalar conservation law
Consider a scalar conservation law on an open and bounded computational domain Ω ⊂ ℝ𝑑 with an initial condition as follows

𝜕𝑢+ ∇ ⋅ 𝐹 (𝑢) = 0,
𝜕𝑡𝑢 (0, 𝑥) = 𝑢0 (𝑥).

(3.1)
3

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

(

)𝑇

(

]

Here, 𝐹 (𝑢) ∶= 𝑓1 (𝑢), … , 𝑓𝑑 (𝑢) is a ﬂux function and 𝑢 ≡ 𝑢 (𝑡, 𝑥) ∶ 𝐼 × Ω ⟶ ℝ is the solution of the problem. Moreover, 𝐼 ∶= 0, 𝑡𝑓
is a time interval with a ﬁnal time 𝑡𝑓 , and 𝑢0 (𝑥) is the initial function.
In orderto discretize{} the problem (3.1) using FVM, we assume that the domain Ω is covered by a non-overlapping set of volumes(cells) 𝑉1 , 𝑉2 , … , 𝑉𝑁 =∶  . The integral form of the conservation law (3.1) on each volume 𝑉 ∈  at time 𝑡 ∈ 𝐼 is obtained, after
applying an integration by parts, as𝑑1𝑢 +𝐹 (𝑢) ⋅ 𝜈 d𝑠 = 0,
𝑑𝑡 𝑉 |𝑉 | ∫for 𝑉 ∈ ,

(3.2)

𝜕𝑉

where

𝑢𝑉 ≡ 𝑢𝑉 (𝑡) ∶=1
𝑢(𝑡, 𝑥) d𝑥|𝑉 | ∫

for 𝑉 ∈ , 𝑡 ∈ 𝐼,

𝑉

is the cell average value of 𝑢 on volume 𝑉 at time 𝑡. We assume that boundary of 𝑉 is denoted by 𝜕𝑉 and composed of distinct𝑞segments Γ1 , … , Γ𝑞 , i.e., 𝜕𝑉 = ∪𝑗=1 Γ𝑗 . Also, 𝜈 is the unit outward normal vector to 𝜕𝑉 . Then, equation (3.2) can be rewritten as
𝑞
𝑑
1 ∑
𝑢𝑉 +
𝐹 (𝑢(𝑡, 𝑠)) ⋅ 𝜈𝑗 d𝑠 = 0,𝑑𝑡
|𝑉 | 𝑗=1 ∫

(3.3)

Γ𝑗

where 𝜈𝑗 are normals to boundary segments Γ𝑗 . The boundary integrals in (3.3) can be approximated by a 𝑛-point Gaussian quadrature (or else) to obtain
𝑞 𝑛
𝑑
1 ∑∑
𝑢𝑉 +
𝜔 𝐹 (𝑢(𝑡, 𝑧𝑗,𝑘 )) ⋅ 𝜈𝑗 ≈ 0,𝑑𝑡
|𝑉 | 𝑗=1 𝑘=1 𝑗,𝑘(3.4)

where 𝜔𝑗,𝑘 and 𝑧𝑗,𝑘 for 𝑘 = 1, … , 𝑛 are integration weights and integration points, respectively, on interface Γ𝑗 .

The next step replaces the exact ﬂux 𝐹 (𝑢) ⋅ 𝜈 by a numerical ﬂux function 𝐹̃(𝑢+ , 𝑢− ; 𝜈) in order to couple the subproblems deﬁned
on each FVM cells, and to ensure that information in the problem travels in the direction of the characteristic curves of the equation.
Here, 𝑢+ and 𝑢− are solutions on the boundary 𝜕𝑉 from inside and outside, respectively. From the list of available numerical ﬂuxes,
this paper employs the Lax-Friedrichs numerical ﬂux function [33,34])1(
𝜎𝐹 (𝑢+ ) + 𝐹 (𝑢− ) ⋅ 𝜈 + (𝑢+ − 𝑢− )𝐹̃(𝑢+ , 𝑢− ; 𝜈) =22
where 𝜎 is the maximum local speed and is obtained by
]
[
𝑑𝑓
𝑑𝑓1
,…, 𝑑 .
𝜎=max
|𝐹 ′ (𝑢) ⋅ 𝜈𝑗 |, 𝐹 ′ (𝑢) =𝑑𝑢𝑑𝑢min(𝑢+ ,𝑢− )≤𝑢≤max(𝑢+ ,𝑢− )When substituted into (3.4), 𝑢+ (⋅, 𝑧𝑗,𝑘 ) represents the solution at integration point 𝑧𝑗,𝑘 within volume 𝑉 itself, whereas 𝑢− (⋅, 𝑧𝑗,𝑘 )represents the solution at the same point, but obtaining from an adjacent volume that shares Γ𝑗 as a common interface with volume
𝑉.
Finally, the cell average values {𝑢𝑉 (𝑡)}𝑉 ∈ are updated by solving the system of ODEs𝑑𝑢 (𝑡) = 𝐿𝑉 (𝑢𝑉 (𝑡)),𝑑𝑡 𝑉𝑉 ∈ ,

(3.5)

where

𝐿𝑉 (𝑢𝑉 (𝑡)) = −

𝑞 𝑛
1 ∑∑
𝜔 𝐹̃(𝑢+ (𝑡, 𝑧𝑗,𝑘 ), 𝑢− (𝑡, 𝑧𝑗,𝑘 ); 𝜈𝑗 ).
|𝑉 | 𝑗=1 𝑘=1 𝑗,𝑘(3.6)

Reconstructing 𝑢+ and 𝑢− from the current cell-averaged values {𝑢𝑉 (𝑡)}𝑉 ∈ is essential for updating the solution from time level 𝑡
to the next time level 𝑡 + Δ𝑡. Detailed explanations of some numerical methods for reconstruction will be provided later.
3.2. Time discretization
The system (3.5) can be integrated in time using an appropriate numerical ODE solver. The strong stability preserving (SSP)
methods were originally developed for time integration in semi-discretizations of hyperbolic conservation laws [11–13,46]. In thispaper, we employ the optimal third-order Shu-Osher SSP Runge-Kutta scheme [13,46]𝑢(1)= 𝑢𝑉 (𝑡𝑗 ) + Δ𝑡𝐿𝑉 (𝑢𝑉 (𝑡𝑗 )),
𝑉
3
1
1
𝑢(2)
= 𝑢𝑉 (𝑡𝑗 ) + 𝑢(1)+ Δ𝑡𝐿𝑉 (𝑢(1)
),
𝑉
𝑉
4
4 𝑉
4
1
2
2
𝑢𝑉 (𝑡𝑗+1 ) = 𝑢𝑉 (𝑡𝑗 ) + 𝑢(2)+ Δ𝑡𝐿𝑉 (𝑢(2)
),
𝑉
3
3 𝑉
3
4

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

to move from data {𝑢𝑉 (𝑡𝑗 )}𝑉 ∈ at time 𝑡𝑗 = 𝑗Δ𝑡 to data {𝑢𝑉 (𝑡𝑗+1 )}𝑉 ∈ at time 𝑡𝑗+1 . The above three-stage third-order scheme isdenoted by SSPRK(3,3). The global order of (Δ𝑡)3 is achieved if the CFL condition𝑟𝑉Δ𝑡 ≤ min(3.7)

𝑉 ∈ 𝑠max

is satisﬁed [30,36]. Here 𝑟𝑉 denotes the radius of the inscribed ball of volume 𝑉 , and 𝑠max = max |𝐹 ′ (𝑢) ⋅ 𝜈|, with the maximum takenover all integration points on the boundary of 𝑉 . In cases where the order of spatial discretization exceeds three, we may choose toemploy SSPRK(3,3) with a very smaller time step, or use the ﬁve-stage fourth-order scheme SSPRK(5,4) [11].
3.3. Reconstruction
In the reconstruction step, our goal is to approximate the density values 𝑢+ and 𝑢− in equation (3.6) (at integration points on cellinterfaces) by reconstructing the function 𝑢 from the cell average values {𝑢𝑉 (𝑡)}𝑉 ∈ at each time step. Reconstruction schemes are
typically based on piecewise polynomials as seen in, for example, [34]. However, in this paper we explore a reconstruction method
based on polyharmonic spline interpolation, as detailed in [2,3,23].
Suppose we have a ﬁxed volume 𝑉 ∈  , and consider its corresponding stencil

 ∶= {𝐶}𝐶∈ ⊂ ,
where 𝑉 ∈  . Each cell 𝐶 ∈  is associated with a linear functional (cell average functional), represented as 𝓁𝐶 , deﬁned as follows:𝓁𝐶 (𝑢) =1
𝑢(𝑥) d𝑥,|𝐶| ∫for all 𝐶 ∈ , 𝑢(𝑥) ≡ 𝑢(𝑡, 𝑥).

(3.8)

𝐶

To reconstruct cell averages {𝓁𝐶 } within the stencil  , we employ a generalized RBF interpolation method [53], which is expressedas𝑠(𝑥) =∑𝐶∈

𝛼𝐶 𝓁𝐶𝑦 𝜙(‖𝑥 − 𝑦‖) +

𝑄
∑
𝑗=1

𝛽𝑘 𝑝𝑗 (𝑥).

(3.9)

𝑦

Here, 𝓁𝐶 denotes that the functional 𝓁𝐶 operates on the variable 𝑦 of 𝜙, i.e.,𝓁𝐶𝑦 𝜙(‖𝑥 − 𝑦‖) =1𝜙(‖𝑥 − 𝑦‖) d𝑦,|𝐶| ∫𝐶 ∈ .

𝐶

The coeﬃcients 𝛼𝐶 and 𝛽𝑗 in (3.9) are determined by imposing the interpolation conditions

𝓁𝐶 ′ (𝑢) = 𝓁𝐶 ′ (𝑠),∀𝐶 ′ ∈ ,

along with the side condition∑
𝐶 ′ ∈

𝛼𝐶 ′ 𝓁𝐶 ′ (𝑝𝑗 ) = 0,

𝑗 = 1, 2, … , 𝑄.

The interpolation conditions give𝓁𝐶 ′ (𝑢) = 𝓁𝐶 ′ (𝑠) =∑𝐶∈

𝛼𝐶 𝓁𝐶𝑥 ′ 𝓁𝐶𝑦 𝜙(‖𝑥 − 𝑦‖) +

𝑄
∑
𝑗=1

𝛽𝑘 𝓁𝐶 ′ (𝑝𝑗 ),𝐶, 𝐶 ′ ∈ which together with the side condition can equivalently be expressed as solving the system of linear equations

[

𝐴𝐿
𝑃𝐿𝑇

𝑃𝐿
0

][ ] [ ]
𝛼
𝓁
=  ,
𝛽
0

(3.10)

where

𝐴𝐿 = (𝓁𝐶𝑥 ′ 𝓁𝐶𝑦 𝜙(‖𝑥 − 𝑦‖)),𝐶, 𝐶 ′ ∈ ,

𝑃𝐿 = (𝓁𝐶 (𝑝𝑗 )),𝐶 ∈ , 1 ⩽ 𝑗 ⩽ 𝑄,𝓁 = (𝓁𝐶 (𝑢)),𝐶 ∈ .

(3.11)

To ensure a unique solution, the set {𝓁𝐶 }𝐶∈ must satisfy the property of being ℙ𝑑𝑚 -unisolvent, meaning that if 𝓁𝐶 (𝑝) = 0 for all𝐶 ∈  , then 𝑝 ≡ 0. In other words, any polynomial 𝑝 ∈ ℙ𝑑𝑚 should be uniquely reconstructed from cell averages {𝓁𝐶 (𝑝)}𝐶∈ .
Nonetheless, a high-order reconstruction can lead to non-physical oscillations near discontinuities and shock fronts, and consequently can destabilize the solution.
5

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

3.4. WENO reconstruction
To suppress nonphysical oscillations near discontinuities and steep fronts, various techniques have been developed to enhance
the reconstruction process [34]. In this section, we review the WENO reconstruction, which represents a weighted adaptation of
the originally formulated ENO schemes. An ENO scheme was initially introduced in [16] for one-dimensional conservation laws.
In the WENO reconstruction, ﬁrst proposed in [35,27], a set of neighboring stencils {𝑘 , 𝑘 = 1, … , 𝐾} is selected for each control
volume 𝑉 , and in each stencil a reconstruction function is computed and weighted using an oscillation indicator that quantiﬁes thesmoothness of each function. Finally, the weighted (convex) combination of all 𝐾 reconstructed functions is used as the ultimate
reconstruction. In contrast, ENO schemes only utilize the smoothest reconstruction and discard the others. Diﬀerent variants of
ENO and WENO reconstructions have been developed and are frequently employed for solving various types of hyperbolic PDEs.
For detailed examples refer to [45]. The authors of [2,3,23], employed the polyharmonic spline interpolation for reconstruction, as
discussed in Section 3.3. They deﬁne the oscillation indicator  by𝑠 ∈ BL𝑘 (ℝ𝑑 ),(𝑠) = |𝑠|BL𝑘 (ℝ𝑑 ) ,as the Beppo-Levi space BL𝑘 (ℝ𝑑 ) is the optimal recovery space for polyharmonic splines [53]. A reconstructed function 𝑠𝑗 (associatedwith stencil 𝑗 ) with a smaller (𝑠𝑗 ) is assumed to be smoother and should therefore receive a higher weight. They introduce
intermediate values𝑤̃ 𝑗 =

1,(𝜖 + (𝑠𝑗 ))𝜌

𝑗 = 1, 2, … , 𝐾

for some 𝜖, 𝜌 > 0, and then compute the ﬁnal weights𝑤̃ 𝑗
𝑤 𝑗 = ∑𝐾

̃𝑖
𝑖=1 𝑤

,

𝑗 = 1, 2, … , 𝐾,

ensuring that 𝑤1 + ⋯ + 𝑤𝐾 = 1. The WENO reconstruction is then deﬁned as𝑠 = 𝑤1 𝑠 1 + ⋯ + 𝑤𝐾 𝑠 𝐾 .
To ﬁnd other RBF-ENO/WENO reconstructions see [4,18,19,56], and to explore other variations involving diﬀerent approximation
techniques, refer to [5,25,26].
To render the WENO process applicable to a wide range of problems, the conﬁguration and size of stencils are crucial andsigniﬁcantly impact the quality and accuracy of the reconstruction. Three types of such stencils are experimentally examined andintroduced in [3].
4. Approximate interpolation and smoothing
Our approach to avoid nonphysical oscillations in the solution (see Section 5 below) is motivated by approximate interpolation
on scattered data points. Let us consider 𝑁 centers 𝑋 = {𝑥1 , … , 𝑥𝑁 } distributed within the bounded region Ω ⊂ ℝ𝑑 , along with
their corresponding values, 𝑦1 , … , 𝑦𝑁 ∈ ℝ. We assume that these values are sampled from an unknown function 𝑓 ∈ 𝐶(Ω) where𝑦𝑗 = 𝑓 (𝑥𝑗 ) for 𝑗 = 1, … , 𝑁 . The aim is to ﬁnd a function 𝑠 ∈ 𝐶(Ω) that approximates the function 𝑓 by utilizing the available valuesof 𝑓 at 𝑋 . In an interpolation process, the error function 𝑓 − 𝑠 vanishes at point set 𝑋 while in an approximation process, such asa curve ﬁtting method, small errors at centers are admissible. When the data values are noisy or inaccurate, the second approachprovides a better and smoother overall approximation. This process is referred to as approximate interpolation if we permit slightdeviations from interpolation constraints to regularize the solution. This concept is related to Tikhonov regularization and smoothing
splines in the context of approximation theory [40].
Suppose  ⊂ 𝐶(Ω) is a normed linear space of continuous functions. Also, assume that  is another normed linear space, and𝑇 ∶  →  is a linear mapping. If 𝐽 ∶  → ℝ is deﬁned as 𝐽 (𝑠) = |𝑇 𝑠|, then for the given data sites 𝑋 = {𝑥1 , … , 𝑥𝑁 } and theircorresponding values 𝑓 (𝑥1 ), … , 𝑓 (𝑥𝑁 ) ∈ ℝ, the smoothing spline 𝑠𝜆 ∈  is deﬁned as the solution of the following minimization
problem:

{𝑁
}
∑
2
min
[𝑓 (𝑥𝑗 ) − 𝑠(𝑥𝑗 )] + 𝜆𝐽 (𝑠) ∶ 𝑠 ∈  ,

(4.1)

𝑗=1

where 𝜆 ⩾ 0 is a smoothing parameter. If 𝜆 gets smaller toward zero, then (4.1) tends to be an interpolation problem. Too large
values of 𝜆, on the other hand, violate the interpolation conditions and result in a solution which might smear out the data too much.
Thus, determining a proper value for this parameter is a challenging task.
Duchon [8] (see also Wendland and Rieger [54]) considered the smoothing spline problem using polyharmonic splines 𝜙𝑑,𝑘 (‖𝑥‖),𝑥 ∈ ℝ𝑑 . The approximation space  is assumed to be the Beppo-Levi space = BL𝑘 (ℝ𝑑 ),𝑘 > 2𝑑,and the functional 𝐽 on this space is deﬁned as
6

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

Fig. 1. Standard PHS interpolation vs. the smoothed PHS interpolation. (For interpretation of the colors in the ﬁgure(s), the reader is referred to the web version ofthis article.)

𝐽 (𝑓 ) ∶= |𝑓 |2BL (ℝ𝑑 ) .
𝑘

The solution of the minimization problem is then given by𝑠(𝑥) =

𝑁∑𝑗=1

𝛼𝑗 𝜙𝑑,𝑘 (‖𝑥 − 𝑥𝑗 ‖) +

𝑄
∑
𝑗=1

𝛽𝑗 𝑝𝑗 (𝑥)(4.2)

where 𝜙𝑑,𝑘 is the PHS kernel (2.4), {𝑝1 , … , 𝑝𝑄 } is a basis for polynomial space ℙ𝑑𝑘−1 , and the coeﬃcients 𝛼𝑗 and 𝛽𝑗 solve the linear
system of equations[

𝐴 + 𝜆𝐼𝑃𝑇𝑃
0][ ] [ ]
𝛼
𝑓
= 𝑋
𝛽
0

(4.3)

provided that 𝑋 is a ℙ𝑑𝑘−1 -unisolvent set. Here, 𝜙𝑑,𝑘 is considered to be conditionally positive deﬁnite of order 𝑘 instead of 𝑘 −⌈𝑑∕2⌉ + 1. The matrix 𝐼 is the identity matrix of size 𝑁 . The proof for a generalized form of this problem is given in Theorem 4.2
below.
The diﬀerence between the approximate interpolation system in (4.3) and the pure interpolation in (2.2) is a diagonal incrementwith 𝜆, causing the interpolant to pass close to interpolation points (if 𝜆 is small) but not exactly through them. For a given set ofnoisy data, this leads to a smooth (non-oscillatory) approximation. See Fig. 1. However, this type of smoothing is not appropriate forsuppressing the oscillations near discontinuities in the FVM setting.
In this paper, we introduce a more general approximate interpolation problem that allows diﬀerent weights for diﬀerent interpolation points:

min

{𝑁
∑ 1
𝑗=1

𝜆𝑗

}
[𝑓 (𝑥𝑗 ) − 𝑠(𝑥𝑗 )]2 + |𝑓 |2BL (ℝ𝑑 ) ∶ 𝑠 ∈ BL𝑘 (ℝ𝑑 )
𝑘

,

(4.4)

where 𝜆 is replaced by 𝑁 given smoothing parameters 𝜆𝑗 > 0, each corresponding to interpolation point 𝑥𝑗 for 𝑗 = 1, … , 𝑁 . To
obtain the form of the solution of this problem we need the following result.
Lemma 4.1. Assume that 𝜙 is the native space corresponding to conditionally positive deﬁnite kernel 𝜙 with respect to polynomial spaceℙ𝑑𝑚 = span{𝑝1 , … , 𝑝𝑄 }. For a given ℙ𝑑𝑚 -unisolvent set 𝑋 = {𝑥1 , … , 𝑥𝑁 } ⊂ Ω, assume that 𝑃 = (𝑝𝑗 (𝑥𝑖 )) ∈ ℝ𝑁×𝑄 , and 𝐻𝜙,𝑋 is a ﬁnite
dimensional space of the form

𝐻𝜙,𝑋 ∶=

{𝑁
∑
𝑗=1

}
𝑁

𝑇

𝛼𝑗 𝜙(‖ ⋅ −𝑥𝑗 ‖) ∶ 𝛼 ∈ ℝ with 𝑃 𝛼 = 0.

Moreover, assume that  is the orthogonal complement of subspace 𝐻𝜙,𝑋 + ℙ𝑑𝑚 of 𝜙 . Then for every function 𝑣 ∈  the vector 𝑣𝑋 =(𝑣(𝑥1 ), … , 𝑣(𝑥𝑁 ))𝑇 is perpendicular to the null space of 𝑃 𝑇 , i.e., 𝑣𝑇𝑋 𝛼 = 0 provided that 𝑃 𝑇 𝛼 = 0.
Proof. According to deﬁnition of space 𝐻𝜙,𝑋 , every function 𝑤 ∈ 𝐻𝜙,𝑋 + ℙ𝑑𝑚 can be written as 𝑤 = 𝜇 𝑥 𝜙(‖ ⋅ −𝑥‖) + 𝑞 with 𝑞 ∈ ℙ𝑑𝑚
and a certain linear functional

𝜇=

𝑁
∑
𝑗=1

𝛼𝑗 𝛿𝑥𝑗7

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

for some coeﬃcients 𝛼𝑗 with 𝜇(𝑝) = 0 for all 𝑝 ∈ ℙ𝑑𝑚 . Here, 𝛿𝑥 is the point evaluation functional; 𝛿𝑥 𝑓 ∶= 𝑓 (𝑥). The superscript 𝑥 on
functional 𝜇 𝑥 means that 𝜇 acts on variable 𝑥 of 𝜙(‖𝑦 − 𝑥‖). Consider a set of ℙ𝑑𝑚 -unisolvent points 𝐸 = {𝜉1 , … , 𝜉𝑄 } ⊂ Ω and assume
that Π ∶ 𝜙 → ℙ𝑑𝑚 is the polynomial interpolation operator on centers 𝐸 . Also assume that 𝐺(⋅, 𝑥) = 𝜙(‖ ⋅ −𝑥‖) − Π𝑥 𝜙(‖ ⋅ −𝑥‖), thepointwise polynomial interpolation error of basis function 𝜙 with respect to the second argument. Now, for a function 𝑣 in  wehave

0 = ⟨𝑣, 𝑤⟩𝜙
= ⟨𝑣, 𝜇 𝑥 𝜙(‖ ⋅ −𝑥‖) + 𝑞⟩𝜙
= ⟨𝑣, 𝜇 𝑥 𝜙(‖ ⋅ −𝑥‖)⟩𝜙= 𝜇(Π𝑣) + ⟨𝑣, 𝜇 𝑥 𝐺(⋅, 𝑥)⟩𝜙= 𝜇(𝑣) =𝑁∑
𝑗=1

𝛼𝑗 𝑣(𝑥𝑗 )for every function 𝑤 ∈ 𝐻𝜙,𝑋 + ℙ𝑑𝑚 . The third equality holds true because ⟨𝑣, 𝑞⟩𝜙 = 0. For the fourth and ﬁfth equalities we use thefact that 𝜇(Π𝑣) = 0, and Theorem 10.17 of [53]. This means that for every vector 𝛼 in the null space of 𝑃 𝑇 we have𝑁∑𝑗=1

𝛼𝑗 𝑣(𝑥𝑗 ) = 0. □Theorem 4.2. The solution of problem (4.4) is given by (4.2) where 𝛼 and 𝛽 satisfy the following system of equations

[

𝐴+Λ 𝑃
0
𝑃𝑇

][ ] [ ]
𝛼
𝑓
= 𝑋
𝛽
0

(4.5)

for Λ = diag{𝜆1 , … , 𝜆𝑁 }. Here, 𝐴, 𝑃 , and 𝑓𝑋 are deﬁned as (2.3) for 𝜙 = 𝜙𝑑,𝑘 and 𝑚 = 𝑘 − 1.
Proof. The native space corresponding to kernel 𝜙𝑑,𝑘 is known to be the Beppo-Levi space BL𝑘 (ℝ𝑑 ), and 𝐻𝜙𝑑,𝑘 ,𝑋 + ℙ𝑑𝑚 is a subspaceof the native space [53, Chapter 10]. Every function 𝑠 ∈ BL𝑘 (ℝ𝑑 ) can be written as 𝑠 = 𝑤 + 𝑣 where 𝑤 ∈ 𝐻𝜙𝑑,𝑘 ,𝑋 + ℙ𝑑𝑚 and 𝑣 is someelement in  ⊂ BL𝑘 (ℝ𝑑 ), the orthogonal complement of subspace 𝐻𝜙𝑑,𝑘 ,𝑋 + ℙ𝑑𝑚 . In fact,

𝑤=

𝑁
∑
𝑗=1

𝛼𝑗 𝜙(‖ ⋅ −𝑥𝑗 ‖) +

𝑄
∑
𝑗=1

𝛽𝑗 𝑝𝑗

for some vectors 𝛼 ∈ ℝ𝑁 and 𝛽 ∈ ℝ𝑄 with 𝑃 𝑇 𝛼 = 0, and 𝑣 ∈ . It is known that the native space semi-norm of 𝑤 is 𝛼 𝑇 𝐴𝛼 as thenorm of polynomials in ℙ𝑑𝑚 is zero [53, Chapter 10]. This means that|𝑠|2BL (ℝ𝑑 ) = 𝛼 𝑇 𝐴𝛼 + |𝑣|2BL (ℝ𝑑 ) .
𝑘𝑘Using these, the target function of minimization problem (4.4) can be reformulated as

[
]𝑇
]
[
𝑇 (𝛼, 𝛽, 𝑣) ∶= 𝑓𝑋 − (𝐴𝛼 + 𝑃 𝛽 + 𝑣𝑋 ) Λ−1 𝑓𝑋 − (𝐴𝛼 + 𝑃 𝛽 + 𝑣𝑋 ) + 𝛼 𝑇 𝐴𝛼 + |𝑣|2BL (ℝ𝑑 )𝑘which must be minimized over all [𝛼 𝑇 𝛽 𝑇 ] ∈ ℝ𝑁+𝑄 with 𝑃 𝑇 𝛼 = 0 and all 𝑣 ∈ . First, for a given 𝑣 ∈  we minimize 𝑇 (𝛼, 𝛽, 𝑣) overall [𝛼 𝑇 𝛽 𝑇 ] ∈ ℝ𝑁+𝑄 with 𝑃 𝑇 𝛼 = 0. This is a quadratic optimization problem as the target function 𝑇 (𝛼, 𝛽, 𝑣) can be rewritten as[][ ]𝐴Λ−1 𝐴 + 𝐴 𝐴Λ−1 𝑃
𝛼
𝑇
−1
𝑇
−1
𝛽
𝑃 Λ 𝐴
𝑃 Λ 𝑃
[
]
−1
𝐴Λ (𝑓𝑋 − 𝑣𝑋 )+ (𝑓𝑋 − 𝑣𝑋 )𝑇 Λ−1 (𝑓𝑋 − 𝑣𝑋 ) + |𝑣|2BL (ℝ𝑑 ) .
− 2[𝛼 𝑇 𝛽 𝑇 ]𝑃 𝑇 Λ−1 (𝑓𝑋 − 𝑣𝑋 )𝑘

𝑇 (𝛼, 𝛽, 𝑣) =[𝛼 𝑇 𝛽 𝑇 ]Since 𝐴 is positive deﬁnite on the null space of 𝑃 𝑇 , the block matrix is positive deﬁnite, and the unique solution of the minimization
problem satisﬁes

[

or

[

𝐴Λ−1 𝐴 + 𝐴 𝐴Λ−1 𝑃
𝑃 𝑇 Λ−1 𝐴
𝑃 𝑇 Λ−1 𝑃𝐴Λ−1𝑃 𝑇 Λ−10−𝐼][][ ] []
𝛼𝐴Λ−1 (𝑓𝑋 − 𝑣𝑋 )=𝑇
−1𝛽𝑃 Λ (𝑓𝑋 − 𝑣𝑋 )

𝐴+Λ 𝑃
0
𝑃𝑇

][ ] [𝛼𝐴Λ−1=𝛽𝑃 𝑇 Λ−10−𝐼][

𝑓𝑋 − 𝑣𝑋0which is equivalent to8]

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

Fig. 2. Standard PHS interpolation (left) vs. the generalized smoothed PHS interpolation (right) for a discontinuous 1D function.

[

𝐴+Λ 𝑃
0
𝑃𝑇

][ ] []
𝛼𝑓 − 𝑣𝑋
= 𝑋
𝛽
0

(4.6)

as the premultiplied matrix is non-singular on the null space of 𝑃 𝑇 . Note that the linear system (4.6) is non-singular because 𝑃 isfull-rank and 𝐴 + Λ is positive deﬁnite on the null space of 𝑃 𝑇 (since 𝐴 has this property and Λ is non-negative). The target function
evaluated at optimal vectors 𝛼 and 𝛽 becomes (as a function of 𝑣)

𝑇 (𝛼, 𝛽, 𝑣) = 𝛼 𝑇 (𝐴 + Λ)𝛼 + |𝑣|2BL (ℝ𝑑 )𝑘= 𝛼 𝑇 (𝑓𝑋 − 𝑣𝑋 ) + |𝑣|2BL (ℝ𝑑 )= 𝛼 𝑇 𝑓𝑋 + |𝑣|2BL (ℝ𝑑 ) .

𝑘𝑘

In the ﬁnal equality we have used 𝛼 𝑇 𝑣𝑋 = 0 from Lemma 4.1. The optimal solution is clearly obtained at 𝑣 = 0, also rendering (4.6)
equivalent to (4.5). □In the setting given above, the polyharmonic kernel 𝜙𝑑,𝑘 can be replaced by other (conditionally) positive deﬁnite functions. In
this case, the space  and functional 𝐽 are replaced by the native space of the new kernel and the (semi-) norm deﬁned on the native
space, respectively.
In Fig. 2, the standard and generalized smoothing interpolations are plotted for a discontinuous 1D function. When we prioritize
the accuracy of the approximation on the left side of the discontinuity, with no emphasis on the right side, the generalized smoothing
approach is an appropriate choice. A comparable scenario arises in transport equations when waves propagate either to the right orleft side, requiring a careful selection of either a backward or forward stencil to guarantee the stability of solution.
In the sequel, we will use an extended and localized version of the above generalized smoothing interpolation problem to reconstruct function values on cell boundaries from cell averages in the course of the FVM algorithm for conservation laws.
5. FVM with weighted smoothed reconstruction
As previously described, the central idea behind WENO reconstruction technique is to employ one-sided stencils in regions of thedomain where the solution function changes rapidly or contains discontinuities and shocks. This process involves creating diﬀerent
stencils for each volume, interpolating within each of these stencils, calculating the oscillation indicator of each interpolant, andﬁnally determining a weighted interpolation. In this section, we propose an alternative technique based on a tricky smoothing
approach that employs only one central stencil, even in regions with high gradients or discontinuities. The new method suppresses
the oscillations by degrading the impact of unﬁtted interpolation points. This approach is built upon the generalized smoothing
technique we introduced in Section 4.
5.1. Reconstruction from cell averages using WSR
The smoothing spline problem (4.4) can be extended to the more general setting

min

{𝑁
∑ 1
𝑗=1

𝜆𝑗

}
[𝓁𝑗 (𝑓 ) − 𝓁𝑗 (𝑠)]2

+ |𝑓 |2BL (ℝ𝑑 ) ∶ 𝑠 ∈ BL𝑘 (ℝ𝑑 )
𝑘

,

(5.1)

where 𝓁𝑗 , 𝑗 = 1, … , 𝑁 are 𝑁 linearly independent functionals from BL∗𝑘 (ℝ𝑑 ), the dual space of BL𝑘 (ℝ𝑑 ), with the property that𝓁𝑗 (𝑝) = 0 for all 𝑗 = 1, … , 𝑁 and 𝑝 ∈ ℙ𝑑𝑚 implies that 𝑝 = 0. The same proof shows that the solution of (5.1) is of the form
9

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

𝑠(𝑥) =

𝑁∑𝑗=1

𝛼𝑗 𝓁𝑗𝑦 𝜙𝑑,𝑘 (‖𝑥 − 𝑦‖) +

𝑄
∑
𝑗=1

𝛽𝑗 𝑝𝑗 (𝑥)where 𝛼 and 𝛽 satisfy the following system of equations[

𝐴 𝐿 + Λ 𝑃𝐿
𝑃𝐿𝑇
0][ ] [
]
𝛼
𝓁(𝑓 )
=
𝛽
0

where Λ = diag{𝜆1 , … , 𝜆𝑁 }, and 𝐴, 𝑃 , and 𝓁(𝑓 ) are deﬁned as

𝐴𝐿 = (𝓁𝑖𝑥 𝓁𝑗𝑦 𝜙(‖𝑥 − 𝑦‖)),𝑃𝐿 = (𝓁𝑗 (𝑝𝑖 )),𝑖, 𝑗 = 1, … , 𝑁,

𝑗 = 1, … , 𝑁, 𝑖 = 1, … , 𝑄,𝓁(𝑓 ) = (𝓁1 (𝑓 ), … , 𝓁𝑁 (𝑓 ))𝑇 .
To apply this approximation on FV stencils, consider a speciﬁc volume 𝑉 and a central stencil  around it. Deﬁne the cell averagefunctionals {𝓁𝐶 ∶ 𝐶 ∈ } as (3.8), and for the purpose of reconstruction from these cell averages, employ the smoothing problem(5.1) with 𝑓 = 𝑢, 𝑁 = || (number of cells in  ), 𝓁𝑗 = 𝓁𝐶 , and 𝜆𝑗 = 𝜆𝐶 for 𝐶 ∈  . The reconstructed function then has the form𝑠(𝑥) =∑𝐶∈

𝛼𝐶 𝓁𝐶𝑦 𝜙𝑑,𝑘 (‖𝑥 − 𝑦‖) +

𝑄
∑
𝑗=1

𝛽𝑗 𝑝𝑗 (𝑥),where the coeﬃcients 𝛼 and 𝛽 satisfy the linear system[

𝐴 𝐿 + Λ 𝑃𝐿
𝑃𝐿𝑇
0][ ] [ ]
𝛼
𝓁
= 
0
𝛽

(5.2)

where

Λ = diag{𝜆𝐶 }𝐶∈ ,(5.3)

and 𝐴𝐿 , 𝑃𝐿 and 𝓁 are deﬁned as (3.11). Furthermore, we assume that

𝜆𝐶 = 𝜓(𝑢; 𝐶; ),𝐶 ∈where 𝜓 is a suitable nonnegative function that depends on the solution function 𝑢 and all cells in the stencil  . This function should
be chosen such that it gives higher weights to more important cells and lower weights to the others. We apply the following criterion:
The closer the values of the average of function 𝑢 on cells in stencil  are to its average value on the central stencil 𝑉 , the smaller thesmoothing parameters associated with these cells become, thereby granting them greater importance. In this paper, we use the function

𝜆𝐶 = 𝜓(𝑢; 𝐶, ) = |𝑢𝑉 − 𝑢𝐶 |𝜂 ,𝐶 ∈ , 𝜂 ≥ 1,

(5.4)

as a model for this criterion. See subsection 5.4 below. This process is referred to as weighted smoothed reconstruction (WSR) andyields a smooth approximation around the cell 𝑉 within the stencil  where the function 𝑓 has a discontinuity or a steep gradient.
The approximate interpolation obtained from WSR can be employed to compute the numerical ﬂux along the boundaries of volume
𝑉 . The remaining steps of the classical FVM are executed as usual. This new scheme is abbreviated as FVM-WSR.
We note that the WSR method can use not only the polyharmonic spline kernels but also other radial basis functions and various
types of approximation functions such as polynomials. However, in this paper, our focus is primarily on polyharmonic splines.
5.2. About conditioning of local systems
The issue of ill-conditioning arises in the PHS reconstruction when dealing with small meshsizes to set up the system (5.2). To
overcome this problem, one can solve the generalized interpolation problem (5.2) on a blown-up stencil ̂ by scaling the original
stencil  . Out of all the well-known RBFs, this technique is only applicable to PHS kernels (2.4) [7,22,37]. In our setting with special
functionals 𝓁𝐶 we use a method similar to the one described in [29]. Consider the weighted smoothed approximation 𝑠(𝑥) on the
stencil  = {𝐶}𝐶∈ corresponding to the volume 𝑉 ∈  . To perform this reconstruction on a scaled stencil, we use the mapping

𝑇𝑉 (𝑥) ∶=𝑥 − 𝑏𝑉∶  ⟶ ̂ℎwhere 𝑏𝑉 is the barycenter of volume 𝑉 , and ℎ is the meshsize. Then we can show that the weighted smoothed approximation on
the new shifted and scaled stencil ̂ = {𝐶̂ = 𝑇𝑉 (𝐶) ∶ 𝐶 ∈ } is equivalent to the original approximation, i.e.,𝑠(𝑥) =∑̂ ̂𝐶∈𝛼̂ 𝐶̂ 𝓁 𝑦̂̂ 𝜙(‖𝑥̂ − 𝑦‖)
̂ +
𝐶

𝑄
∑
𝑗=1

𝛽̂𝑗 𝑝𝑗 (𝑥),̂where 𝑥̂ = 𝑇𝑉 (𝑥) and the coeﬃcients 𝛼̂ and 𝛽̂ satisfy the linear system
10

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

Fig. 3. The blue triangle stands for central triangle 𝑉 , and the yellow triangles form a central stencil for 𝑉 . From left to right: the 𝑁 (𝑉 ) stencil, the 𝑀 (𝑉 ) stencil,
a subset of 𝑀 (𝑉 ) of size 7, a subset of 𝑀 (𝑉 ) of size 9.

[

̂ 𝑃̂𝐿
𝐴̂ 𝐿 + Λ𝑃̂𝐿𝑇0][ ] [ ]
𝛼̂
𝓁
= 
0
𝛽̂

for

̂ = ℎ−2𝑘+𝑑 Λ,Λwhere 2𝑘 − 𝑑 is the exponent of the polyharmonic kernel (2.4), and𝐴̂ 𝐿 = (𝓁 𝑥̂ 𝓁 𝑦̂ ′ 𝜙(‖𝑥 − 𝑦‖)),𝐶 𝐶̂̂ 𝐶̂ ′ ∈ ,𝐶,𝑃̂𝐿 = (𝓁𝐶̂ (𝑝𝑗 )),̂ 1 ⩽ 𝑗 ⩽ 𝑄,𝐶̂ ∈ ,𝓁 = (𝓁𝐶 (𝑢)),𝐶 ∈ .

The proof closely follows that given in [29, Chapter 3] with a slight modiﬁcation for additional matrix Λ in the coeﬃcient matrix
̂ in the scaled version. The matrices 𝐴̂ 𝐿 and 𝑃̂𝐿 are now independentof linear system (5.2). This matrix needs to be replaced by Λof the mesh parameter ℎ, and if we appropriately scale the matrix Λ, the condition number of the scaled system remains constant,̃ ⋅ (ℎ𝜂 ) where Λ̃ depends on the
independent of ℎ. For example, if the smoothing parameters 𝜆𝐶 are chosen as (5.4) then Λ = Λ̂ independent of ℎ and results in an ℎ-independent
gradient of 𝑢 and is independent of ℎ. Therefore, scaling Λ by ℎ2𝑘−𝑑−𝜂 renders Λ
condition number for the entire system. However, this scaling of Λ may not always yield the most accurate results. If the user opts
̂ may lead to some worse condition numbers for very small values of ℎ. In such cases, thefor a diﬀerent scaling, the contribution of Λ
adaptation of preconditioning techniques for saddle point systems would be helpful, although this aspect has not been addressed inthis work.
5.3. Central stencils
In this paper, we use a conforming triangulation to form the volumes. A requirement for applying polyharmonic spline RBFs in
the reconstruction step of FVM-WSR is that the stencil size matches or exceeds the associated polynomial space dimension. In theexperimental section we employ polyharmonic splines 𝜙2,2 (𝑟) = 𝑟2 log(𝑟), 𝜙2,3 (𝑟) = 𝑟4 log(𝑟) and 𝜙2,4 (𝑟) = 𝑟6 log(𝑟). The polynomial
space dimensions associated with these kernels in the two-dimensional space are 3, 6 and 10, respectively. We address how a central
stencil can be formed to fulﬁll the solvability requirement and yield the expected accuracy.
For a given volume (triangle) 𝑉 , the Neumann neighborhood of 𝑉 is denoted by 𝑁 (𝑉 ) and includes all triangles that sharean edge with 𝑉 . On the other hand, the Moore neighborhood of 𝑉 is denoted by 𝑀 (𝑉 ) and includes all triangles that share anedge or a vertex with 𝑉 . It is clear that 𝑁 (𝑉 ) ⊂ 𝑀 (𝑉 ). The Moore neighborhood is usually well distributed in all directions
around the central volume 𝑉 . Experimental results showed that the optimal stencil size for WSR depends on the smoothness of thesolution: smaller stencils are preferable for smooth solutions, while larger stencils are needed for solutions with discontinuities orsteep gradients. We may use a subset of 𝑀 (𝑉 ) with certain number of triangles for diﬀerent types of RBFs. In Fig. 3 some stencils
are depicted for a sample triangle 𝑉 (the blue triangle).
5.4. Smoothing parameters
The choice of smoothing parameters in (5.3) requires some consideration to achieve an approximation with minimal oscillationand optimal rate of convergence. This choice depends on the behavior of the solution within the stencil where we intend to reconstruct the function. As previously explained, the value of 𝜆𝐶 inversely correlates with the inﬂuence of volume 𝐶 on the weighted
approximation. For instance, volume 𝑉 holds the highest inﬂuence as 𝜆𝑉 = 0. When a discontinuity traverses through a stencil, sharp
diﬀerences between cell averages within that stencil can arise. To address this, the smoothing function should assign a small 𝜆𝐶 (resulting in a large weight) to a cell 𝐶 whose cell average closely aligns with that of volume 𝑉 . Conversely, cells with cell averages
signiﬁcantly diﬀerent from that of 𝑉 should be assigned a large 𝜆𝐶 (resulting in a small weight). To model such a behavior, we
employ the expression (5.4) for 𝜆𝐶 values, although the user can design other possible formulations. See Fig. 4 for diﬀerent powers𝜂.
Our experiments show that using a power value of 𝜂 = 3 yields accurate and less oscillatory solutions for problems with bothsmooth and discontinuities (and shock front) solutions. However, for problems with smooth solutions speciﬁcally, using a power
11

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

Fig. 4. A model for the smoothing function.

value of 𝜂 between 1 and 2 results in convergence orders that are closer to the expected theoretical orders. To introduce adaptability
into this approach for problems with discontinuous and shock solutions, we follow this procedure:
1. Calculate 𝑢𝑚𝑎𝑥 = max{|𝓁𝐶 (𝑢)|} and 𝑢𝑚𝑖𝑛 = min {|𝓁𝐶 (𝑢)|},
𝐶∈

√

𝐶∈

2. Compute 𝑝 = (𝑢𝑚𝑎𝑥 − 𝑢𝑚𝑖𝑛 )∕(ℎ ||),3. If 𝑝 < 𝐾 , set 𝜂 = 2; otherwise, set 𝜂 = 3.
Here, || represents the size of stencil, ℎ denotes the meshsize, and 𝐾 is a user-deﬁned constant. For our experiments, we havechosen 𝐾 = 2.
6. Experimental results
In this section, we apply FVM-WSR to both linear and nonlinear conservation law equations with smooth and discontinuous
initial conditions. We compare errors with those obtained using FVM-WENO. For temporal discretization, we employ SSPRK(3,3) inall examples except for one case, in which we use SSPRK(5,4). The code is implemented in MATLAB and executed on a machine
equipped with a Core(TM) i7-6700k CPU running at 4 GHz, and 16 GB of RAM.
Since all experiments are provided for two-dimensional hyperbolic PDEs, we employ the notation 𝑥 = (𝑥1 , 𝑥2 ) to denote a point
in ℝ2 . We have also implemented the method for various 1D conservation laws and obtained the expected results. However, in theinterest of brevity and to control the length of paper, we do not include them here. Instead, we focus on presenting 2D problems
exclusively.
In the reconstruction step, we utilize the PHS kernels 𝜙2,2 (𝑟) = 𝑟2 log(𝑟) (2nd order), 𝜙2,3 (𝑟) = 𝑟4 log(𝑟) (3rd order), and 𝜙2,4 (𝑟) =𝑟6 log(𝑟) (4th order), each with central stencils of sizes 7, 9, and 15, respectively. For the highest-order RBF, we employ SSPRK(5,4)
for time discretization while we use SSPRK(3,3) for the others. In all cases, the steplength Δ𝑡 is chosen small enough to satisfy the
CFL condition (3.7) on all stencils.
6.1. Linear advection with a smooth initial condition
Consider the linear advection equation

𝑢𝑡 + 𝑢𝑥1 + 𝑢𝑥2 = 0,
where 𝑢 ∶= 𝑢(𝑡, 𝑥) ∶ [0, 𝑡𝑓(6.1)] × [−0.5, 0.5]2 ⟶ ℝ, with periodic boundary conditions. We apply the initial condition:

))
))
( (
( (
1
1
sin2 𝜋 𝑥2 +,𝑢(0, 𝑥) = sin2 𝜋 𝑥1 +
2
2

(6.2)

as depicted in the top-left of Fig. 5. As time advances, the exact solution of (6.1) propagates in the direction [1, 1] and retains its
initial shape in an extended domain at 𝑡 = 1.
Table 2 presents the norm-1 and norm-inﬁnity errors at the ﬁnal time 𝑡 = 1 for both WSR and WENO methods. The results are
1 1 1 1
1
1obtained on a sequence of triangulations with meshsizes ℎ = , , , , and. The steplength Δ𝑡 is initially set to 0.02 for ℎ =8 16 32 64
128
8
and is halved successively as the values of ℎ are divided by 2. As shown, FVM with WSR achieves (approximately) the expected
convergence orders in norm-1 and norm-inﬁnity (also in norm 2, though not reported in the table), except for the norm inﬁnity of
the 4th order scheme, where the order appears to stay about 3.2. However, the experiment demonstrates a superior convergencecompared to the WENO scheme. The results of the WENO method depend on the selection of forward and backward stencils where
diﬀerent choices may yield more or less accurate solutions. In this study, we follow the criteria outlined in [3] for our selection
12

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

Table 2The 𝓁1 -norm and 𝓁∞ errors (L1E and InfE) and orders of convergence for the numerical solutions of the advection equation (6.1) with initialcondition (6.2) at 𝑡 = 1 using FVM-WSR and FVM-WENO, with polyharmonic spline RBFs 𝑟2 log 𝑟 (2nd order), 𝑟4 log 𝑟 (3rd order), and 𝑟6 log 𝑟 (4thorder).

ℎ

W
S
R

W
E
N
O

2nd order

3rd order

4th order

L1E

𝑞

InfE

𝑞

L1E

𝑞

InfE

𝑞

L1E

𝑞

InfE

𝑞

1∕8
1∕16
1∕32
1∕64
1∕128

3.05e − 2
7.11e − 3
1.94e − 3
5.37e − 4
1.40e − 4

−
2.10
1.87
1.85
1.94

1.17e − 1
3.57e − 2
8.59e − 3
3.17e − 3
7.46e − 4

−
1.71
2.06
1.44
2.09

2.14e − 2
3.28e − 3
4.53e − 4
5.82e − 57.39e − 6

−
2.70
2.86
2.96
2.98

6.80e − 2
1.11e − 2
1.43e − 3
1.85e − 42.55e − 5

−
2.61
2.96
2.94
2.86

4.21e − 3
4.24e − 4
2.80e − 5
2.01e − 6
1.37e − 7

−
3.31
3.92
3.79
3.88

2.17e − 22.74e − 3
3.04e − 4
3.17e − 5
3.42e − 6

−
2.98
3.17
3.26
3.21

1∕8
1∕16
1∕32
1∕64
1∕128

3.02e − 2
8.79e − 3
2.33e − 3
6.96e − 4
2.18e − 4

−
1.78
1.91
1.74
1.68

9.41e − 2
3.12e − 2
9.69e − 3
4.02e − 3
1.25e − 3

−
1.59
1.69
1.27
1.69

2.05e − 2
6.22e − 3
1.15e − 3
2.14e − 44.45e − 5

−
1.72
2.43
2.42
2.27

6.73e − 2
2.11e − 2
6.02e − 3
1.21e − 3
3.29e − 4

−
1.67
1.81
2.31
1.88

−
−
−
−
−

−
−
−
−
−

−
−
−
−
−

−
−
−
−
−

Fig. 5. The initial condition (6.2) and numerical solutions using FVM-WSR for the advection equation at time levels 𝑡 = 0.25, 𝑡 = 0.5, and 𝑡 = 1.

Table 3
Runtime comparison: FVM-WSR vs. FVM-WENO.

ℎ

1
8

1
16

1
32

1
64

1
128

WENO
WSR

2.25

2.24

2.12

1.97

1.84

process. However, as there is no reported selection criterion for a 4th order RBF-based method in the mentioned reference or any
other sources, we have chosen not to include the corresponding results in the table.
1
The numerical solutions, surface plots and heatmaps, employing FVM-WSR with RBF 𝑟4 log 𝑟 and meshsize ℎ =are illustrated
128
in Fig. 5 at various time levels.
FVM-WSR oﬀers the advantage of a lower computational cost in comparison to FVM-WENO, as it employs only one stencil for eachreference volume. Table 3 presents a comparison of running times for both methods. Our experiments demonstrate (approximately)a 2x speedup with the new approach. There could still be some opportunities to speed up the WSR scheme when it comes to solving
local linear systems of the form (5.2) for a varying diagonal matrix Λ. We will pursue these possibilities in our future studies.
13

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

Fig. 6. The initial condition (6.3) (left plots) and numerical solutions using FVM-WSR for the advection equation at time levels 𝑡 = 0.25, 𝑡 = 0.5 and 𝑡 = 1.

6.2. Linear advection with steep gradient or discontinuous initial conditions
Consider equation (6.1) deﬁned on the same domain as before but with diﬀerent initial conditions

{

(
exp

𝑢(0, 𝑥) =‖𝑥−𝑐‖2‖𝑥−𝑐‖2 −𝑅20,), ‖𝑥 − 𝑐‖ < 𝑅,

(6.3)

otherwise,and{𝑢(0, 𝑥) =1, ‖𝑥 − 𝑐‖ < 𝑅,0, otherwise,(6.4)

where 𝑐 = (−0.25, −0.25) and 𝑅 = 0.1. These functions are depicted in the top-left of Figs. 6 and 7, respectively. As 𝑡 increases,
the initial proﬁles move in the direction [1, 1], but preserve their initial forms. Numerical solutions using FVM-WSR with RBF
1
𝜙2,3 (𝑟) = 𝑟4 log(𝑟) and meshsize ℎ = 128
are presented in Figs. 6 and 7 at various time levels. In the reconstruction step, central
stencils of 𝑀 (𝑉 ) type are employed. Both cases exhibit minimal damping and oscillations, with slight smearing observed in the
discontinuous solution. For a closer examination, we plot a cross-section of the solution along the direction [1, 1] (on the plane𝑥1 = 𝑥2 ) at the ﬁnal time in Fig. 8, and compare them with the exact and WENO solutions. We observe that the new method excels in
capturing discontinuities. The same meshsize and temporal steplength are used for both methods. We followed the approach given
in [3] to select the (forward and backward) stencils in the RBF-WENO method. However, one may opt for an alternative stencil
selection or use a diﬀerent oscillation indicator to enhance the result obtained by the WENO reconstruction.
6.3. Burgers’ equation with a smooth initial condition
Consider the nonlinear Burgers’ equation

𝑢𝑡 + 𝑢𝑢𝑥1 + 𝑢𝑢𝑥2 = 0,

(6.5)

where 𝑢 ∶= 𝑢(𝑡, 𝑥) ∶ [0, 𝑡𝑓 ] × [−0.5, 0.5]2 ⟶ ℝ, with periodic boundary conditions and with initial condition (6.3) for 𝑐 =(−0.25, −0.25) and 𝑅 = 0.1. As the time advances, a rarefaction develops behind the wave, while a shock forms in front of it.
1
Numerical solutions using FVM-WSR with ℎ =at time levels 𝑡 = 0.25, 𝑡 = 1 and 𝑡 = 2.3 are shown in Figs. 9. Cross-section so128
lutions on plane 𝑥1 = 𝑥2 at diﬀerent time levels are also shown in Fig. 10. As we observe from the ﬁgures, the shock is perfectly
captured and moves in the northeast direction without any signiﬁcant oscillations. The polyharmonic RBF 𝜙2,3 (𝑟) = 𝑟4 log(𝑟) andcentral stencils of type 𝑀 (𝑉 ) are used in the reconstruction step.
14

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

Fig. 7. The initial condition (6.4) (left plots) and numerical solutions using FVM-WSR for the advection equation at time levels 𝑡 = 0.25, 𝑡 = 0.5 and 𝑡 = 1.

Fig. 8. The cross-section solutions in plane 𝑥1 = 𝑥2 . Numerical solutions using WSR and WENO schemes with ℎ = 1∕128 are compared with exact solutions at ﬁnal
time 𝑡 = 1 with initial condition (6.3) (ﬁrst row) and initial condition (6.3) (second row).
15

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

Fig. 9. The initial condition (upper left) and numerical solutions using FVM-WSR for the Burgers’ equation at time levels 𝑡 = 0.25, 𝑡 = 1 and 𝑡 = 2.3.

6.4. Burgers’ equation with discontinuous initial condition
Consider the nonlinear Burgers’ equation (6.5) for 𝑢 ∶= 𝑢(𝑡, 𝑥) ∶ [0, 0.5] × [0, 1]2 ⟶ ℝ, and with initial condition⎧ −1.0, if 𝑥1 > 0.5, 𝑥2 > 0.5,
⎪ −0.2, if 𝑥1 < 0.5, 𝑥2 > 0.5,
𝑢(0, 𝑥) = ⎨⎪ 0.5, if 𝑥1 < 0.5, 𝑥2 < 0.5,
⎩ 0.8, if 𝑥1 > 0.5, 𝑥2 < 0.5.

(6.6)

The color contour plot of the (discontinuous) initial condition is shown in Fig. 11 (left panel). This is a challenging problem due tothe emergence of diﬀerent shocks and rarefaction waves as time increases. Following [19] we extend the problem domain to [−1, 2]2
and apply a periodic boundary condition. The results are ﬁnally presented in the original domain [0, 1]2 . The plot of the FVM-WSR
1solution with meshsize ℎ =at ﬁnal time 𝑡 = 0.5 is shown in Fig. 11. We have employed the polyharmonic RBF 𝜙2,3 (𝑟) and 𝑀 (𝑉 )128
stencils.
6.5. Kurganov-Petrova-Popov (KKP) rotating wave
Consider the rotating wave equation

𝑢𝑡 + sin(𝑢)𝑢𝑥1 + cos(𝑢)𝑢𝑥2 = 0,

(6.7)

where 𝑢 ∶= 𝑢(𝑡, 𝑥) ∶ [0, 1] × [−2, 2] × [−2.5, 1.5] ⟶ ℝ, with periodic boundary conditions and with initial condition
16

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

Fig. 10. Cross-section solutions at time levels 𝑡 = 0.25, 𝑡 = 1 and 𝑡 = 2.3 for the Burgers’ equation.

Fig. 11. The initial condition (6.6) (left) and FVM-WSR solution at ﬁnal time 𝑡 = 0.5.

{
𝑢(0, 𝑥) =3.5𝜋, 𝑥2 + 𝑦2 < 1,0.25𝜋, otherwise.

(6.8)

The KKP problem, which involves nonconvex advective ﬂuxes, was initially introduced by Kurganov [32] (see also [17,49]). Thisequation is speciﬁcally designed to serve as a test problem for evaluating diﬀerent schemes in handling entropy-violating solutions.
1
In Fig. 12, we present the FVM-WSR solution at the ﬁnal time 𝑡 = 1, computed with meshsize ℎ = . The reconstruction step employs
64
the polyharmonic RBF 𝜙2,3 and 𝑀 (𝑉 ) stencils.
17

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

Fig. 12. The initial proﬁle (left) and FVM-WSR solution at time 𝑡 = 1 (right) for the KKP rotating equation.

7. Conclusions
In this paper, we introduced the novel FVM-WSR method for solving scalar conservation laws on unstructured meshes. This
scheme incorporates a weighted smoothed reconstruction (WSR) and eﬀectively suppresses nonphysical oscillations in the presenceof discontinuities, shocks, and sharp fronts. Our experimental results show that this method achieves accuracy on par with (andsometimes superior to) other numerical approaches, such as RBF-WENO, while using only a single central stencil in its reconstruction
step. This property contributes to the method’s brevity, simplicity, and reduced computational cost. Additionally, this method canbe easily adapted to deal with more complex problems such as 3D hyperbolic problems and systems of conservation laws. In an
upcoming research paper, we will further extend the application of this technique for solving Euler’s and Navier-Stokes equations.
CRediT authorship contribution statement
Davoud Mirzaei: Conceptualization, Formal analysis, Investigation, Methodology, Resources, Software, Supervision, Validation,
Visualization, Writing – original draft, Writing – review & editing. Navid Soodbakhsh: Formal analysis, Investigation, Methodology,
Resources, Software, Validation, Visualization, Writing – original draft, Writing – review & editing.
Declaration of competing interest
The authors declare the following ﬁnancial interests/personal relationships which may be considered as potential competinginterests: Navid Soodbakhsh reports ﬁnancial support was provided by Institute for Research in Fundamental Sciences.
18

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh

Data availability
Data will be made available on request.
References[1] R. Abgrall, On ENO schemes on unstructured meshes: analysis and implementation, J. Comput. Phys. 114 (1994) 45–58.
[2] T. Aboiyar, E.H. Georgoulis, A. Iske, High order WENO ﬁnite volume schemes using polyharmonic spline reconstruction, in: Proceedings of the International
Conference on Numerical Analysis and Approximation Theory, Cluj-Napoca (Romania), 2006, pp. 1–14.
[3] T. Aboiyar, E.H. Georgoulis, A. Iske, Adaptive ADER methods using kernel-based polyharmonic spline WENO reconstruction, SIAM J. Sci. Comput. 32 (2010)
3251–3277.
[4] T. Arbogast, C.S. Huang, M.H. Kuo, RBF WENO reconstructions with adaptive order and applications to conservation laws, J. Sci. Comput. 91 (2) (2022) 51.
[5] D.S. Balsara, S. Garain, V. Florinski, W. Boscheri, An eﬃcient class of WENO schemes with adaptive order for unstructured meshes, J. Comput. Phys. 404 (2020)
109062.
[6] T.J. Barth, H. Deconinck, High Order Methods for Computational Physics, Springer Verlag, 1999.
[7] O. Davydov, R. Schaback, Optimal stencils in Sobolev spaces, IMA J. Numer. Anal. 39 (2019) 398–422.
[8] J. Duchon, Splines minimizing rotation-invariant semi-norms in Sobolev spaces, in: W. Schempp, K. Zeller (Eds.), Constructive Theory of Functions of Several
Variables, Springer, Berlin, 1977, pp. 85–100.
[9] G.E. Fasshauer, Meshfree Approximations Methods with Matlab, World Scientiﬁc, Singapore, 2007.
[10] O. Friedrich, Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids, J. Comput. Phys. 144 (1998) 194–212.
[11] S. Gottlieb, On high order strong stability preserving Runge-Kutta and multi step time discretization, J. Sci. Comput. 25 (2005) 105–127.
[12] S. Gottlieb, C.W. Shu, Total variation diminishing Runge-Kutta schemes, Math. Comput. 67 (1998) 73–85.
[13] S. Gottlieb, C.W. Shu, E. Tadmor, Strong stability preserving high order time discretization methods, SIAM Rev. 43 (2001) 89–112.
[14] J.L. Guermond, R. Pasquetti, B. Popov, Entropy viscosity method for nonlinear conservation laws, J. Comput. Phys. 230 (2011) 4248–4267.
[15] A. Harten, S.R. Chakravarthy, Multidimensional ENO schemes for general geometries, Technical Report, ICASE, 221:91–76, 1991.
[16] A. Harten, B. Engquist, S. Osher, S.R. Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes III, J. Comput. Phys. 71 (1987) 231–303.
[17] J.S. Hesthaven, Numerical Methods for Conservation Laws: From Analysis to Algorithms, SIAM, 2017.
[18] J.S. Hesthaven, F. Monkeberg, Entropy stable essentially nonoscillatory methods based on RBF reconstruction, ESAIM, M2AN 53 (3) (2019) 925–958.
[19] J.S. Hesthaven, F. Monkeberg, Two-dimensional RBF-ENO method on unstructured grids, J. Sci. Comput. 82 (3) (2020) 1–24.
[20] C. Hu, C.W. Shu, Weighted essentially non-oscillatory schemes on triangular meshes, J. Comput. Phys. 150 (1999) 97–127.
[21] T.J.R. Hughes, T.E. Tezduyar, Finite element methods for ﬁrst-order hyperbolic systems with particular emphasis on the compressible Euler equations, Comput.
Methods Appl. Mech. Eng. 45 (1984) 217–284.
[22] A. Iske, On the approximation order and numerical stability of local Lagrange interpolation by polyharmonic splines, in: International Series of Numerical
Mathematics, vol. 145, Birkhäuser Verlag, Basel, 2003, pp. 153–165.
[23] A. Iske, On the construction of kernel-based adaptive particle methods in numerical ﬂow simulation, in: Notes on Numerical Fluid Mechanics and Multidisciplinary Design (NNFM), Springer, Berlin, 2013, pp. 197–221.
[24] L. Ivan, C.P.T. Groth, High-order solution-adaptive central essentially non-oscillatory (CENO) method, J. Comput. Phys. 275 (2014) 830–862.
[25] Z. Ji, T. Liang, L. Fu, A class of new high-order ﬁnite-volume TENO schemes for hyperbolic conservation laws with unstructured meshes, J. Sci. Comput. 92 (2)(2022) 61.
[26] Z. Ji, T. Liang, L. Fu, High-order ﬁnite-volume TENO schemes with dual ENO-like stencil selection for unstructured meshes, J. Sci. Comput. 95 (3) (2023) 76.
[27] G. Jiang, C.W. Shu, Eﬃcient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228.
[28] C. Johnson, A. Szepessy, Shock-capturing streamline diﬀusion ﬁnite element methods for nonlinear conservation laws, in: T.J.R. Hughes, T.E. Tezduyar (Eds.),
Recent Developments in Computatioonal Fluid Dynamics, AMD, New York, 1988, pp. 75–92.
[29] L. Kadrnka, The Finite Volume Particle Method: a Meshfree Method of Second Order for the Numerical Solution of Hyperbolic Conservation Laws, PhD thesis,
Universität Hamburg, Hamburg, 2014.
[30] M. Kaser, A. Iske, ADER schemes on adaptive triangular meshes for scalar conservation laws, J. Comput. Phys. 205 (2005) 486–508.
[31] D. Kroner, Numerical Schemes for Conservation Laws, John Wiley and Sons, 1997.
[32] A. Kurganov, G. Petrova, B. Popov, Adaptive semidiscrete central-upwind schemes for nonconvex hyperbolic conservation laws, SIAM J. Sci. Comput. 29 (2007)
2381–2401.
[33] P.D. Lax, Weak solutions of nonlinear hyperbolic equations and their numerical approximation, Commun. Pure Appl. Math. 7 (1954) 159–193.
[34] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, Cambridge, 2002.
[35] X.D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994) 200–212.
[36] A. Meister, J. Struckmeier, Hyperbolic Partial Diﬀerential Equations: Theory, Numerics and Applications, Springer, 2002.
[37] D. Mirzaei, The direct radial basis function partition of unity (D-RBF-PU) method for solving PDEs, SIAM J. Sci. Comput. 43 (2021) A54–A83.
[38] K.W. Morton, T. Sonar, Finite volume methods for hyperbolic conservation laws, Acta Numer. 16 (2007) 155–238.
[39] M. Nazarov, J. Hoﬀman, Residual based artiﬁcial viscosity for simulation of turbulent compressible ﬂow using adaptive ﬁnite element methods, Int. J. Numer.
Methods Fluids 71 (2013) 339–357.
[40] L.L. Schumaker, Numerical Methods for Conservation Laws: From Analysis to Algorithms, Wiley-Interscience Publication, 1981.
[41] V. Shankar, G.B. Wright, A. Narayan, A robust hyperviscosity formulation for stable RBF-FD discretizations of advection-diﬀusion-reaction equations on manifolds, SIAM J. Sci. Comput. 42 (2020) A2371–A2401.
[42] J. Shi, Y.T. Zhang, C.W. Shu, Resolution of high order WENO schemes for complicated ﬂow structures, J. Comput. Phys. 186 (2003) 690–696.
[43] C.W. Shu, High order ENO and WENO schemes for computational ﬂuid dynamics, in: High Order Methods for Computational Physics, Springer, Berlin, 1991,
pp. 439–852.
[44] C.W. Shu, High order weighted essentially nonoscillatory schemes for convection dominated problems, SIAM Rev. 51 (2009) 82–126.
[45] C.W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes, Acta Numer. 29 (2020) 701–762.
[46] C.W. Shu, S. Osher, Eﬃcient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys. 77 (1988) 439–471.
[47] C.W. Shu, S. Osher, Eﬃcient implementation of essentially non-oscillatory shock capturing schemes, II, J. Comput. Phys. 83 (1989) 32–78.
[48] T. Sonar, On the construction of essentially non-oscillatory ﬁnite volume approximations to hyperbolic conservation laws on general triangulations: polynomial
recovery, accuracy and stencil selection, Comput. Methods Appl. Mech. Eng. 140 (1997) 157–181.
[49] V. Stiernström, L. Lundgren, M. Nazarov, K. Mattsson, A residual-based artiﬁcial viscosity ﬁnite diﬀerence method for scalar conservation laws, J. Comput. Phys.
430 (2021) 110100.
[50] I. Tominec, M. Nazarov, Residual viscosity stabilized RBF-FD methods for solving nonlinear conservation laws, J. Sci. Comput. 94 (14) (2023).
[51] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, Springer, 2013.
19

Journal of Computational Physics 508 (2024) 112981

D. Mirzaei and N. Soodbakhsh[52] J. Von Neumann, R.D. Richtmyer, A method for the numerical calculation of hydrodynamic shocks, J. Appl. Phys. 21 (1950) 232–237.
[53] H. Wendland, Scattered Data Approximation, Cambridge University Press, 2005.
[54] H. Wendland, C. Rieger, Approximate interpolation with applications to selecting smoothing parameters, Numer. Math. 101 (2005) 729–748.
[55] W.R. Wolf, J.L.F. Azevedo, High-order ENO and WENO schemes for unstructured grids, Int. J. Numer. Methods Fluids 55 (2007) 917–943.
[56] N. Zhan, R. Chen, Y. You, Three-dimensional high-order ﬁnite-volume method based on compact WENO reconstruction with hybrid unstructured grids, J.
Comput. Phys. 49 (2023) 112300.
[57] Y.T. Zhang, J. Shi, C.W. Shu, Y. Zhou, Numerical viscosity and resolution of high-order weighted essentially nonoscillatory schemes for compressible ﬂows withhigh Reynolds numbers, Phys. Rev. E 68 (2003) 046709.

20