Conservation Properties in Discretization: Ensuring Mass and Energy Balance

TL;DR

• Local conservation is guaranteed by FVM flux consistency; FEM and DG must design fluxes carefully.
• First-order upwind adds artificial viscosity — at CFL = 1, it vanishes, a surprising result most practitioners miss.
• MOOSE’s RDG(P0P1) with flux limiting reduces numerical diffusion by over 80% compared to full upwinding.
• Skew-symmetric formulations conserve kinetic energy by preventing aliasing errors from creating energy spuriously.
• New 2026 relaxation methods can enforce mass and energy conservation to machine precision (10⁻¹⁴).

Why Discretization Changes What Your Simulation Conserves

When you discretize a partial differential equation, you’re not just choosing an accuracy order — you’re choosing what physical quantities your solver actually conserves. Mass, energy, momentum: these properties that are analytically guaranteed by the continuous equations can quietly disappear (or mysteriously appear) at the discrete level depending on your numerical flux design.

This is one of the most consequential trade-offs in scientific simulation, yet it’s often overlooked during method selection. You might spend hours tuning parameters only to discover that the discretization scheme itself is introducing artificial mass loss or spurious energy generation.

The core principle is simple: flux leaving one cell must equal flux entering its neighbor. When this is true by construction, you have local conservation. When it requires careful design, you need to verify it. Understanding this distinction — and the tools to verify it — is the difference between trusting your simulation and auditing every result.

The Flux Consistency Principle: Local vs Global Conservation

The conservation property of a numerical scheme is determined by whether numerical fluxes are consistent across cell interfaces. This is known in the literature as the Flux Consistency Principle, and it explains why some methods naturally conserve mass while others don’t.

Finite Volume Method: Conservation by Construction

The finite volume method (FVM) guarantees local conservation because it balances fluxes across cell faces by definition. When mass flows out of one control volume, the exact same amount flows into its neighbor. No numerical flux can be created or destroyed — the algebraic structure of the discretization enforces global mass conservation automatically.

FiPy implements this through its conservation equation structure, where every convection term represents a flux balance across cell faces. The FiPy conservation equation documentation describes how the general conservation equation enforces exactly this flux consistency.

# FiPy conservation equation — fluxes are balanced by construction
equation = TransientTerm(var=phi) == UpwindConvectionTerm(coeff=velocity)

Finite Element Method: Conservation Requires Design

Standard continuous finite element formulations do not automatically conserve mass locally. The weak form integrates over elements, and fluxes across internal boundaries can cancel globally without being conserved locally. This is why standard FEM formulations can lose mass conservation in advection-dominated problems.

However, FEM can achieve local conservation through specially designed numerical fluxes. The discontinuous Galerkin (DG) method, for example, achieves strong local conservation through its interface fluxes — see the DG method tutorial for a detailed derivation. DG combines high-order polynomial approximation with flux-based conservation that mirrors FVM’s flux consistency.

What This Means for Method Selection

The When to Use FEM, FVM, or FDM guide on matforge.org explains these trade-offs in detail, including how conservation requirements should drive method selection before geometry or accuracy are considered.

Truncation Error and Artificial Viscosity

Every numerical scheme introduces truncation error — a difference between the exact differential operator and its discrete approximation. For upwind schemes, this truncation error has a physical interpretation: it behaves like artificial diffusion (also called numerical viscosity).

The Modified Equation Derivation

The CFD University modified equation analysis provides the complete derivation. Here’s the essential result.

For first-order upwind advection with constant velocity $a$:

$$\frac{\partial \phi}{\partial t} + a \frac{\partial \phi}{\partial x} = \underbrace{\frac{a \Delta x}{2}}_{\text{Artificial viscosity coefficient}} \frac{\partial^2 \phi}{\partial x^2} + O(\Delta x^2)$$

The leading truncation error acts like a diffusion term with coefficient $\nu = a \Delta x / 2$. This explains both why upwind is stable (it adds diffusion) and why it smears fronts (that diffusion spreads profiles).

The CFL Number Changes Everything

