Finite volume¶

Overview¶

The finite_volume module provides composable building blocks for solving 1D PDEs with the cell-centered finite volume method. The primitives batch over arbitrary leading dimensions and plug into TransientDriver, ImplicitUpdate, and ComposedModel like any other NEML2 model, so finite-volume transport models compose with the rest of the framework through the same interfaces a constitutive model uses.

The primitives generalize to a range of cell-centered finite-volume problems; the provided examples and regression tests focus on 1D transport (advection-diffusion-reaction). The module ships:

Discretization operators — LinearlyInterpolateToCellEdges, FiniteVolumeGradient, FiniteVolumeUpwindedAdvectiveFlux.
Boundary-condition utilities — FiniteVolumeAppendBoundaryCondition.
Coordinate and source helpers — SemiInfiniteCoordinateTransform, SmearedDeltaSource, DumpInSmallestBin.

Composition (via ComposedModel) wires these together with generic algebraic models like ScalarLinearCombination and time-integration primitives like ScalarBackwardEulerTimeIntegration to produce a full discretized transport equation that ImplicitUpdate then advances in time.

Math¶

Consider a conserved quantity \(u(x,t)\) with total flux \(J\) and reaction source \(R\) on a domain \(\Omega = [a,b]\), governed by

\[ u_t + J_x = R. \]

Partition \([a,b]\) into \(N\) cells

\[ I_i = \left[x_{i-\tfrac{1}{2}}, x_{i+\tfrac{1}{2}}\right], \quad i = 1, \dots, N, \]

with widths \(\Delta x_i = x_{i+\tfrac{1}{2}} - x_{i-\tfrac{1}{2}}\) and cell centers \(x_i\). The cell average is

\[ \bar{u}_i = \frac{1}{\Delta x_i}\int_{I_i} u\, dx. \]

Integrating the PDE over \(I_i\) and dividing by \(\Delta x_i\) yields the canonical semi-discrete finite-volume update

\[ \frac{d \bar{u}_i}{dt} + \frac{1}{\Delta x_i}\left(J\big\rvert_{x_{i+\tfrac{1}{2}}} - J\big\rvert_{x_{i-\tfrac{1}{2}}}\right) = \bar{R}_i, \]

which the catalog reconstructs in three pieces — edge values, gradients, and advective fluxes — each provided by its own primitive.

Cell-edge interpolation¶

For nonuniform grids, cell-edge values are obtained by linear interpolation between adjacent cell centers,

\[ q_{i+\tfrac{1}{2}} = \frac{x_{i+1} - x_{i+\tfrac{1}{2}}}{x_{i+1} - x_i}\, q_i + \frac{x_{i+\tfrac{1}{2}} - x_i}{x_{i+1} - x_i}\, q_{i+1}. \]

This is implemented by LinearlyInterpolateToCellEdges, and is used to lift cell-centered diffusivities or velocities onto edges where the fluxes live.

Gradients¶

Gradient (and, in 1D, divergence) terms use the first-order expression

\[ \mathrm{grad}_{u, i+\tfrac{1}{2}} = -p_{i+\tfrac{1}{2}}\, \frac{\bar{u}_{i+1} - \bar{u}_i}{x_{i+1} - x_i}, \]

where \(p\) is an optional cell-edge prefactor (e.g. a diffusivity for Fickian diffusion). FiniteVolumeGradient computes this for both diffusive fluxes (edge prefactor = \(D\), spacing = cell-center spacing) and flux divergences (no prefactor, spacing = cell widths).

Advective flux¶

Advective fluxes are stabilized with first-order upwinding,

(5)¶\[\begin{align} \hat{J}_{\mathrm{adv}, i+\tfrac{1}{2}} &= v^{+}_{i+\tfrac{1}{2}}\, \bar{u}_i + v^{-}_{i+\tfrac{1}{2}}\, \bar{u}_{i+1}, \\ v^{\pm}_{i+\tfrac{1}{2}} &= \tfrac{1}{2}\left(v_{i+\tfrac{1}{2}} \pm \left|v_{i+\tfrac{1}{2}}\right|\right), \end{align}\]

