Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems

Vincenzo Angelo Pennati; Antoine Celestin Kengni Jotsa; Jacques Tagoudjeu

doi:doi:10.11648/j.pamj.20251404.11

Research Article |

| Peer-Reviewed

Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems

Vincenzo Angelo Pennati^*

, Antoine Celestin Kengni Jotsa

, Jacques Tagoudjeu

Published in Pure and Applied Mathematics Journal (Volume 14, Issue 4)

Received: 8 June 2025 Accepted: 30 June 2025 Published: 15 August 2025

Views: Downloads:

Download PDF

Share This Article

Twitter
Linked In
Facebook

Abstract

Several recently published studies regarding flow problems propose schemes of high order of accuracy designed as evolution of traditional methods. A drawback common to these new schemes is the necessity to adopt uniform mesh refinement for solving sharp problems, by increasing the computational cost. Even the so called essentially non-oscillatory and weight essentially non-oscillatory methods suffer of the same drawback and are not suitable to cope with h-adaptive methods due to their definition on finite volumes necessarily of equal diameter. Therefore, in order to overcome the above drawback, the formulation of dynamically locally self h-adaptive processes is designed to achieve the dual purpose to increase the accuracy and to keep as small as possible the number of finite volumes. To define a locally h-adaptive finite volume (FV) scheme need two simple but important tools, namely a particular FV named Bridge FV positioned between two adjacent subdomains and the definition of suitable profiles approximating the fluxes on the FV faces. In this article a new FV method for the numerical solution of convective-diffusive 1D problems is developed. It is conservative, second order in time and space for equal FV, and allows the partitioning of the domain by equal or unequal finite volumes, thus dynamically locally self h-adaptive. The definition of the monotonic profiles is accomplished by means of cubic weighted ν-splines and Taylor expansions. The profile analysis respect to the numerical properties is conducted in the normalized plane with the velocity varying in time and space and gives the flux value on the FV faces. Moreover the flux is assigned by Upwind or by second order back-ward Characteristics if the estimated flux is outside of the unit square or the transformation into the normalized plane is not possible, respectively. The initial-boundary stability and convergence properties of the new method are examined in detail, also in presence of h-adaptivity. In addition, a generalization of the new scheme to 2D and 3D problems is presented. Finally, some numerical test are carried out to verify the properties of the new method, including two CFD problems.

Published in	Pure and Applied Mathematics Journal (Volume 14, Issue 4)
DOI	10.11648/j.pamj.20251404.11
Page(s)	69-92
Creative Commons	This is an Open Access article, distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution and reproduction in any medium or format, provided the original work is properly cited.
Copyright	Copyright © The Author(s), 2025. Published by Science Publishing Group

Keywords

FVM for Convective-diffusive Problems, High Order Monotonic Schemes, Convergence for Initial-boundary Values Problems, Dynamically Local Self h-Adaptivity

1. Introduction

Convection-diffusion equations are of fundamental importance in computational fluid dynamics (CFD) applications to solve 1D, 2D and 3D problems. Numerical solutions of convection-diffusion equations often use finite elements methods (FEM)

[1-7]

, finite differences methods (FDM)

[8, 9]

and finite volumes methods (FVM), both in cell-centered and cell-vertex schemes

[10, 11]

and their relevant generalizations

[12-17]

. In order to achieve better accuracy and lower numerical dissipation, high order FVM have been developed and intensively studied

[18-21]

. However, these methods can be very costly for solving problems containing violent change in states, due to uniform mesh refinement. Adaptive strategies have been introduced to overcome this serious drawback, gaining increasing scientific attention

[22-27]

Generally, the recently published studies regarding flow problems proposed schemes designed as evolutions of traditional methods or that belong to the class of essentially non-oscillatory (ENO) and weight-ENO (WENO) methods. For example, Ayalew et al. developed an Upwind FV scheme to solve the 1D convection-diffusion equation by reformulating the well-known Quick scheme in the method of deferred corrections

[28]

. Although alleviating stability problems, problematic negative coefficients needed to be introduced into the source term in order to retain the positive main coefficients. The scheme was conditionally stable, third order accurate in space and first-order in time.

Similarly, Peng proposed a new positivity-preserving finite volume scheme was proposed to solve the convection-diffusion equation on 2D or 3D distorted meshes

[29]

. The nonlinear two-point flux approximation was applied to the diffusion flux while the discretization of convection flux was based on the second-order Upwind method with a slope limiter. Time approximation was not given.

Kumar et al. presented the numerical solution of the incompressible Navier–Stokes equations by a finite-volume discretization over staggered grids. The velocity components were calculated by solving local boundary value problems, with the pressure gradient and the gradient of the transverse flux as source terms. The convection-diffusion operator was a weighted average of the Central and Upwind approximations, thus a second order convergence was provided without introducing any spurious oscillations. The time integration was pursued with a Runge–Kutta method

[30]

A novel difference scheme by Sukhinov et al. based on a linear weighted combination of the Upwind Leapfrog and the Standard Leapfrog schemes, with the weight coefficients determined by minimizing the approximation error, allowed to solve transport problems effectively, even when the Péclet number fell in the range from 2 to 20

[31]

Kurganov solved the 2D Saint-Venant system of non-linear hyperbolic partial differential equations in

[32]

. The system admits non-smooth solutions that may contain shock and rarefaction waves and that may break down even when the initial data are smooth. Therefore, a numerical scheme capable of preserving the steady-state solutions by ensuring the well balance between the fluxes and the source terms was required. Since the convective flux approximation should take advantage of up-winding feature, the space approximation was estimated by conservative-central-Upwind well-balanced second order scheme and the time approximation by a Runge-Kutta third order method.

Neelan et al. described a WENO scheme for solving Euler equations based on a weighting procedure that did not directly rely on the smoothness indicator. Once the smooth polynomials were determined, a new switch function managed the further steps. Contrary to other approaches, this switch function did not involve any conditional statements. The presented scheme used three three-points sub-stencils and achieved fourth-order accuracy in space, while the time approximation was accomplished by the method of lines

[33]

Alternatively, Kossaczká et al. improved the well-known fifth order WENO shock-capturing scheme was improved by using learning techniques

[34]

. The smoothness indicators of the WENO scheme were the products of the original smooth indicators, Jiang et al.

[35]

and of the multipliers outputs of a properly fine-tuned neural network algorithm. The fluxes were approximated with a five points stencil and the time approximation was accomplished by the method of lines, that is by a system of ordinary differential equations integrated by a third order TVD Runge-Kutta method. Some numerical tests, including the 1D Euler equations, were solved.

It is important to underline that the properties of conservation, up-winding, local monotonicity and small size of the stencils are required to all the methods for transport equations. The proposed schemes should provide an optimal compromise between high order accuracy in time and space and computational cost. However, the study for the approximation of boundary conditions and the application of local h-adaptivity are rarely performed. For example, according to Xu, high order ENO and WENO methods require usually rather large stencils which lead to an increase of the computational cost and the difficulty to deal with boundary conditions

[36]

Therefore, it is of great importance to develop numerical schemes with the same size stencil as that of the classical FVM on one hand and that are accurate and robust on the other hand for simulating complex flow phenomena.

In a previous article by Pennati et al., the method named FAST (Flow Approximation by Spline Techniques) for the solution of convective-diffusive 1D problems was defined

[37]

. The profiles were based on the union of two cubic weighted

ν

-spline polynomials (i.e. two third-degree polynomials satisfying suitable conditions relevant to their continuity and their first and second derivatives), so that the local monotonicity property was guaranteed by the definition itself and the construction of h-adaptive techniques was feasible.

Download: Download full-size image

Figure 1. Grids with halving (a) and doubling (b) the net size and Bridge FV.

In the current article, we propose a new scheme for solving 1D convective-diffusive problems named FASTEST (Fast with Estimated Second derivative Time expansion), which uses equal or unequal FV and is defined by applying the technique given by Roache

[38]

. The method is conservative, second order accurate in time and space (with equal FV), locally monotonic even with velocity u variable in time and space, and appropriate to develop dynamically locally self h-adaptive schemes.

The h-adaptive schemes are characterized by the generation of suitable partitions of the domain

Ω

into subdomains

Ω_{j},

each formed by equal FV (at least four), and by the presence of a single particular finite volume, Bridge finite volume, positioned between each couple of subdomains

Ω_{j}

and

Ω_{j + 1}

, see Figure 1(a) and 1(b).

Given our interest to solve problems respecting the restriction on the grid spacing

h \leq \frac{(2 Γ)}{u}, Γ

diffusion coefficient and u velocity, the focus of the paper is on the approximation of the convective term while the diffusive term is approximated by second order finite difference formulae. The article is organized as follows: after this Introduction, Section 2 gives details about the Fastest method. Section 3 examines the stability and convergence properties of the novel method. Section 4 presents the dynamically locally self h-adaptive method and the Pick & Roll technique for the generation of new partitions. Section 5 describes shortly the interpolation process used for the Characteristics and for the h-adaptive schemes. Section 6 generalizes Fastest to 2D and 3D problems. The results obtained by applying Fastest in test cases, including two CFD problems, are discussed in Section 7. Conclusions are drawn in Section 8. Finally, Appendix provides additional data relevant for the numerical tests.

2. Fastest Method

2.1. Scheme Definition

Let consider the unsteady convection-diffusion equation, defined on a bounded domain

Ω

\frac{\partial ξ}{\partial t} + u \frac{\partial ξ}{\partial x} - Γ \frac{\partial^{2} ξ}{\partial x^{2}} = S

(1)

and the appropriate initial and boundary conditions, where

ξ (x, t)

is the dependent variable,

u

the constant convective velocity agreed with

x

axis and never changing its direction,

Γ

the diffusion coefficient,

S

the source term and

\partial Ω

the boundary of

Ω

. Suppose that

Ω

is partitioned with FV of equal diameter

D

respecting the condition

D \leq (2 Γ) / u

, with

uL / Γ = R

cell Reynolds number and L characteristic length, thus the numerical solution of the transport equation (1) is qualitatively similar to the one of a parabolic differential equation dominated by the convective term.

To obtain a high order in time and space profile, we consider the Taylor expansion in time

ξ^{n + 1} = ξ^{n} + Δt \frac{\partial ξ}{\partial t} + \frac{1}{2} Δ t^{2} \frac{\partial^{2} ξ}{\partial t^{2}} + O (Δ t^{3})

(2)

and result first- and second- time derivatives from the homogeneous advection equation

\frac{\partial ξ}{\partial t} = - u \frac{\partial ξ}{\partial x} and \frac{\partial^{2} ξ}{{\partial t}^{2}} = u^{2} \frac{\partial^{2} ξ}{{\partial x}^{2}}

(3)

by replacing both the expressions (3) into (2), we have (see

[38]

)

ξ^{n + 1} = ξ^{n} - uΔt \frac{\partial ξ}{\partial x} + \frac{1}{2} {u^{2} Δt}^{2} \frac{\partial^{2} ξ}{{\partial x}^{2}} + O ({Δt}^{3}) .

(4)

By adopting second-order FV approximations

{[\frac{\partial ξ}{\partial x}]}_{i}^{n}

and

{[\frac{\partial^{2} ξ}{{\partial x}^{2}}]}_{i}^{n}

for

\frac{\partial ξ}{\partial x}

and

\frac{\partial^{2} ξ}{\partial x^{2}}

respectively, we obtain

ξ_{i}^{n + 1} = ξ_{i}^{n} - uΔt {[\frac{\partial ξ}{\partial x}]}_{i}^{n} + \frac{1}{2} {u^{2} Δt}^{2} {[\frac{\partial^{2} ξ}{{\partial x}^{2}}]}_{i}^{n} + O ({Δt}^{3}, {Δx}^{3}) .

(5)

Integrating this approximation we obtain

\int_{FV} ξ_{i}^{n + 1} dx = \int_{FV} ξ_{i}^{n} dx - u_{i}^{n} Δt \int_{FV} {[\frac{\partial ξ}{\partial x}]}_{i}^{n} dx + \frac{1}{2} {(u_{i}^{n})}^{2} {Δt}^{2} \int_{FV} {[\frac{\partial^{2} ξ}{{\partial x}^{2}}]}_{i}^{n} dx .

(6)

Let apply the Mean Value Integral theorem and indicate with

D_{i}

the diameter of a generic FV

ξ_{i}^{n + 1} = ξ_{i}^{n} - \frac{u_{i}^{n} Δt}{D_{i}} [F_{r}^{n} - F_{l}^{n}]

(7)

where

F_{r}^{n} = ξ_{r}^{n} - \frac{u_{r}^{n} Δt}{2} \frac{ξ_{i + 1}^{n} - ξ_{i}^{n}}{x_{i + 1} - x_{i}} and F_{l}^{n} = ξ_{l}^{n} - \frac{u_{l}^{n} Δt}{2} \frac{ξ_{i}^{n} - ξ_{i - 1}^{n}}{x_{i} - x_{i - 1}} .

(8)

The numerical fluxes

F_{f}

, with

f = r, l

denoting respectively the right and left face, are the sum of the value

ξ_{f}

estimated by third-degree

ν

-spline polynomials (see

[37]

) and the second term representing the

ξ

gradient calculated on the same face

f

. Omitting the time index for simplicity, the expression of

ξ_{f = r} = \frac{- σ^{2}}{4 (σ + 1)} ξ_{i - 1} + \frac{(σ + 2)}{4} ξ_{i} + \frac{(σ + 2)}{4 (σ + 1)} ξ_{i + 1} .

(9)

The parameter

σ = h / k

with

h = x_{i} - x_{i - 1}

and

k = x_{i + 1} - x_{i}

, determines the profile to use considering the variations of the grid. Let choose

σ = 1, 1 / 2, 2

, which define the profiles for equally spaced nodes (the most frequent situation), half of the distance

k

with respect to

h

or double of the distance

k

with respect to

h

(in these last cases

i

results to be the node of a Bridge FV). The three obtained profiles are

ξ_{r} = \frac{- 1}{8} ξ_{i - 1} + \frac{3}{4} ξ_{i} + \frac{3}{8} ξ_{i + 1} for σ = 1

(10)

ξ_{r} = \frac{- 1}{24} ξ_{i - 1} + \frac{5}{8} ξ_{i} + \frac{5}{12} ξ_{i + 1} for σ = \frac{1}{2}

(11)

ξ_{r} = \frac{- 1}{3} ξ_{i - 1} + ξ_{i} + \frac{1}{3} ξ_{i + 1} for σ = 2

(12)

By the substitution of (10), (11) and (12) into (8) the right fluxes are

F_{r} = \frac{- 1}{8} ξ_{i - 1} + (\frac{3}{4} + \frac{c_{f}}{2}) ξ_{i} + (\frac{3}{8} - \frac{c_{f}}{2}) ξ_{i + 1} for σ = 1

(13)

F_{r} = \frac{- 1}{24} ξ_{i - 1} + (\frac{5}{8} + \frac{c_{f}}{2}) ξ_{i} + (\frac{5}{12} - \frac{c_{f}}{2}) ξ_{i + 1} for σ = \frac{1}{2}

(14)

F_{r} = \frac{- 1}{3} ξ_{i - 1} + (1 + \frac{c_{f}}{2}) ξ_{i} + (\frac{1}{3} - \frac{c_{f}}{2}) ξ_{i + 1} for σ = 2

(15)

where

c_{f} = \frac{u_{f} Δt}{x_{i + 1} - x_{i}}

with

u_{f}

denoting the velocity at the considered face. The expressions of

ξ_{l}

and

F_{l}

at the left face are similar the above relations (10) – (15) after decreasing of a unit the indexes of the nodes.

2.2. Normalized Variable Diagram

To control the behaviour of the numerical solution, we analyze the diagram of the normalized variable

{\bar{ξ}}_{i}

, i.e. the variable

ξ

transformed by

\bar{ξ} = \frac{(ξ - ξ_{i - 1})}{(ξ_{i + 1} - ξ_{i - 1})}

(16)

and inversely

ξ = (ξ_{i + 1} - ξ_{i - 1}) \bar{ξ} + ξ_{i - 1} .

(17)

It is easy to verify that the transformed values

{\bar{ξ}}_{i + 1}

ξ_{i + 1}

and

{\bar{ξ}}_{i - 1}

ξ_{i - 1}

are

1

and

0

, respectively. This means that while

F_{r}

depends on both

ξ_{i - 1}

ξ_{i}

and

ξ_{i + 1}

, the normalized variable

{\bar{F}}_{r}

depends only on

{\bar{ξ}}_{i}

. Therefore, a Normalized Variable Diagram (NVD) allows a graphic representation of the functional relationship between the values

{\bar{ξ}}_{i}

and

{\bar{F}}_{r}

[16]

. Since Fastest has to be transformed into unit square by piecewise profiles, three pairs of abscissas

E_{i}

and

