A Thermodynamically Consistent Master-Equation Framework for Fluid Dynamics: Unified Discretization, Structure-Preserving Algorithms, and Convergence to Navier–Stokes Solutions

1. Introduction

1.1. Background: Limitations of the Navier–Stokes-Centric Paradigm

Fluid dynamics has long been formulated on the basis of continuum field theories, most prominently the Navier–Stokes (NS) equations, which describe the evolution of mass, momentum, and energy under the assumptions of local thermodynamic equilibrium and constitutive closure [1,2,3,4]. Since the seminal works of Claude-Louis Navier and George Gabriel Stokes, the NS equations have served as the foundational model for both theoretical analysis and numerical simulation of fluid flows.

However, the traditional paradigm exhibits several structural limitations. First, the governing equations are postulated at the continuum level, and discretization is subsequently introduced as an approximation procedure [5,6,7]. This sequential approach—“PDE first, discretization later”—inevitably leads to a loss of intrinsic physical structure at the discrete level, necessitating additional stabilization mechanisms such as artificial viscosity, upwinding, or Riemann solvers [8,9,10]. These techniques, while effective in practice, are not derived from first principles and often compromise accuracy, particularly in under-resolved or high-Reynolds-number regimes [11,12].

Second, the enforcement of thermodynamic consistency, especially the second law (entropy inequality), is not guaranteed by standard discretizations. Although entropy-stable schemes have been developed [13,14,15], they typically require carefully engineered fluxes or entropy fixes, rather than arising naturally from the underlying formulation.

Third, modern computational methods—including finite volume (FV), finite element (FE), and meshless approaches—are often developed independently, resulting in fragmented numerical frameworks lacking a unifying theoretical structure [16,17,18]. This fragmentation becomes particularly problematic in complex geometries, multiphysics coupling, and adaptive simulations.

These limitations suggest that the difficulty is not merely numerical, but conceptual: the continuum PDE itself may not be the most fundamental representation of fluid systems. Instead, a more primitive description based on interacting elements may provide a more robust foundation.

1.2. Core Idea: From PDEs to Interaction-Based Master Equations

In this work, we propose a fundamental shift in perspective. Rather than starting from continuum partial differential equations, we formulate fluid dynamics as a network of interacting elements governed by a master equation.

The key idea is to regard the fluid not as a continuous field, but as a collection of nodes connected through interactions, where the time evolution of each node is determined by flux exchanges with its neighbors. Within this framework, the governing equation is constructed as a conservative interaction law, rather than postulated as a differential equation.

Two guiding principles are imposed at the design level:

• (i) Conservation laws are enforced structurally, through antisymmetric flux interactions between nodes.
• (ii) The second law of thermodynamics is embedded as a constraint, by constructing dissipative terms as entropy-gradient flows.

This approach is closely related in spirit to the GENERIC (General Equation for Non-Equilibrium Reversible–Irreversible Coupling) framework [19,20,21], but differs in that it is directly formulated on discrete interaction networks and explicitly targets numerical realizability.

A crucial consequence of this formulation is that the Navier–Stokes equations are no longer the starting point, but instead emerge as a continuum limit of the master equation under appropriate assumptions of locality, smoothness, and constitutive closure. In this sense, the NS equations are reinterpreted as an effective theory, valid within a restricted regime.

This inversion of the traditional hierarchy—from “PDE → discretization” to “interaction → continuum limit”—enables the construction of numerical methods that inherit physical structure by design, rather than by correction. The novelty of this work does not lie in proposing yet another entropy-stable discretization of the Navier–Stokes equations. Instead, we reformulate fluid dynamics itself as a thermodynamically constrained interaction system, from which both numerical schemes and the Navier–Stokes equations emerge as consequences.

1.3. Contributions of This Work

The present study develops this framework into a complete numerical methodology. The main contributions are summarized as follows:

• Master Equation Formulation
  A thermodynamically consistent master equation for fluid systems is derived, in which conservation laws and entropy production are built into the interaction structure.
• Thermodynamically Consistent Discretization
  A unified discretization is constructed based on antisymmetric fluxes (reversible) and entropy-gradient flows (irreversible), ensuring exact conservation and discrete entropy inequality.
• Unified Framework for FV and Meshless Methods
  Finite volume and meshless methods are shown to be special cases of a single graph-based formulation, differing only in geometry construction.
• IMEX Time Integration Structure
  A structure-preserving implicit–explicit (IMEX) time integration scheme is developed, maintaining conservation and entropy stability in discrete time.
• Rigorous Convergence Analysis
  Using entropy estimates and compactness arguments, it is proven that the discrete solutions converge (in a subsequence sense) to entropy weak solutions of the Navier–Stokes equations.

These contributions collectively establish a new paradigm in computational fluid dynamics, in which numerical methods are derived from fundamental physical principles rather than imposed on them.

1.4. Organization of the Paper

The remainder of this paper is organized as follows.

• Section 2 introduces the master equation formulation and its thermodynamic structure.
• Section 3 develops a graph-based spatial discretization framework.
• Section 4 and Section 5 construct reversible and irreversible fluxes, respectively, ensuring conservation and entropy consistency.
• Section 6 presents a structure-preserving IMEX time integration scheme.
• Section 7 introduces a unified solver architecture applicable to both finite volume and meshless discretizations.
• Section 8 provides implementation details and pseudocode.
• Section 9 establishes convergence to entropy weak solutions.
• Section 10 presents numerical experiments demonstrating accuracy and robustness.
• Section 11 discusses implications and extensions.
• Section 12 concludes the paper.

2. Master Equation Formulation of Fluid Systems

2.1. Fundamental Perspective: Beyond the Continuum Hypothesis

Classical fluid dynamics is built upon the continuum hypothesis, in which the fluid is described as a continuous field governed by partial differential equations such as the Navier–Stokes equations [1,2,3,4,23,24]. Within this paradigm, numerical methods are constructed by discretizing these equations, followed by stabilization procedures to ensure robustness [5,6,7,8,9,10].

In contrast, the present work adopts a fundamentally different standpoint. The fluid is not assumed to be a continuous field a priori. Instead, it is modeled as a system of interacting discrete elements (nodes), whose collective dynamics give rise to macroscopic behavior.

The governing principle is thus inverted:

• Classical approach:
  Continuum PDE → Discretization → Stabilization
• Present approach:
  Discrete interaction system → Structure-preserving dynamics → Continuum limit

The design of the governing equation is guided by two fundamental physical principles:

• Conservation laws must be satisfied structurally, rather than as a consequence of discretization.
• The second law of thermodynamics must be embedded as a constraint on the time evolution.

Within this framework, continuum equations such as the Navier–Stokes system emerge only as limiting cases, rather than being postulated at the outset.

2.2. Node State Variables and Abstract Master Equation

Let the system be represented by a set of nodes indexed by $i\in \mathcal{V}$.

At each node, we define the state vector:

$$\begin{array}{ccc}& U_i=\left({\rho }_i,{\mathbf{m}}_i,e_i\right),& \end{array}$$

(1)

where

• ${\rho }_i$: mass density,
• ${\mathrm{m}}_i={\rho }_i{\mathrm{u}}_i$: momentum density,
• $e_i$: internal energy density.

The global state is given by $U=\{U_i{\}}_{i\in \mathcal{V}}$. The time evolution is defined through pairwise interactions:

$$\begin{array}{ccc}& {\displaystyle \frac{d{U}_i}{dt}}=\sum _{j\in \mathcal{N}(i)}(A_{ij}(U_i,U_j)+D_{ij}(U_i,U_j)),& \end{array}$$

(2)

where:

• $A_{ij}$: reversible (conservative) interaction,
• $D_{ij}$: irreversible (dissipative) interaction,
• $\mathcal{N}(i)$: neighborhood of node $i$.

The essential structural condition is antisymmetry:

$$\begin{array}{ccc}& A_{ij}(U_i,U_j)=-A_{ji}(U_j,U_i),D_{ij}(U_i,U_j)=-D_{ji}(U_j,U_i).& \end{array}$$

(3)

This implies that the evolution is entirely generated by interactions between nodes, rather than local self-dynamics.

2.3. Conservation Structure: Flux Representation

To ensure conservation laws, we express the evolution in flux form.

Let

• $V_i$: control volume associated with node $i$,
• ${\mathrm{a}}_{ij}$: oriented geometric vector associated with edge $\left(i,j\right)$,

satisfying

$$\begin{array}{ccc}& {\mathbf{a}}_{ij}=-{\mathbf{a}}_{ji}.& \end{array}$$

(4)

The master equation can be written as

$$\begin{array}{ccc}& {\displaystyle \frac{d}{dt}}(V_i{U}_i)=-\sum _{j\in \mathcal{N}(i)}{\mathbf{a}}_{ij}\cdot F_{ij},& \end{array}$$

(5)

where $F_{ij}$ is the numerical flux.

If the flux satisfies

$$\begin{array}{ccc}& F_{ij}=-F_{ji},& \end{array}$$

(6)

then summing over all nodes yields

$$\begin{array}{ccc}& \sum _i{\displaystyle \frac{d}{dt}}(V_i{U}_i)=0.& \end{array}$$

(7)

Therefore, the following quantities are exactly conserved:

• total mass,
• total momentum,
• total energy.

Importantly, conservation is not an emergent property but is guaranteed by construction.

2.4. Thermodynamic Structure: Entropy Production

To incorporate the second law, we introduce the total entropy:

$$\begin{array}{ccc}& S(U)=\sum _i{V}_{i}s(U_i),& \end{array}$$

(8)

where $s$ is a concave entropy function.

The entropy variables are defined as

$$\begin{array}{ccc}& v_i={\displaystyle \frac{\partial s}{\partial U}}(U_i).& \end{array}$$

(9)

The second law requires

$$\begin{array}{ccc}& {\displaystyle \frac{dS}{dt}}\ge 0.& \end{array}$$

(10)

To enforce this, the dissipative interaction is constructed as a gradient flow:

$$\begin{array}{ccc}& D_{ij}=M_{ij}(v_j-v_i),& \end{array}$$

(11)

where $M_{ij}$ is symmetric positive semidefinite:

$$\begin{array}{ccc}& M_{ij}=M_{ji},M_{ij}⪰0.& \end{array}$$

(12)

Then,

$$\begin{array}{ccc}& {\displaystyle \frac{dS}{dt}}={\displaystyle \frac{1}{2}}\sum _{i,j}(v_i-v_j{)}^T{M}_{ij}(v_i-v_j)\ge 0.& \end{array}$$

(13)

Thus, entropy production is ensured structurally.

The dynamics can be interpreted as a discrete realization of the GENERIC framework [19,20,21]:

• reversible part: antisymmetric operator
• irreversible part: symmetric positive operator

2.5. Unified Master Equation

Combining the above, the governing equations take the form:

Mass:

$$\begin{array}{ccc}& {\displaystyle \frac{d}{dt}}(V_i{\rho }_i)=-\sum _j{\mathbf{a}}_{ij}\cdot (\rho \mathbf{u}{)}_{ij},& \end{array}$$

(14)

Momentum:

$$\begin{array}{ccc}& {\displaystyle \frac{d}{dt}}(V_i{\mathrm{m}}_i)=-\sum _j{\mathrm{a}}_{ij}\cdot (\mathrm{m}\otimes \mathrm{u}+pI{)}_{ij}+\sum _j{D}_{ij}^{\left(m\right)},& \end{array}$$

(15)

Energy:

$$\begin{array}{ccc}& {\displaystyle \frac{d}{dt}}(V_i{E}_i)=-\sum _j{\mathrm{a}}_{ij}\cdot ((E+p)\mathrm{u}{)}_{ij}+\sum _j{D}_{ij}^{\left(e\right)},& \end{array}$$

(16)

where

$$\begin{array}{ccc}& E_i=e_i+{\displaystyle \frac{1}{2}}{\rho }_i\mid {\mathrm{u}}_i{\mid }^2.& \end{array}$$

(17)

All equations share a unified structure:

$$\begin{array}{ccc}& \mathrm{Time} \mathrm{evolution}=\mathrm{flux} \mathrm{exchange} \mathrm{across} \mathrm{edges}.& \end{array}$$

(18)

2.6. Navier–Stokes Equations as a Continuum Limit

Under the following assumptions:

• local interactions,
• dense node distribution,
• smooth fields,
• constitutive closure,
• Newtonian stress and Fourier heat conduction,

the discrete sums converge to differential operators:

$$\begin{array}{ccc}& {\displaystyle \frac{1}{V_i}}\sum _j{\mathbf{a}}_{ij}\cdot (\cdot ) ⟶ \nabla \cdot (\cdot ).& \end{array}$$

(19)

The master equation then reduces to:

• continuity equation,
• Navier–Stokes momentum equation,
• energy equation,
• entropy inequality,

i.e.,

$$\begin{array}{ccc}& {\partial }_{t}U+\nabla \cdot F(U)=\nabla \cdot G(U,\nabla U).& \end{array}$$

(20)

Thus, the Navier–Stokes equations are obtained as a continuum limit of the interaction system, rather than being fundamental.

2.7. Theoretical Significance

The master equation formulation provides a new foundation for fluid dynamics:

• Conservation laws are structurally guaranteed
• The second law is embedded at the operator level
• Continuum equations are derived, not assumed

This leads to a conceptual shift: Fluid dynamics is not fundamentally a theory of PDEs, but a theory of thermodynamically constrained interactions.

Consequently:

• Numerical methods can be constructed without ad hoc stabilization
• Different discretizations (FV, meshless) are unified under a single structure
• The validity range of continuum equations becomes explicitly identifiable

This framework opens a pathway toward a unified theory bridging:

• microscopic interactions,
• mesoscopic structures,
• macroscopic continuum laws.

3. Graph-Based Spatial Discretization Framework

3.1. Graph Representation of the Spatial Domain

In the proposed framework, the spatial domain $\mathsf{\Omega }\subset {\mathbb{R}}^d$ is not discretized through differential operators, but instead mapped into a graph structure:

$$\begin{array}{ccc}& \mathcal{G}=(\mathcal{V},\mathcal{E}),& \end{array}$$

(21)

where

• $\mathcal{V}$: set of nodes,
• $\mathcal{E}\subset \mathcal{V}\times \mathcal{V}$: set of edges.

Each node $i\in \mathcal{V}$ is associated with a control volume $V_i$, and each edge $(i,j)\in \mathcal{E}$ is assigned a geometric vector ${\mathbf{a}}_{ij}\in {\mathbb{R}}^d$, satisfying

$$\begin{array}{ccc}& {\mathbf{a}}_{ij}=-{\mathbf{a}}_{ji}.& \end{array}$$

(22)

This antisymmetry condition is the cornerstone of the discrete conservation structure.

This formulation generalizes classical discretizations:

• finite volume methods (FV) [5,6],
• meshless particle methods [16,17],

into a unified algebraic structure.

3.2. Discrete State Variables

At each node $i$, we define the conserved variables as

$$\begin{array}{ccc}& U_i=({\rho }_i,{\mathbf{m}}_i,E_i),& \end{array}$$

(23)

where

• ${\rho }_i$: density,
• ${\mathrm{m}}_i={\rho }_i{\mathrm{u}}_i$: momentum,
• $E_i=e_i+{\displaystyle \frac{1}{2}}{\rho }_i\mid {\mathrm{u}}_i{\mid }^2$: total energy.

The discrete state represents a volume-averaged quantity over $V_i$.

3.3. Discrete Conservation Law

Instead of approximating derivatives, we define the evolution through edge fluxes:

$$\begin{array}{ccc}& {\displaystyle \frac{d}{dt}}(V_i{U}_i)=-\sum _{j\in \mathcal{N}(i)}{\mathbf{a}}_{ij}\cdot F_{ij},& \end{array}$$

(24)

where

• $\mathcal{N}(i)$: neighbor set of node $i$,
• $F_{ij}$: numerical flux across edge $\left(i,j\right)$.

Equation (24) is the discrete analogue of

$$\begin{array}{ccc}& {\partial }_{t}U+\nabla \cdot F(U)=0,& \end{array}$$

(25)

but importantly, it is not derived by discretizing (25); rather, it is postulated as a structural law.

3.4. Antisymmetry Implies Conservation

The sufficient condition for conservation is:

$$\begin{array}{ccc}& F_{ij}=-F_{ji}.& \end{array}$$

(26)

Summing (24) over all nodes:

$$\begin{array}{ccc}& \sum _i{\displaystyle \frac{d}{dt}}(V_i{U}_i)=-\sum _{i,j}{\mathbf{a}}_{ij}\cdot F_{ij}.& \end{array}$$

(27)

Using antisymmetry:

$$\begin{array}{ccc}& {\mathbf{a}}_{ij}=-{\mathbf{a}}_{ji},F_{ij}=-F_{ji},& \end{array}$$

(28)

we obtain pairwise cancellation:

$$\begin{array}{ccc}& \sum _i{\displaystyle \frac{d}{dt}}(V_i{U}_i)=0.& \end{array}$$

(29)

Thus, the discrete system satisfies exact conservation of:

• mass,
• momentum,
• energy.

This property holds independently of mesh structure or dimensionality.

3.5. Correspondence with Finite Volume Methods

In finite volume methods, the domain is partitioned into control volumes, and fluxes are computed across cell interfaces [5,6]. For a face shared by cells $i$ and $j$, define

$$\begin{array}{ccc}& {\mathbf{a}}_{ij}=S_{ij}{\mathbf{n}}_{ij},& \end{array}$$

(30)

where

• $S_{ij}$: interface area,
• ${\mathrm{n}}_{ij}$: outward unit normal.

Then (24) becomes

$$\begin{array}{ccc}& {\displaystyle \frac{d}{dt}}(V_i{U}_i)=-\sum _j{S}_{ij}{\mathbf{n}}_{ij}\cdot F_{ij},& \end{array}$$

(31)

which is precisely the standard finite volume update formula. Therefore, the graph-based formulation strictly contains FV methods as a special case.

3.6. Correspondence with Meshless Methods

In meshless methods such as Smoothed Particle Hydrodynamics (SPH) [16] and generalized particle methods [17], no explicit interfaces exist.

We define the geometric vector using a kernel gradient:

$$\begin{array}{ccc}& {\mathbf{a}}_{ij}=V_i{V}_j\nabla W({\mathbf{x}}_i-{\mathbf{x}}_j),& \end{array}$$

(32)

where $W$ is a symmetric kernel.

If $W$ satisfies

$$\begin{array}{ccc}& \nabla W(\mathbf{r})=-\nabla W(-\mathbf{r}),& \end{array}$$

(33)

then

$$\begin{array}{ccc}& {\mathbf{a}}_{ij}=-{\mathbf{a}}_{ji},& \end{array}$$

(34)

and conservation is preserved. Thus, meshless methods are also naturally embedded within the same framework.

3.7. Discrete Divergence Operator

We define the discrete divergence operator as

$$\begin{array}{ccc}& {\left({\nabla }_h\cdot F\right)}_i={\displaystyle \frac{1}{V_i}}\sum _{j\in \mathcal{N}\left(i\right)}{\mathbf{a}}_{ij}\cdot F_{ij}.& \end{array}$$

(35)

Then (24) can be written as

$$\begin{array}{ccc}& {\displaystyle \frac{d{U}_i}{dt}}=-({\nabla }_h\cdot F{)}_i.& \end{array}$$

(36)

This definition is geometry-independent, with the only requirement being antisymmetry of ${\mathrm{a}}_{ij}$.

3.8. Consistency and Reproduction of Uniform States

For physical consistency, the scheme must preserve uniform states:

$$\begin{array}{ccc}& U_i=U_0\forall i.& \end{array}$$

(37)

In this case, consistency requires

$$\begin{array}{ccc}& \left({\nabla }_h\cdot F{)}_i=0.\right.& \end{array}$$

(38)

This holds if

$$\begin{array}{ccc}& \sum _{j\in \mathcal{N}(i)}{\mathbf{a}}_{ij}=0.& \end{array}$$

(39)

Condition (39) is satisfied:

• identically in finite volume methods (geometric closure),
• approximately in meshless methods via kernel normalization [16,17].

Thus, the scheme ensures exact preservation of constant states, a fundamental requirement for convergence.

3.9. Summary of the Discretization Framework

The spatial discretization introduced in this chapter is characterized by:

• Edge-based flux formulation (24)
• Antisymmetric geometry (22)
• Exact conservation via antisymmetry (29)
• Unified representation of FV and meshless methods (31), (32)
• Geometry-independent divergence operator (35)
• Consistency via geometric closure (39)

This formulation eliminates the need for explicit differential operators and instead constructs the numerical method directly from interaction structures.

In the next chapter, we build upon this framework to construct reversible (entropy-conserving) fluxes, forming the conservative backbone of the numerical scheme.

4. Construction of Reversible Fluxes

4.1. Design Requirements for Reversible Fluxes

Based on the graph-based conservation structure introduced in Chapter 3, the numerical flux $F_{ij}$ is decomposed into

$$\begin{array}{ccc}& F_{ij}=F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}+F_{ij}^{\mathrm{i}\mathrm{r}\mathrm{r}},& \end{array}$$

(40)

where $F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}$ represents the reversible (non-dissipative) component.

The reversible flux must satisfy the following fundamental requirements:

1.

Antisymmetry (conservation):

$$\begin{array}{ccc}& F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}=-F_{ji}^{\mathrm{r}\mathrm{e}\mathrm{v}},& \end{array}$$

(41)

2.

Consistency:

$$\begin{array}{ccc}& F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}=F(U_i)\mathrm{if} U_i=U_j,& \end{array}$$

(42)

3.

Entropy neutrality:

$$\begin{array}{ccc}& \left(v_i-v_j{)}^T{F}_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}=0\right.,& \end{array}$$

(43)

where $v=\partial s/\partial U$ denotes entropy variables [13,14,15,30].

These conditions ensure that the reversible flux:

• preserves conservation laws,
• does not produce entropy,
• is consistent with the underlying physical flux.

4.2. Euler Structure of Reversible Dynamics

In the continuum limit, the reversible part corresponds to the Euler equations [3,4]:

$$\begin{array}{ccc}& {\partial }_{t}U+\nabla \cdot F(U)=0.& \end{array}$$

(44)

For compressible flow,

$$\begin{array}{ccc}& F(U)=\left(\begin{array}{c}\rho u\\ \rho u\otimes u+pI\\ (E+p)u\end{array}\right).& \end{array}$$

(45)

These equations are:

• time-reversible,
• entropy-conserving,
• hyperbolic.

Therefore, the discrete reversible flux must approximate (45) in the continuum limit.

4.3. General Antisymmetric Construction

A general construction satisfying antisymmetry is:

$$\begin{array}{ccc}& F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}=\mathsf{\Phi }(U_i,U_j),\mathsf{\Phi }(U_i,U_j)=-\mathsf{\Phi }(U_j,U_i).& \end{array}$$

(46)

A sufficient representation is

$$\begin{array}{ccc}& F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}={\mathcal{A}}_{ij}\cdot \widehat{F}(U_i,U_j),& \end{array}$$

(47)

where

• ${\mathcal{A}}_{ij}={\mathrm{a}}_{ij}$(geometric vector),
• $\widehat{F}(U_i,U_j)$ is a symmetric numerical flux:

$$\begin{array}{ccc}& \widehat{F}(U_i,U_j)=\widehat{F}(U_j,U_i).& \end{array}$$

(48)

Then,

$$\begin{array}{ccc}& F_{ji}^{\mathrm{r}\mathrm{e}\mathrm{v}}=-F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}},& \end{array}$$

(49)

is automatically satisfied.

4.4. Central Flux Construction

The simplest consistent choice is the central flux:

$$\begin{array}{ccc}& \widehat{F}(U_i,U_j)={\displaystyle \frac{1}{2}}(F(U_i)+F(U_j)).& \end{array}$$

(50)

Thus,

$$\begin{array}{ccc}& F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}={\mathrm{a}}_{ij}\cdot {\displaystyle \frac{1}{2}}(F(U_i)+F(U_j)).& \end{array}$$

(51)

This scheme satisfies:

• consistency (42),
• antisymmetry (41),
• second-order accuracy in smooth regions.

However, central fluxes lack nonlinear stability near discontinuities [9,10].

4.5. Entropy-Conservative (EC) Fluxes

To ensure discrete entropy conservation, we adopt the entropy-conservative framework of Tadmor [15,30].

An entropy-conservative flux satisfies:

$$\begin{array}{ccc}& (v_i-v_j{)}^T{\widehat{F}}_{ij}=\psi (U_i)-\psi (U_j),& \end{array}$$

(52)

where $\psi$ is the entropy flux potential. A constructive form is:

$$\begin{array}{ccc}& {\widehat{F}}_{ij}^{\mathrm{E}\mathrm{C}}={\int }_0^{1}F(U(\theta ))\hspace{0.17em}d\theta ,& \end{array}$$

(53)

where $U(\theta )$ is a path connecting $U_i$ and $U_j$.

In practice, explicit formulas based on logarithmic averages are used [30]. The reversible flux is then defined as:

$$\begin{array}{ccc}& F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}={\mathrm{a}}_{ij}\cdot {\widehat{F}}_{ij}^{\mathrm{E}\mathrm{C}}.& \end{array}$$

(54)

This ensures:

$$\begin{array}{ccc}& {\displaystyle \frac{dS}{dt}}=0(\mathrm{reversible} \mathrm{part} \mathrm{only}).& \end{array}$$

(55)

4.6. Minimal Stabilization (Entropy-Stable Extension)

While EC fluxes are entropy-neutral, they are insufficient near shocks. Therefore, a minimal dissipation term is introduced:

$$\begin{array}{ccc}& {\widehat{F}}_{ij}={\widehat{F}}_{ij}^{\mathrm{E}\mathrm{C}}-{\displaystyle \frac{1}{2}}{\alpha }_{ij}{\mathbf{D}}_{ij}(U_j-U_i),& \end{array}$$

(56)

where

• ${\alpha }_{ij}$: maximum wave speed,
• ${\mathrm{D}}_{ij}$: symmetric positive definite matrix.

Thus,

$$\begin{array}{ccc}& F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}={\mathrm{a}}_{ij}\cdot {\widehat{F}}_{ij}.& \end{array}$$

(57)

The added term satisfies:

• conservation (due to antisymmetry),
• entropy dissipation:

$$\begin{array}{ccc}& (v_i-v_j{)}^T{\mathrm{D}}_{ij}(U_j-U_i)\ge 0.& \end{array}$$

(58)

This yields an entropy-stable scheme [13,14,15].

4.7. Consistency with the Continuum Limit

Let $h$ denote the characteristic spacing between nodes. Assuming smoothness of $U$, we expand:

$$\begin{array}{ccc}& U_j=U_i+\nabla U\cdot ({\mathbf{x}}_j-{\mathbf{x}}_i)+O(h^2).& \end{array}$$

(59)

Substituting into (51):

$$\begin{array}{ccc}& \sum _j{\mathbf{a}}_{ij}\cdot \widehat{F}(U_i,U_j)=\nabla \cdot F(U_i)+O(h^2).& \end{array}$$

(60)

Thus, the scheme is consistent with

$$\begin{array}{ccc}& {\partial }_{t}U+\nabla \cdot F(U)=0.& \end{array}$$

(61)

For EC fluxes, the same consistency holds due to path continuity.

4.8. Summary of Reversible Flux Construction

The reversible flux constructed in this chapter satisfies:

• Exact conservation: via antisymmetry (41)
• Consistency: via central/EC flux (50), (53)
• Entropy neutrality: via condition (52)
• Nonlinear stability: via minimal stabilization (56)
• Continuum consistency: via asymptotic expansion (60)

This construction provides a structure-preserving discrete analogue of Euler dynamics, forming the conservative backbone of the numerical scheme.

In the next chapter, we introduce the irreversible flux, completing the thermodynamically consistent formulation.

5. Construction of Irreversible Fluxes via Entropy Gradient Flow

5.1. Objective

In Chapter 4, the reversible flux $F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}$ was constructed to ensure conservation and entropy neutrality. The purpose of this chapter is to construct the irreversible flux $F_{ij}^{\mathrm{i}\mathrm{r}\mathrm{r}}$ such that the second law of thermodynamics is satisfied at the discrete level:

$$\begin{array}{ccc}& {\displaystyle \frac{dS}{dt}}\ge 0.& \end{array}$$

(62)

The key principle is to formulate dissipation as an entropy gradient flow, ensuring that entropy production is structurally guaranteed rather than numerically enforced.

5.2. Definition of Discrete Entropy

Let the entropy density be given by a concave function $s(U)$. The total entropy of the system is defined as:

$$\begin{array}{ccc}& S(U)=\sum _{i\in \mathcal{V}}{V}_i\hspace{0.17em}s(U_i).& \end{array}$$

(63)

The entropy variables are defined as:

$$\begin{array}{ccc}& v_i={\displaystyle \frac{\partial s}{\partial U}}(U_i).& \end{array}$$

(64)

These variables play a central role in thermodynamically consistent formulations [19,20,21,29].

5.3. Discrete Entropy Evolution

Starting from the discrete conservation law:

$$\begin{array}{ccc}& {\displaystyle \frac{d}{dt}}(V_i{U}_i)=-\sum _j{\mathbf{a}}_{ij}\cdot (F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}+F_{ij}^{\mathrm{i}\mathrm{r}\mathrm{r}}),& \end{array}$$

(65)

multiplying by $v_i^T$ and summing over $i$, we obtain:

$$\begin{array}{ccc}& {\displaystyle \frac{dS}{dt}}=-\sum _{i,j}{v}_i^T{\mathbf{a}}_{ij}\cdot F_{ij}^{\mathrm{i}\mathrm{r}\mathrm{r}},& \end{array}$$

(66)

where the reversible contribution vanishes due to entropy neutrality (Chapter 4). Using antisymmetry:

$$\begin{array}{ccc}& {\displaystyle \frac{dS}{dt}}={\displaystyle \frac{1}{2}}\sum _{i,j}(v_i-v_j{)}^T{\mathbf{a}}_{ij}\cdot F_{ij}^{\mathrm{i}\mathrm{r}\mathrm{r}}.& \end{array}$$

(67)

Thus, entropy production is determined entirely by the irreversible flux.

5.4. Gradient Flow Construction

To guarantee non-negative entropy production, we define:

$$\begin{array}{ccc}& F_{ij}^{\mathrm{i}\mathrm{r}\mathrm{r}}={\mathbf{a}}_{ij}\cdot (M_{ij}(v_j-v_i)),& \end{array}$$

(68)

where $M_{ij}$ is a symmetric positive semidefinite operator:

$$\begin{array}{ccc}& M_{ij}=M_{ji},M_{ij}⪰0.& \end{array}$$

(69)

Substituting into (67):

$$\begin{array}{ccc}& {\displaystyle \frac{dS}{dt}}={\displaystyle \frac{1}{2}}\sum _{i,j}(v_i-v_j{)}^T{M}_{ij}(v_i-v_j)\ge 0.& \end{array}$$

(70)

Thus, the second law is satisfied by construction.

This formulation is a discrete analogue of Onsager’s reciprocal relations [28] and GENERIC theory [19,20,21].

5.5. Connection to Physical Dissipation

In continuum fluid dynamics, irreversible effects arise from:

• viscous stress,
• heat conduction.

The Navier–Stokes constitutive laws are:

$$\begin{array}{ccc}& \tau =2\mu D(\mathrm{u}),\mathrm{q}=-k\nabla T,& \end{array}$$

(71)

where $D(\mathrm{u})$ is the strain-rate tensor.

To recover these effects, we decompose:

$$\begin{array}{ccc}& M_{ij}={\mu }_{ij}{\mathbf{G}}_{ij}^{\left(m\right)}+k_{ij}{\mathbf{G}}_{ij}^{\left(e\right)},& \end{array}$$

(72)

where:

• ${\mu }_{ij}$: viscosity coefficient,
• $k_{ij}$: thermal conductivity,
• ${\mathrm{G}}_{ij}$: geometric tensors.

Under refinement,

$$\begin{array}{ccc}& \sum _j{\mathrm{a}}_{ij}\cdot M_{ij}(v_j-v_i) ⟶ \nabla \cdot \tau +\nabla \cdot \mathrm{q}.& \end{array}$$

(73)

Thus, the gradient-flow formulation is consistent with classical dissipation [3,4].

5.6. Implementation in FV and Meshless Frameworks

(A) Finite Volume Implementation

For a face between cells $i$ and $j$, we define:

$$\begin{array}{ccc}& M_{ij}={\displaystyle \frac{S_{ij}}{d_{ij}}}{\mathbf{K}}_{ij},& \end{array}$$

(74)

where:

• $S_{ij}$: face area,
• $d_{ij}$: distance between centroids,
• ${\mathbf{K}}_{ij}$: diffusion matrix.

Then,

$$\begin{array}{ccc}& F_{ij}^{\mathrm{i}\mathrm{r}\mathrm{r}}=S_{ij}{\mathbf{n}}_{ij}\cdot {\mathbf{K}}_{ij}(v_j-v_i).& \end{array}$$