and computed by FiniteVolumeUpwindedAdvectiveFlux from cell-centered \(\bar{u}\) and an edge velocity \(v_{\mathrm{edge}}\).

Boundary conditions¶

Dirichlet and Neumann boundary conditions are imposed by appending a prescribed value to the left or right end of an intermediate tensor with FiniteVolumeAppendBoundaryCondition. Applied to the flux array \(J\), this realizes a Neumann condition; applied to \(u\) on a ghost edge, it realizes a Dirichlet condition. The two ends are configured by chaining the primitive twice, once per side.

Example: combined advection-diffusion-reaction¶

The following input file discretizes

\[ \partial_t u + \partial_x\!\left(v u - D\, \partial_x u\right) = -k\, u, \qquad u(x, 0) = e^{-\frac{(x - 0.25)^2}{2\,(0.05)^2}}, \]

on \(x \in [0, 1]\) with \(v = 0.4\), \(D = 0.5\), \(k = 0.05\), zero-flux boundary conditions, and backward-Euler time integration:

# neml2
[Tensors]
  [edges]
    # Coarsened for regression-test runtime: 201 -> 51 (50 cells). Full
    # resolution lives under tests/regression/; here we only detect drift.
    type = Python
    expr = 'Scalar.linspace(0.0, 1.0, 51).sub_batch.retag(1)'
  []
  [centers]
    type = Python
    expr = 'Scalar(0.5 * (edges.data[..., :-1] + edges.data[..., 1:]), sub_batch_ndim=1)'
  []
  [dx_centers]
    type = Python
    expr = 'Scalar(centers.data[..., 1:] - centers.data[..., :-1], sub_batch_ndim=1)'
  []
  [dx]
    type = Python
    expr = 'Scalar(edges.data[..., 1:] - edges.data[..., :-1], sub_batch_ndim=1)'
  []
  [ic]
    type = Python
    expr = 'Scalar(1.0 * torch.exp(-0.5 * ((centers.data - 0.25) / 0.05) ** 2), sub_batch_ndim=1)'
  []
  [time]
    type = Python
    expr = 'Scalar.linspace(0.0, 1.0, 25)'
  []
  [D_cells]
    type = Python
    expr = 'Scalar(0.5 * torch.ones_like(centers.data), sub_batch_ndim=1)'
  []
  [v_cells]
    type = Python
    expr = 'Scalar(0.4 * torch.ones_like(centers.data), sub_batch_ndim=1)'
  []
[]

[Drivers]
  [driver]
    type = TransientDriver
    model = 'model'
    prescribed_time = 'time'
    ic_Scalar_names = 'concentration'
    ic_Scalar_values = 'ic'
    save_as = 'result.pt'
  []
  [regression]
    type = TransientRegression
    driver = 'driver'
    reference = 'gold/result.pt'
  []
[]