G_{i} (i = 1, 2, 3)

, projections of three pairs of points

B_{i}

and

A_{i}

on the

{\overset{̀}{ξ}}_{i}

axis, need to be calculated see Figure 2 as an example for

σ = 1

After some simplifications, the transformed profiles are

{\overset{̀}{F}}_{r} = K_{1} {\overset{̀}{ξ}}_{i} + K_{2} with K_{1} = \frac{3}{4} + \frac{c_{r}}{2} K_{2} = \frac{3}{8} - \frac{c_{r}}{2} for σ = 1

(18)

{\overset{̀}{F}}_{r} = K_{3} {\overset{̀}{ξ}}_{i} + K_{4} with K_{3} = \frac{5}{8} + \frac{c_{r}}{2} K_{4} = \frac{5}{12} - \frac{c_{r}}{2} for σ = \frac{1}{2}

(19)

{\overset{̀}{F}}_{r} = K_{5} {\overset{̀}{ξ}}_{i} + K_{6} with K_{5} = 1 + \frac{c_{r}}{2} K_{6} = \frac{1}{3} - \frac{c_{r}}{2} for σ = 2 .

(20)

The solutions of the six systems

A_{1} = \{\begin{matrix} {\overset{̀}{F}}_{r} = K_{1} {\overset{̀}{ξ}}_{i} + K_{2} \\ {\overset{̀}{F}}_{r} = 1.0 \end{matrix}

and

B_{1} = \{\begin{matrix} {\overset{̀}{F}}_{r} = K_{1} {\overset{̀}{ξ}}_{i} + K_{2} \\ {\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i} / c_{r} \end{matrix} for σ = 1

(21)

A_{2} = \{\begin{matrix} {\overset{̀}{F}}_{r} = K_{3} {\overset{̀}{ξ}}_{i} + K_{4} \\ {\overset{̀}{F}}_{r} = 1.0 \end{matrix}

and

B_{2} = \{\begin{matrix} {\overset{̀}{F}}_{r} = K_{3} {\overset{̀}{ξ}}_{i} + K_{4} \\ {\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i} / c_{r} \end{matrix} for σ = 1 / 2

(22)

A_{3} = \{\begin{matrix} {\overset{̀}{F}}_{r} = K_{5} {\overset{̀}{ξ}}_{i} + K_{6} \\ {\overset{̀}{F}}_{r} = 1.0 \end{matrix}

and

B_{3} = \{\begin{matrix} {\overset{̀}{F}}_{r} = K_{5} {\overset{̀}{ξ}}_{i} + K_{6} \\ {\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i} / c_{r} \end{matrix} for σ = 2

(23)

result in three couples of values

E_{1} = (\frac{K_{2} c_{r}}{1 - K_{1} * c_{r}}, 0.0) G_{1} = (\frac{1 - K_{2}}{K_{1}}, 0.0) for σ = 1

(24)

E_{2} = (\frac{K_{4} c_{r}}{1 - K_{3} * c_{r}}, 0.0) G_{2} = (\frac{1 - K_{4}}{K_{3}}, 0.0) for σ = \frac{1}{2}

(25)

E_{3} = (\frac{K_{6} c_{r}}{1 - K_{5} * c_{r}}, 0.0) G_{3} = (\frac{1 - K_{6}}{K_{5}}, 0.0) for σ = 2

(26)

The analysis of the profiles

{\overset{̀}{F}}_{r}

into the unit square when

{\overset{̀}{F}}_{r} = {\overset{̀}{F}}_{r} ({\overset{̀}{ξ}}_{i}, c_{r}, σ)

, gives for

0 < c_{r} \leq 1

, a set of straight non-parallel lines. In particular, for

σ = 1

and

c_{r} = 0

the transformed flux (18) is reduced into

{\overset{̀}{F}}_{r} = \frac{3}{4} {\overset{̀}{ξ}}_{i} + \frac{3}{8}

which for

{\overset{̀}{ξ}}_{i} = \frac{1}{2}

gives

{\overset{̀}{F}}_{r} = \frac{3}{4}

, that is the point

Q = (\frac{1}{2}, \frac{3}{4})

. Instead for

σ = 1

and

c_{r} = 1

the transformed flux (18) is reduced into

{\overset{̀}{F}}_{r} = \frac{5}{4} {\overset{̀}{ξ}}_{i} - \frac{1}{8}

which for

{\overset{̀}{ξ}}_{i} = 0

and

{\overset{̀}{ξ}}_{i} = 1

Download: Download full-size image

Figure 2. Normalized Variable Diagram of Fastest with

σ = 1

and Courant numbers

c_{r} = 0, 1 / 2, 1 .

results

{\overset{̀}{F}}_{r} = \frac{- 1}{8}

and

{\overset{̀}{F}}_{r} = \frac{9}{8}

respectively, thus by definition values

{\overset{̀}{F}}_{r} = 0

and

{\overset{̀}{F}}_{r} = 1

. Lastly for

{\overset{̀}{ξ}}_{i} = \frac{1}{2}

we obtain

{\overset{̀}{F}}_{r} = \frac{1}{2}

, hence the Upwind profile see Figure 2.

The graphic of Fastest, assigned

σ = 1

and

c_{r} = 0.5

, is shown in Figure 2. It is composed of three segments:

\bar{OB}

delimited by the points O=(0,0) and

B

intersection of the two straight lines

{\overset{̀}{F}}_{r} = \frac{{\overset{̀}{ξ}}_{i}}{c_{r}}

and

{\overset{̀}{F}}_{r} = K_{1} {\overset{̀}{ξ}}_{i} + K_{2}, \bar{BA}

with

A =

intersection of the two straight lines

{\overset{̀}{F}}_{r} = 1

and

{\overset{̀}{F}}_{r} = K_{1} {\overset{̀}{ξ}}_{i} + K_{2}

and finally

\bar{AD}

with D=(1.0,1.0). To better understand what the segment

\bar{OB}

represents (see Section 2.3 and

[16])

Given

ξ_{i}^{n}

the corresponding

F_{r}^{n}

value is obtained by Fastest applying the Algorithm 2.1, both with velocity u constant and with the velocity u varying in time and space, see section 2.3.

Algorithm 2.1.

After fixing the parameters

σ

and

j

(

j = 1

for

σ = 1

j = 2

for

σ = \frac{1}{2}

and

j = 3

for

σ = 2

), five cases are possible:

|ξ_{i + 1}^{n} - ξ_{i - 1}^{n}| \leq ε

and

0 < ε ≪ 1

, the transformation in the normalised plane is not possible therefore

ξ_{i}^{n + 1}

is estimated by means of second-order back-ward Characteristic starting from

x_{i}

, abscissa of the

ith

node. The value on foot

x_{p}

is interpolated with a third-degree polynomial, which is function of four

ξ^{n}

values calculated on the partition

P^{m - 1}

P^{m}

(for further details see Section 5)

{\overset{̀}{F}}_{r} = \frac{{\overset{̀}{ξ}}_{i}}{c_{r}}

for

0.0 \leq {\overset{̀}{ξ}}_{i} < E_{j}

{\overset{̀}{F}}_{r} = K_{l} {\overset{̀}{ξ}}_{i} + K_{l + 1}

for

E_{j} \leq {\overset{̀}{ξ}}_{i} < G_{j}

(with

l = 1

for

j = 1

l = 3

for

j = 2

and

l = 5

for

j = 3

)

{\overset{̀}{F}}_{r} = 1.0

for

G_{j} \leq {\overset{̀}{ξ}}_{i} \leq 1.0

{\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i}

for

{\overset{̀}{ξ}}_{i} <

0 or

{\overset{̀}{ξ}}_{i}

1.0

After computing

{\overset{̀}{F}}_{r}

in the cases 2, 3, 4 and 5 we obtain with inverse transformation

F_{r}^{n} = (ξ_{i + 1}^{n} - ξ_{i - 1}^{n}) {\overset{̀}{F}}_{r} + ξ_{i - 1}^{n} .

(27)

Remark 2.1 In the case 5,

ξ_{i}

can be a local extremum and the solution might present sharp changes due to the geometry of the domain. The use of high order numerical schemes may thus generate spurious oscillations and the no convergence of the solution. In order to comply with the condition established by the Godunov’s theorem in

[39]

, one may reformulate Fastest similarly as in some high order schemes

[40]

. However, a more effective strategy may be adopted by implementing h-adaptivity. Generated a suitable subdomain

\bar{Ω}

in the region where the wiggles manifest

,

it is assigned Upwind to all the faces of the FV belonging to

\bar{Ω} .

Thus, Fastest would accordingly turn into a Total Variation Diminishing method that notoriously avoids spurious oscillations near a jump.

2.3. Velocity Varying in Time and Space

In case the Courant, Friedrichs and Lewy number c=c(x, t) varies in time and space, the local monotonicity is guaranteed in the unit square as long as certain conditions are respected. Since the monotonicity can increase or decrease, it is necessary to examine these conditions with respect to the right and left face in both the cases:

1) Increasing monotonicity:

ξ_{i}^{n + 1} \geq ξ_{i - 1}^{n} and ξ_{i}^{n + 1} \leq ξ_{i + 1}^{n}

(28)

at the right and the left face respectively

2) Decreasing monotonicity:

ξ_{i}^{n + 1} \leq ξ_{i - 1}^{n} and ξ_{i}^{n + 1} \geq ξ_{i + 1}^{n}

(29)

at the right and the left face respectively.

Firstly, consider the condition (28) at the right face

ξ_{i}^{n + 1} \geq ξ_{i - 1}^{n}

and the inequality chain

ξ_{i - 1}^{n} \leq F_{l}^{n} \leq ξ_{i}^{n} \leq F_{r}^{n} \leq ξ_{i + 1}^{n}

. By applying the expression (7) to

ξ_{i}^{n + 1}

, we obtain

ξ_{i}^{n} - {c_{r} F}_{r}^{n} + c_{l} F_{l}^{n} \geq ξ_{i - 1}^{n}

. Then by reordering

c_{r} F_{r}^{n} \leq ξ_{i}^{n} - ξ_{i - 1}^{n} + c_{l} F_{l}^{n}

and applying to

F_{l}^{n}

the worst estimation

ξ_{i - 1}^{n}

, we obtain

c_{r} F_{r}^{n} \leq ξ_{i}^{n} - ξ_{i - 1}^{n} + c_{l} ξ_{i - 1}^{n}

. Thus, considering

F_{r} \leq \frac{ξ_{i}^{n} + ξ_{i - 1}^{n} (c_{l} - 1)}{c_{r}}

, the condition required to be satisfied in the NVD plane is

{\overset{̀}{F}}_{r} \leq \frac{{\overset{̀}{ξ}}_{i}}{c_{r}} .

(30)

For the increasing monotonicity at the left face,

ξ_{i}^{n + 1} \leq ξ_{i + 1}^{n}

is the condition to consider. By applying (7) to

ξ_{i}^{n + 1}

, we have

ξ_{i}^{n} - {c_{r} F}_{r}^{n} + {c_{l} F}_{l}^{n} \leq ξ_{i + 1}^{n}

. Reordering

ξ_{i}^{n} - ξ_{i + 1}^{n} + {c_{l} F}_{l}^{n} \leq {c_{r} F}_{r}^{n}

and applying the worst estimation

ξ_{i}^{n}

F_{l}^{n}

, we obtain

ξ_{i}^{n} - ξ_{i + 1}^{n} + {c_{l} ξ}_{i}^{n} \leq {c_{r} F}_{r}^{n}

. The inequality in NVD plane becomes:

{\overset{̀}{F}}_{r} \geq \frac{{\overset{̀}{ξ}}_{i} - 1}{c_{r}} + \frac{c_{l} {\overset{̀}{ξ}}_{i}}{c_{r}}

. Since

\frac{{\overset{̀}{ξ}}_{i} - 1}{c_{r}} \leq 0

for

0 \leq {\overset{̀}{ξ}}_{i} \leq 1

, the condition ultimately achieved is

{\overset{̀}{F}}_{r} \geq {\overset{̀}{ξ}}_{i} \frac{c_{l}}{c_{r}} .

(31)

Similarly, for decreasing monotonicity at the right and left face, the conditions to be satisfied in the NVD plane are, respectively

{\overset{̀}{F}}_{r} \geq \frac{{\overset{̀}{ξ}}_{i}}{c_{r}}

(32)

and

{\overset{̀}{F}}_{r} \leq {\overset{̀}{ξ}}_{i} \frac{c_{l}}{c_{r}}

. (33)

Thus, the four limitations are: for increasing monotonicity

{\overset{̀}{F}}_{r} \leq \frac{{\overset{̀}{ξ}}_{i}}{c_{r}}

and

{\overset{̀}{F}}_{r} \geq {\overset{̀}{ξ}}_{i} \frac{c_{l}}{c_{r}}

at the right and left face, respectively; for decreasing monotonicity

{\overset{̀}{F}}_{r} \geq \frac{{\overset{̀}{ξ}}_{i}}{c_{r}}

and

{\overset{̀}{F}}_{r} \leq {\overset{̀}{ξ}}_{i} \frac{c_{l}}{c_{r}}

at the right and left face, respectively. Therefore, if

{\overset{̀}{F}}_{r} = \frac{{\overset{̀}{ξ}}_{i}}{c_{r}}

and

{\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i} \frac{c_{l}}{c_{r}}

all the above conditions are satisfied. Since the increasing and decreasing monotonicity have the same representation in the unit square, a single analysis is required satisfying the aforementioned necessary conditions.

The above condition

{\overset{̀}{F}}_{r} = \frac{{\overset{̀}{ξ}}_{i}}{c_{r}}

is used to define the flux with constant velocity u in the case 2 of Algorithm 2.1 (see

[16])

. Additionally, an example of flux

{\overset{̀}{F}}_{r}

is shown in a diagrammatic form by the segment

\bar{OB}

in Figure 2.

The second condition

{\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i} \frac{c_{l}}{c_{r}}

splits in three cases

{\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i} \frac{c_{l}}{c_{r}}

and

\frac{c_{l}}{c_{r}} < 1

. Since

{\overset{̀}{F}}_{r}

is less than

{\overset{̀}{ξ}}_{i}

{\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i}

is assigned by Upwind

{\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i} \frac{c_{l}}{c_{r}}

and

\frac{c_{l}}{c_{r}} = 1

. Since

{\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i}

, the value is again assigned by Upwind

{\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i} \frac{c_{l}}{c_{r}}

and

\frac{c_{l}}{c_{r}} > 1

with

\frac{c_{l}}{c_{r}} \leq \frac{1}{{\overset{̀}{ξ}}_{i}}

. The condition

{\overset{̀}{F}}_{r} = {\overset{̀}{ξ}}_{i} \frac{c_{l}}{c_{r}}

is used, similarly to case 2 of Algorithm 2.1, to define the flux with velocity u varying in time and space. To note that

{\overset{̀}{ξ}}_{i} \frac{c_{l}}{c_{r}} < \frac{{\overset{̀}{ξ}}_{i}}{c_{r}}

. Even the flux

{\overset{̀}{F}}_{r}

can be graphically shown similarly to the segment

\bar{OB}

in Figure 2.

3. Consistency, Stability and Convergence of Fastest

3.1. Consistency

Given the one-way wave equation:

\frac{\partial ξ}{\partial t} + u \frac{\partial ξ}{\partial x} = 0

(34)

The corresponding partial differential operator

P

P = \frac{\partial}{\partial t} + u \frac{\partial}{\partial x} .

(35)

For the Fastest scheme, the corresponding approximation operator

P (Δt, Δx, σ)

P (Δt, Δx, σ) ξ = ξ_{i}^{n + 1} - ξ_{i}^{n} + \frac{u_{i}^{n} Δt}{D_{i}} [F_{r}^{n} - F_{l}^{n}],

(36)

where

ξ_{i}^{n} = ξ (nΔt, iΔx)

and index

i

identifies the FV where (36) is applied.

Definition 3.1

Let suppose the function ξ(t,x) suitably differentiable with respect to time and space, then the scheme (36) is consistent with the partial differential equation (34) if the difference

Pξ - P (Δt, h, σ) ξ \to 0 as (Δt, h) \to 0

Theorem 3.1

The consistency of the Fastest method with the equation (34) is demonstrated by using the definition 3.1.

Proof

Since the expressions of

F_{r}^{n}

and

F_{l}^{n}

depend on parameter

σ

(

σ = 1, 1 / 2, 2

), when

σ = 1 / 2

or σ =

2

both the Bridge

{FV}_{i}

and its next

{FV}_{i + 1}

need to define a suitable flux

F_{f}

. Therefore, the cases to examine are five. For simplicity, the time index of the fluxes is neglected.

σ = 1

(FV with equal diameter)

The difference

F_{r}^{n} - F_{l}^{n} = (\frac{- 1}{8} ξ_{i - 1} + \frac{3}{4} ξ_{i} + \frac{3}{8} ξ_{i + 1} - \frac{uΔt}{2} \frac{ξ_{i + 1} - ξ_{i}}{x_{i + 1} - x_{i}})

(\frac{- 1}{8} ξ_{i - 2} + \frac{3}{4} ξ_{i - 1} + \frac{3}{8} ξ_{i} - \frac{uΔt}{2} \frac{ξ_{i} - ξ_{i - 1}}{x_{i} - x_{i - 1}})

(37)

can be reordered as

F_{r}^{n} - F_{l}^{n} = \frac{1}{8} ξ_{i - 2} + (\frac{- 7}{8} - \frac{uΔt}{2 h}) ξ_{i - 1} + (\frac{3}{8} + \frac{uΔt}{h}) ξ_{i} + (\frac{3}{8} - \frac{uΔt}{2 h}) ξ_{i + 1},

(38)

By replacing suitable Taylor series in the scheme (36) with

D_{i} = h

and

(x_{i} - x_{i - 1}) = (x_{i + 1} - x_{i}) = h

, the difference

Pξ - P (Δt, Δx, σ) ξ = \frac{- {Δt}^{2}}{6} ξ_{ttt} - \frac{h^{2}}{24} u ξ_{xxx} + O ({Δt}^{3}) + O (h^{3}) \to 0 as (Δt, Δh) \to 0 .

(39)

2 . σ = 1 / 2

(Bridge FV with grid of the next subdomain

Ω_{j + 1}

halved)

The fluxes difference reordered is

F_{r}^{n} - F_{l}^{n} = \frac{1}{8} ξ_{i - 2} + (\frac{- 19}{24} - \frac{uΔt}{2 h}) ξ_{i - 1} + (\frac{1}{4} + \frac{3 uΔt}{2 h}) ξ_{i} + (\frac{5}{12} - \frac{uΔt}{h}) ξ_{i + 1}

(40)

By replacing suitable Taylor series in the scheme (36) with

D_{i} = \frac{3}{4} h

x_{i + 1} - x_{i} = \frac{h}{2}

and

x_{i} - x_{i - 1} = h

, the difference

Pξ - P (Δt, Δx, σ) ξ = \frac{h}{8} {uξ}_{xx} - \frac{{Δt}^{2}}{6} ξ_{ttt}

\frac{{5 h}^{2}}{144} u ξ_{xxx} - \frac{h}{12} u^{2} Δt ξ_{xxx} + O ({Δt}^{3}) + O (h^{3})

→0

as (Δt, Δh) \to 0

(41)

3 . σ = 1 / 2

(First FV of the next subdomain

Ω_{j + 1}

with halfed grid spacing)

The fluxes difference reordered is

F_{r}^{n} - F_{l}^{n} = \frac{1}{24} ξ_{i - 2} + (\frac{- 3}{4} - \frac{uΔt}{h}) ξ_{i - 1} + (\frac{1}{3} + \frac{2 uΔt}{h}) ξ_{i} + (\frac{3}{8} - \frac{uΔt}{h}) ξ_{i + 1}

(42)

By replacing suitable Taylor series in the scheme (36) with

D_{i} = \frac{h}{2}

x_{i + 1} - x_{i} = \frac{h}{2}

and

x_{i} - x_{i - 1} = \frac{h}{2}

, the difference

Pξ - P (Δt, Δx, σ) ξ = \frac{- {Δt}^{2}}{6} ξ_{ttt} - \frac{h^{3}}{64} u ξ_{xxxx} + O ({Δt}^{3}) + O (h^{4}) \to 0 as (Δt, h) \to 0 .

(43)

4 . σ = 2

(Bridge FV and next subdomain

Ω_{j + 1}

with doubled grid spacing)

The fluxes difference reordered is

F_{r}^{n} - F_{l}^{n} = \frac{1}{8} ξ_{i - 2} + (\frac{- 13}{12} - \frac{uΔt}{2 h}) ξ_{i - 1} + (\frac{5}{8} + \frac{3 uΔt}{4 h}) ξ_{i} + (\frac{1}{3} - \frac{uΔt}{4 h}) ξ_{i + 1}

(44)

By replacing suitable Taylor series in the scheme (36) with

D_{i} = \frac{3}{2} h

x_{i + 1} - x_{i} = 2 h

and

x_{i} - x_{i - 1} = h

, the difference

Pξ - P (Δt, Δx, σ) ξ = \frac{- {Δt}^{2}}{6} ξ_{ttt} - \frac{h}{4} u ξ_{xx} + \frac{h}{6} u^{2} Δt ξ_{xxx} - \frac{{11 h}^{2}}{36} u ξ_{xxx} + O ({Δt}^{3}) + O (h^{3}) \to 0 as (Δt, h) \to 0 .

(45)

5 . σ = 2

(First FV of the next subdomain

Ω_{j + 1}

with grid double spacing)

The fluxes difference reordered is

F_{r}^{n} - F_{l}^{n} = \frac{1}{3} ξ_{i - 2} + (\frac{- 9}{8} - \frac{uΔt}{4 h}) ξ_{i - 1} + (\frac{5}{12} + \frac{uΔt}{2 h}) ξ_{i} + (\frac{3}{8} - \frac{uΔt}{4 h}) ξ_{i + 1}

(46)

By replacing suitable Taylor series in the scheme (36) with

D_{i} = 2 h

x_{i + 1} - x_{i} = 2 h

and

x_{i} - x_{i - 1} = 2 h

the difference

Pξ - P (Δt, Δx, σ) ξ = \frac{- {Δt}^{2}}{6} ξ_{ttt} - \frac{h^{2}}{4} u ξ_{xxx} + O ({Δt}^{3}) + O (h^{3}) \to 0 as (Δt, h) \to 0 .

(47)

In conclusion, the method Fastest is consistent with the partial differential equation (34) with

σ = 1

and

σ = \frac{1}{2}, 2

for Bridge FV and first FV of

Ω_{j + 1},

respectively.

3.2. Order of Accuracy of the Schemes

Given the discretized problem

P_{Δt, h} ξ = R_{Δt, h} f

, where the expressions

P_{Δt, h} ξ

and

R_{Δt, h} f

are evaluated at the grid point

(t^{n}, x_{i})

, a consistent scheme for the partial differential equation

Pξ = f

is accurate of order

p

in time and order

q

in space if for any smooth function

φ (t, x)

P_{Δt, h} φ - R_{Δt, h} Pφ = O ({Δt}^{p}) + O (h^{q})

. (48)

Definition 3.2

The difference

P_{Δt, h} φ - R_{Δt, h} Pφ

defines the truncation error of the scheme

[41]

Thus, the accuracy of the above schemes are respectively:

1 . O ({Δt}^{2}) + O (h^{2})

, 2.

O ({Δt}^{2}) + O (h)

, 3.

O ({Δt}^{2}) + O (h^{3}),

O ({Δt}^{2}) + O (h),

O ({Δt}^{2}) + O (h^{2})

(49)

3.3. Stability for Initial Value Problem and σ=1

Consider the one-way wave equation (33) on the interval

[0, 1]

, where the solution satisfies the boundary conditions

ξ (t, 0) = ξ (t, 1)

for all nonnegative values of

t

and the initial condition given by periodic data. The application of Fourier analysis to the equation (33) shows that the initial value problem is well-posed

[41]

, therefore both the continuity condition

{‖ξ (t, •)‖}_{2} \leq C * {‖ξ (0, •)‖}_{2}

(50)

with

C

= constant independent from

ξ,

and the robustness property are satisfied. The definitions of

L^{2}

norm for functions and for grid functions are, respectively

{‖φ‖}_{2} = {(\int_{- \infty}^{+ \infty} {|φ (x)|}^{2} dx)}^{\frac{1}{2}}

(51)

{‖ψ‖}_{h} = {(h \sum_{m = - \infty}^{+ \infty} {|ψ_{m}|}^{2})}^{\frac{1}{2}} .

(52)

Theorem 3.2

In the case

σ = 1

, Fastest is conditionally stable and the stability region is the interval

0 < c \leq \bar{c}, with \bar{c} = 0.78

Proof

Let the interval

[0, 1]

be partitioned in

M

equal sub-intervals and consider the discrete Fourier transform of the numerical solution

{\overset{̀}{ξ}}_{l}^{n + 1} = \sum_{i = 0}^{M - 1} e^{- 2 \frac{πjli}{M}} ξ_{i}^{n + 1}

(53)

with

i

and

l = 0, \dots, M - 1

, where

j

is the imaginary unit. Its inverse is

ξ_{i}^{n + 1} = \frac{1}{M} \sum_{l = 0}^{M - 1} e^{2 \frac{πjil}{M}} {\overset{̀}{ξ}}_{l}^{n + 1} .

(54)

Let apply the transform (54) to scheme (36), which for

σ = 1

ξ_{i}^{n + 1} = \frac{- c}{8} ξ_{i - 2}^{n} + (\frac{7 c}{8} + \frac{c^{2}}{2}) ξ_{i - 1}^{n} + (1 - \frac{3 c}{8} - c^{2}) ξ_{i}^{n} + (\frac{- 3 c}{8} + \frac{c^{3}}{2}) ξ_{i + 1}^{n}

(55)

Download: Download full-size image

Figure 3. Graphic of surface

\sqrt{{|G (Θ, c)|}^{2}}

for

Θ \in (0, 2 π]

and

c \in (0, 1] .

Download: Download full-size image

Figure 4. Graphical representation of the section

S = S (\bar{c}, hs)

where

\bar{c}

is the value c whose corresponding section

S

has the maximum height

hs \leq 1

obtaining

ξ_{i}^{n + 1}

\frac{1}{M} [\frac{- c}{8} \sum_{l = 0}^{M - 1} e^{\frac{2 πj (i - 2) l}{M}} {\overset{̀}{ξ}}_{l}^{n} + (\frac{7 c}{8} + \frac{c^{2}}{2}) \sum_{l = 0}^{M - 1} e^{\frac{2 πj (i - 1) l}{M}} {\overset{̀}{ξ}}_{l}^{n} + (1 - \frac{3 c}{8} - c^{2}) \sum_{l = 0}^{M - 1} e^{\frac{2 πjil}{M}} {\overset{̀}{ξ}}_{l}^{n} + (\frac{- 3 c}{8} + \frac{c^{2}}{2}) \sum_{l = 0}^{M - 1} e^{\frac{2 πj (i + 1) l}{M}} {\overset{̀}{ξ}}_{l}^{n} ⦌

(56)

thus from (54) and (56), the equality

\frac{1}{M} \sum_{l = 0}^{M - 1} e^{\frac{2 πjil}{M}} {\overset{̀}{ξ}}_{l}^{n + 1} = \frac{1}{M} \sum_{l = 0}^{M - 1} e^{\frac{2 πjil}{M}} {\overset{̀}{ξ}}_{l}^{n} ⦋ \frac{- c}{8} e^{\frac{- 4 πjl}{M}} + (\frac{7 c}{8} + \frac{c^{2}}{2}) e^{\frac{- 2 πjl}{M}} + (1 - \frac{3 c}{8} - c^{2}) + (\frac{- 3 c}{8} + \frac{c^{2}}{2}) e^{\frac{2 πjl}{M}} ⦌

(57)

By applying a term-by-term equality and defining the phase angle

θ = 2 π l / M

, we obtain

{\overset{̀}{ξ}}_{l}^{n + 1} = ⦋ \frac{- c}{8} e^{- 2 jΘ} + (\frac{7 c}{8} + \frac{c^{2}}{2}) e^{- jΘ}

(1 - \frac{3 c}{8} - c^{2}) + (\frac{- 3 c}{8} + \frac{c^{3}}{2}) e^{jΘ}

⦌

{\overset{̀}{ξ}}_{l}^{n}

(58)

Thereby

{\overset{̀}{ξ}}_{l}^{n + 1} = {G (Θ, c) \overset{̀}{ξ}}_{l}^{n}

(59)

where

G (Θ, c)

is the amplification factor.

By iteratively applying the (59),

{\overset{̀}{ξ}}_{l}^{n + 1} = {[G (Θ, c)]}^{n} {\overset{̀}{ξ}}_{l}^{0}

, where

{\overset{̀}{ξ}}_{l}^{0}

is the transformed initial condition, consequently the behaviour in time of the solution

\overset{̀}{ξ}

is controlled by the initial condition and the powers of the amplification factor.

Therefore

\sqrt{{|G (Θ, c)|}^{2}} \leq 1

with

0 < Θ \leq 2 π

and

0 < c \leq 1

(60)

is the condition to respect in order to ensure the stability. The above restriction is the well-known CFL condition. The amplification factor, for

σ = 1,

can be rewritten as

\sqrt{{|G (Θ, c)|}^{2}}

= {

[{1 - \frac{c}{4} {(1 - cosΘ)}^{2} - c^{2} (1 - cosΘ)]}^{2}

{[\frac{c}{4} sinΘ (cosΘ - 5)]}^{2}

}^1/2(61)

In order to determine the interval of

c

values that respect the condition (60), one can analyze the graphic of surface (61) and determine the section

S = S (\bar{c}, hs)

, where

\bar{c}

is the value whose corresponding section

S

has the maximum height

hs \leq 1

, Figure 3 and Figure 4. It is possible to conclude that, for

σ = 1

the stability region is the interval

0 < c \leq

0.78 and that Fastest is conditionally stable.

Furthermore, with varying velocity u, the stability condition is

Δt \leq \frac{\bar{c} h}{\sup_{x \in R, t > 0} |u (x, t)|}

and in the case that even the net size varies, the stability condition is

Δt \leq \min_{l} \frac{\bar{c} h_{l}}{\sup_{x \in (x_{l + 1} - x_{l}), t > 0} |u (x, t)|}

(62)

with

h_{l} = x_{l + 1} - x_{l} .

3.4. Stability for Initial-boundary Value Problem with σ=1

To obtain a solution of the initial-boundary value problems such as the one-way wave equation (34), the inflow boundary condition required by the partial differential equation need to be used in addition to the initial data. Moreover, to avoid the solution is overdetermined no data have to be assigned at the outflow. However many schemes, including Fastest, require additional numerical boundary conditions, in order to determine the solution uniquely. Specifically, the boundary conditions assigned to the problem (36) are the solution

ξ

x = 0

and a quasi-characteristic extrapolation at outflow, respectively.

Theorem 3.3

Let be formulated the following assumptions (see

[41]

1) The differential equation and the initial data are extended to the entire

x

axis and the initial value problem results to be well-posed and stable in case of a difference scheme.

2) The solution

ξ

can be expressed as

φ + ξ^{'}

allowing to obtain an initial-boundary value problem for

ξ^{'}

Ω

similar to the original problem for

ξ

except for

f

and

ξ_{0}

equal to zero.

3) The time interval is extended from

(0, + \infty)

(- \infty, + \infty)

so that the Laplace transform can be used in the analysis of the boundary conditions.

4) Well-posedness of boundary conditions of a partial differential equation is essentially a local property and consequently, only the differential equation and boundary condition on

\partial Ω

need to be considered.

Then, the Fastest solution to the boundary value problem is bounded by the norm

|β|

of the data on the boundary multiplied for a suitable constant

\overset{̀}{C}

{‖ξ‖}_{η, h, k}^{2} \leq \overset{̀}{C} * {|β|}_{η, h, k}^{2}

Proof

Let define the Laplace transform of a discrete

ξ_{i}^{n}

function

\overset{̀}{ξ} (s) = \frac{1}{\sqrt{2 π}} k \sum_{n = - \infty}^{+ \infty} e^{- snk} ξ^{n} or \overset{̀}{ξ} (z) = \frac{1}{\sqrt{2 π}} \sum_{n = - \infty}^{+ \infty} z^{- n} ξ^{n}

(63)

with

k = Δt

n =

time index,

i =

spatial index,

s = η + jτ

with

j = \sqrt{- 1}

τ =

dual variable,

η > 0

specifies that we are considering the initial-boundary problem for

t

in positive direction,

z = e^{sk}

. Since the functions

\overset{̀}{ξ} (s)

and

\overset{̀}{ξ} (z)

are respectively analytic functions of the complex variables

s

and

z

, the inversion formula is

ξ^{n} = \frac{1}{\sqrt{2 π}} \int_{\frac{- π}{k}}^{\frac{π}{k}} e^{snk} \overset{̀}{ξ} (s) dτ

ξ^{n} = \frac{1}{\sqrt{2 π}} \int_{\frac{- π}{k}}^{\frac{π}{k}} z^{n} \overset{̀}{ξ} (s) dτ .

(64)

The Fastest integration formula for

σ = 1

ξ_{i}^{n + 1} = {αξ}_{i - 2}^{n} + β ξ_{i - 1}^{n} + γ ξ_{i}^{n} + δ ξ_{i + 1}^{n},

(65)

Where

α = \frac{- c}{8}

β = (\frac{7 c}{8} + \frac{c^{2}}{2})

γ = (1 - \frac{3 c}{8} {- c}^{2})

δ = (\frac{- 3 c}{8} + \frac{c^{2}}{2})

Applying the inverse transformation to (65), one gets

\frac{1}{\sqrt{2 π}} \int_{\frac{- π}{k}}^{\frac{π}{k}} z^{n + 1} {\overset{̀}{ξ}}_{i} (s) dτ = \frac{1}{\sqrt{2 π}} \int_{\frac{- π}{k}}^{\frac{π}{k}} z^{n} (α {\overset{̀}{ξ}}_{i - 2} (s) + β {\overset{̀}{ξ}}_{i - 1} (s) + γ {\overset{̀}{ξ}}_{i} (s) + δ {\overset{̀}{ξ}}_{i + 1} (s)) dτ .

(66)

The resolvent equation obtained by executing a term-by-term equality and simplifying for

z^{n},

{\overset{̀}{ξ}}_{i}

α {\overset{̀}{ξ}}_{i - 2} + β {\overset{̀}{ξ}}_{i - 1} + γ {\overset{̀}{ξ}}_{i} + δ {\overset{̀}{ξ}}_{i + 1}

.(67)

In order to find solutions to the resolvent equation that are function of

z

L^{2}

{\overset{̀}{ξ}}_{i}

is replaced by

K^{i} (z)

, with

K (z)

analytic and continuous function of

z

, for

i \geq 0

z K^{i} = α K^{i - 2} + β K^{i - 1} + γ K^{i} + δ K^{i + 1} .

(68)

Lemma 3.1

The general integral of the resolvent equation (67) function of z in

L^{2}

has the form

{\overset{̀}{ξ}}_{i} = A (z) K^{i} (z),

with

A (z)

real function and

|K (z)| < 1

for

|z| > 1

Proof

The right term of equation (68) represents the third degree characteristic polynomial in K(z) with four real coefficients, when posed equal to zero it has three roots and consequently the general solution of the resolvent equation (67) is the linear combination of the three roots

K_{1} (z)

K_{2} (z)

K_{3} (z)

with three coefficient functions

A_{1} (z)

A_{2} (z)

A_{3} (z)

\overset{̀}{ξ} = A_{1} (z) K_{1} (z) + A_{2} (z) K_{2} (z) + A_{3} (z) K_{3} (z)

(69)

The analysis of the three roots by applying the algorithm of Jury reported in

[42]

and the polynomial zeros theorem of Marden given in

[43]

, here omitted for sake of brevity, leads to the following conclusion: since just the solutions of the resolvent equation that belonging to

L^{2}

when

|z| > 1

are taken into account, there is only a real root

K (z)

satisfying

|K (z)| < 1

for

0 < c < 1

. Hence, the general integral reduce to the particular form

{\overset{̀}{ξ}}_{i} = A (z) K^{i} (z),

where

|K (z)| < 1

for

|z| > 1

.(70)

An extrapolation formula approximates the boundary condition at inflow

ξ_{i - 2}^{n} = 2 ξ_{i - 1}^{n} - ξ_{i}^{n} .

(71)

Replacing the expression (71) in the generic left flux

F_{l} = \frac{- 1}{8} ξ_{i - 2} + (\frac{3}{4} + \frac{c}{2}) ξ_{i - 1} + (\frac{3}{8} - \frac{c}{2}) ξ_{i}

(72)

one obtains for i =1

F_{l} = \frac{1}{2} (1 + c) ξ_{0} + \frac{1}{2} (1 - c) ξ_{1}

(73)

where

ξ_{0}^{n} = D^{n}

is the Dirichlet value.

By integrating the convection equation in the first FV

ξ_{1}^{n + 1} = (\frac{5 c}{8} + \frac{c^{2}}{2}) ξ_{0}^{n} + (1 - \frac{c}{4} - c^{2}) ξ_{1}^{n} + (\frac{- 3}{8} c + \frac{c^{2}}{2}) ξ_{2}^{n} .

(74)

Finally the boundary condition to analyze is

ξ_{0}^{n + 1} = a ξ_{0}^{n} + b ξ_{1}^{n} + {β (z)}^{n + 1},

(75)

where

a = \frac{5 c}{8} + \frac{c^{2}}{2}

b = 1 - \frac{c}{4} - c^{2}

and

β {(z)}^{n + 1}

results by the subtracting solutions so that the initial function is zero.

The coefficient

A (z)

is determined by transforming the boundary function

{\overset{̀}{β} (z)}^{n + 1}

Considering the equation (75)

ξ_{0}^{n + 1} = \frac{1}{\sqrt{2 π}} \int_{\frac{- π}{k}}^{\frac{π}{k}} z^{n + 1} ({\overset{̀}{ξ}}_{0} (z) - \overset{̀}{β} (z)) dτ = \frac{1}{\sqrt{2 π}} \int_{\frac{- π}{k}}^{\frac{π}{k}} z^{n + 1} \frac{(a {\overset{̀}{ξ}}_{0} (z) - b {\overset{̀}{ξ}}_{1} (z))}{z} dτ

(76)

that is the equation

{\overset{̀}{ξ}}_{0} (z) = z^{- 1} (a {\overset{̀}{ξ}}_{0} (z) - b {\overset{̀}{ξ}}_{1} (z)) + \overset{̀}{β} (z) .

(77)

Recalling (70)

, {\overset{̀}{ξ}}_{0} = A (z)

and

{\overset{̀}{ξ}}_{1} = A (z) K^{1} (z),

then from (77)

zA (z) = aA (z) + bA (z) K (z) + z \overset{̀}{β} (z)

(78)

A (z) [z - a - bK (z)] = z \overset{̀}{β} (z) .

(79)

The norm of the solution

{\overset{̀}{ξ}}_{i} \in L^{2}

{‖\overset{̀}{ξ} (z)‖}^{2} = h \sum_{i = 0}^{\infty} {|{\overset{̀}{ξ}}_{i} (z)|}^{2} = h {|A (z)|}^{2} \sum_{i = 0}^{\infty} {K (z)}^{2 i} = {|A (z)|}^{2} \frac{h}{1 - {|K (z)|}^{2}}

.(80)

With regard to the function

ξ_{i}^{n}

, by the Parseval’s relations, the norm is

{‖ξ‖}_{η, h, k}^{2} = \int_{\frac{- π}{h}}^{\frac{π}{h}} {|A (e^{sk})|}^{2} \frac{h}{1 - {|K (e^{sk})|}^{2}} dτ,

(81)

where

s = η + jτ

. Then, by substitution of

A (e^{sk})

with (79)

{‖ξ‖}_{η, h, k}^{2} = \int_{\frac{- π}{h}}^{\frac{π}{h}} \frac{{{|z|}^{2} |β (e^{sk})|}^{2}}{{|e^{sk} - a - bK (e^{sk})|}^{2}} \frac{h}{1 - {|K (e^{sk})|}^{2}} dτ .

(82)

Note that lower bound is necessary on

|e^{sk} - a - bK (e^{sk})|

Now, the analysis consists to verify that nontrivial solutions of form (78) solve the homogeneous boundary condition, therefore one needs to check whether there is a

K (z)

satisfying (67) such that

A (z) [z - a - bK (z)] = 0 .

(83)

First, let consider the case

|z| > 1

and

|K| < 1 .

By writing

|z| = 1 + ε

and

|K| = 1 - ρ

, where both

ε

and

ρ

are greater than 0 and less than 1, set a constant

\bar{c}

, with

0 < \bar{c} ≪ 1

, for the values

c

with

\bar{c} < c <

0.78, the first factor of expression (82) gives

\frac{{|z|}^{2}}{{|z - a - bK (z)|}^{2}} = \frac{{(1 + ε)}^{2}}{{|(1 + ε) - (\frac{5}{8} c + \frac{c^{2}}{2}) - (1 - \frac{c}{4} - c^{2}) (1 - ρ)|}^{2}} = {\overset{̀}{c}}_{1},

where the ratio

{\overset{̀}{c}}_{1}

> 0. (84)

Then, the second factor of the expression (82)

\frac{{|β (e^{sk})|}^{2} h}{1 - {|K (e^{sk})|}^{2}} = \frac{{|β (z)|}^{2} h}{1 - {|K (z)|}^{2}} and \frac{h}{1 - {|K (z)|}^{2}} \leq \frac{h}{1 - |K (z)|}

with

|K| < 1

.(85)

Since

K (z)

is a root of the equation (68) with

|K| = 1

for

|z| = 1

, according to Strikwerda (see Lemma 11.3.2, in

[41]

) there exists a constant

\overset{̀}{c}

such that

||K (z)| - 1| > \overset{̀}{c} (|z| - 1) .

(86)

Let set the expansion

z = e^{sk} = e^{jτ} (1 + ηk + O {(ηk)}^{2})

, since

k / h

is constant, it results

1 - |K (z)| \geq \overset{̀}{c} ηk

(87)

and finally

\frac{h}{1 - |K (z)|} \leq \frac{{\overset{̀}{c}}_{2}}{η} .

(88)

To conclude, the dependence of the solution on the data, for

|z| > 1

and

|K| < 1

, is given by

{‖ξ‖}_{η, h, k}^{2} \leq C_{1} * {|β|}_{η, h, k}^{2},

(89)

where the constant

C_{1}

takes into account

{\overset{̀}{c}}_{1}

and

\frac{{\overset{̀}{c}}_{2}}{η}

Second, let consider the case

|z| = 1

and

|K| = 1

Set a constant

\bar{c}

, with 0 <

\bar{c} ≪

1, for the values

c

with

\bar{c} < c <

0.78 the first factor of expression (82) gives

\frac{{|z|}^{2}}{{|z - a - bK (z)|}^{2}} = \frac{1}{{|\frac{- 3}{8} c + \frac{c^{2}}{2}|}^{2}} = {\overset{̀}{c}}_{1} .

(90)

Since the relation (88) is once again applicable to the second factor of (82), it is possible to conclude that the dependence of the solution on the data is given by

{‖ξ‖}_{η, h, k}^{2} \leq \overset{̀}{C} * {|β|}_{η, h, k}^{2},

(91)

where

\overset{̀}{C}

is a constant that takes into account

C_{1}

and

{\overset{̀}{c}}_{1}

The stability of the initial-boundary value problem for the one-way wave equation, here solved with the Fastest method and inflow Dirichlet condition, is demonstrated.

However, in order to determine uniquely the numerical solution, Fastest requires the

ξ^{n}

value at outflow. The extrapolation formula

ξ_{outflow}^{n} = 2 ξ_{N}^{n} - ξ_{N - 1}^{n}

is therefore proposed.

3.5. Stability of Fastest for Initial-boundary Value Problem with h-Adaptivity

The convergence of the numerical solution of transport problems in the context of Domain Decomposition methods is shown in

[9].

A theorem based on the previous study is formulated to prove the stability of Fastest in presence of h-adaptivity.

Theorem 3.4

Let consider the differential equation (34) on a domain

Ω

where the solution

ξ

satisfies the initial-boundary data and let partition

Ω

into N subdomains

Ω_{k}

(k=1,…, N) which FV diameter D_k have different values. In addition, Bridge

{FV}_{k}

(k=1,…, N-1) are placed between each couple of subdomains. Suppose that the stability property of Fastest is enjoyed in each subdomain

Ω_{k}

with regard to the initial or initial-boundary values depending on the

Ω_{k}

in exam.

Therefore, a condition for Fastest stability in

Ω

in presence of h-adaptivity is that the continuity condition

uξ

is accomplished on each interfaces

Γ_{k} .

Note that the continuity condition is used to solve the linear system arising from the discretization of the partial differential equation and not for its solution.

Proof

Theorem 3.4 is demonstrated by a constructive proof, with the assumption

σ = 1 / 2

. It is useful considering the Bridge

{FV}_{k}

divided into two parts, each of them added to the adjacent subdomain.

Let consider that:

1) the first subdomain

Ω_{1}

is partitioned with

n_{1}

equal

{FV}_{i}

i = 1, \dots, n_{1}

and let define the subdomain

Ω_{1}^{˄}

by including the inflow boundary and the left part of the first Bridge FV;

2) the generic subdomain

Ω_{j}

with 1< j < N, (2<N), is partitioned with

n_{j}

equal

{FV}_{i}

i = 1, \dots, n_{j}

, and let define the subdomain

Ω_{j}^{˄}

by including the right part of (j-1)th Bridge FV and the left part of jth Bridge FV;

3) the last subdomain

Ω_{N}

is partitioned with

n_{N}

equal

{FV}_{i}

i = 1, \dots, n_{N}

and let define the subdomain

Ω_{N}^{˄}

by including the right part of (N-1)th Bridge FV and the outflow boundary;

4) the node of each Bridge

