Operator Assembly¶

This page documents how the finite-volume operators gaussGreenGrad, gaussGreenDiv, and gaussGreenLaplacian assemble contributions into field vectors and linear systems. It is a companion to Matrix Addressing which explains the CSR/COO indexing used throughout.

Source files:

src/finiteVolume/cellCentred/operators/gaussGreenGrad.cpp
src/finiteVolume/cellCentred/operators/gaussGreenDiv.cpp
src/finiteVolume/cellCentred/operators/gaussGreenLaplacian.cpp

Conventions¶

Face orientation

The face-area vector \(\mathbf{S}_f\) points from the owner cell to the neighbour cell. By construction, owner[f] < neighbour[f] for all internal faces.

Therefore \(F_f > 0\) means flux leaves the owner cell and enters the neighbour cell. This sign convention is only valid because \(\mathbf{S}_f\) points from owner to neighbour.

Interpolation weight

weightsV[f] (alias w) is the weight for the owner cell, so the face-centred value is

\[\varphi_f = w \, \varphi_\text{own} + (1-w) \, \varphi_\text{nei}.\]

Weights are produced by SurfaceInterpolation::weight(); for linear interpolation this is the geometric inverse-distance weight biased toward the owner.

Matrix entry naming

Following the LDU convention (see Matrix Addressing):

upperIdx(own, f) = rowOffs[own] + ownerOffset[f] → position of \(A[\text{own},\text{nei}]\) (upper-triangular, own < nei).
lowerIdx(nei, f) = rowOffs[nei] + neighbourOffset[f] → position of \(A[\text{nei},\text{own}]\) (lower-triangular, nei > own).

Operator kernels use the direct offset form for performance:

// A[own, nei] — upper triangular
values[rowOwnStart + ownOffs[facei]] ...
// values[matIt->upperIdx(own, facei)] ...   // equivalent via matrixIterator

// A[nei, own] — lower triangular
values[rowNeiStart + neiOffs[facei]] ...
// values[matIt->lowerIdx(nei, facei)] ...   // equivalent via matrixIterator

Assembly calling order

For implicit operators: internal faces → physical boundary faces. The physical boundary loop runs over facei ∈ [nInternalFaces, sGamma.size()).

gaussGreenGrad — Explicit Only¶

GaussGreenGrad provides only an explicit (field-to-field) gradient. There is no implicit (matrix-assembling) overload; calling dsl::solve with a grad term will raise an error at runtime.

Algorithm¶

Interpolate phi to faces: phif = surfInterp.interpolate(phi).

Loop over internal faces i ∈ [0, nInternalFaces):

Vec3 flux = faceAreaS[i] * phif[i];    // S_f * φ_f (S_f points owner→nei by construction)
grad[owner[i]]    += flux;              // +S_f * φ_f (S_f is outward from owner)
grad[neighbour[i]] -= flux;             // −S_f * φ_f (S_f is inward to neighbour)

Loop over boundary faces i ∈ [nInternalFaces, phif.size()):

Vec3 flux = faceAreaS[i] * phif[i];
grad[own] += flux;                      // +S_f * φ_f (S_f outward from owner; no neighbour)

Normalise: grad[celli] *= operatorScaling[celli] / vol[celli].

The result is the standard Green–Gauss gradient:

\[(\nabla \varphi)_C \approx \frac{1}{V_C} \sum_f \mathbf{S}_f \, \varphi_f.\]

S_f points from owner to neighbour, so it is the outward area vector for the owner and the inward area vector for the neighbour. The subtract for the neighbour follows directly.

gaussGreenDiv — Explicit¶

computeDivExp interpolates phi to faces and calls the shared computeDiv helper. The result is the positive mathematical divergence \(+\nabla \cdot (F\varphi)\).

Note

Consistent sign convention. Both the explicit computeDiv and the implicit computeDivImp assemble \(+\nabla \cdot (F\varphi)\) — the positive divergence operator — matching OpenFOAM’s fvm::div convention. The DSL equation form mirrors OpenFOAM:

dsl::solve(ddt(phi) + div(faceFlux, phi) - laplacian(gamma, phi) == source, ...)

Internal faces¶

ValueType flux = faceFlux[i] * phif[i];    // F_f * φ_f
div[owner[i]]     += flux;                  // F_f is outward from owner (S_f points owner→nei)
div[neighbour[i]] -= flux;                  // F_f is inward to neighbour → subtract

Sign rationale:

F_f > 0 means flux leaves the owner cell (S_f points from owner to neighbour by construction).
The divergence at owner is the net outward flux: a positive outward \(F_f\) gives a positive divergence contribution at the owner cell. This is why div[own] += flux is correct — it does not contradict the fact that the owner is losing \(\varphi\). The conservation equation uses \(\partial\varphi/\partial t = -\nabla\cdot(F\varphi)\), where the sign flip is applied externally (or in the implicit operator).
For the neighbour, \(F_f\) is inward (S_f points away from neighbour), so the outward flux from the neighbour’s perspective is \(-F_f\), hence div[nei] -= flux.