[Models]
  [diffusivity]
    type = LinearlyInterpolateToCellEdges
    cell_values = 'D_cells'
    cell_centers = 'centers'
    cell_edges = 'edges'
    edge_values = 'D'
  []
  [advection_velocity]
    type = LinearlyInterpolateToCellEdges
    cell_values = 'v_cells'
    cell_centers = 'centers'
    cell_edges = 'edges'
    edge_values = 'v_edge'
  []
  [diffusive_flux]
    type = FiniteVolumeGradient
    u = 'concentration'
    prefactor = 'D'
    dx = 'dx_centers'
  []
  [advective_flux]
    type = FiniteVolumeUpwindedAdvectiveFlux
    u = 'concentration'
    v_edge = 'v_edge'
  []
  [reaction]
    type = ScalarLinearCombination
    from = 'concentration'
    to = 'R'
    weights = '-0.05'
  []
  [total_flux]
    type = ScalarLinearCombination
    from = 'grad_u flux'
    to = 'J'
    weights = '1 1'
  []
  [left_bc]
    type = FiniteVolumeAppendBoundaryCondition
    input = 'J'
    bc_value = 0.0
    side = 'left'
  []
  [right_bc]
    type = FiniteVolumeAppendBoundaryCondition
    input = 'J_with_bc_left'
    bc_value = 0.0
    side = 'right'
  []
  [flux_divergence]
    type = FiniteVolumeGradient
    u = 'J_with_bc_left_with_bc_right'
    dx = 'dx'
    grad_u = 'flux_div'
  []
  [rate_of_change]
    type = ScalarLinearCombination
    from = 'R flux_div'
    to = 'concentration_rate'
    weights = '1 1'
  []
  [integrate_u]
    type = ScalarBackwardEulerTimeIntegration
    variable = 'concentration'
  []
  [implicit_rate]
    type = ComposedModel
    models = 'diffusivity advection_velocity diffusive_flux advective_flux reaction total_flux left_bc right_bc flux_divergence rate_of_change integrate_u'
  []
[]

[EquationSystems]
  [eq_sys]
    type = NonlinearSystem
    model = 'implicit_rate'
    unknowns = 'concentration'
  []
[]

[Solvers]
  [newton]
    type = Newton
    linear_solver = 'lu'
  []
  [lu]
    type = DenseLU
  []
[]

[Models]
  [predictor]
    type = ConstantExtrapolationPredictor
    unknowns_Scalar = 'concentration'
  []
  [model]
    type = ImplicitUpdate
    equation_system = 'eq_sys'
    solver = 'newton'
    predictor = 'predictor'
  []
[]

Explanation¶

The example builds a single implicit-update model out of three layers.

Mesh and field tensors ([Tensors]) — edges is a 1D Scalar carrying the cell edges with sub_batch_ndim=1. That puts the spatial axis in sub-batch, so any leading batch dimension (parameter sweep, ensemble of simulations) composes on top without any code change here. centers and dx_centers derive cell-center positions and center-to-center spacings; dx is the cell-width vector. The initial condition ic is a Gaussian centered at \(x = 0.25\), and D_cells, v_cells populate cell-centered diffusivity and velocity.

Edge lift — two LinearlyInterpolateToCellEdges instances map D_cells and v_cells from cell centers onto cell edges, producing the edge tensors D and v_edge that the flux operators expect.

Flux assembly — the diffusive flux is computed by FiniteVolumeGradient with prefactor = D and dx = dx_centers, giving \(-D\,\partial_x u\) at each interior edge. The advective flux \(v u\) is computed by FiniteVolumeUpwindedAdvectiveFlux. Both are scalars on the \((N-1)\)-sized edge axis. A ScalarLinearCombination sums them into the total flux J.

Boundary conditions — J lives on the \(N-1\) interior edges, so two FiniteVolumeAppendBoundaryCondition calls append zero-flux values at the left and right ends, producing J_with_bc_left_with_bc_right on the full \(N+1\)-edge axis.

Reaction and divergence — a ScalarLinearCombination with weight \(-k\) gives the reaction term \(R = -k u\). A second FiniteVolumeGradient (no prefactor; dx = dx) computes the flux divergence \(J_x\) per cell. Their sum is the cell-centered rate of change.

Time integration and solve — ScalarBackwardEulerTimeIntegration turns the rate-of-change expression into a residual whose root in \(u\) is the backward-Euler update. The [EquationSystems] block wraps it as a NonlinearSystem in the unknown concentration, solved each step by Newton with the DenseLU linear solver and a ConstantExtrapolationPredictor initial guess. The outer ImplicitUpdate is what TransientDriver actually evaluates at every step of the prescribed time grid.

The data flow is therefore: concentration → {diffusive, advective} flux → J → bc → flux_div and concentration → R, summed into concentration_rate, then implicitly integrated to advance concentration to the next step.