{FV}_{k}

has distances from the left and right faces equal to D_k/2 and D_k+1/2, respectively (see Figure 1a and 1b);

5) the Bridge FV nodes coincide with the interfaces

Γ_{k} .

Furthermore considering that:

in each subdomain

Ω_{k}^{˄}

the approximations of the numerical fluxes of the kth Bridge FV can be assimilated to the outflow and inflow conditions related to the

Γ_{k}

interface;

by hypothesis the schemes are stable in each subdomain

Ω_{k}

;

the integration formula relevant to the Bridge FV assigns the

ξ_{h}

value on its node (that is the continuity condition u

ξ_{h}

is fulfilled on the interface

Γ_{k}

by construction).

One can conclude that Fastest is stable even in presence of h-adaptivity.

Note that the partition with two subdomains is not included in the previous definitions, however the condition N=2 is easily intelligible. Moreover, the theorem 3.4 with the assumption

σ = 2

can be demonstrated similarly.

Remark 3.1 The assimilation of the numerical right flux of the kth Bridge FV to a numerical inflow condition for

Ω_{k + 1}^{˄}

subdomain would require a stability study similar to that developed in section 3.4. Anyway, an empirical explanation of the good functioning (i.e. stability) of the method consists in the presence of two continuity conditions, of

uξ