Here’s the non-obvious result that most practitioners miss: when a time discretization is included, the artificial viscosity coefficient becomes CFL-dependent. The complete modified equation shows:

$$\nu = \frac{a \Delta x}{2} (1 – \text{CFL})$$

At CFL = 1, the artificial diffusion vanishes. This explains when upwinding is “free” — no numerical diffusion at all. Many simulation practitioners unknowingly operate at CFL = 1 and wonder why their first-order upwind scheme produces unexpectedly clean results. The CFL-dependent formula explains why.

This result has significant practical implications:

• At CFL close to 0 (very small time steps), upwinding introduces maximum artificial diffusion
• At CFL = 1 (time step equals one cell transit time), first-order upwind adds no artificial diffusion
• At CFL > 1, the formula suggests potential instability — the modified equation changes sign

Note that temporal discretization also contributes to numerical diffusion. The time integrator choice (Backward Euler, Crank–Nicolson, Runge-Kutta) affects damping characteristics independently of the spatial scheme. The Time Integration Methods for PDE Solvers: Explicit vs Implicit Schemes guide discusses how temporal error interacts with spatial numerical diffusion and how to separate them in verification.

Grid Misalignment Creates “False Diffusion”

Even with zero physical diffusion coefficient and perfectly uniform flow, unstructured meshes introduce numerical diffusion proportional to the angle between the flow direction and cell face normals. This is known as “false diffusion” or “crosswind diffusion”.

On structured Cartesian grids, flow aligned with a coordinate axis is represented cleanly. On unstructured meshes (triangles, tetrahedra), even simple uniform flow introduces artificial diffusion when the flow crosses the mesh obliquely. The diffusion magnitude scales with $\sin^2(\theta)$, where $\theta$ is the misalignment angle.

This explains why simulations on unstructured meshes often require finer resolution than structured meshes for the same advection problem — the artificial diffusion isn’t a property of the scheme, but of the mesh orientation.

Upwind Schemes and the Diffusion-Stability Trade-Off

Upwind schemes add numerical diffusion to stabilize advection calculations. But how much diffusion is introduced? And what alternatives exist?

What Upwind Actually Does

First-order upwind replaces the symmetric central difference:

$$\frac{\partial \phi}{\partial x} \approx \frac{\phi_{i+1} – \phi_{i-1}}{2\Delta x}$$

with a direction-aware approximation:

$$\frac{\partial \phi}{\partial x} \approx \frac{\phi_i – \phi_{i-1}}{\Delta x} \quad \text{(for } a > 0\text{)}$$

As the CFD University analysis shows, this upwind approximation can be written exactly as:

$$\text{Upwind} = \text{Central} – \frac{\Delta x}{2} \frac{\partial^2}{\partial x^2}$$

The second-derivative term is the artificial diffusion. It suppresses short waves and sharp gradients — which is why upwind is stable but also why it smears fronts.

MOOSE Numerical Diffusion Comparison

The MOOSE numerical diffusion documentation provides the most comprehensive numerical diffusion comparison available in any open-source solver documentation. Their tracer advection test compares seven discretization strategies:

The RDG(P0P1) scheme with the superbee flux limiter reduces numerical diffusion by over 80% compared to full upwinding while maintaining monotonicity. This is concrete evidence for method selection: linear reconstruction with flux limiting is vastly superior to first-order upwind for transport problems where resolution matters.

When to Use Full Upwind vs Limited Schemes

The MOOSE upwinding documentation provides practical guidance:

• Use full upwind when: robustness is the priority, solution boundedness matters more than resolution, or you’re solving strong nonlinear problems where oscillations would break the solver
• Use limited higher-order schemes when: transport rates are being measured quantitatively, sharp fronts must be resolved, or the problem has weak nonlinearities that don’t require full upwinding

The practical tradeoff: Full upwind is robust but over-diffusive. Limited schemes preserve resolution but add nonlinear complexity. For conservation-critical simulations (reactive transport, phase-field models), the extra complexity is usually worth the resolution gain.

Skew-Symmetric Formulations for Energy Preservation

