← Back to Mesh-based solvers Module · scimba_jax

Finite volumes

Finite volumes are the simplest of the three mesh-based methods. They use one value per cell, namely the cell average, and no basis functions. When a PDE is written as a conservation law, the evolution of each cell average only depends on the fluxes through the faces of the cell. A finite volume scheme discretizes these fluxes.

1. A Poisson problem as a conservation law

\(-\nabla\cdot(\nabla u) = f\)

This equation is already in conservation form: the flux is grad u, and its divergence balances the source. AbstractConservativePDE requires the same decomposition: a diffusive flux A, an optional advective flux F, a reaction term R and a source f. On a Cartesian mesh, TPFADiffusionFlux (two-point flux approximation) is the simplest consistent approximation of the diffusive flux. The flux through a face is proportional to the difference between the values of the two neighbouring cells, divided by the distance between them. 

import jax
import jax.numpy as jnp
from scimba_jax.linear_approximation.finite_volume import (
    EllipticFVScheme,
    TPFADiffusionFlux,
)
from scimba_jax.linear_approximation.meshes.cartesian_mesh import cartesian_mesh
from scimba_jax.linear_approximation.variables.variables_fv import VariablesFV
from scimba_jax.nonlinear_approximation.model_class.funcparam_matrix import (
    ParamMatrixFunction,
)
from scimba_jax.nonlinear_approximation.model_class.funcparam_scalar import (
    ParamScalarFunction,
)
from scimba_jax.physical_models.abstract_conservative_pde import (
    AbstractConservativePDE,
)
from scimba_jax.physical_models.abstract_physical_conservative_model import (
    AbstractPhysicalConservativeModel,
)
class Poisson2D(AbstractConservativePDE):
    """``-div(grad u) = f``, written as a conservation law: div(flux) = source."""
    def __init__(self):
        super().__init__(dim=2)
    def construct_A(self, u):  # the diffusion tensor multiplying grad(u)
        return ParamMatrixFunction(u.dims, lambda _sp, _x: jnp.eye(2), f_type=u.f_type)
    def construct_f(self):  # the source, one value per cell
        return ParamScalarFunction(
            {"x": 2, "mu": 0},
            lambda _sp, x: 2.0 * (x[0] * (1.0 - x[0]) + x[1] * (1.0 - x[1])),
            f_type="x",
        )
    def construct_F(self, u):
        return None  # no first-order (advective) flux here
    def construct_R(self, u):
        return None  # no reaction term here
mesh = cartesian_mesh([32, 32], quad_order=3)
model = AbstractPhysicalConservativeModel.from_pde(
    Poisson2D(), dirichlet=lambda x: jnp.zeros(1)
)
# TPFA: flux through a face = (u_right - u_left) / distance, times A.
scheme = EllipticFVScheme(
    model,
    VariablesFV(mesh, projection="average"),
    TPFADiffusionFlux(),
    assemble_first_order=False,
    assemble_reaction=False,
)
scheme = EllipticFVScheme.solve(scheme)
# RMS error against x(1-x)y(1-y): 3.7e-05 at 32x32 cells.

2. Flux classes

Every finite volume flux is an AbstractFVFlux, which has up to four methods, one per combination of interior/boundary face and first-order (advective)/second-order (diffusive) term. A scheme only calls the methods required by its problem. For instance, TPFADiffusionFlux above only overrides the two diffusive methods, since the Poisson problem has no advection. Conversely, for a hyperbolic conservation law, RusanovFlux only overrides the two advective methods. It only requires a bound on the local wave speed to be stable and monotone. 

class AbstractFVFlux(ScimbaPytree):
    """A finite volume flux implements up to four methods on a face."""
    def first_order(self, model, u_left, u_right, context):
        """Advective flux density leaving the left cell (a Riemann solver)."""
        return jnp.zeros_like(u_left)
    def second_order(self, model, u_left, u_right, context):
        """Diffusive flux density leaving the left cell (e.g. two-point flux)."""
        return jnp.zeros_like(u_left)
    def boundary_first_order(self, model, u_inside, boundary_condition, context):
        """Same as first_order on a boundary face, with context.point instead of u_right."""
        return jnp.zeros_like(u_inside)
    def boundary_second_order(self, model, u_inside, boundary_condition, context):
        """Diffusive counterpart of boundary_first_order."""
        return jnp.zeros_like(u_inside)