Γ_{k}

and of the flux

F_{l}^{n}

on the left face of first FV of

Ω_{k + 1}^{˄}

The stability analysis of a parabolic initial-boundary value problem is very similar to that of a hyperbolic one, indeed there is no need to check for admissible solutions of the resolvent equation that satisfy the homogeneous boundary condition for

|z| = 1,

see section 3.4.

3.6. Convergence

The convergence of the Fastest method is established on the Lax-Richtmeyer equivalence theorem

[41]

. Fastest is a one-step scheme, consistent for well-posed initial value problem (34), conditionally stable and therefore convergent. Moreover, for smooth initial condition

φ_{0}

, the order of accuracy of the solution is equal to the order of accuracy of the scheme.

Remark 3.2 From (49) the Fastest scheme is second order accurate in time while the accuracy with respect to the space is third order in the aforementioned case 3, second order in cases 1 and 5 and first order in cases 2 and 4 respectively. To achieve a higher space accuracy, the second and third spatial derivatives of the truncation errors can be approximate by Generalized Finite Difference (GFD) formulae (see

[44]

) and then, multiplied for the suitable factors, to be added to the

ξ_{i}^{n}

value. Since the stencils of the GFD formulae are with four nodes, the derivative approximations result to be of second and first order respectively. Therefore, it is possible to obtain a second order in time and third order in space method.

4. Dynamically Locally Self h-Adaptive Method and Pick & Roll Technique

4.1. h-Adaptive FV Methods

An extensive literature exists regarding the h-adaptive FD or FV methods, with or without domain transformations, enabling to improve the accuracy of solutions of fluid-dynamic problems

[45-48]

. However, the numerical solution of partial differential equation problems by high order FV methods, defined on irregular grids, remains still challenging

[23]

. A self h-adaptive method requires the definition of two factors: the so-called error indicators and the techniques that efficiently modify the discretization of the domain Ω. In our adaptive technique the error indicators are calculated on three suitable nodes, respectively two on the left (left flags LF1 and LF2) and one on the right (right flag RF) of each critical region. The differences between the initial and the adjourned numerical solutions are calculated on these nodes and, whenever predetermined values are not respected, a new partition is generated.

4.2. Pick & Roll Technique

In order to save computational cost, the novel partition enables to dynamically halves or doubles the FV diameters according to the belonging to critical or favorable regions. The process, named Pick & Roll, may be considered a movement of FV in accordance with an adequate number of FV moves ideally from the right to the left hand.

Algorithm 4.1.

Suppose to consider the subdomain

Ω_{j}

with FV of diameter D_j and the three flag nodes LF1, LF2 and RF. Define the numerical solution calculated on partition

P^{m}

at time step

n

, as

ξ^{n, m} .

A numerical solution is considered "acceptable" if each of the three flags comply with certain conditions described further on.

ξ^{n, m}

meets the acceptability conditions, then the solution

ξ^{n + 1, m}