When simulating compressible flow or reacting systems, energy conservation becomes critical. Standard conservative formulations of nonlinear convection can spuriously generate or destroy kinetic energy through aliasing errors. The solution, identified in a widely cited review, is to reformulate the convective operator.

Coppola 2019: The Skew-Symmetric Split

The seminal 2019 ASME review by Coppola (95 citations) on discrete energy conservation demonstrates that non-linear convective terms conserve kinetic energy only when formulated in split form. The key insight is that the convective operator must be split into:

1. A conservative form (divergence form): $\nabla \cdot (\mathbf{u} \phi)$
2. A skew-symmetric form: $\mathbf{u} \cdot \nabla \phi – \frac{1}{2} (\nabla \cdot \mathbf{u}) \phi$

When combined, the discrete convective operator becomes skew-adjoint, which means it preserves quadratic invariants like kinetic energy without artificial dissipation. This prevents aliasing errors from spuriously creating (or destroying) energy.

For incompressible flow ($\nabla \cdot \mathbf{u} = 0$), the skew-symmetric form simplifies to:

$$\mathbf{u} \cdot \nabla \phi$$

For compressible flow, the divergence term matters, and the skew-symmetric split becomes essential for energy conservation.

Why This Matters for DG Methods

The discontinuous Galerkin energy-stable method (DGSEM), as developed by Gassner and others, uses this skew-symmetric splitting to ensure discrete energy preservation. See the DG method tutorial for DG flux design details. DG methods with entropy-stable fluxes and skew-symmetric convective splitting are the state-of-the-art for energy-preserving simulations of compressible flow.

When Skew-Symmetric Formulations Are Essential

Verifying Conservation in Practice

Theory tells us which schemes conserve mass and energy. In practice, you need to verify it. Here’s a step-by-step workflow.

Step 1: Define the Conservation Quantity

Identify which physical quantities must be conserved. Common choices:

• Mass (integral of concentration or density over the domain)
• Energy (thermal energy, kinetic energy, or total energy)
• Species mass (individual chemical species)

Step 2: Compute the Discrete Integral

At each time step, compute the domain integral:

# Mass conservation check in a scalar transport problem
import numpy as np
from fipy import Grid2D, CellVariable

mesh = Grid2D(nx=100, ny=100, dx=0.01, dy=0.01)
phi = CellVariable(mesh=mesh, value=0.0)

# Initial mass
initial_mass = np.sum(phi.value * mesh.cellVolumes)
print(f"Initial mass: {initial_mass:.6f}")

# After simulation
final_mass = np.sum(phi.value * mesh.cellVolumes)
print(f"Final mass: {final_mass:.6f}")
print(f"Mass loss: {(initial_mass - final_mass) / initial_mass * 100:.4f}%")

Step 3: Check Boundary Fluxes

Conservation violations usually originate at boundaries. Verify that:

• Inflow fluxes match prescribed boundary values
• Outflow boundaries don’t over-constrain the problem
• No-flux boundaries truly have zero flux

Step 4: Compare Against Reference

The Verification vs Validation guide explains that verification requires comparing against analytical or manufactured solutions. For conservation verification:

1. Run a transport problem with known initial mass
2. Verify the numerical mass remains constant (within tolerance) over time
3. Test on multiple mesh resolutions — conservation should improve or remain constant as mesh refines, not degrade

Step 5: Document the Tolerance

Set a conservation tolerance appropriate for your problem:

• Scientific publication: Mass error < 0.1%
• Engineering simulation: Mass error < 1%
• Research prototyping: Mass error < 5% (but document it)

Advanced Conservation Enforcement: Relaxation-Based Methods

When standard discretization doesn’t achieve sufficient conservation, recent methods offer post-processing corrections.

The arXiv 2604.05940 Result

A recent preprint by Rajvanshi and Ketcheson (arXiv 2604.05940, April 2026) introduces relaxation-based conservation enforcement for mass and energy balance. The method works with any time-stepping scheme by adding a post-time-step correction.

Key results:

• Achieves machine-precision conservation up to 10⁻¹⁴ error
• Two variants: multiple-relaxation and projection-relaxation
• Projection-relaxation is more robust and only requires solving a scalar root-finding problem (2-3× computational cost)
• Validated on 2D and 3D cosmological test cases