# TPFADiffusionFlux only overrides the diffusive methods. RusanovFlux only overrides
# first_order (local Lax-Friedrichs flux, with a bound on the wave speed).
from scimba_jax.linear_approximation.finite_volume.hyperbolic import RusanovFlux
flux = RusanovFlux(wave_speed=lambda _model, u, _context: jnp.abs(u[0]))
# {F(uL)+F(uR)}/2 - alpha/2 (uR - uL), alpha = max wave speed on the face.

3. Time-dependent problems

\(\partial_t u + \partial_x\!\left(\dfrac{u^2}{2}\right) = 0\)

TimeDiscreteFVscheme wraps a FiniteVolumeScheme in the same way that TimeDiscreteFEscheme wraps an elliptic FEM problem: the spatial scheme gives the right-hand side, and a Butcher tableau defines the time stepping. For hyperbolic problems with shocks, such as the Burgers equation, the usual choice is an explicit tableau under a CFL condition. Indeed, these problems are not stiff, so an implicit scheme would not allow a much larger time step. 

from scimba_jax.linear_approximation.finite_volume import (
    FiniteVolumeScheme,
    TimeDiscreteFVscheme,
)
from scimba_jax.linear_approximation.finite_volume.hyperbolic import RusanovFlux
from scimba_jax.linear_approximation.meshes.cartesian_mesh import cartesian_mesh
from scimba_jax.linear_approximation.variables.variables_fv import VariablesFV
from scimba_jax.nonlinear_approximation.model_class.funcparam_scalar import (
    ParamScalarFunction,
)
from scimba_jax.physical_models.abstract_conservative_pde import AbstractConservativePDE
from scimba_jax.physical_models.abstract_physical_conservative_model import (
    AbstractPhysicalConservativeModel,
)
from scimba_jax.physical_models.weak_boundary_conditions import Dirichlet
from scimba_jax.time_discrete.butcher_tableau import build_explicit_euler_tableau
# Burgers: du/dt + d(u^2/2)/dx = 0, with a Riemann problem (u_L=1, u_R=-0.25) that
# forms a shock. FiniteVolumeScheme gives the right-hand side of the time scheme.
class BurgersPDE(AbstractConservativePDE):
    def __init__(self):
        super().__init__(dim=1)
    def construct_F(self, u):
        return ParamScalarFunction(u.dims, lambda sp, x: 0.5 * u(sp, x) ** 2, f_type=u.f_type)
    def construct_A(self, u):
        return None
    def construct_R(self, u):
        return None
model = AbstractPhysicalConservativeModel(BurgersPDE())
model.add_boundary_condition("west", Dirichlet(lambda _x: jnp.array([1.0])))
model.add_boundary_condition("east", Dirichlet(lambda _x: jnp.array([-0.25])))
flux = RusanovFlux(wave_speed=lambda _model, u, _context: jnp.abs(u[0]))
spatial_scheme = FiniteVolumeScheme(
    model,
    VariablesFV(cartesian_mesh([40], quad_order=2)),
    flux,
    assemble_second_order=False,
    assemble_reaction=False,
    assemble_source=False,
)
scheme = TimeDiscreteFVscheme(spatial_scheme, build_explicit_euler_tableau(), dt=0.4 / 40)
initial = scheme.initialize(
    lambda x: jnp.where(x[0] < 0.5, jnp.array([1.0]), jnp.array([-0.25]))
)
final, history = scheme.solve(initial, t0=0.0, nt=20)
# Explicit Euler with dt <= CFL h / max wave speed. 40 cells, CFL=0.4, t=0.2:
# the mass is exactly conserved (0.46875), without overshoot.

Since finite volumes do not use basis functions, the Basis functions page does not apply to them: the mesh alone defines the discretization. A nonlinear diffusion coefficient is handled in the same way as for FEM and DG (see Nonlinear solvers).