at the next time step

n + 1

can be calculated by the same partition

m

used to calculate the accepted solution

ξ^{n, m}

The conditions of unacceptability are:

1 . RF

is external to

Ω_{j}

or, if internal, there are less than twelve FV between it and the boundary of

Ω_{j}

. A new partition

m + 1

is generated by adding four FV with diameter D_j to the right boundary of

Ω_{j}

and by interpolation of the missing values, then a new solution

ξ^{n, m + 1}

is calculated. This second solution is in turn verified with respect to the acceptability and consequently, two events may occur:

the second solution is acceptable and thus the transient can proceed to calculate the solution

ξ^{n + 1, m + 1}

;

the second solution is not acceptable. Then a new partition

m +

2 is created by adding other four FV with diameter D_jto the right boundary of

Ω_{j}

and by interpolation of the missing values. Thus, a new solution

ξ^{n, m + 2}

is calculated. The transient can proceed to calculate the solution

ξ^{n + 1, m + 2}

2 . LF 2

is external to

Ω_{j}

and the solution

ξ^{n, m}

is not acceptable. A new partition

m + 1

is created by adding four FV with diameter D_j to the left boundary of

Ω_{j}

, the missing values are interpolated and a new solution

ξ^{n, m + 1}

calculated. Thus, the transient can proceed to calculate the solution

ξ^{n + 1, m + 1}

to next time step

n + 1

by means of the partition

m + 1

. Since the velocity direction is from the left to the right, this event is not frequent.

3 . LF 1

is internal to

Ω_{j}

and more than twelve FV are present between LF1 and the left boundary of

Ω_{j}

. The solution

ξ^{n, m}

is not acceptable and four FV with diameter D_j are subtracted from the left boundary of

Ω_{j}

, the missing values are interpolated and a new solution

ξ^{n, m + 1}

calculated. The transient can proceed to calculate the solution

ξ^{n + 1, m + 1}

. Considering the direction of the velocity, this event is frequent.

Two 1D hyperbolic problems (i.e. the advancing step and the isolated wave tests) are solved in Section 7 to test the efficiency of new method. Since it is known the location of the singularities as well as the larger variations of the analytical solutions by the initial conditions, it is possible to determine where the domains

Ω

have FV with the smallest diameter. For the two analyzed tests, the domain is partitioned into five subdomains with the

Ω_{3}

central subdomains having FV with the smallest diameter.

Table 1. Diameters. Di of the FV belonging to the subdomains Ω_i, i=1, …, 5 and of Bridge FV (BFVj), j=1,2,3,4.

D1=0.04	D2=0.02	D3=0.01	D4=0.02	D5=0.04
BFV1=0.03	BFV2=0.015	BFV3=0.015	BFV4=0.03

In order to optimize the adaptive process, some of the FV belonging to the previous partition are transformed to limit the increase of the total number of FV.

For example, let’s suppose that after a time interval from the beginning of the transient the error indicator RF lies externally to

Ω_{3}

, consequently the h-adaptive process adds four FV to the right boundary of

Ω_{3} .

The four FV are obtained from the following process: the first FV of

Ω_{5}

is subtracted from

Ω_{5}

(i.e. picked) and added (i.e. rolled) to

Ω_{4}

like two FV of diameter

0.02 .

Then the two first FV of diameter

0.02

Ω_{4}

are picked and rolled to

Ω_{3}

like four FV of diameter

0.01 .

Similarly, if

LF 2

is external to

Ω_{3}

, the last FV of

Ω_{1}

with diameter

0.04

is halved, picked and rolled to

Ω_{2}

like two FV of diameter

0.02

, then the final two FV of

Ω_{2}

are halved, picked and rolled like four FV of diameter

0.01

at the left boundary of

Ω_{3}

. Again, if

LF 1

is internal to

Ω_{3}

and more than twelve FV are present between LF1 and the left boundary of

Ω_{3}

, the first four FV of

Ω_{3}

are picked and rolled to

Ω_{2}

like two FV of diameter

0.02

, then first two FV of

Ω_{2}

are then picked and rolled to

Ω_{1}

like a FV of diameter

0.04

The efficiency of the Pick & Roll process lies in the fact that the number of FV belonging to the subdomains

Ω_{2}

and

Ω_{4}

is constant, the number of FV belonging to

Ω_{5}

decreases, the number of FV belonging to

Ω_{1}

increases while the number of FV of

Ω_{3}

changes in relation to the flags

RF

LF 1

and

LF 2

. In conclusion, a significant increase of the total number of FV derived from each partition

m

is not expected. Numerical tests corroborate this prediction.

4.3. Data Organization

The application of the above mentioned techniques is efficient if an appropriate organization of the data is respected. In particular, the data determined by the initial partition are stored in two matrices as follows: 1) the topological matrix

NOD

contains the first and last node as well as the number of the nodes of each subdomain, the progressive node number of the Bridge FV and lastly the code regarding the grid modification (reduction or enlargement) of the next subdomain. 2) The geometrical matrix

ABSNF

contains, for each subdomain, the abscissa of first and last node as well as the abscissa of the left face of the first FV and of the right face of the last FV, respectively. Moreover,

ABSNF

contains the node abscissa of the Bridge FV, the grid modification factor of the next subdomain and the abscissas of the left and right face of Bridge FV.

The matrices NOD and ABSNF allow to generate the partitions defined by the h-adaptive processes. The data relating to the events examined above are codified in the matrices:

ITNOD 1

when four FV are added to the right of

Ω_{3}

boundary,

ITNOD 2

when four FV are subtracted to the left of

Ω_{3}

boundary,

ITNOD 3

when four FV are added to the left of

Ω_{3}

boundary. The Appendix reports the matrices

NOD

ABSNF

and

ITNOD 1

used in the problems 1 and 2.

5. Interpolations

Third degree polynomials are used in order to calculate the value of ξ on the foot of characteristic as well as ξ and u on the new nodes and faces generate from new partitions. The polynomial coefficients are usually calculated by numerically solving the system

AB = C

where

A

is the matrix of the powers of the abscissas

a_{i}

of the four nodes belonging to the polynomial support,

B

is the vector of the coefficients

b_{i}

of the third-degree polynomial and

C

is the vector of known values c_{i =}

ξ_{i}

u_{i}

. Since the matrices

A

are ill conditioned (i.e. Vandermonde matrices) and in order to optimize the computational time one can solve formally the generic system

AB = C

obtaining sixteen algebraic expressions

φ_{1}, \dots, φ_{16}

, which are functions of the four generic

a_{i}

. By the functions

φ_{j}

it is possible to compute the coefficients

γ_{l, j}

φ_{j}

(

a_{i}

) where l = 1,…,nc (with nc = number of all the possible combinations l of the four nodes belonging to the polynomial support), j=1,…,16 and i=1,…,4. To conclude, the four coefficients

b_{i}

are obtained by linear combinations of c_iand the adequate quadruples of

γ_{l, j}

. For example, let define a local reference coordinate system whose origin lies on the left face of a FV where the foot of a characteristic is placed. The four FV forming the support of the third-degree polynomial are identified by the topological NOD matrix, consequently the local abscissas a_i of the related four nodes by the

COLOCH

matrix. Therefore, by the functions

φ_{j}

it is possible to compute the scalars

γ_{l, j}

of the

COEFCH

(19,16) matrix, in which the first four values refer to first coefficient b₁, the second four values to b₂, etc.

To note that the calculation of the four polynomial coefficients

b_{i}

is reduced to the determination of the suitable row l of

COEFCH

matrix generated before the start of the transient, and to the four linear combinations

b_{1} = {ξ_{1} * γ}_{l, 1} + {ξ_{2} * γ}_{l, 2} + {ξ_{3} * γ}_{l, 3} + {ξ_{4} * γ}_{l, 4}

(92)

b_{2} = {ξ_{1} * γ}_{l, 5} + {ξ_{2} * γ}_{l, 6} + {ξ_{3} * γ}_{l, 7} + {ξ_{4} * γ}_{l, 8}

(93)

b_{3} = {ξ_{1} * γ}_{l, 9} + {ξ_{2} * γ}_{l, 10} + {ξ_{3} * γ}_{l, 11} + {ξ_{4} * γ}_{l, 12}

(94)

b_{4} = {ξ_{1} * γ}_{l, 13} + {ξ_{2} * γ}_{l, 14} + {ξ_{3} * γ}_{l, 15} + {ξ_{4} * γ}_{l, 16}

(95)

The interpolation of missing values

ξ

and u due to an adaptive partition follows the same approach. After generating the appropriate matrices

COLOC

and

COEF (13, 16)

with scalars

δ_{l, j}

l = 1, \dots, 13

and

j = 1, \dots, 16,

the four linear combinations between

δ_{l, j}

and

ξ_{i}

u_{i}

give the four

b_{i}

. Appendix shows

COLOCH

matrix used to solve problems 1 and 2.

6. Generalization of Fastest to 2D and 3D Problems

The extension of Fastest to two and three dimensions, when the FV are generated by structured grids with constant grid spacing, is a straightforward procedure.

6.1. 2D Generalization

The analogous of the formula (7) in two dimensions is

ξ_{i, j}^{n + 1} = ξ_{i, j}^{n} - \frac{u_{i, j}^{n} Δt}{h} [(F_{r}^{n} - F_{l}^{n}) + (F_{t}^{n} - F_{b}^{n})]

(96)

where

h = x_{r} - x_{l} = y_{t} - y_{b}

and with the top (

t

) and bottom (

b

) faces added to the right and left faces of the one dimension case. For example, let consider two dimensional case where the convection is cross the left face, and the velocity components

u

and

v

are arranged as shown in Figure 5.

Download: Download full-size image

Figure 5. Two-dimensional five-node stencil to estimate the left face value using Fastest (velocity component u positive).

The two dimensional flux, in the direction normal to the

middle point of the face l, is approximated by Fastest in

F_{l} = (\frac{- 1}{8} ξ_{i - 2, j} + \frac{3}{4} ξ_{i - 1, j} + \frac{3}{8} ξ_{i, j} - \frac{u_{l} Δt}{2} \frac{ξ_{i, j} - ξ_{i - 1, j}}{x_{i, j} - x_{i - 1, j}}) + \frac{1}{24} (ξ_{i - 1, j + 1} - 2 ξ_{i - 1, j} + ξ_{i - 1, j - 1})

(97)

where the first term is the usual one dimensional profile Fastest and the second term is the transverse curvature (see

[11]

) representing the effects of transverse transport in a two dimensional flow.

Let reorder the five nodes profile (97) in

F_{l} = \frac{- 1}{8} ξ_{i - 2, j} + (\frac{3}{4} + \frac{c_{l}}{2} - \frac{1}{12}) ξ_{i - 1, j} + (\frac{3}{8} - \frac{c_{l}}{2}) ξ_{i, j} + \frac{1}{24} ξ_{i - 1, j + 1} + \frac{1}{24} ξ_{i - 1, j - 1}

(98)

where

c_{l} = \frac{u_{l} Δt}{x_{i, j} - x_{i - 1, j}}

If the velocity components are in opposite directions, a five nodes stencil specular to the previous one with respect to the face itself has to be used. This leads ultimately to use thirteen nodes for all the velocity directions on all four faces.

6.2. 3D Generalization

The three dimensional profile to calculate the convective term on the left control volume face, involves two transverse curvature terms, see Figure 6 for ingoing flux

F_{l} = (\frac{- 1}{8} ξ_{i - 2, j, z} + \frac{3}{4} ξ_{i - 1, j, z} + \frac{3}{8} ξ_{i, j, z} - \frac{u_{l} Δt}{2} \frac{ξ_{i, j, z} - ξ_{i - 1, j, z}}{x_{i, j, z} - x_{i - 1, j, z}}) + \frac{1}{24} (ξ_{i - 1, j + 1, z} - 2 ξ_{i - 1, j, z} + ξ_{i - 1, j - 1, z}) + \frac{1}{24} (ξ_{i - 1, j, z + 1} - 2 ξ_{i - 1, j, z} + ξ_{i - 1, j, z - 1})

(99)

With outgoing velocity

u

from the face, seven nodes stencil specular to that above with respect to the face itself, has to be used. Appropriate rotations and translations of the stencil allow its application to the other faces of the FV and ensue the 3D integration formula

ξ_{i, j, z}^{n + 1} = ξ_{i, j, z}^{n} - \frac{u_{i, j, z}^{n} Δt}{h} [(F_{r}^{n} - F_{l}^{n}) + (F_{t}^{n} - F_{bo}^{n}) + (F_{f}^{n} - F_{ba}^{n})]

(100)

Download: Download full-size image

Figure 6. Three-dimensional seven-node stencil to estimate the left face value using Fastest (velocity component

u

positive).

where

h = x_{r} - x_{l} = y_{t} - y_{bo} = z_{f} - z_{ba}

and where

r

= right,

l

= left,

t

= top,

bo

= bottom,

f

= forth,

ba

= back faces respectively.

As previously stated by Leonard et al., the high approximation of the convective flux normal to the FV face is highly important, even if the accuracy may be pursued adding the curvature terms

[17]

. Thus five and seven nodes are sufficient for a good approximation of the flux in

2 D

and

3 D

problems, respectively. A study developed by Deponti corroborates the same conclusions

[49]

. One-dimensional studies based on

NVD

are not directly applicable to two- and three-dimensional flows, nevertheless numerical experiences have shown how the stabilization of the profile along the main-stream direction prevents numerical oscillations. Therefore, the one-dimensional NVD studies can be appropriately applied even in two- and three-dimensional cases.

Remark 6.1 The generalization of h-adaptive Fastest solutions for 2D and 3D problems may appear to be challenging. However, taking in consideration structured grids, the adaptive partitions for domains both 2D and 3D can be obtained fairly easily using topologic, geometric and translation matrices

[12, 13]

. Finally, the missing values of ξ and v on new nodes and faces can be estimated using parametric transformations of the FV clusters in a unit reference square and defining the related polynomial basis functions.

7. Numerical Tests

Five problems, with known analytical solution, are here solved in order to verify the properties of the new method.

7.1. Advection Equation (Advancing Step)

The first proposed problem concerns the 1D advection equation

\frac{\partial ξ}{\partial t} + u \frac{\partial ξ}{\partial x} = 0

(101)

with velocity

u = 1

, boundary condition of Dirichlet at inflow and numerical boundary condition at outflow. The analytical solution is represented by the "jump" function and the initial condition is obtained from