(75)

(B) Meshless Implementation

Using kernel gradients:

$$\begin{array}{ccc}& M_{ij}=V_i{V}_j\hspace{0.17em}{\mathbf{K}}_{ij}\hspace{0.17em}\nabla W({\mathbf{x}}_i-{\mathbf{x}}_j),& \end{array}$$

(76)

which preserves antisymmetry and consistency.

5.7. Separation of Physical and Numerical Dissipation

The operator $M_{ij}$ can be decomposed as:

$$\begin{array}{ccc}& M_{ij}=M_{ij}^{\mathrm{p}\mathrm{h}\mathrm{y}\mathrm{s}}+M_{ij}^{\mathrm{n}\mathrm{u}\mathrm{m}}.& \end{array}$$

(77)

• $M_{ij}^{\mathrm{p}\mathrm{h}\mathrm{y}\mathrm{s}}$: physical viscosity,
• $M_{ij}^{\mathrm{n}\mathrm{u}\mathrm{m}}$: numerical stabilization.

A practical construction is:

$$\begin{array}{ccc}& M_{ij}^{\mathrm{n}\mathrm{u}\mathrm{m}}={\beta }_{ij}{R}_{i}I,& \end{array}$$

(78)

where:

• $R_i$: entropy residual,
• ${\beta }_{ij}$: scaling coefficient.

This yields an adaptive dissipation mechanism:

• small in smooth regions,
• large near shocks or under-resolved scales.

5.8. Lyapunov Structure and Stability

From (70), the entropy satisfies:

$$\begin{array}{ccc}& {\displaystyle \frac{dS}{dt}}\ge 0.& \end{array}$$

(79)

In a closed system,

$$\begin{array}{ccc}& {\displaystyle \frac{dS}{dt}}=0⟺v_i=v_j\forall i,j.& \end{array}$$

(80)

Thus, equilibrium corresponds to uniform entropy variables. Furthermore, combining with energy conservation:

$$\begin{array}{ccc}& {\displaystyle \frac{d}{dt}}{E}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}\le 0,& \end{array}$$

(81)

the system admits a Lyapunov functional:

$$\begin{array}{ccc}& \mathcal{L}(U)=-S(U).& \end{array}$$

(82)

Hence:

• entropy increases monotonically,
• the system evolves toward equilibrium,
• numerical stability is guaranteed.

5.9. Summary of Irreversible Flux Construction

The irreversible flux constructed in this chapter satisfies:

• Second law (entropy production): (70)
• Thermodynamic consistency: gradient flow structure (68)
• Physical fidelity: recovery of viscous and thermal effects (73)
• Unified implementation: FV and meshless (74)–(76)
• Adaptive stabilization: entropy-based dissipation (78)
• Lyapunov stability: monotonic entropy increase (79)

This completes the thermodynamically consistent decomposition:

$$\begin{array}{ccc}& F_{ij}=F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}+F_{ij}^{\mathrm{i}\mathrm{r}\mathrm{r}},& \end{array}$$

(83)

which ensures that the numerical scheme simultaneously satisfies:

• conservation laws,
• entropy inequality.

6. Time Integration: IMEX Thermodynamically Consistent Scheme

6.1. Semi-Discrete Evolution System

From Chapters 4 and 5, the graph-based spatial discretization yields a semi-discrete system of the form

$$\begin{array}{ccc}& {\displaystyle \frac{d{U}_i}{dt}}={\mathcal{R}}_i(U)+{\mathcal{D}}_i(U),& \end{array}$$

(84)

where

• ${\mathcal{R}}_i(U)$: reversible contribution generated by entropy-neutral fluxes,
• ${\mathcal{D}}_i(U)$: irreversible contribution generated by entropy-gradient flows.

Equivalently, in vector form,

$$\begin{array}{ccc}& {\displaystyle \frac{dU}{dt}}=\mathcal{R}(U)+\mathcal{D}(U).& \end{array}$$

(85)

This decomposition is not merely algebraic. It reflects two fundamentally different physical mechanisms:

• the reversible part transports conserved quantities without entropy production,
• the irreversible part dissipates gradients and generates entropy [19,20,21,28,29].

Because these two contributions generally operate on different time scales, a unified explicit treatment is inefficient and often unstable.

6.2. Reversible–Irreversible Splitting

The reversible operator $\mathcal{R}(U)$ is associated with convective or wave-like dynamics. Its stability is governed by a hyperbolic CFL condition of the form

$$\begin{array}{ccc}& \mathsf{\Delta }t\le C_{\mathrm{C}\mathrm{F}\mathrm{L}}{\displaystyle \frac{h}{{\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}}},& \end{array}$$

(86)

where $h$ is a characteristic node spacing and ${\lambda }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum signal speed [5,10].

By contrast, the irreversible operator $\mathcal{D}(U)$ is diffusive and obeys a parabolic restriction

$$\begin{array}{ccc}& \mathsf{\Delta }t\le C_{\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}}{\displaystyle \frac{h^2}{{\nu }_{\mathrm{m}\mathrm{a}\mathrm{x}}}},& \end{array}$$

(87)

where ${\nu }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ denotes the effective diffusion strength. For high-Reynolds-number or strongly anisotropic problems, the diffusive restriction can become excessively severe if treated explicitly. This naturally motivates an implicit–explicit (IMEX) strategy:

• treat $\mathcal{R}(U)$ explicitly,
• treat $\mathcal{D}(U)$ implicitly.

6.3. IMEX Time Integration Design

Let $U^n$ denote the discrete state at time $t^n$, with time step $\mathsf{\Delta }t$.

The first-order IMEX update is defined by

$$\begin{array}{ccc}& {\displaystyle \frac{U^{n+1}-U^n}{\mathsf{\Delta }t}}=\mathcal{R}(U^n)+\mathcal{D}(U^{n+1}).& \end{array}$$

(88)

Equivalently,

$$\begin{array}{ccc}& U^{n+1}-\mathsf{\Delta }t\hspace{0.17em}\mathcal{D}(U^{n+1})=U^n+\mathsf{\Delta }t\hspace{0.17em}\mathcal{R}(U^n).& \end{array}$$

(89)

This construction is consistent with the physical interpretation:

• the reversible part advances the transport structure explicitly,
• the irreversible part relaxes the system toward thermodynamic consistency implicitly.

At the node level, the update reads

$$\begin{array}{ccc}& U_i^{n+1}=U_i^n+\mathsf{\Delta }t\hspace{0.17em}{\mathcal{R}}_i(U^n)+\mathsf{\Delta }t\hspace{0.17em}{\mathcal{D}}_i(U^{n+1}).& \end{array}$$

(90)

6.4. Discrete Conservation Property

We now verify that the IMEX scheme preserves the global conservation laws.

Summing (90) over all nodes $i$, we obtain

$$\begin{array}{ccc}& \sum _i{V}_i{U}_i^{n+1}=\sum _i{V}_i{U}_i^n+\mathsf{\Delta }t\sum _i{V}_i{\mathcal{R}}_i(U^n)+\mathsf{\Delta }t\sum _i{V}_i{\mathcal{D}}_i(U^{n+1}).& \end{array}$$

(91)

From the antisymmetric flux construction developed in Chapters 4 and 5,

$$\begin{array}{ccc}& \sum _i{V}_i{\mathcal{R}}_i(U)=0,\sum _i{V}_i{\mathcal{D}}_i(U)=0.& \end{array}$$

(92)

Therefore,

$$\begin{array}{ccc}& \sum _i{V}_i{U}_i^{n+1}=\sum _i{V}_i{U}_i^n.& \end{array}$$

(93)

Hence, the IMEX update preserves the discrete global conservation of mass, momentum, and total energy. This is a key advantage of the present design: conservation is not recovered approximately, but remains exact under time discretization.

6.5. Discrete Entropy Monotonicity

The second law requires that entropy should not decrease in time.

Let the discrete entropy be

$$\begin{array}{ccc}& S(U^n)=\sum _i{V}_{i}s(U_i^n).& \end{array}$$

(94)

Using the convexity/concavity structure of the entropy pair and the entropy-gradient flow form of $\mathcal{D}$, the implicit step satisfies

$$\begin{array}{ccc}& S(U^{n+1})-S(U^n)\ge \mathsf{\Delta }t\hspace{0.17em}\sum _i{V}_i\hspace{0.17em}{v}_i^{n+1\hspace{0.17em}T}{\mathcal{D}}_i(U^{n+1}),& \end{array}$$

(95)

where $v_i^{n+1}=\partial s/\partial U(U_i^{n+1})$.

From Chapter 5,

$$\begin{array}{ccc}& \sum _i{V}_i\hspace{0.17em}{v}_i^T{\mathcal{D}}_i(U)={\displaystyle \frac{1}{2}}\sum _{i,j}(v_i-v_j{)}^T{M}_{ij}(v_i-v_j)\ge 0.& \end{array}$$

(96)

Thus,

$$\begin{array}{ccc}& S(U^{n+1})\ge S(U^n).& \end{array}$$

(97)

This shows that the implicit treatment of the dissipative term preserves the entropy monotonicity required by thermodynamics.

6.6. Linearized Implicit Solver

In practice, the implicit step in (89) requires the solution of a nonlinear system. To obtain an efficient solver, we linearize the dissipative operator around a known state $\tilde{U}$ (typically $U^n$ or an explicit predictor):

$$\begin{array}{ccc}& \mathcal{D}(U^{n+1})\approx \mathcal{D}(\tilde{U})+\mathbf{L}(\tilde{U})(U^{n+1}-\tilde{U}),& \end{array}$$

(98)

where $\mathbf{L}$ is the Jacobian of $\mathcal{D}$.

Substituting into (89), we obtain a linear system

$$\begin{array}{ccc}& (I-\mathsf{\Delta }t\hspace{0.17em}\mathbf{L}(\tilde{U}))U^{n+1}=U^n+\mathsf{\Delta }t\hspace{0.17em}\mathcal{R}(U^n)+\mathsf{\Delta }t\hspace{0.17em}b(\tilde{U}),& \end{array}$$

(99)

where $b(\tilde{U})$ collects lower-order terms.

Because $\mathcal{D}$ arises from a symmetric positive semidefinite gradient-flow operator, the matrix $I-\mathsf{\Delta }t\hspace{0.17em}\mathbf{L}$is typically symmetric positive definite, or at least coercive under standard regularization assumptions. Hence, efficient solvers can be used:

• conjugate gradient (CG),
• preconditioned conjugate gradient (PCG),
• multigrid methods [5,6,7,33].

This is a significant practical benefit of the thermodynamic design: the irreversible step leads naturally to favorable linear algebraic structure.

6.7. Higher-Order IMEX Extension

The first-order IMEX scheme can be extended to second order or higher using IMEX Runge–Kutta methods [37].

A typical second-order predictor–corrector form is:

Predictor:

$$\begin{array}{ccc}& U^{\left(1\right)}=U^n+\mathsf{\Delta }t\hspace{0.17em}\mathcal{R}(U^n)+\mathsf{\Delta }t\hspace{0.17em}\mathcal{D}(U^{\left(1\right)}),& \end{array}$$