Boundary faces¶

ValueType flux = faceFlux[i] * phif[i];
div[own] += flux;           // F_f outward from owner (no neighbour on this rank)

After both loops: div[celli] *= operatorScaling[celli] / vol[celli].

gaussGreenDiv — Implicit¶

computeDivImp assembles \(+\nabla \cdot (F\varphi)\) into a LinearSystem, identical in sign to computeDivExp and to OpenFOAM’s fvm::div.

Divergence convention¶

With \(F_f > 0\) meaning flux from owner P to neighbour N (valid because \(\mathbf{S}_f\) points from owner to neighbour):

Owner row P: total contribution +F_f (\(\mathbf{S}_f\) is outward from P → positive divergence at P).
Neighbour row N: total contribution −F_f (\(\mathbf{S}_f\) is inward to N from this face → negative contribution at N).

Global verification: \(+F_f - F_f = 0\) per face ✓

Internal faces¶

The face-interpolated transported quantity is:

\[\varphi_f = w \, \varphi_\text{own} + (1-w) \, \varphi_\text{nei}\]

The face flux \(F_f \varphi_f\) decomposes linearly into owner and neighbour contributions:

\[F_f \, \varphi_f = w \, F_f \, \varphi_\text{own} + (1-w) F_f \, \varphi_\text{nei}\]

In the code the variables ownFluxContrib and neiFluxContrib carry a pre-applied minus sign so that uniform +=/-= patterns in the assembly loops still produce the correct positive divergence:

ownFluxContrib = -w * F_f           // note: NEGATIVE; code uses this in += and atomic_sub
neiFluxContrib = -(1-w) * F_f       // note: NEGATIVE; code uses this in -= and atomic_add

Implicit divergence — internal face entries¶
Matrix entry	Value	Sign note	Code
A[nei, own] (lower tri.)	\(-w \cdot F_f \cdot \text{scalingNei}\)	own value → negative contribution to nei’s divergence (inward face for N)	`values[rowNeiStart + neiOffs[f]] += ownFluxContrib * scalingNei`
diag[own]	\(+w \cdot F_f \cdot \text{scalingOwn}\)	own value → positive contribution to P’s divergence (outward face for P)	`atomic_sub(&values[diagOwn], ownFluxContrib * scalingOwn)`
A[own, nei] (upper tri.)	\(+(1-w) \cdot F_f \cdot \text{scalingOwn}\)	nei value → positive contribution to P’s divergence	`values[rowOwnStart + ownOffs[f]] -= neiFluxContrib * scalingOwn`
diag[nei]	\(-(1-w) \cdot F_f \cdot \text{scalingNei}\)	nei value → negative contribution to N’s divergence (inward face for N)	`atomic_add(&values[diagNei], neiFluxContrib * scalingNei)`

Row-sum verification per face:

Owner row: \(+w F_f + (1-w) F_f = +F_f\) ✓
Neighbour row: \(-w F_f - (1-w) F_f = -F_f\) ✓

Note

The diagonal is assembled explicitly — there is no separate negSumDiag pass. Each face contributes \(+w F_f\) to the owner diagonal and \(-(1-w) F_f\) to the neighbour diagonal, so the full sums converge to \(+F_f\) and \(-F_f\) respectively.

Note

Matches OpenFOAM. OpenFOAM’s gaussConvectionScheme sets lower[f] = -w F_f, upper[f] = (1-w) F_f, then negSumDiag() gives diag[own] += w F_f and diag[nei] -= (1-w) F_f — the same structure as NeoN.

Physical boundary¶

The boundary loop runs over faces facei ∈ [nInternalFaces, faceFluxV.size()).

Boundary weight: flux = bweights[bcfacei] * faceFlux[facei], where bweights is the boundary interpolation weight for the owner cell.

Mixed-BC decomposition uses valueFraction[bcfacei] (\(f_D\)):

\(f_D = 1\) → pure Dirichlet (fixedValue): boundary value \(\varphi_b = \text{refValue}\) is fully prescribed, no implicit diagonal coupling.
\(f_D = 0\) → pure Neumann: gradient is prescribed (\(\nabla_n \varphi = \text{refGradient}\)), full diagonal coupling.

Implicit divergence — physical BC entries¶
Contribution	Value	Code
diag[own]	\(+\text{flux} \cdot \text{scalingOwn} \cdot (1-f_D)\)	`atomic_add(&values[diagOffs[own]], valueMat)`
rhs[own]	\(-\text{flux} \cdot \text{scalingOwn} \cdot (f_D \cdot \text{refValue} + (1-f_D) \cdot \text{refGrad} / \delta)\)	`atomic_sub(&rhs[own], valueRhs)`

where \(\delta\) is deltaCoeffs[bcfacei] = inverse distance from cell centre to boundary face centre.

Note