ξ (x, t) = \{\begin{matrix} 0.5 for 0 \leq x \leq 1.15 + t \\ 0 for x > 1.15 + t \end{matrix}

(102)

Download: Download full-size image

Figure 7. Solution of an advection problem regarding a "jump" function at 4000th step with velocity

u =

1 and

Δt =

0.0005. The data relevant to the subdomains at the final step are Ω₁=[0.02;1.58], Ω₂=[1.61;1.82], Ω₃=[1.835;3.355], Ω₄=[3.37;3.55], Ω₅=[3.58;9.98]. Each subdomain Ω_k=[XL_k; XR_k] has XL_k= abscissa of left face of the firs FV of

Ω_{k}

, XR_k=abscissa of right face of the last FV of Ω_k. Solution of an advection problem regarding a "jump" function at 4000th step with velocity

u =

1 and

Δt =

Ω_{k}

, XR_k=abscissa of right face of the last FV of Ω_k.

Download: Download full-size image

Figure 8. A zoom of the Figure 7 with FV relating to the subdomains

Ω_{3}

Ω_{4}

and the beginning of

Ω_{5}

This classic 1D advection problem is solved to quantify the amount of numerical viscosity of the method in presence of singularities. Let consider the domain

Ω =

[0,10] divided in five non-overlapping subdomains

Ω_{1}, Ω_{2}, Ω_{3}, Ω_{4}, Ω_{5}

interspersed by four Bridge FV. Initially, the subdomain

Ω_{1}

has 18 volumes with diameter 0.04,

Ω_{2}

10 volumes with diameter 0.02,

Ω_{3}

35 volumes with diameter 0.01,

Ω_{4}

10 volumes with diameter 0.02 and

Ω_{5}

210 volumes with diameter 0.04, respectively. Figure 7 and Figure 8 shown the results obtained for a transient of 4000 steps with

Δt =

0.0005, which exhibit a very low numerical viscosity. To note, the new method does not generate wiggles since the local monotonicity is guaranteed by definition itself and dynamically locally self h-adaptive partitions may be applied.

Download: Download full-size image

Figure 9. Solution of an advection problem regarding an "isolated wave" function at first step and 2000th step with velocity u=1 and Δt= 0.0005. The data relevant to the subdomains at the final step are: Ω₁=[0.02; 2.74],Ω₂=[2.77; 297], Ω₃=[2.985; 3.395], Ω₄=[3.41; 3.59], Ω₅=[3.62; 9.98].Each subdomain Ω_k=[XL_k; XR_k] has XL_k= abscissa of left face of the firs FV of

Ω_{k}

{XR}_{k} =

abscissa of right face of the last FV of Ω_k.

7.2. Advection Equation (Isolated Wave)

A numerical experience still concerning the 1D advection equation carried out by Leonard

[11]

, Pennati et al.

[37]

and Cordero et al.

[50]

is here presented. It involves analyzing the accuracy of the numerical solution by considering the analytical solution with a continuously turning gradient and with a single local maximum. The analytical solution is represented by the "isolated wave" function with velocity

u = 1,

boundary conditions of Dirichlet at inflow and numerical boundary condition at outflow. The initial condition is deduced from the function (103)

ξ (x, t) = \{\begin{matrix} \sin^{2} \frac{π (x - 1.06 - ut)}{0.2} for 1.06 + ut \leq x \leq 1.26 + ut \\ 0 elsewhere \end{matrix}

(103)

Let consider the domain

Ω =

[0,10] initially partitioned as in problem 1. The solutions obtained at the first and the 2000th time step, respectively, executing a transient of 2000 steps with

∆ t

= 0.0005, are shown in Figure 9. A zoom of the solution to the FV of subdomains

Ω_{3}

Ω_{4}

and the beginning of

Ω_{5}

, is illustrated in Figure 10. It is thus possible to verify the accuracy of the new method and in particular the small reduction of the wave amplitude (e.g. the maximum error at step 2000 is 0.073), even in presence of a solution with continuous changes in gradient. As a result of the efficiency of the Pick & Roll technique, the numbers of FV belonging to

Ω_{3}

and

Ω_{5}

at the 2000th time step, are respectively 109 and 302.

Download: Download full-size image

Figure 10. A zoom of the Figure 9 with the FV relating to the subdomains

Ω_{2}

Ω_{3}

Ω_{4}

and the beginning of

Ω_{5}

7.3. Steady State Boundary Values Problems

It is considered the numerical solution of a steady state boundary values problem previously presented by Leonard et al.

[17]

, Pennati et al.

[37]

, Cordero et al.

[50]

, expressed by the equation (104)

u \frac{\partial ξ}{\partial x} = Γ \frac{\partial^{2} ξ}{\partial x^{2}} + S (x)

(104)

where the velocity is

u = 1

, the diffusive coefficient

Γ

can assume either the values 0.002 or 0.001, and the source term is

S (x) = \{\begin{matrix} 10 - 50 x, for 0 \leq x < 0.3 \\ 50 x - 20, for 0.3 \leq x < 0.4 \\ 0 for 0, 0.4 \leq x \leq 1 \end{matrix}

(105)

with Dirichlet boundary conditions:

ξ (0) = 0, ξ (1) = 0.1

. It is worth mentioning that the test takes into account the presence of the source term, which is of particular interest when dealing with convection-diffusion problems. It can represent either a real 1D source term or the effects of transverse transport in 2D or 3D flows. To note that most of the classical schemes give accurate solutions only if no source term is present. In presence of the source term (105), the solution assumes three different behaviors on the domain: quadratic, constant and exponential. The first problem is solved with Peclet number Pe = 5 and

Γ =

0.002 while the second with

Pe

=10 and

Γ =

0.001. Let consider the domain

Ω =

[0,1] partitioned in three non-overlapping subdomains

Ω_{1}

Ω_{2}

and

Ω_{3}

, interspersed by two Bridge FV. The subdomain

Ω_{1}

has 18 volumes with diameter 0.04,

Ω_{2}

has 4 volumes with diameter 0.02 and

Ω_{3}

has 13 volumes with diameter 0.01. The numerical solutions are obtained by iterative processes, simulating transients with

Δt =

0.002 parameter that allows for achieving the steady states, i.e. the convergence to the analytical solutions.

The

L^{2}

error norm

Err

assumes, after 1000 steps, the values

Err =

0.000578 and

Err =

0.000501 in the first and second problem, respectively. Figure 11 shows the solution for 1000 steps with Pe=10 and

Γ

=0.001. It is notable that the numerical and analytical solutions are very close (in both the tests) and the absolute absence of oscillations. For comparisons with the results obtained from other methods see

[17, 37, 50].

Download: Download full-size image

Figure 11. Graphical solution of the steady state problem with Pe = 10 and

Γ

= 0.001

L^{2}

error norm 5.01E-4.

7.4. Distribution of a Pollutant in a River

A study of the distribution of a pollutant in a river, developed by Kengni Jotsa in

[4]

, concerns the solution of the convection–diffusion-reaction equation (106)

\frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} - Γ \frac{\partial^{2} C}{\partial x^{2}} + ac = S

(106)

where

c (x, t)

is the concentration of the pollutant and the velocity

u

is:

u (x, t) = \frac{[0.9 - 0.2 \cos (2 Fx) \sin (2 Ft)]}{[3 + 0.2 \sin (2 Fx) \cos (2 Ft)]} .

(107)

Download: Download full-size image

Figure 12. Solution of a convection-diffusion-reaction problem for the distribution of a pollutant c in a river after 10800s.

The source function is

S (x, t) = - Fcos (Fx) \cos (Ft) + uFsin (Fx) \sin (Ft) - Γ F^{2} \cos (Fx) \sin (Ft) + a (1 - \cos (Fx) \sin (Ft))

(108)

with F=

\frac{4 π}{3800}

, the diffusion coefficient

Γ

= 0.001and the reaction coefficient

a =

0.0115.

The analytical solution of the problem (106) is

c (x, t) = 1 - \cos (Fx) \sin (Ft) .

(109)

The initial and boundary conditions of Dirichlet at inflow and Neumann at outflow are deduced from the analytical solution. The approximation of the diffusive term is by classical FD schemes. Let partition the domain

Ω =

[0,3800] with 949 equal FV and be

Δt =

4s for a transient of 2700 steps. Despite the fact that the velocity varies both in space and time and the h-adaptivity strategy is not applied, the graphical results illustrated in Figure 12, appears to be very accurate with the maximum error equal to 1.13E-2.

7.5. 1D Model for Shallow Water Flows

The fifth problem concerns the incompressible fluid-dynamic theory of the long waves

[51]

. The partial differential equation system consists of a momentum equation, where the unit width discharge

q

is the dependent variable, and of the continuity equation (or mass conservation), where the water height

h

relative to a reference level is the dependent variable (

h

is here replaced by the elevation

ξ

of the free surface)

\{\begin{matrix} \frac{\partial q}{\partial t} + \frac{\partial}{\partial x} (\frac{qq}{h}) - μ \frac{\partial^{2} q}{\partial x^{2}} + gh \frac{\partial ξ}{\partial x} = S \\ \frac{\partial ξ}{\partial t} + \frac{\partial q}{\partial x} = 0 \end{matrix}

(110)

where

S = - gh \frac{\partial z_{b}}{\partial x} + g S_{f}

μ

is the diffusion coefficient,

g

is the gravity acceleration,

\frac{\partial z_{b}}{\partial x}

is the bottom slope and

S_{f}

is the bottom friction effects depending on the Strickler parameter

K

S_{f} = \frac{qq}{{(K)}^{2} {(h)}^{\frac{7}{3}}}

(111)

The solution of the above system by Ambrosi et al.

[52]

and Kengni Jotsa et al.

[53]

is based on a fractional step procedure for time advancing, which allows that the physical contributions are decoupled. The non linear convective term presented in the first equation of (110) is linearized at first order "freezing" velocity at the previous time step. The spatial discretization of the equations is by equal FV, so that the numerical solution is second order in space. In detail, the fractional step consists in

\frac{q^{n + \frac{1}{3}} - q^{n}}{Δt} + u^{n} \frac{{\partial q}^{n}}{\partial x} = 0

with

u^{n} = \frac{q^{n}}{h^{n}}

(112)

\frac{q^{n + \frac{2}{3}} - q^{n + \frac{1}{3}}}{Δt} = - g h^{n} \frac{\partial z_{b}}{\partial x} + g S_{f} + μ \frac{\partial^{2} q}{{\partial x}^{2}}

with

S_{f} = \frac{q^{n + \frac{1}{3}} q^{n + \frac{2}{3}}}{{(K)}^{2} {(h^{n})}^{\frac{7}{3}}}

(113)

ξ^{n + 1} {- (Δt)}^{2} g h^{n} \frac{\partial^{2} ξ^{n + 1}}{\partial x^{2}} + \frac{Δt}{h^{n}} q^{n + \frac{2}{3}} \frac{\partial ξ^{n + 1}}{\partial x} = ξ^{n} + \frac{Δt}{h^{n}} q^{n + \frac{2}{3}} \frac{\partial ξ^{n}}{\partial x} - Δt \frac{{\partial q}^{n + \frac{2}{3}}}{\partial x}

(114)

\frac{q^{n + 1} - q^{n + \frac{2}{3}}}{Δt} = - g h^{n} \frac{\partial ξ^{n + 1}}{\partial x} + \frac{q^{n + \frac{2}{3}}}{h^{n}} \frac{{\partial q}^{n + 1}}{\partial x}

with

ξ^{n} = h^{n} - h_{b}

(115)

and with

h_{b}

bottom depth respect to a reference level. To obtain the equation (114), of Helmotz type (see Kengni Jotsa et al.

[53])

Algorithm 7.1.

The four numerical steps are:

1) let obtain

q^{n + \frac{1}{3}} :

calculated

u^{n} = \frac{q^{n}}{h^{n}}

, solve explicitly the advection equation (112) by Fastest, assigning Dirichlet boundary condition at inflow and a numerical boundary extrapolation at outflow

2) let obtain

q^{n + \frac{2}{3}} :

solve the equation (113), assigned the bottom slope

α = \frac{\partial z_{b}}{\partial x}

, the diffusion parameter

μ

and the Strickler friction parameter

K

3) let obtain

ξ^{n + 1} :

solve the implicit equation (114), assigning homogeneous Neumann boundary condition at inflow and homogeneous Dirichlet boundary condition at outflow, by the Thomas solver (to note, the matrix of the system has the diagonally dominant property and is well-conditioned)

4) let obtain

q^{n + 1} :

solve the equation (115) with the updated elevation

ξ^{n + 1}

and finally, update the values of

h^{n}

A shallow water problem with known analytical solution has been solved by means of the above scheme by Kengni Jotsa et al.

[53]

and Vignoli et al.

[54]

. Let defined A=2

π / 3800

B = \cos (Ax)

C = \cos (At)

D = \sin (Ax)

E = \sin (At)

, the expression of the elevation, the discharge and the water height, respectively are

ξ(x,t)=0.2*C*D, q(x,t) = 0.9 - 0.2*B*E, h(x,t) = 3 +ξ(x,t)(116)

giving rise to the source function

f(x,t)=5.886*A*B*C+0.3924*A*B*

C^{2}

*D+

\frac{0.36 * A * D * E}{3 + 0.2 * D * C} - \frac{0.08 * A * D * E^{2}}{3 + 0.2 * D * C} \frac{- 0.2 * A * B * C * {(0.9 - 0.2 * B * E)}^{2}}{{(3 + 0.2 * D * C)}^{2}}

. (117)

Download: Download full-size image

Figure 13. Solution of a 1D Shallow Water problem elevation

ξ

after 10800

s

Download: Download full-size image

Figure 14. Solution of a 1D Shallow Water problem unit discharge

q

after 10800s.

The physical parameters are:

μ =

0.00001 and

g =

9.81, and the slope and friction effects are considered null. Let partition the domain

Ω =

[3800] in 1899 equal FV and be

Δt =

1s. The initial and boundary conditions are obtained by the analytical solution (116). Figure 13 and Figure 14 show graphically the elevation

ξ

and the unit discharge

q

, at the end of a transient of 10800s. The values of error norms

L^{2}

and

L_{\infty}

are respectively 0.2405 and 0.0054 for

ξ,

0.2864 and 0.0095 for

q

. A large number of numerical tests have been carried out with optimal results (here omitted for sake of brevity).

The system of equations (110) assumes importance for a generalization to shallow water 2D and 3D problems, as well as to study the propagation of moderately long waves while considering their dispersion (e.g. Boussinesq equations), and finally for the simulation of flows in Open Channel problems

[53]

8. Conclusions

In this article a new conservative scheme named Fastest, is developed to solve 1D convection-diffusion problems. The scheme achieves the second order of accuracy in time and space for equal FV and second order in time and first order in space for not equal FV. A technique based on cubic

ν

-spline polynomials and on Taylor expansions in time is applied to define Fastest. A significant property of the new method consists of being locally monotonic so that, in presence of continuous velocity gradients, spurious oscillations are avoided by the definition itself. The accuracy is preserved even when the transformation into the normalized plane is not possible by using second order back-ward Characteristics. The analysis of Fastest is completed in the normalized plane in contexts in which the velocity is dependent of both time and space.

The dynamically locally self h-adaptive strategy is achieved by the partition of the domain into subdomains which FV (at least four) are equal, and by the respect the rule whereby the ratio

σ

of the diameters of FV belonging to adjacent subdomains is

σ = 2

σ = 1 / 2

. To guarantee the consistency and stability properties of the method, a particular FV denoted as the Bridge FV is positioned between adjacent subdomains. It is characterized to have its node no longer positioned on the middle point (as usual for other FV) and its diameter determined by the average of the diameters of the FV belonging to the adjacent subdomains. Appropriate profiles and integration formulae are developed and applied if a change in subdomains occurs.

The h-adaptive strategy is efficiently implemented by the technique named Pick & Roll that generates new partitions by conveniently updating the matrices referred to topological and geometrical data. Any missing data relevant to the velocity and the dependent function, both on the nodes and the faces, are interpolated by third-degree polynomials. Since additional interpolations are required to determine the functional value on the foot of characteristic, the polynomial coefficients are calculated by linear combinations of functional or velocity values and suitable scalars calculated before the start of transient.

GKSO theory is applied to theoretically study the consistency, stability and convergence properties of Fastest by analysing the initial-boundary value of hyperbolic and parabolic problems. It is therefore possible to define the accuracy of the method and of the solution also in presence of h-adaptivity.

Generalizations of Fastest to 2D and 3D problems are presented, in particular by resorting the so called transverse curvatures which represent the effects of transverse transport in two- and three-dimensional flows.

The advantages of the Fastest method are demonstrated by solving problems with known analytical solution such as classical advection problems including advancing step and solitary wave, a stationary parabolic problem with source function and Dirichlet boundary conditions, a transient convection-diffusion-reaction problem regarding the distribution of a pollutant in a river and lastly a 1D transient shallow water problem.

To conclude, it should be noted how the application of h-adaptivity has been highly profitable by reducing significantly the computational cost if compared to the use of fine grid over the entire domain.

Abbreviations

CFD	Computational Fluid Dynamics
FEM	Finite Element Methods
FDM	Finite Difference Methods
FVM	Finite Volume Methods
ENO	Essentially Non-Oscillatory Methods
WENO	Weight Essentially Non-Oscillatory Methods
TVD	Total Variation Diminishing Methods
NVD	Normalized Variable Diagram
GFD	Generalized Finite Difference Formulae

Acknowledgments

Authors thank M. D., Ph. D., G. V. Pennati for the review of the manuscript and thank the Reviewer for his careful work and the useful comments on the manuscript.

Author Contributions

Vincenzo Angelo Pennati: Conceptualization, Data curation, Formal Analysis, Methodology, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

Antoine Celestin Kengni Jotsa: Conceptualization, Data curation, Formal Analysis, Methodology, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

Jacques Tagoudjeu: Conceptualization, Data curation, Formal Analysis, Methodology, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

Data Availability

The manuscript has not associated data. Sources codes used in this study are available upon request.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial or not-for-profit sectors.

Conflict of Interest

The authors declare no conflicts of interest.

Appendix

Table 2. Matrix NOD (9,3) relevant to the initial topological data concerning problems 1 and 2.

IIN1 = 1	IFI1 = 18	NFV1 = 18
INB1 = 19	NU1 = 42	NFB1 = 1
IIN2 = 20	IFI2 = 29	NFV2 = 10
INB2 = 30	NU2 = 21	NFB2 = 1
IIN3 = 31	IFI3=65	NFV3=35
INB3 = 66	NU3=12	NFB3=1
IIN4 = 67	IFI4=76	NFV4=10
INB4 = 77	NU4=24	NFB4=1
IIN5 = 78	IFI5=287	NFV5=210

IINi, IFIi and NFVi refer to the first and last node and the number of nodes of each subdomain

Ω_{i} i = 1, \dots, 5

, respectively. INBi, NUi, NFBi refer to the number of the node, the code of the diameters of the adjacent FV and the code 1 of each Bridge FV,

i = 1, \dots, 4,

respectively.

Table 3. Matrix ITNOD1 (9,3) for the transformation of the matrix NOD concerning the right flag RF external to

Ω_{3}