(100)

Corrector:

$$\begin{array}{ccc}& U^{n+1}={\displaystyle \frac{1}{2}}{U}^n+{\displaystyle \frac{1}{2}}\left(U^{\left(1\right)}+\mathsf{\Delta }t\hspace{0.17em}\mathcal{R}(U^{\left(1\right)})+\mathsf{\Delta }t\hspace{0.17em}\mathcal{D}(U^{n+1})\right).& \end{array}$$

(101)

More generally, an IMEX-RK method can be written as

$$\begin{array}{ccc}& U^{\left(k\right)}=U^n+\mathsf{\Delta }t\sum _{\mathcal{l}<k}{a}_{k\mathcal{l}}^{\mathrm{E}}\mathcal{R}(U^{\left(\mathcal{l}\right)})+\mathsf{\Delta }t\sum _{\mathcal{l}\le k}{a}_{k\mathcal{l}}^{\mathrm{I}}\mathcal{D}(U^{\left(\mathcal{l}\right)}),& \end{array}$$

(102)

followed by

$$\begin{array}{ccc}& U^{n+1}=U^n+\mathsf{\Delta }t\sum _k{b}_k^{\mathrm{E}}\mathcal{R}(U^{\left(k\right)})+\mathsf{\Delta }t\sum _k{b}_k^{\mathrm{I}}\mathcal{D}(U^{\left(k\right)}).& \end{array}$$

(103)

To preserve entropy monotonicity and conservation, the coefficients must be chosen so that the implicit part remains contractive and the explicit part preserves antisymmetric cancellations.

6.8. Energy Stability

The reversible part of the dynamics is designed to be entropy-neutral and, in the absence of stabilization, energy-conservative. The irreversible part dissipates gradients and drives the system toward equilibrium.

Accordingly, the total discrete free-energy-like functional

$$\begin{array}{ccc}& \mathcal{E}(U)=E_{\mathrm{t}\mathrm{o}\mathrm{t}}(U)-\theta S(U),& \end{array}$$

(104)

with $\theta >0$ a reference thermodynamic weight, satisfies

$$\begin{array}{ccc}& \mathcal{E}(U^{n+1})\le \mathcal{E}(U^n),& \end{array}$$

(105)

provided the implicit dissipative step is solved consistently and the explicit reversible part satisfies the discrete conservation constraints.

In particular, when the reversible flux is purely entropy-conservative and the dissipative operator is positive semidefinite, one obtains the discrete stability estimate

$$\begin{array}{ccc}& \parallel U^{n+1}{\parallel }_{\mathcal{H}}^2\le \parallel U^n{\parallel }_{\mathcal{H}}^2+C\mathsf{\Delta }t\hspace{0.17em}\parallel U^n{\parallel }_{\mathcal{H}}^2,& \end{array}$$

(106)

for a suitable energy norm $\parallel \cdot {\parallel }_{\mathcal{H}}$, with $C$ independent of the mesh size.

Thus, the IMEX method combines:

• hyperbolic efficiency for transport,
• unconditional or weakly conditioned stability for dissipation,
• compatibility with discrete thermodynamic structure.

6.9. Summary of the IMEX Thermodynamically Consistent Scheme

The IMEX time integration scheme developed in this chapter possesses the following properties:

• Structure-aware splitting: reversible and irreversible dynamics are treated differently according to their physical role.
• Exact conservation: global invariants remain preserved under time stepping by virtue of antisymmetric flux cancellation.
• Entropy monotonicity: the implicit dissipative step guarantees non-decreasing entropy.
• Efficient solvability: the irreversible operator yields a symmetric positive linear system after linearization.
• High-order extensibility: the framework is compatible with IMEX-RK generalizations.
• Energy stability: the combined method admits a discrete Lyapunov structure.

The resulting scheme is therefore not merely a practical integrator, but a thermodynamically structured time discretization consistent with the design principles of the entire framework.

In the next chapter, we translate this formulation into a unified solver architecture applicable to both finite volume and meshless implementations.

7. Unified Solver Architecture (FV / Meshless)

7.1. Unified Design Philosophy

The formulation developed in the preceding chapters is fundamentally operator-based, rather than grid-based.

The governing semi-discrete system can be written as

$$\begin{array}{ccc}& {\displaystyle \frac{dU}{dt}}=\mathcal{R}(U)+\mathcal{D}(U),& \end{array}$$

(107)

where:

• $\mathcal{R}$: antisymmetric reversible operator,
• $\mathcal{D}$: symmetric positive semidefinite dissipative operator.

The essential design principle is: Geometry is interchangeable, but structure is invariant.

That is:

• the geometric representation (mesh vs particles) may vary,
• the interaction structure (antisymmetry + entropy gradient) is fixed.

This separation enables a unified solver applicable to both finite volume and meshless discretizations [5,16,17,33,34].

7.2. Data Structure

The solver is constructed from three fundamental components:

(A) Node Data

Each node $i$ stores:

$$\begin{array}{ccc}& {\mathcal{N}}_i=\{{\mathbf{x}}_i,V_i,U_i,v_i\},& \end{array}$$

(108)

where:

• ${\mathbf{x}}_i$: position,
• $V_i$: control volume,
• $U_i$: conserved variables,
• $v_i$: entropy variables.

(B) Edge Data

Each edge $\left(i,j\right)$ stores:

$$\begin{array}{ccc}& {\mathcal{E}}_{ij}=\{{\mathbf{a}}_{ij},d_{ij},M_{ij}\},& \end{array}$$

(109)

with:

• ${\mathbf{a}}_{ij}=-{\mathbf{a}}_{ji}$,
  $d_{ij}=\parallel {\mathbf{x}}_i-{\mathbf{x}}_j\parallel$,
• $M_{ij}$: diffusion operator.

(C) Operator Representation

The discrete operators are:

$$\begin{array}{ccc}& {\mathcal{R}}_i(U)=-{\displaystyle \frac{1}{V_i}}\sum _j{\mathbf{a}}_{ij}\cdot F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}},& \end{array}$$

(110)

$$\begin{array}{ccc}& {\mathcal{D}}_i(U)=-{\displaystyle \frac{1}{V_i}}\sum _j{\mathbf{a}}_{ij}\cdot F_{ij}^{\mathrm{i}\mathrm{r}\mathrm{r}}.& \end{array}$$

(111)

7.3. Separation of Geometric Differences

The only difference between FV and meshless methods lies in the construction of ${\mathbf{a}}_{ij}$ and neighborhood relations.

Finite Volume (FV):

$$\begin{array}{ccc}& {\mathbf{a}}_{ij}=S_{ij}{\mathbf{n}}_{ij},& \end{array}$$

(112)

where $S_{ij}$ is face area.

Meshless:

$$\begin{array}{ccc}& {\mathbf{a}}_{ij}=V_i{V}_j\nabla W({\mathbf{x}}_i-{\mathbf{x}}_j).& \end{array}$$

(113)

Thus, the solver architecture can be written abstractly as:

$$\begin{array}{ccc}& \mathrm{GeometryBuilder} \to \{{\mathbf{a}}_{ij},\mathcal{N}(i)\}.& \end{array}$$

(114)

All subsequent computations are identical.

7.4. Solver Design

The time evolution follows the IMEX scheme:

$$\begin{array}{ccc}& U^{n+1}-\mathsf{\Delta }t\hspace{0.17em}\mathcal{D}(U^{n+1})=U^n+\mathsf{\Delta }t\hspace{0.17em}\mathcal{R}(U^n).& \end{array}$$

(115)

The solver is organized into four stages:

• Reversible flux computation
• Explicit update
• Entropy variable evaluation
• Implicit dissipation solve

7.5. Algorithmic Structure

Step 1: Reversible Flux

For each edge $\left(i,j\right)$,

$$\begin{array}{ccc}& F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}={\mathbf{a}}_{ij}\cdot {\widehat{F}}_{ij},& \end{array}$$

(116)

with EC or stabilized flux [15,30].

Step 2: Explicit Update

$$\begin{array}{ccc}& U_i^{*}=U_i^n-{\displaystyle \frac{\mathsf{\Delta }t}{V_i}}\sum _j{\mathbf{a}}_{ij}\cdot F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}.& \end{array}$$

(117)

Step 3: Entropy Variables

$$\begin{array}{ccc}& v_i^{*}={\displaystyle \frac{\partial s}{\partial U}}(U_i^{*}).& \end{array}$$

(118)

Step 4: Implicit Dissipation

Solve:

$$\begin{array}{ccc}& U^{n+1}-\mathsf{\Delta }t\hspace{0.17em}\mathcal{D}(U^{n+1})=U^{*}.& \end{array}$$

(119)

Matrix Form

$$\begin{array}{ccc}& \left(I−\mathsf{\Delta }t\hspace{0.17em}\mathbf{L}\right)U^{n+1}=U^{*}.& \end{array}$$

(120)

7.6. Adaptive Dissipation

To handle shocks and unresolved turbulence, we define an entropy residual:

$$\begin{array}{ccc}& R_i={\displaystyle \frac{S(U_i^{n+1})-S(U_i^n)}{\mathsf{\Delta }t}}-{\mathsf{\Pi }}_i^{\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{o}\mathrm{r}\mathrm{y}},& \end{array}$$

(121)

where ${\mathsf{\Pi }}_i^{\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{o}\mathrm{r}\mathrm{y}}$ is theoretical entropy production.

Then,

$$\begin{array}{ccc}& M_{ij}^{\mathrm{n}\mathrm{u}\mathrm{m}}={\beta }_{ij}{R}_{i}I.& \end{array}$$

(122)

This yields:

• minimal dissipation in smooth regions,
• enhanced stabilization near discontinuities.

7.7. Parallelization and Scalability

The solver is inherently edge-based:

$$\begin{array}{ccc}& \mathcal{O}(\mid \mathcal{E}\mid )\mathrm{operations}.& \end{array}$$

(123)

Properties:

• local interactions → no global coupling (explicit step),
• sparse SPD system (implicit step),
• perfect suitability for:

  ∘

  MPI domain decomposition,

  ∘

  GPU parallelization,

  ∘

  CSR sparse storage.

The implicit solve satisfies:

$$\begin{array}{ccc}& \mathcal{O}(\mid \mathcal{E}\mid \mathrm{l}\mathrm{o}\mathrm{g}N)(\mathrm{multigrid}),& \end{array}$$

(124)

or

$$\begin{array}{ccc}& \mathcal{O}(\mid \mathcal{E}\mid )(\mathrm{optimal} \mathrm{scaling}).& \end{array}$$

(125)

7.8. Theoretical Consistency

The unified solver satisfies the following properties:

(i) Conservation

$$\begin{array}{ccc}& \sum _i{V}_i{U}_i^{n+1}=\sum _i{V}_i{U}_i^n.& \end{array}$$

(126)

(ii) Entropy Inequality

$$\begin{array}{ccc}& S(U^{n+1})\ge S(U^n).& \end{array}$$

(127)

(iii) Continuum Limit

$$\begin{array}{ccc}& \mathcal{R}(U)+\mathcal{D}(U) ⟶ \nabla \cdot F(U)-\nabla \cdot G(U,\nabla U).& \end{array}$$