The diagonal term uses the Neumann fraction \((1-f_D)\). For pure Dirichlet the diagonal is unchanged and the full contribution is in the RHS; for pure Neumann the diagonal receives the flux and the RHS receives the gradient source.

gaussGreenLaplacian — Explicit¶

computeLaplacianExp delegates face-normal gradients to a FaceNormalGradient helper and accumulates the positive Laplacian \(+\nabla \cdot (\gamma \nabla \varphi)\).

Internal faces¶

fnGrad[f] is computed by FaceNormalGradient as:

fnGrad[f] = nonOrthDeltaCoeffs[f] * (phi[nei] − phi[own])

fnGrad > 0 when \(\varphi_N > \varphi_P\) (gradient points outward from owner, i.e. S_f direction). In this case, diffusion brings \(\varphi\) into the owner cell and out of the neighbour cell, which is the correct physical interpretation (diffusion flows from high to low concentration).

ValueType flux = faceArea[i] * fnGrad[i];   // |S_f| * (∂φ/∂n)_f   (outward from owner)
laplacian[owner[i]]     += flux;             // fnGrad>0: φ_N>φ_P, owner gains φ → positive
laplacian[neighbour[i]] -= flux;             // neighbour loses φ → negative

Boundary faces¶

ValueType flux = faceArea[i] * fnGrad[i];
laplacian[own] += flux;                      // S_f outward from owner; no neighbour on this rank

After both loops: laplacian[celli] *= operatorScaling[celli] / vol[celli].

gaussGreenLaplacian — Implicit¶

computeLaplacianImpl assembles the diffusion operator \(+\nabla \cdot (\gamma \nabla \varphi)\) into the linear system.

Internal faces¶

Face diffusion coefficient: flux = deltaCoeffs[facei] * sGamma[facei] * magFaceArea[facei], i.e. \(\delta_f \cdot \gamma_f \cdot |S_f|\).

The Laplacian is symmetric — the same flux value enters both the owner and neighbour rows with appropriate signs. \(\mathbf{S}_f\) points from owner to neighbour, so diffusion out of one cell equals diffusion into the other.

Implicit Laplacian — internal face entries¶
Matrix entry	Value	Sign note	Code
A[nei, own] (lower tri.)	\(+\text{flux} \cdot \text{scalingNei}\)	enters neighbour → positive	`values[rowNeiStart + neiOffs[facei]]`
diag[own]	\(-\text{flux} \cdot \text{scalingOwn}\)	leaves owner → negative	`atomic_sub(&values[diagOffs[own]], flux)`
A[own, nei] (upper tri.)	\(+\text{flux} \cdot \text{scalingOwn}\)	enters owner → positive	`values[rowOwnStart + ownOffs[facei]]`
diag[nei]	\(-\text{flux} \cdot \text{scalingNei}\)	leaves neighbour → negative	`atomic_sub(&values[diagOffs[nei]], flux)`

The assembled matrix is symmetric and negative semi-definite for positive \(\gamma\) and consistent operatorScaling, representing \(-\nabla \cdot (\gamma \nabla \varphi)\).

Note

Matches OpenFOAM. OpenFOAM’s gaussLaplacianScheme sets upper() = deltaCoeffs * gammaMagSf (positive) and calls negSumDiag() (diagonal = −sum of off-diagonals), yielding the same structure.

Physical boundary¶

The boundary loop covers facei ∈ [nInternalFaces, sGamma.size()).

Face diffusion coefficient (no deltaCoeffs here, unlike internal faces):

flux = sGamma[facei] * magFaceArea[facei]     // γ · |S_f|

Mixed-BC decomposition with valueFraction (\(f_D\)):

Implicit Laplacian — physical BC entries¶
Contribution	Value	Code
diag[own]	\(-\text{flux} \cdot \text{scalingOwn} \cdot f_D \cdot \delta_f\)	`atomic_sub(&values[diagOffs[own]], valueMat)`
rhs[own]	\(-\text{flux} \cdot \text{scalingOwn} \cdot (f_D \cdot \delta_f \cdot \text{refValue} + (1-f_D) \cdot \text{refGradient})\)	`atomic_sub(&rhs[own], valueRhs)`

where \(\delta_f\) = deltaCoeffs[facei] (all-faces array index, consistent throughout the boundary kernel).

The diagonal term uses the Dirichlet fraction \(f_D\):

Pure Dirichlet (\(f_D = 1\)): diag[own] -= flux * δ, providing full implicit coupling to \(\varphi_C\).
Pure Neumann (\(f_D = 0\)): no diagonal change; gradient source enters via the RHS only.

Note

The boundary flux = γ · |S_f| lacks the deltaCoeffs factor present in the internal face flux (δ · γ · |S_f|). The deltaCoeffs factor is applied explicitly in valueMat and valueRhs instead, so the total implicit coupling is still δ · γ · |S_f| · f_D. This differs from the internal face assembly where flux already absorbs δ.