, through the Pick & Roll technique.

$IIN 1^{m + 1} = IIN 1^{m}$	$IFI 1^{m + 1} = IFI 1^{m}$	$NFV 1^{m + 1} = NFV 1^{m}$
$INB 1^{m + 1} = INB 1^{m}$	$NU 1^{m + 1} = NU 1^{m}$	$NFB 1^{m + 1} = NFB 1^{m}$
$IIN 2^{m + 1} = IIN 2^{m}$	$IFI 2^{m + 1} = IFI 2^{m}$	$NFV 2^{m + 1} = NFV 2^{m}$
$INB 2^{m + 1} = INB 2^{m}$	$NU 2^{m + 1} = NU 2^{m}$	$NFB 2^{m + 1} = NFB 2^{m}$
$IIN 3^{m + 1} = IIN 3^{m}$	$IFI 3^{m + 1} = IFI 3^{m} + 4$	$NFV 3^{m + 1} = NFV 3^{m} + 4$
$INB 3^{m + 1} = INB 3^{m} + 4$	$NU 3^{m + 1} = NU 3^{m}$	$NFB 3^{m + 1} = NFB 3^{m}$
$IIN 4^{m + 1} = IIN 4^{m} + 4$	$I FI 4^{m + 1} = IFI 4^{m} + 4$	$NFV 4^{m + 1} = NFV 4^{m}$
$INB 4^{m + 1} = INB 4^{m} + 4$	$NU 4^{m + 1} = NU 4^{m}$	$NFB 4^{m + 1} = NFB 4^{m}$
$IIN 5^{m + 1} = IIN 5^{m} + 4$	$IFI 5^{m + 1} = IFI 5^{m} + 3$	$NFV 5^{m + 1} = NFV 5^{m} - 1$

Table 4. Matrix ABSNF (9,4) relevant to the initial geometrical data concerning the problems 1 and 2.

$ANIN 1 = 0.04$	$ANFI 1 = 0.72$	$AFSI 1 = 0.02$	$AFDF 1 = 0.74$
$AINB 1 = 0.76$	$DVB 1 = 0.03$	$AFSB 1 = 0.74$	$AFDB 1 = 0.77$
$ANIN 2 = 0.78$	$ANFI 2 = 0.96$	$AFSI 2 = 0.77$	$AFDF 2 = 0.97$
$AINB 2 = 0.98$	$DVB 2 = 0.015$	$AFSB 2 = 0.97$	$AFDB 2 = 0.985$
$ANIN 3 = 0.99$	$ANFI 3 = 1.33$	$AFSI 3 = 0.985$	$AFDF 3 = 1.335$
$AINB 3 = 1.34$	$DVB 3 = 0.015$	$AFSB 3 = 1.33$	$AFDB 3 = 1.35$
$ANIN 4 = 1.36$	$ANFI 4 = 1.54$	$AFSI 4 = 1.35$	$AFDF 4 = 1.55$
$AINB 4 = 1.56$	$DVB 4 = 0.03$	$AFSB 4 = 1.55$	$AFDB 4 = 1.58$
$ANIN 5 = 1.60$	$ANFI 5 = 9.96$	$AFSI 5 = 1.58$	$AFDF 5 = 9.98$

ANINi and ANFIi refer to the abscissa of first and last node, AFSIi and AFDFi refer to the abscissa of the left face of the first FV and of the right face of the last FV of each subdomain

Ω_{i} i = 1, \dots, 5

, respectively. AINBi, AFSBi and AFDBi refer to the abscissas of the node, the left and right face of each Bridge FV

i = 1, \dots, 4,

respectively. DVBi refer to the diameter value of each Bridge FV,

i = 1, \dots, 4

Table 5. Matrix COLOCH (19,4) and column ICOLH for the interpolation of ξ on the foot of the characteristic for the problems 1 and 2.

$ICOLH$	$COLOCH$				$ICOLH$	$COLOCH$
1111	$- 1.5$	$- 0.5$	$+ 0.5$	$+ 1.5$	4444	$- 6.0$	$- 2.0$	$+ 2.0$	$+ 6.0$
5111	$- 1.75$	$- 0.5$	$+ 0.5$	$+ 1.5$	4443	$- 6.0$	$- 2.0$	$+ 2.0$	$+ 5.5$
2511	$- 2.5$	$- 0.75$	$+ 0.5$	$+ 1.5$	4432	$- 6.0$	$- 2.0$	$+ 1.5$	$+ 4.0$
2251	$- 3.0$	$- 1.0$	$+ 0.75$	$+ 2.0$	4322	$- 5.0$	$- 1.5$	$+ 1.0$	$+ 3.0$
2225	$- 3.0$	$- 1.0$	$+ 1.0$	$+ 2.75$	3222	$- 3.5$	$- 1.0$	$+ 1.0$	$+ 3.0$
1115	$- 1.5$	$- 0.5$	$+ 0.5$	$+ 1.75$	2222	$- 3.0$	$- 1.0$	$+ 1.0$	$+ 3.0$
1152	$- 1.5$	$- 0.5$	$+ 0.75$	$+ 2.5$	2223	$- 3.0$	$- 1.0$	$+ 1.0$	$+ 3.5$
1522	$- 2.0$	$- 0.75$	$+ 1.0$	$+ 3.0$	2234	$- 3.0$	$- 1.0$	$+ 1.5$	$+ 5.0$
5222	$- 2.75$	$- 1.0$	$+ 1.0$	$+ 3.0$	2344	$- 4.0$	$- 1.5$	$+ 2.0$	$+ 6.0$
					3444	$- 5.5$	$- 2.0$	$+ 2.0$	$+ 6.0$

ICOLH: codes relevant to all the possible combinations of the four FV composing the support of the third-degree polynomial. COLOCH: local abscissa a_i of the nodes belonging to the support, multiplied by 100

References

[1]	Franca, L. P., Frey, S. L., Hughes, T. J. R.: Stabilized finite element methods: I. Application to the advective-diffusive model. Comput. Methods Appl. Mech. Engrg. 95, 253-276 (1992). https://doi.org/10.1016/0045-7825(92)90143-8
[2]	Hughes, T. J. R., Franca, L. P., Hulbert, G. M.: A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations. Comput. Methods Appl. Mech. Engrg. 73, 173-189 (1989). https://doi.org/10.1016/0045-7825(89)90111-4
[3]	Kengni Jotsa, A. C., Pennati, V. A.: Solution of the 2D Navier-Stokes equations by a new FE fractional step method. IMACS/ISGG WORKSHOP, 141-150, Roma, Italy, 2011.
[4]	Kengni Jotsa, A. C.: Solution of the 2D Navier-Stokes equations by a new FE fractional step method. Doctoral dissertation, Università degli Studi dell’Insubria, Como, Italy. (2012). Available: http://hdl.handle.net/11383/2090232 Accessed 30 June 2024
[5]	Kengni Jotsa, A. C., Pennati, V. A.: A cost-effective FE method for 2D Navier–Stokes equations. Engineering Applications of Computational Fluid Mechanics. 9, 66-83 (2015). https://doi.org/10.1080/19942060.2015.1004811
[6]	Khelifa, A., Robert, J.-L., Ouellet, Y.: A Douglas-Wang finite element approach for transient advection-diffusion problems. Comput. Methods Appl. Mech. Engrg. 10, 113-129 (1993). https://doi.org/10.1016/0045-7825(93)90023-Q
[7]	Siegel, P., Mosé, R., Ackerer, P., Jaffre, J.: Solution of the Advection-Diffusion Equation Using a Combination of Discontinuous and Mixed Finite Elements. Int. J. Numer. Meth. Fluids. 24, 595-613 (1997). https://doi.org/10.1002/(SICI)1097-0363(19970330)24:6<595::AID-FLD512>3.0.CO;2-I
[8]	Harten, A., Engquist, B., Osher, S., Chakravarthy, S. R.: Uniformly high order accurate essentially non-oscillatory schemes, III. J. Comput. Phys. 71, 231-303 (1987). https://doi.org/10.1016/0021-9991(87)90031-3
[9]	Quarteroni, A.: Modellistica numerica per problemi differenziali. Springer-Verlag, Milano (2008).
[10]	Leonard, B. P.: A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comput. Methods Appl. Mech. Engrg. 19, 59-98 (1979). https://doi.org/10.1016/0045-7825(79)90034-3
[11]	Leonard, B. P.: Simple high-accuracy resolution program for convective modelling of discontinuities. Int. J. Numer. Meth. Fluids. 8, 1291-1318 (1988). https://doi.org/10.1002/fld.1650081013
[12]	De Biase, L., Feraudi, F., Pennati, V. A.: A finite volume method for the solution of convection-diffusion 2D problems by a quadratic profile with smoothing. Int. J. Num. Meth. Heat Fluid Flow. 6, 3-24 (1996). https://doi.org/10.1108/09615539610123414
[13]	Deponti, A., Pennati, V. A., De Biase, L.: A fully 3D finite volume method for incompressible Navier–Stokes equations. Int. J. Num. Meth. Fluids. 52, 617-638 (2006). https://doi.org/10.1002/fld.1190
[14]	Deponti, A., Pennati, V. A., De Biase, L., Maggi, V., Berta, F.: A new fully three-dimensional numerical model for ice dynamics. Journal of Glaciology. 52, 365-376 (2006). https://doi.org/10.3189/172756506781828629
[15]	Feraudi, F., Pennati, V. A.: Trasporto del Calore in Fluidi Incomprimibili: un Nuovo Approccio Numerico per Problemi Bidimensionali Non Stazionari. L’Energia Elettrica. 74, 238-251 (1993).
[16]	Leonard, B. P.: The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection. Comput. Methods Appl. Mech. Engrg. 88, 17-74 (1991). https://doi.org/10.1016/0045-7825(91)90232-U
[17]	Leonard, B. P., Mokhtari, S.: Beyond first-order upwinding: The ultra-sharp alternative for non-oscillatory steady-state simulation of convection. Int. J. Numer. Meth. Eng., 30, 729-766 (1990). https://doi.org/10.1002/nme.1620300412
[18]	Babu, G., Bansal, K.: A high order robust numerical scheme for singularly perturbed delay parabolic convection diffusion problems. Journal of Applied Mathematics and Computing. 68, 363-389 (2022). https://doi.org/10.1007/s12190-021-01512-1
[19]	Du, Y., Wang, Y., Yuan, L.: A high-order modified finite-volume method on Cartesian grids for nonlinear convection–diffusion problems. Comp. Appl. Math. 39, 214 (2020). https://doi.org/10.1007/s40314-020-01253-0
[20]	Xu, M.: A high-order finite volume scheme for unsteady convection-dominated convection–diffusion equations. Numerical Heat Transfer, Part B: Fundamentals. 76, 253-272 (2019). https://doi.org/10.1080/10407790.2019.1665421
[21]	Zhu, X., Rui, H.: High-order compact difference scheme of 1D nonlinear degenerate convection–reaction–diffusion equation with adaptive algorithm. Numerical Heat Transfer, Part B: Fundamentals. 75, 43-66 (2019). https://doi.org/10.1080/10407790.2019.1591858
[22]	Huang, L., Jiang, Z., Lou, S., Zhang, X., Yan, C.: Simple and robust h-adaptive shock-capturing method for flux reconstruction framework. Chinese Journal of Aeronautics. 36, 348-365 (2023). https://doi.org/10.1016/j.cja.2023.04.010
[23]	Nishikawa, H.: The QUICK Scheme is a Third-Order Finite-Volume Scheme with Point-Valued Numerical Solutions. arXiv preprint arXiv:2006.15143 [math. NA]. (2020).
[24]	Oh, J.: Adaptive Mesh-free Methods for Partial Differential Equations. Doctoral dissertation. Department of Mathematics of the University of Southern Mississippi, USA. (2018). Available: https://aquila.usm.edu/dissertations/1544 Accessed 30 June 2024.
[25]	Sabbagh-Yazdi, S.-R., Amiri, T., Asil Gharebaghi, S.: A Proposed Damping Coefficient of Quick Adaptive Galerkin Finite Volume Solver for Elasticity Problems. Journal of the Serbian Society for Computational Mechanics. 13, 56-79 (2019). https://doi.org/10.24874/jsscm.2019.13.01.04
[26]	Zheng, Q., Liu, Z.: Uniform convergence analysis of a new adaptive upwind finite difference method for singularly perturbed convection reaction - diffusion boundary value problem. Journal of Applied and Computing. 70, 601-618 (2024). https://doi.org/10.1007/s12190-023-01977-2
[27]	Kumar, S., Das, P.: A uniformly convergent analysis for multiple scale parabolic singularly perturbed convection diffusion coupled systems: Optimal accuracy with less computational time. Appl. Numer. Math. 207, 534-557 (2025). https://doi.org/10.1016/j.apnum.2024.09.020
[28]	Ayalew, M., Aychluh, M., Suthar, D. L., Purohit, S. D.: Quadratic upwind differencing scheme in the finite volume method for solving the convection-diffusion equation. Mathematical and Computer Modelling of Dynamical Systems. 29, 265-285 (2023). https://doi.org/10.1080/13873954.2023.2282974
[29]	Peng, G.: A positivity-preserving finite volume scheme for convection–diffusion equation on general meshes. International Journal of Computer Mathematics. 99, 355-369 (2022). https://doi.org/10.1080/00207160.2021.1910943
[30]	Kumar, N., Thije Boonkkamp, J. H. M., Koren, B.: Approximation of the convective flux in the incompressible Navier–Stokes equations using local boundary-value problems. Journal of Computational and Applied Mathematics. 340, 523-536 (2018). https://doi.org/10.1016/j.cam.2018.01.030
[31]	Sukhinov, A., Chistyakov, A., Kuznetsova, I., Belova, Y., Rahimbaeva, E.: Development and Research of a Modified Upwind Leapfrog Scheme for Solving Transport Problems. Mathematics. 10, 3564 (2022).
[32]	Kurganov, A.: Finite-volume schemes for shallow-water equations. Cambridge University Press. Acta Numerica. 27, 289-351 (2018). https://doi.org/10.1017/S0962492918000028
[33]	Neelan, A. A. G., Nair, M. T., Bürger, R.: Three-level order-adaptive weighted essentially non-oscillatory schemes. Results in Applied Mathematics. 12, 100217 (2021). https://doi.org/10.1016/j.rinam.2021.100217
[34]	Kossaczká, T., Ehrhardt, M., Günther, M.: Enhanced fifth order WENO shock-capturing schemes with deep learning. Results in Applied Mathematics. 12, 100201 (2021). https://doi.org/10.1016/j.rinam.2021.100201
[35]	Jiang, G.-S., Shu, C.-W.: Efficient Implementation of Weighted ENO Schemes. J. Comput. Phys. 126, 202-228 (1996). https://doi.org/10.1006/jcph.1996.0130
[36]	Xu, M.: A finite volume scheme for unsteady linear and nonlinear convection-diffusion-reaction problems. International Communications in Heat and Mass Transfer. 139, 106417 (2022). https://doi.org/10.1016/j.icheatmasstransfer.2022.106417
[37]	Pennati, V. A., Marelli, M., De Biase, L. M.: FAST - A generalized FV solution of convection-diffusion problems by means of monotonic v-splines profiles. Int. J. Num. Meth. Heat Fluid Flow. 6, 3-30 (1996). https://doi.org/10.1108/eb017549
[38]	Roache, P. J.: Computational fluid dynamics. Hermosa Publishers, Albuquerque, NM (1972).
[39]	Godunov, S. K.: A difference scheme for numerical solution of discontinuous solution of hydrodynamic equations. Mathematics of the USSR-Sbornik. 47, 271-306 (1959). Translated US Joint Publ. Ross. Service, JPRS7226, 1969.
[40]	Dullemond, C. P.: Lecture Numerical Fluid Dynamics. Chapter 4. University of Heidelberg: Heidelberg University Press (2011). https://www.ita.uni-heidelberg.de/~dullemond/lectures/num_fluid_2011. Accessed 30 June 2024
[41]	Strikwerda, J. C.: Finite difference schemes and partial differential equations. Wadsworth & Brooks/Cole Advanced Books & Software, Pacific Grove, California, USA (1989).
[42]	Jury, E. I.: On the roots of a real polynomial inside the unit circle and a stability criterion for linear discrete systems. IFAC Proceedings Volumes. 1, 142-153 (1963). https://doi.org/10.1016/S1474-6670(17)69648-4
[43]	Marden, M.: The Geometry of the Zeros of a Polynomial in a Complex Variable. American Mathematical Society, New York (1949).
[44]	Reali, M., Dassie, G., Pennati, V. A.: Direct general finite difference techniques for elliptic problems defined on bounded or unbounded 2-D domains. Lewis, R. W. et al. (ed.). Numerical Methods for Transient and Coupled Problems, pp. 43-58. Wiley. (1987).
[45]	Corti, S., De Biase, L., Pennati, V. A.: A multistep three dimensional generalized finite difference technique. Sixth International Conference on Numerical methods in Thermal problems. Pineridge Press, Swansea, U. K. (1989).
[46]	Kelly, D. W., Mills, R. J., Reizes, J. A., Miller, A. D.: A posteriori error estimates in finite difference techniques. J. Comput. Phys. 74, 214-232 (1988). https://doi.org/10.1016/0021-9991(88)90077-0
[47]	Pennati, V., Corti, S.: Generalized finite-differences solution of 3D elliptical problems involving Neumann boundary conditions. Commun. Numer. Meth. Engng. 10, 43-58 (1994). https://doi.org/10.1002/cnm.1640100106
[48]	Thompson, J. F., Soni, B. K., Weatherill, N. P.: Handbook of grid generation. CRC press (1998).
[49]	Deponti, A.: Mass and thermal flows in Alpine glaciers. Application to Lys Glacier (Monte Rosa, Italian Alps). Doctoral dissertation, University of Milano–Bicocca. Italy (2003).
[50]	Cordero, E., De Biase, L., Pennati, V. A.: A new finite volume method for the solution of convection–diffusion equations: analysis of stability and convergence. Commun. Numer. Meth. Engng. 13, 923-940 (1997).
[51]	Cunge, J. A., Holly, F. M., Verwey, A.: Practical Aspects of Computational River Hydraulics. Pitman, London (1980).
[52]	Ambrosi, D., Corti, S., Pennati, V., Saleri, F.: Numerical Simulation of Unsteady Flow at Po River Delta. Journal of Hydraulic Engineering. 122, 735-743 (1996). https://doi.org/10.1061/(ASCE)0733-9429(1996)122:12(735)
[53]	Kengni Jotsa, A. C., Pennati, V. A., Di Guardo, A., Morselli, M.: Shallow Water 1D Model for Pollution River Study. Pure and Applied Mathematics Journal. 6, 76-88 (2017). https://doi.org/10.11648/j.pamj.20170602.13
[54]	Vignoli, G., Titarev, V. A., Toro, E. F.: ADER schemes for the shallow water equations in channel with irregular bottom elevation. J. Comput. Phys. 227, 2463-2480 (2008). https://doi.org/10.1016/j.jcp.2007.11.006