The method applies to implicit-explicit (ImEx) Runge-Kutta schemes and can enforce both mass conservation and energy balance simultaneously. This is particularly valuable for simulations where conservation violations accumulate over long time spans and degrade physical fidelity.

When to use relaxation-based methods:

• Long-duration simulations where small conservation errors accumulate
• Multi-physics coupling where conservation violations cause instability
• High-order schemes that sacrifice strict conservation for resolution
• Production simulations where mass/energy balance is non-negotiable

Practical Decision Framework

Choosing a discretization for conservation-critical simulations requires considering multiple factors simultaneously.

Method Selection Matrix

Decision Rules

1. If conservation is non-negotiable (mass/energy balance must be exact): Start with FVM or DG. The flux consistency principle guarantees it by design.
2. If sharp fronts must be resolved without spurious oscillations: Use limited higher-order schemes (MOOSE RDG P0P1, Kuzmin-Turek with limiter). Avoid full upwind — it reduces diffusion by less than what limited schemes achieve.
3. If compressible flow energy conservation matters: Use skew-symmetric split formulations (Coppola 2019, Gassner DGSEM). Standard conservative forms can generate energy through aliasing.
4. If you need strict mass conservation on long simulations: Consider relaxation-based enforcement (arXiv 2604.05940) as a post-processing step at 2-3× computational cost.
5. If mesh resolution is limited: Be aware of false diffusion on unstructured meshes. Structured grids avoid this. On unstructured meshes, expect higher artificial diffusion proportional to the flow-mesh misalignment angle.

Summary: What to Avoid

Summary and Next Steps

Conservation properties in discretization are determined by flux consistency across cell interfaces. FVM achieves it by construction. DG achieves it through numerical fluxes. Standard FEM requires careful design.

The key insight is that numerical diffusion is not an unavoidable property of advection schemes — it depends on CFL number, mesh orientation, and scheme choice. First-order upwind adds artificial viscosity that vanishes at CFL = 1. Limited higher-order schemes reduce numerical diffusion by over 80% compared to full upwinding. Skew-symmetric formulations prevent spurious energy generation in compressible flow.

For conservation-critical simulations, always verify mass and energy balance quantitatively. Don’t rely on method choice alone — the flux consistency principle is a theoretical guarantee, but your specific implementation must be checked.

Recommended next steps:

1. Run a mass conservation test on your current simulation (compute domain integrals at each time step)
2. If conservation error exceeds your tolerance, consider upgrading to limited higher-order schemes
3. For compressible flow, verify your convective operator uses skew-symmetric splitting
4. Consider relaxation-based enforcement for long-duration simulations where small errors accumulate

The Advection-Dominated Problems guide on matforge.org provides complementary guidance on stabilization strategies (SUPG, flux limiters, Peclet analysis) that work alongside the conservation properties covered here.

Further Reading

• Coppola, S. (2019). On the conservation of quadratic invariants in numerical schemes for the compressible Euler and Navier–Stokes equations. ASME Applied Mechanics Reviews, 71(1). 95 citations.
• Rajvanshi, A., & Ketcheson, D. I. (2026). Efficient high-order mass-conserving and energy-balancing schemes. arXiv:2604.05940.
• MOOSE numerical diffusion documentation: https://mooseframework.inl.gov/modules/porous_flow/numerical_diffusion.html
• CFD University modified equation analysis: https://cfd.university/blog/what-is-numerical-dissipation-in-cfd-and-why-do-we-need-it/

Related Guides

• Advection-Dominated Problems: Upwind Schemes and SUPG Stabilization — Cross-reference for stabilization strategies (SUPG, flux limiters, Peclet analysis)
• Discontinuous Galerkin Method Tutorial — DG flux design and local conservation details
• Verification vs Validation in Scientific Simulations — MMS and verification methodology
• When to Use FEM, FVM, or FDM — Method selection including conservation requirements
• Time Integration Methods for PDE Solvers — Temporal error discussion and CFL restrictions