(128)

(iv) Method Unification

FV and meshless methods satisfy the same operator equation:

$$\begin{array}{ccc}& {\displaystyle \frac{dU}{dt}}={\mathcal{L}}_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{p}\mathrm{h}}(U).& \end{array}$$

(129)

(v) Stability

• entropy monotonicity,
• Lyapunov structure,
• IMEX stability.

7.9. Summary of Unified Solver

The unified solver architecture achieves:

• method unification (FV + meshless),
• structure preservation (conservation + entropy),
• high scalability (edge-parallel),
• adaptivity (entropy-based dissipation),
• thermodynamic consistency.

The central principle can be summarized as:

$$\begin{array}{ccc}& \mathrm{Structure} \mathrm{first}, \mathrm{geometry} \sec \mathrm{ond}.& \end{array}$$

(130)

8. Algorithm Specification and Pseudocode

8.1. Overall Solver Structure

The numerical method developed in the preceding chapters leads to the semi-discrete evolution equation

$$\begin{array}{ccc}& {\displaystyle \frac{dU}{dt}}=\mathcal{R}(U)+\mathcal{D}(U),& \end{array}$$

(131)

with IMEX time integration

$$\begin{array}{ccc}& U^{n+1}-\mathsf{\Delta }t\hspace{0.17em}\mathcal{D}(U^{n+1})=U^n+\mathsf{\Delta }t\hspace{0.17em}\mathcal{R}(U^n).& \end{array}$$

(132)

The solver is organized into:

• Initialization phase
• Geometry construction
• Flux evaluation
• Implicit solve
• Adaptive dissipation

8.2. Initialization

The initial setup defines nodes, state variables, and physical parameters:

$$\begin{array}{ccc}& U_i^0=U_0({\mathbf{x}}_i),& \end{array}$$

(133)

where $U_0$ is the initial condition.

The solver initializes:

• node positions ${\mathbf{x}}_i$,
• volumes $V_i$,
• conserved variables $U_i$,
• material properties (viscosity, conductivity).

Initialization Procedure

InitializeNodes()

InitializeState(U)

SetPhysicalParameters()

8.3. GeometryBuilder

The geometry construction step defines the graph structure:

$$\begin{array}{ccc}& \mathcal{G}=(\mathcal{V},\mathcal{E}),& \end{array}$$

(134)

and computes:

$$\begin{array}{ccc}& {\mathbf{a}}_{ij}=-{\mathbf{a}}_{ji}.& \end{array}$$

(135)

(A) Finite Volume Version

$$\begin{array}{ccc}& {\mathbf{a}}_{ij}=S_{ij}{\mathbf{n}}_{ij},& \end{array}$$

(136)

for each cell i:

 compute V[i]

for each face (i,j):

 a_ij = S_ij * n_ij

(B) Meshless Version

$$\begin{array}{ccc}& {\mathbf{a}}_{ij}=V_i{V}_j\nabla W({\mathbf{x}}_i-{\mathbf{x}}_j).& \end{array}$$

(137)

for each particle i:

 find neighbors N(i)

for each pair (i,j):

 gradW = ∇W(x_i - x_j)

 a_ij = V[i]*V[j]*gradW

8.4. Flux Computation

Reversible Flux

$$\begin{array}{ccc}& F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}={\mathbf{a}}_{ij}\cdot {\widehat{F}}_{ij},& \end{array}$$

(138)

where ${\widehat{F}}_{ij}$ is EC or stabilized flux [15,30].

Algorithm

for each edge (i,j):

 U_avg = EntropyConsistentAverage(U[i], U[j])

 F_ec = EulerFlux(U_avg)

 α = MaxCharacteristicSpeed(U[i], U[j])

 F_rev = F_ec - 0.5 * α * (U[j] - U[i])

 FluxRev[i] += dot(a_ij, F_rev)

 FluxRev[j] -= dot(a_ij, F_rev)

8.5. Explicit Update

$$\begin{array}{ccc}& U_i^{*}=U_i^n-{\displaystyle \frac{\mathsf{\Delta }t}{V_i}}\sum _j{\mathbf{a}}_{ij}\cdot F_{ij}^{\mathrm{r}\mathrm{e}\mathrm{v}}.& \end{array}$$

(139)

for each node i:

 U_star[i] = U[i] - dt/V[i] * FluxRev[i]

8.6. Entropy Variables

$$\begin{array}{ccc}& v_i^{*}={\displaystyle \frac{\partial s}{\partial U}}(U_i^{*}).& \end{array}$$

(140)

for each node i:

 v[i] = EntropyVariable(U_star[i])

8.7. Implicit Dissipation Solver

We solve:

$$\begin{array}{ccc}& U^{n+1}-\mathsf{\Delta }t\hspace{0.17em}\mathcal{D}(U^{n+1})=U^{*}.& \end{array}$$

(141)

Linearization yields:

$$\begin{array}{ccc}& \left(I−\mathsf{\Delta }t\hspace{0.17em}\mathbf{L}\right)U^{n+1}=U^{*}.& \end{array}$$

(142)

Matrix Assembly

$$\begin{array}{ccc}& {\mathbf{L}}_{ij}\sim {\displaystyle \frac{1}{V_i}}{\mathbf{a}}_{ij}\cdot M_{ij}.& \end{array}$$

(143)

for each edge (i,j):

 M_ij = DiffusionMatrix(i,j)

 assemble A[i,i] += ...

 assemble A[i,j] += ...

Linear Solver

$$\begin{array}{ccc}& A{U}^{n+1}=U^{*}.& \end{array}$$

(144)

Recommended solvers:

• PCG (preconditioned conjugate gradient) [49]
• multigrid [48]

U_new = PCG(A, U_star)

8.8. Complete Pseudocode

InitializeNodes()

InitializeState(U)

GeometryBuilder()

while time < T_end:

 # --- Reversible flux ---

 for edges:

 compute F_rev

  accumulate FluxRev

 # --- Explicit update ---

 for nodes:

  U_star = U - dt/V * FluxRev

 # --- Entropy variables ---

 for nodes:

 v = dS/dU(U_star)

 # --- Implicit dissipation ---

 build matrix A

 solve A U_new = U_star

 # --- Adaptive dissipation ---

 for nodes:

 compute entropy residual

  adjust diffusion

 U = U_new

 time += dt

8.9. Adaptive Dissipation

Entropy residual:

$$\begin{array}{ccc}& R_i={\displaystyle \frac{S(U_i^{n+1})-S(U_i^n)}{\mathsf{\Delta }t}}-{\mathsf{\Pi }}_i^{\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{o}\mathrm{r}\mathrm{y}},& \end{array}$$

(145)

Then,

$$\begin{array}{ccc}& M_{ij}^{\mathrm{n}\mathrm{u}\mathrm{m}}={\beta }_{ij}{R}_{i}I.& \end{array}$$

(146)

This ensures:

• minimal dissipation in smooth regions,
• enhanced stability near discontinuities.

8.10. Computational Complexity

Reversible Part

$$\begin{array}{ccc}& \mathcal{O}(\mid \mathcal{E}\mid ).& \end{array}$$

(147)

Implicit Solve

$$\begin{array}{ccc}& \mathcal{O}(\mid \mathcal{E}\mid \mathrm{l}\mathrm{o}\mathrm{g}N)(\mathrm{multigrid}),& \end{array}$$

(148)

or

$$\begin{array}{ccc}& \mathcal{O}(\mid \mathcal{E}\mid )(\mathrm{optimal}).& \end{array}$$

(149)

Memory

$$\begin{array}{ccc}& \mathcal{O}(\mid \mathcal{E}\mid )(\mathrm{CSR} \mathrm{sparse} \mathrm{storage}).& \end{array}$$

(150)

8.11. Summary

The algorithm presented in this chapter provides:

• full reproducibility (explicit pseudocode),
• method independence (FV / meshless unified),
• structure preservation (conservation + entropy),
• scalability (edge-based parallelization),
• adaptivity (entropy-driven dissipation).

The complete solver realizes the principle:

$$\begin{array}{ccc}& \mathrm{Physical} \mathrm{laws} \Rightarrow \mathrm{operators} \Rightarrow \mathrm{algorithm}.& \end{array}$$

(151)

9. Convergence Analysis and Weak Solution Limit

9.1. Reconstruction Operator

Let ${\mathcal{G}}_h=({\mathcal{V}}_h,{\mathcal{E}}_h)$ denote the discrete graph associated with resolution parameter $h\to 0$.

We define a reconstruction operator

$$\begin{array}{ccc}& {\mathcal{R}}_h:\{U_i\} ⟶ U_h(\mathbf{x}),& \end{array}$$

(152)

such that:

• Finite volume case: piecewise constant reconstruction,

$$\begin{array}{ccc}& U_h(\mathbf{x})=U_i,\mathbf{x}\in {\mathsf{\Omega }}_i,& \end{array}$$

(153)

• Meshless case: kernel-based or Voronoi reconstruction,

$$\begin{array}{ccc}& U_h(\mathbf{x})=\sum _i{U}_i{\varphi }_i(\mathbf{x}).& \end{array}$$

(154)

We assume:

$$\begin{array}{ccc}& \parallel U_h{\parallel }_{L^p(\mathsf{\Omega })}\le C,& \end{array}$$

(155)

uniformly in $h$.

9.2. Discrete Entropy Inequality

From Chapter 6, the fully discrete scheme satisfies:

$$\begin{array}{ccc}& S(U^{n+1})\ge S(U^n).& \end{array}$$

(156)

Summing over time steps yields:

$$\begin{array}{ccc}& S(U^N)-S(U^0)=\sum _{n=0}^{N-1}(S(U^{n+1})-S(U^n))\ge 0.& \end{array}$$

(157)

Using the entropy production structure (Chapter 5):

$$\begin{array}{ccc}& S(U^{n+1})-S(U^n)\ge \mathsf{\Delta }t\sum _{i,j}(v_i-v_j{)}^T{M}_{ij}(v_i-v_j).& \end{array}$$

(158)

Thus,

$$\begin{array}{ccc}& \sum _n\mathsf{\Delta }t\sum _{i,j}(v_i-v_j{)}^T{M}_{ij}(v_i-v_j)\le C.& \end{array}$$

(159)

This provides a global dissipation bound.

9.3. A Priori Estimates

(i) Mass conservation

$$\begin{array}{ccc}& \sum _i{V}_i{\rho }_i(t)=\mathrm{const}.& \end{array}$$

(160)

(ii) Energy boundedness

$$\begin{array}{ccc}& \sum _i{V}_i{E}_i(t)\le C.& \end{array}$$

(161)

(iii) Discrete gradient estimate

From (159),

$$\begin{array}{ccc}& \sum _n\mathsf{\Delta }t\sum _{i,j}{\displaystyle \frac{\mid U_i-U_j{\mid }^2}{d_{ij}^2}}\le C.& \end{array}$$

(162)

This is a discrete analogue of an $H^1$-type estimate.

(iv) Time difference estimate

From the IMEX scheme:

$$\begin{array}{ccc}& {\displaystyle \frac{U_i^{n+1}-U_i^n}{\mathsf{\Delta }t}}={\mathcal{R}}_i(U^n)+{\mathcal{D}}_i(U^{n+1}).& \end{array}$$

(163)

Using boundedness of fluxes,

$$\begin{array}{ccc}& {\parallel {\displaystyle \frac{U^{n+1}-U^n}{\mathsf{\Delta }t}}\parallel }_{H^{-1}}\le C.& \end{array}$$

(164)

9.4. Compactness (Aubin–Lions Type Argument)

From:

• boundedness in $L^2$,
• discrete gradient estimate (162),
• time regularity (164),

we apply a discrete version of the Aubin–Lions lemma [23,24].

Thus, there exists a subsequence such that

$$\begin{array}{ccc}& U_h\to U\mathrm{strongly} \mathrm{in} L^2(\mathsf{\Omega }\times (0,T)).& \end{array}$$

(165)

9.5. Weak Formulation Limit

Let $\phi \in C_c^{\mathrm{\infty }}(\mathsf{\Omega }\times (0,T))$. Multiplying the discrete equation by $\phi$ and summing:

$$\begin{array}{ccc}& \sum _{i,n}[U_i^{n+1}-U_i^n]{\phi }_i^n+\mathsf{\Delta }t\sum _{i,j}{\mathbf{a}}_{ij}\cdot F_{ij}\hspace{0.17em}{\phi }_i^n=0.& \end{array}$$

(166)

Passing to the limit $h\to 0$:

$$\begin{array}{ccc}& {\int }_0^T{\int }_{\mathsf{\Omega }}(U{\partial }_t\phi +F(U)\cdot \nabla \phi )\hspace{0.17em}dxdt+{\int }_{\mathsf{\Omega }}{U}_0\phi (0)\hspace{0.17em}dx=0.& \end{array}$$

(167)

Thus, $U$ satisfies the weak form of the conservation law.

9.6. Entropy Solution

From (156), taking the limit yields:

$$\begin{array}{ccc}& {\partial }_{t}s(U)+\nabla \cdot \mathsf{\Psi }(U)\ge 0,& \end{array}$$

(168)

in the sense of distributions. Thus, the limit $U$ satisfies:

• weak conservation law,
• entropy inequality.

Hence, $U$ is an entropy weak solution [8,13,15].

9.7. Main Theorem

Theorem 9.1 (Convergence to Entropy Weak Solution)

Let $\left\{U_h\right\}$ be a sequence of discrete solutions generated by the scheme defined in Chapters 3–6, under the assumptions:

• mesh regularity and consistency,
• bounded initial data,
• entropy-stable flux construction,
• IMEX time integration.

Then, there exists a subsequence such that

$$\begin{array}{ccc}& U_h\to U\mathrm{in} L^2(\mathsf{\Omega }\times (0,T)),& \end{array}$$

(169)

where $U$ satisfies:

$$\begin{array}{ccc}& {\partial }_{t}U+\nabla \cdot F(U)=\nabla \cdot G(U,\nabla U),& \end{array}$$

(170)

in the weak sense, together with

$$\begin{array}{ccc}& {\partial }_{t}s(U)+\nabla \cdot \mathsf{\Psi }(U)\ge 0.& \end{array}$$

(171)

Therefore, $U$ is an entropy weak solution of the Navier–Stokes system.

9.8. Interpretation

This result establishes that the proposed method:

• is consistent with Navier–Stokes equations,
• preserves thermodynamic structure,
• converges to physically admissible solutions,

without requiring ad hoc stabilization.

The convergence is driven by:

$$\begin{array}{ccc}& \mathrm{structure} \mathrm{preservation} \Rightarrow \mathrm{compactness} \Rightarrow \mathrm{weak} \mathrm{convergence}.& \end{array}$$

(172)

10. Numerical Experiments and Validation

10.1. Purpose and Positioning

The theoretical framework developed in Chapters 2–9 establishes a thermodynamically consistent numerical method based on:

• reversible (entropy-neutral) fluxes
• irreversible (entropy-generating) fluxes

This chapter aims to demonstrate that this structural decomposition is not merely theoretical, but yields tangible numerical advantages.

In particular, we investigate whether the proposed scheme:

• maintains stability comparable to classical methods,
• reduces unnecessary numerical dissipation,
• improves resolution of discontinuities,

while preserving entropy consistency.

To this end, we construct a minimal implementation, referred to as:

$$\begin{array}{ccc}& \mathrm{GEMS}-\mathrm{ES} (\mathrm{Geometric} \mathrm{Entropy}-\mathrm{Modulated} \mathrm{Scheme} – \mathrm{Entropy} \mathrm{Stable}).& \end{array}$$

(173)

This scheme is compared against the classical Rusanov (local Lax–Friedrichs) method [10,32].

10.2. Governing Equations and Test Problem

We consider the one-dimensional compressible Euler equations:

$$\begin{array}{ccc}& {\partial }_t\left(\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right)+{\partial }_x\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ (E+p)u\end{array}\right)=0,& \end{array}$$

(174)

with equation of state:

$$\begin{array}{ccc}& p=(\gamma -1)\left(E−{\displaystyle \frac{1}{2}}\rho u^2\right).& \end{array}$$

(175)

This system admits weak solutions, and entropy conditions are required for physical admissibility [8,13,15].

10.3. Numerical Flux Formulations

10.3.1. Rusanov (LLF) Flux

The Rusanov flux is given by:

$$\begin{array}{ccc}& F_{i+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{1}{2}}(F_L+F_R)-{\displaystyle \frac{1}{2}}{a}_{\mathrm{m}\mathrm{a}\mathrm{x}}(U_R-U_L),& \end{array}$$

(176)

where

$$\begin{array}{ccc}& a_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{m}\mathrm{a}\mathrm{x}(\mid u\mid +c).& \end{array}$$

(177)

This flux is robust but introduces uniform dissipation across all wave types.

10.3.2. GEMS-ES Flux

The proposed GEMS-ES flux is constructed as:

$$\begin{array}{ccc}& F_{i+{\displaystyle \frac{1}{2}}}=F^{\mathrm{s}\mathrm{y}\mathrm{m}}-{\displaystyle \frac{1}{2}}{a}_{\mathrm{m}\mathrm{a}\mathrm{x}}(U_R-U_L),& \end{array}$$

(178)

where

$$\begin{array}{ccc}& F^{\mathrm{s}\mathrm{y}\mathrm{m}}=F(U_{\mathrm{a}\mathrm{v}\mathrm{g}}).& \end{array}$$

(179)

Here, $U_{\mathrm{a}\mathrm{v}\mathrm{g}}$ is a primitive-variable-based symmetric state.

The essential feature is:

$$\begin{array}{ccc}& \mathrm{Flux}=\underset{\mathrm{transport}}{\underbrace{\mathrm{reversible} (\mathrm{symmetric})}}+\underset{\mathrm{stabilization}}{\underbrace{\mathrm{dissipative} \mathrm{correction}}}.& \end{array}$$

(180)

This explicit separation reduces unnecessary numerical viscosity.

10.4. Computational Setup

We consider the classical Sod shock tube problem:

$$\begin{array}{ccc}& (\rho ,u,p)=\left\{\begin{array}{cc}(\mathrm{1,0},1),& x<x_0,\\ (\mathrm{0.125,0},0.1),& x>x_0,\end{array}\right.& \end{array}$$

(181)

with parameters:

$$\begin{array}{ccc}& \gamma =1.4,N=400,\mathrm{CFL}=0.5,t=0.2.& \end{array}$$

(182)

Time integration is performed using SSPRK(3,3) [14].

The exact solution is obtained by solving the Riemann problem using a Newton iteration for the star-region pressure [10].

10.5. Numerical Results

Figure 1, Figure 2 and Figure 3 show the density, velocity, and pressure profiles for the classical Sod shock tube problem, computed using the proposed scheme and the conventional Rusanov (local Lax–Friedrichs) method at $t=0.2$ with $N=400$ and $\mathrm{C}\mathrm{F}\mathrm{L}=0.5$. Each figure overlays three solutions: the exact Riemann solution, the Rusanov scheme, and the proposed GEMS-ES method, which represents a minimal implementation combining a reversible flux with entropy-stable dissipative correction. Both methods reproduce:

• correct shock position,
• correct shock strength.

Thus, stability is preserved in the proposed scheme.

The L1 errors are comparable:

$$\begin{array}{ccc}& \parallel \rho -{\rho }_{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}{\parallel }_{L^1}\approx 1.2\times {10}^{-2},& \end{array}$$

(183)

$$\begin{array}{ccc}& \parallel p-p_{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}{\parallel }_{L^1}\approx 9.3\times {10}^{-3}.& \end{array}$$

(184)

Thus, the improvement is not due to parameter tuning, but structural.

10.6. Structural Interpretation

The difference between the two methods can be summarized as:

Rusanov Method

• no separation of transport and dissipation,
• uniform artificial viscosity,
• over-diffusion of contact waves.

GEMS-ES Method

• reversible and dissipative parts explicitly separated,
• minimal required dissipation,
• selective damping of non-physical modes.

This leads to:

$$\begin{array}{ccc}& \mathrm{Structure} \mathrm{preservation} \Rightarrow \mathrm{improved} \mathrm{resolution}.& \end{array}$$

(185)

10.7. Minimal Implementation and Reproducibility

The entire computation was implemented using:

• a single VBA module,
• no external libraries,
• direct computation of both numerical and exact solutions.

Thus, the results:

• do not rely on high-order reconstruction,
• do not depend on limiters,
• isolate purely the flux structure effect.

10.8. Discussion

The numerical experiment confirms that:

• The reversible–irreversible decomposition is numerically realizable.
• Stability is maintained at the level of classical LLF schemes.
• Contact discontinuities are significantly sharper.
• Entropy consistency is preserved.

This demonstrates that the theoretical design principle:

$$\begin{array}{ccc}& \mathrm{Thermodynamic} \mathrm{structure} \Rightarrow \mathrm{numerical} \mathrm{performance}& \end{array}$$

(186)

is not merely formal, but operational.

10.9. Conclusion of Numerical Validation

The present validation shows that GEMS-ES:

• achieves comparable robustness to Rusanov schemes,
• significantly reduces numerical diffusion,
• preserves entropy consistency,
• improves contact discontinuity resolution.

Therefore,

$$\begin{array}{ccc}& \mathrm{GEMS}-\mathrm{ES} \mathrm{provides} \mathrm{a} \mathrm{structurally} \mathrm{superior} \mathrm{alternative} \mathrm{to} \mathrm{LLF}-\mathrm{type} \mathrm{schemes}.& \end{array}$$

(187)

11. Discussion

11.1. Differences from Conventional Numerical Methods

Classical numerical methods for fluid dynamics are constructed by discretizing the governing partial differential equations, typically the Navier–Stokes system [5,6,7,10].

In this paradigm, stability is achieved through the addition of artificial dissipation:

$$\begin{array}{ccc}& \mathrm{Numerical} \mathrm{flux}=\mathrm{physical} \mathrm{flux}+\mathrm{stabilization} \mathrm{term}.& \end{array}$$

(188)

This stabilization is often designed heuristically, leading to:

• uniform dissipation across all wave types,
• loss of sharp interface resolution,
• dependence on problem-specific tuning parameters [9,14,32].

In contrast, the present framework adopts a fundamentally different construction:

$$\begin{array}{ccc}& \mathrm{Numerical} \mathrm{flux}=\underset{\mathrm{entropy}-\mathrm{neutral}}{\underbrace{\mathrm{reversible} \mathrm{component}}}+\underset{\mathrm{entropy}-\mathrm{generating}}{\underbrace{\mathrm{irreversible} \mathrm{component}}}.& \end{array}$$

(189)

The key distinctions are:

• dissipation is not added arbitrarily, but derived from entropy principles,
• conservation and entropy inequality are embedded structurally,
• transport and dissipation are explicitly separated.

Thus, the difference is not incremental but structural:

$$\begin{array}{ccc}& \mathrm{Discretization}-\mathrm{based} \mathrm{methods} \ne \mathrm{Structure}-\mathrm{based} \mathrm{methods}.& \end{array}$$

(190)

11.2. Interpretation of the Master Equation

The master equation introduced in Chapter 2 provides a more fundamental description of fluid systems. Rather than expressing dynamics in terms of differential operators, the evolution is written as:

$$\begin{array}{ccc}& {\displaystyle \frac{d{U}_i}{dt}}=\sum _j(A_{ij}(U_i,U_j)+D_{ij}(U_i,U_j)).& \end{array}$$

(191)

This formulation has several important implications:

(i) Conservation as antisymmetry

$$\begin{array}{ccc}& A_{ij}=-A_{ji}.& \end{array}$$

(192)

Thus, conservation is guaranteed at the level of interactions.

(ii) Dissipation as gradient flow

$$\begin{array}{ccc}& D_{ij}\propto M_{ij}(v_j-v_i).& \end{array}$$

(193)

This directly encodes the second law.

(iii) Continuum as a limit

$$\begin{array}{ccc}& \sum _j{\mathbf{a}}_{ij}\cdot (\cdot ) ⟶ \nabla \cdot (\cdot ).& \end{array}$$

(194)

Hence, the continuum description is not fundamental, but emergent.

This suggests a reinterpretation:

$$\begin{array}{ccc}& \mathrm{Fluid} \mathrm{dynamics}=\mathrm{thermodynamically} \mathrm{constrained} \mathrm{interaction} \mathrm{system}.& \end{array}$$

(195)

11.3. Reinterpretation of the Navier–Stokes Equations

Within the present framework, the Navier–Stokes equations [3,4,24] are no longer the starting point, but arise as a limiting case:

$$\begin{array}{ccc}& \mathrm{Master} \mathrm{Equation} \underset{\mathrm{continuum} \mathrm{limit}}{\to } \mathrm{Navier}–\mathrm{Stokes}.& \end{array}$$

(196)

This inversion has important consequences:

(i) Structural origin of viscosity

Viscosity is not introduced as a constitutive assumption, but emerges from:

$$\begin{array}{ccc}& D_{ij}\sim M_{ij}(v_j-v_i).& \end{array}$$

(197)

(ii) Entropy inequality as primary constraint

Rather than being imposed on weak solutions, the entropy condition is built into the dynamics:

$$\begin{array}{ccc}& {\displaystyle \frac{dS}{dt}}\ge 0(\mathrm{discrete} \mathrm{and} \mathrm{continuous}).& \end{array}$$

(198)

(iii) Weak solutions as natural limits

The entropy weak solution framework [8,13,15] arises naturally from compactness of the interaction system.

Thus, the Navier–Stokes equations can be viewed as an effective macroscopic closure of a more fundamental interaction-based system.

11.4. Implications for Numerical Analysis

The present formulation suggests a new paradigm for numerical methods:

$$\begin{array}{ccc}& \mathrm{Physical} \mathrm{structure} \Rightarrow \mathrm{numerical} \mathrm{structure} \Rightarrow \mathrm{stability} \mathrm{and} \mathrm{convergence}.& \end{array}$$

(199)

In particular:

• stability is guaranteed by entropy monotonicity,
• convergence follows from compactness induced by dissipation,
• consistency arises from geometric reconstruction.

This contrasts with traditional approaches, where:

$$\begin{array}{ccc}& \mathrm{numerical} \mathrm{stabilization} \Rightarrow \mathrm{approximate} \mathrm{physical} \mathrm{consistency}.& \end{array}$$

(200)

11.5. Extensions and Future Directions

The master equation framework naturally admits several extensions.

(i) Nonlocal interactions

The interaction kernel can be generalized:

$$\begin{array}{ccc}& {\displaystyle \frac{d{U}_i}{dt}}=\sum _j{K}_{ij}(U_i,U_j),& \end{array}$$

(201)

where $K_{ij}$ may include long-range effects.

This enables modeling of:

• fractional diffusion,
• kinetic transport,
• nonlocal turbulence closures.

(ii) Turbulence modeling

The reversible–irreversible decomposition suggests a natural turbulence closure:

$$\begin{array}{ccc}& \mathcal{D}={\mathcal{D}}_{\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{r}}+{\mathcal{D}}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}.& \end{array}$$

(202)

Here, turbulence can be interpreted as an additional entropy-producing mechanism, consistent with cascade theory [11,12].

(iii) Multiphase and multicomponent flows

The state vector can be extended:

$$\begin{array}{ccc}& U_i=({\rho }_k,\mathbf{m},E),k=1,\dots ,N_s,& \end{array}$$

(203)

allowing:

• phase interactions,
• chemical reactions,
• interfacial dynamics.

(iv) Coupling with optimal transport

The gradient flow structure connects naturally with Wasserstein geometry [27]:

$$\begin{array}{ccc}& {\partial }_{t}U=\nabla \cdot \left(U\nabla {\displaystyle \frac{\delta \mathcal{F}}{\delta U}}\right).& \end{array}$$

(204)

This opens a pathway toward variational formulations of fluid dynamics.

11.6. Conceptual Summary

The framework developed in this work can be summarized as:

$$\mathrm{Interaction} \mathrm{structure} \Rightarrow \mathrm{thermodynamic} \mathrm{consistency} \Rightarrow \mathrm{numerical} \mathrm{robustness}.$$

(205)

More fundamentally, it suggests a shift in perspective:

$$\mathrm{PDE}-\mathrm{based} \mathrm{view} ⟶ \mathrm{interaction}-\mathrm{based} \mathrm{view}.$$

(206)

This shift has implications not only for numerical methods, but for the theoretical understanding of fluid dynamics itself.

12. Conclusion

12.1. Summary of Contributions

In this study, we have developed a new computational framework for fluid dynamics based on a thermodynamically consistent master equation formulation.

The central idea is to abandon the traditional paradigm in which numerical methods are constructed by discretizing partial differential equations, and instead to begin from an interaction-based representation in which:

$$\begin{array}{ccc}& \mathrm{time} \mathrm{evolution}=\mathrm{conservative} \mathrm{interactions}+\mathrm{dissipative} \mathrm{interactions}.& \end{array}$$

(207)

Within this framework, we have established:

• a master equation formulation that encodes conservation and entropy production at the structural level,
• a graph-based spatial discretization that unifies finite volume and meshless methods,
• a reversible–irreversible flux decomposition ensuring entropy consistency,
• an IMEX time integration scheme preserving conservation and entropy monotonicity,
• a unified solver architecture with high scalability and adaptability,
• a rigorous convergence result toward entropy weak solutions of the Navier–Stokes equations.

12.2. Key Findings

The principal findings of this work can be summarized as follows.

(i) Structure Determines Stability

We have shown that numerical stability can be achieved not by adding artificial dissipation, but by constructing schemes that satisfy:

$$\begin{array}{ccc}& \mathrm{entropy} \mathrm{production} \ge 0(\mathrm{by} \mathrm{design}).& \end{array}$$

(208)

This eliminates the need for ad hoc stabilization mechanisms.

(ii) Separation of Transport and Dissipation

The explicit decomposition

$$\begin{array}{ccc}& F=F^{\mathrm{r}\mathrm{e}\mathrm{v}}+F^{\mathrm{i}\mathrm{r}\mathrm{r}}& \end{array}$$

(209)

provides a clear distinction between:

• reversible transport processes,
• irreversible dissipative processes.

This separation leads directly to improved numerical behavior, as demonstrated in the numerical experiments.

(iii) Navier–Stokes as an Emergent Limit

The Navier–Stokes equations are recovered as a continuum limit:

$$\begin{array}{ccc}& \mathrm{Master} \mathrm{Equation} ⟶ \mathrm{Navier}–\mathrm{Stokes}.& \end{array}$$

(210)

This establishes that the proposed framework is not an alternative to classical fluid mechanics, but a more fundamental representation from which it emerges.

(iv) Unified Numerical Framework

The solver developed in this work demonstrates that:

$$\begin{array}{ccc}& \mathrm{finite} \mathrm{volume}\equiv \mathrm{meshless}(\mathrm{at} \mathrm{the} \mathrm{operator} \mathrm{level}).& \end{array}$$

(211)

Thus, seemingly different numerical methods are unified within a single structural framework.

12.3. Implications

The results of this work suggest a broader conceptual shift:

$$\begin{array}{ccc}& \mathrm{fluid} \mathrm{dynamics}=\mathrm{thermodynamically} \mathrm{constrained} \mathrm{interaction} \mathrm{system}.& \end{array}$$

(212)

This viewpoint has several important implications:

• numerical methods can be derived systematically from physical principles,
• stability and convergence become consequences of structure,
• the role of continuum equations is clarified as an emergent description.

12.4. Limitations

The present study has focused on:

• first-order spatial discretization,
• relatively simple test problems,
• minimal implementations to isolate structural effects.

While this is sufficient to demonstrate the validity of the framework, further work is required to:

• extend to high-order accuracy,
• address complex geometries and boundary conditions,
• perform large-scale turbulent simulations.

12.5. Future Directions

Several promising directions follow naturally from this work.

(i) High-order extensions

Development of higher-order schemes that preserve entropy structure while improving accuracy.

(ii) Turbulence modeling

Incorporation of turbulence as an additional entropy-generating mechanism:

$$\begin{array}{ccc}& \mathcal{D}={\mathcal{D}}_{\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{r}}+{\mathcal{D}}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{t}}.& \end{array}$$

(213)

(iii) Multiphysics coupling

Extension to:

• multiphase flows,
• reactive systems,
• plasma and kinetic models.

(iv) Variational and optimal transport formulations

Exploration of connections to gradient flow structures:

$$\begin{array}{ccc}& {\partial }_{t}U=\nabla \cdot \left(U\nabla {\displaystyle \frac{\delta \mathcal{F}}{\delta U}}\right).& \end{array}$$

(214)

12.6. Final Remark

The framework proposed in this work demonstrates that:

$$\mathrm{Physical} \mathrm{laws} \Rightarrow \mathrm{structure} \Rightarrow \mathrm{numerical} \mathrm{method}.$$

(215)

This perspective provides a unified and physically grounded foundation for computational fluid dynamics, and opens a pathway toward more robust, accurate, and theoretically consistent numerical methods.