Cite This Article

Plain Text BibTeX RIS

APA Style

Pennati, V. A., Jotsa, A. C. K., Tagoudjeu, J. (2025). Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems. Pure and Applied Mathematics Journal, 14(4), 69-92. https://doi.org/10.11648/j.pamj.20251404.11

Copy | Download

ACS Style

Pennati, V. A.; Jotsa, A. C. K.; Tagoudjeu, J. Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems. Pure Appl. Math. J. 2025, 14(4), 69-92. doi: 10.11648/j.pamj.20251404.11

Copy | Download

AMA Style

Pennati VA, Jotsa ACK, Tagoudjeu J. Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems. Pure Appl Math J. 2025;14(4):69-92. doi: 10.11648/j.pamj.20251404.11

Copy | Download

@article{10.11648/j.pamj.20251404.11,
  author = {Vincenzo Angelo Pennati and Antoine Celestin Kengni Jotsa and Jacques Tagoudjeu},
  title = {Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems
},
  journal = {Pure and Applied Mathematics Journal},
  volume = {14},
  number = {4},
  pages = {69-92},
  doi = {10.11648/j.pamj.20251404.11},
  url = {https://doi.org/10.11648/j.pamj.20251404.11},
  eprint = {https://article.sciencepublishinggroup.com/pdf/10.11648.j.pamj.20251404.11},
  abstract = {Several recently published studies regarding flow problems propose schemes of high order of accuracy designed as evolution of traditional methods. A drawback common to these new schemes is the necessity to adopt uniform mesh refinement for solving sharp problems, by increasing the computational cost. Even the so called essentially non-oscillatory and weight essentially non-oscillatory methods suffer of the same drawback and are not suitable to cope with h-adaptive methods due to their definition on finite volumes necessarily of equal diameter. Therefore, in order to overcome the above drawback, the formulation of dynamically locally self h-adaptive processes is designed to achieve the dual purpose to increase the accuracy and to keep as small as possible the number of finite volumes. To define a locally h-adaptive finite volume (FV) scheme need two simple but important tools, namely a particular FV named Bridge FV positioned between two adjacent subdomains and the definition of suitable profiles approximating the fluxes on the FV faces. In this article a new FV method for the numerical solution of convective-diffusive 1D problems is developed. It is conservative, second order in time and space for equal FV, and allows the partitioning of the domain by equal or unequal finite volumes, thus dynamically locally self h-adaptive. The definition of the monotonic profiles is accomplished by means of cubic weighted ν-splines and Taylor expansions. The profile analysis respect to the numerical properties is conducted in the normalized plane with the velocity varying in time and space and gives the flux value on the FV faces. Moreover the flux is assigned by Upwind or by second order back-ward Characteristics if the estimated flux is outside of the unit square or the transformation into the normalized plane is not possible, respectively. The initial-boundary stability and convergence properties of the new method are examined in detail, also in presence of h-adaptivity. In addition, a generalization of the new scheme to 2D and 3D problems is presented. Finally, some numerical test are carried out to verify the properties of the new method, including two CFD problems.},
 year = {2025}
}

Copy | Download

TY  - JOUR
T1  - Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems

AU  - Vincenzo Angelo Pennati
AU  - Antoine Celestin Kengni Jotsa
AU  - Jacques Tagoudjeu
Y1  - 2025/08/15
PY  - 2025
N1  - https://doi.org/10.11648/j.pamj.20251404.11
DO  - 10.11648/j.pamj.20251404.11
T2  - Pure and Applied Mathematics Journal
JF  - Pure and Applied Mathematics Journal
JO  - Pure and Applied Mathematics Journal
SP  - 69
EP  - 92
PB  - Science Publishing Group
SN  - 2326-9812
UR  - https://doi.org/10.11648/j.pamj.20251404.11
AB  - Several recently published studies regarding flow problems propose schemes of high order of accuracy designed as evolution of traditional methods. A drawback common to these new schemes is the necessity to adopt uniform mesh refinement for solving sharp problems, by increasing the computational cost. Even the so called essentially non-oscillatory and weight essentially non-oscillatory methods suffer of the same drawback and are not suitable to cope with h-adaptive methods due to their definition on finite volumes necessarily of equal diameter. Therefore, in order to overcome the above drawback, the formulation of dynamically locally self h-adaptive processes is designed to achieve the dual purpose to increase the accuracy and to keep as small as possible the number of finite volumes. To define a locally h-adaptive finite volume (FV) scheme need two simple but important tools, namely a particular FV named Bridge FV positioned between two adjacent subdomains and the definition of suitable profiles approximating the fluxes on the FV faces. In this article a new FV method for the numerical solution of convective-diffusive 1D problems is developed. It is conservative, second order in time and space for equal FV, and allows the partitioning of the domain by equal or unequal finite volumes, thus dynamically locally self h-adaptive. The definition of the monotonic profiles is accomplished by means of cubic weighted ν-splines and Taylor expansions. The profile analysis respect to the numerical properties is conducted in the normalized plane with the velocity varying in time and space and gives the flux value on the FV faces. Moreover the flux is assigned by Upwind or by second order back-ward Characteristics if the estimated flux is outside of the unit square or the transformation into the normalized plane is not possible, respectively. The initial-boundary stability and convergence properties of the new method are examined in detail, also in presence of h-adaptivity. In addition, a generalization of the new scheme to 2D and 3D problems is presented. Finally, some numerical test are carried out to verify the properties of the new method, including two CFD problems.
VL  - 14
IS  - 4
ER  -

Copy | Download

Author Information

Vincenzo Angelo Pennati

Theoretical Analysis Service, Hydraulic & Structural Research Center, ENEL, Milano, Italy

Contact Email

http://orcid.org/0000-0002-0770-7645
Antoine Celestin Kengni Jotsa

Department of Fundamental Sciences, National Advanced School of Mines and Petroleum Industries at Kaele, University of Maroua, Maroua, Cameroon

Contact Email

http://orcid.org/0000-0003-3395-2561
Jacques Tagoudjeu

Department of Mathematics and Physical Sciences, National Advanced School of Engineering of Yaounde, University of Yaounde I, Yaounde, Cameroon

Contact Email

http://orcid.org/0000-0002-9379-3171

Download PDF

Submit an Article

Figure 1

Figure 1. Grids with halving (a) and doubling (b) the net size and Bridge FV.
Figure 2

Figure 2. Normalized Variable Diagram of Fastest with $σ = 1$ and Courant numbers $c_{r} = 0, 1 / 2, 1 .$
Figure 3

Figure 3. Graphic of surface $\sqrt{{|G (Θ, c)|}^{2}}$ for $Θ \in (0, 2 π]$ and $c \in (0, 1] .$
Figure 4

Figure 4. Graphical representation of the section $S = S (\bar{c}, hs)$ where $\bar{c}$ is the value c whose corresponding section $S$ has the maximum height $hs \leq 1$ .
Figure 5

Figure 5. Two-dimensional five-node stencil to estimate the left face value using Fastest (velocity component u positive).
Figure 6

Figure 6. Three-dimensional seven-node stencil to estimate the left face value using Fastest (velocity component $u$ positive).
Figure 7

Figure 7. Solution of an advection problem regarding a "jump" function at 4000th step with velocity $u =$ 1 and $Δt =$ 0.0005. The data relevant to the subdomains at the final step are Ω₁=[0.02;1.58], Ω₂=[1.61;1.82], Ω₃=[1.835;3.355], Ω₄=[3.37;3.55], Ω₅=[3.58;9.98]. Each subdomain Ω_k=[XL_k; XR_k] has XL_k= abscissa of left face of the firs FV of $Ω_{k}$ , XR_k=abscissa of right face of the last FV of Ω_k.
Figure 8

Figure 8. A zoom of the Figure 7 with FV relating to the subdomains $Ω_{3}$ , $Ω_{4}$ and the beginning of $Ω_{5}$ .
Figure 9

Figure 9. Solution of an advection problem regarding an "isolated wave" function at first step and 2000th step with velocity u=1 and Δt= 0.0005. The data relevant to the subdomains at the final step are: Ω₁=[0.02; 2.74],Ω₂=[2.77; 297], Ω₃=[2.985; 3.395], Ω₄=[3.41; 3.59], Ω₅=[3.62; 9.98].Each subdomain Ω_k=[XL_k; XR_k] has XL_k= abscissa of left face of the firs FV of $Ω_{k}$ , ${XR}_{k} =$ abscissa of right face of the last FV of Ω_k.
Figure 10

Figure 10. A zoom of the Figure 9 with the FV relating to the subdomains $Ω_{2}$ , $Ω_{3}$ , $Ω_{4}$ and the beginning of $Ω_{5}$ .
Figure 11

Figure 11. Graphical solution of the steady state problem with Pe = 10 and $Γ$ = 0.001 $L^{2}$ error norm 5.01E-4.
Figure 12

Figure 12. Solution of a convection-diffusion-reaction problem for the distribution of a pollutant c in a river after 10800s.
Figure 13

Figure 13. Solution of a 1D Shallow Water problem elevation $ξ$ after 10800 $s$ .
Figure 14

Figure 14. Solution of a 1D Shallow Water problem unit discharge $q$ after 10800s.

Plain Text BibTeX RIS

APA Style

Pennati, V. A., Jotsa, A. C. K., Tagoudjeu, J. (2025). Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems. Pure and Applied Mathematics Journal, 14(4), 69-92. https://doi.org/10.11648/j.pamj.20251404.11

Copy | Download

ACS Style

Pennati, V. A.; Jotsa, A. C. K.; Tagoudjeu, J. Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems. Pure Appl. Math. J. 2025, 14(4), 69-92. doi: 10.11648/j.pamj.20251404.11

Copy | Download

AMA Style

Pennati VA, Jotsa ACK, Tagoudjeu J. Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems. Pure Appl Math J. 2025;14(4):69-92. doi: 10.11648/j.pamj.20251404.11

Copy | Download

@article{10.11648/j.pamj.20251404.11,
  author = {Vincenzo Angelo Pennati and Antoine Celestin Kengni Jotsa and Jacques Tagoudjeu},
  title = {Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems
},
  journal = {Pure and Applied Mathematics Journal},
  volume = {14},
  number = {4},
  pages = {69-92},
  doi = {10.11648/j.pamj.20251404.11},
  url = {https://doi.org/10.11648/j.pamj.20251404.11},
  eprint = {https://article.sciencepublishinggroup.com/pdf/10.11648.j.pamj.20251404.11},
  abstract = {Several recently published studies regarding flow problems propose schemes of high order of accuracy designed as evolution of traditional methods. A drawback common to these new schemes is the necessity to adopt uniform mesh refinement for solving sharp problems, by increasing the computational cost. Even the so called essentially non-oscillatory and weight essentially non-oscillatory methods suffer of the same drawback and are not suitable to cope with h-adaptive methods due to their definition on finite volumes necessarily of equal diameter. Therefore, in order to overcome the above drawback, the formulation of dynamically locally self h-adaptive processes is designed to achieve the dual purpose to increase the accuracy and to keep as small as possible the number of finite volumes. To define a locally h-adaptive finite volume (FV) scheme need two simple but important tools, namely a particular FV named Bridge FV positioned between two adjacent subdomains and the definition of suitable profiles approximating the fluxes on the FV faces. In this article a new FV method for the numerical solution of convective-diffusive 1D problems is developed. It is conservative, second order in time and space for equal FV, and allows the partitioning of the domain by equal or unequal finite volumes, thus dynamically locally self h-adaptive. The definition of the monotonic profiles is accomplished by means of cubic weighted ν-splines and Taylor expansions. The profile analysis respect to the numerical properties is conducted in the normalized plane with the velocity varying in time and space and gives the flux value on the FV faces. Moreover the flux is assigned by Upwind or by second order back-ward Characteristics if the estimated flux is outside of the unit square or the transformation into the normalized plane is not possible, respectively. The initial-boundary stability and convergence properties of the new method are examined in detail, also in presence of h-adaptivity. In addition, a generalization of the new scheme to 2D and 3D problems is presented. Finally, some numerical test are carried out to verify the properties of the new method, including two CFD problems.},
 year = {2025}
}

Copy | Download

TY  - JOUR
T1  - Fastest - A New High Order FV Method Dynamically Locally Self h-Adaptive for Convective-diffusive Problems

AU  - Vincenzo Angelo Pennati
AU  - Antoine Celestin Kengni Jotsa
AU  - Jacques Tagoudjeu
Y1  - 2025/08/15
PY  - 2025
N1  - https://doi.org/10.11648/j.pamj.20251404.11
DO  - 10.11648/j.pamj.20251404.11
T2  - Pure and Applied Mathematics Journal
JF  - Pure and Applied Mathematics Journal
JO  - Pure and Applied Mathematics Journal
SP  - 69
EP  - 92
PB  - Science Publishing Group
SN  - 2326-9812
UR  - https://doi.org/10.11648/j.pamj.20251404.11
AB  - Several recently published studies regarding flow problems propose schemes of high order of accuracy designed as evolution of traditional methods. A drawback common to these new schemes is the necessity to adopt uniform mesh refinement for solving sharp problems, by increasing the computational cost. Even the so called essentially non-oscillatory and weight essentially non-oscillatory methods suffer of the same drawback and are not suitable to cope with h-adaptive methods due to their definition on finite volumes necessarily of equal diameter. Therefore, in order to overcome the above drawback, the formulation of dynamically locally self h-adaptive processes is designed to achieve the dual purpose to increase the accuracy and to keep as small as possible the number of finite volumes. To define a locally h-adaptive finite volume (FV) scheme need two simple but important tools, namely a particular FV named Bridge FV positioned between two adjacent subdomains and the definition of suitable profiles approximating the fluxes on the FV faces. In this article a new FV method for the numerical solution of convective-diffusive 1D problems is developed. It is conservative, second order in time and space for equal FV, and allows the partitioning of the domain by equal or unequal finite volumes, thus dynamically locally self h-adaptive. The definition of the monotonic profiles is accomplished by means of cubic weighted ν-splines and Taylor expansions. The profile analysis respect to the numerical properties is conducted in the normalized plane with the velocity varying in time and space and gives the flux value on the FV faces. Moreover the flux is assigned by Upwind or by second order back-ward Characteristics if the estimated flux is outside of the unit square or the transformation into the normalized plane is not possible, respectively. The initial-boundary stability and convergence properties of the new method are examined in detail, also in presence of h-adaptivity. In addition, a generalization of the new scheme to 2D and 3D problems is presented. Finally, some numerical test are carried out to verify the properties of the new method, including two CFD problems.
VL  - 14
IS  - 4
ER  -

Copy | Download