3 Implementation Conventions and Coupled Jacobian

This chapter records conventions that are easy to misread when moving between the mathematical equations and the C++ implementation. They are common to the active compressible and incompressible coupled solvers unless stated otherwise.

3.1 Primitive columns and conservative equation rows

Let \(d\) be the spatial dimension. The core aerodynamic primitive vector is \[\mathbf{q}_i=[u_1,\ldots,u_d,p,T]^T_i. \](3.1) The AD directions used to construct the local Jacobian are seeded in these primitive variables. Therefore the matrix columns are primitive-variable derivatives; a post-assembly conservative-to-primitive chain transformation is not required for this path.

The row slots use the same indices for storage convenience, but their physical meaning is different:

When turbulence, particles, or species are active, their equations and primitive corrections are appended according to the compile-time equation layout selected by the matrix driver.

3.2 Residual and matrix sign

Define the finite-volume residual by Eq. (2.8). The current assembly accumulates a face contribution \(\mathbf{F}_fA_f\) as \[\mathbf{B}_P \mathrel{-}= \mathbf{F}_f A_f, \qquad \mathbf{A}_{Pj} \mathrel{+}= \frac{\partial(\mathbf{F}_fA_f)}{\partial\mathbf{q}_j}.\] A physical source \(\mathbf{S}_P\Omega_P\) has the opposite residual sign and is therefore accumulated as \[\mathbf{B}_P \mathrel{+}= \mathbf{S}_P\Omega_P, \qquad \mathbf{A}_{Pj} \mathrel{-}= \frac{\partial(\mathbf{S}_P\Omega_P)}{\partial\mathbf{q}_j}.\] Thus \[\boxed{ \mathbf{B}=-\mathbf{R}, \qquad \mathbf{A}=-\frac{\partial\mathbf{B}}{\partial\mathbf{q}} =\frac{\partial\mathbf{R}}{\partial\mathbf{q}}} \](3.2) and the linear correction satisfies \[\boxed{\mathbf{A}\,\delta\mathbf{q}=\mathbf{B}.} \](3.3) After relaxation/limiting, the primitive variables are updated by addition of the solved increment.

3.3 Face-local AD and sparse blocks

For an interior face shared by cells \(L\) and \(R\), the numerical flux depends on both reconstructed states, \[\mathbf{F}_f=\mathbf{F}_f(\mathbf{q}_L,\mathbf{q}_R;\mathbf{g}_f),\] where \(\mathbf{g}_f\) denotes geometry and any prepared coefficients. Forward-mode AD evaluates the local derivatives \[\mathbf{J}_{f,L}=\frac{\partial(\mathbf{F}_fA_f)}{\partial\mathbf{q}_L}, \qquad \mathbf{J}_{f,R}=\frac{\partial(\mathbf{F}_fA_f)}{\partial\mathbf{q}_R}.\] These become the owner and neighbour blocks in the global sparse matrix. Contributions from every face and every active source associated with the same owner cell are accumulated before the owner block is inserted.

This distinction matters for nonlinear methods: “exact Jacobian” in this guide means exact differentiation of the AD-coded residual with the prepared outer coefficients held at the values used by that assembly. Quantities deliberately refreshed outside the AD expression are not silently differentiated through an additional outer iteration.

3.4 Stored state and pseudo-time state

CMPS maintains primitive variables and conservative or conservative-like histories. The time derivative need not be the identity in primitive variables. Let \[\mathbf{W}=\mathbf{W}(\mathbf{q}), \qquad \boldsymbol{\Gamma}=\frac{\partial\mathbf{W}}{\partial\mathbf{q}}.\] A local pseudo-time contribution has the generic form \[\frac{\Omega_P}{\Delta\tau_P}\boldsymbol{\Gamma}_P\,\delta\mathbf{q}_P.\] For the compressible solver, \(\boldsymbol{\Gamma}\) may be modified by low-Mach time-derivative preconditioning. For the constant-density incompressible solver, the pressure row contains the artificial-compressibility factor \(1/c_{ac}^{2}\) while the stored pressure component remains \(p\).

3.5 Spatial versus physical-time residuals

In steady calculations, pseudo time is a convergence device: it stabilizes the Newton-like iteration but does not add a physical unsteady residual. In transient calculations, a physical BDF residual is added to the spatial residual and dual-time iterations converge that physical step.

For a stored physical state \(\mathbf{W}\), \[\mathbf{R}^{\mathrm{time}}_P =\frac{\Omega_P}{\Delta t} \left(e_0\mathbf{W}^{n+1}-e_1\mathbf{W}^n+e_2\mathbf{W}^{n-1}\right), \](3.4) with \[\begin{aligned} (e_0,e_1,e_2)&=(1,1,0) &&\text{BDF1},\\ (e_0,e_1,e_2)&=\left(\frac32,2,\frac12\right) &&\text{BDF2}. \end{aligned}\] For artificial-compressibility incompressible flow, the pressure equation has no physical pressure-storage term; its pressure time derivative remains a pseudo-time coupling device.

3.6 Why the conventions matter

Several common implementation errors can be avoided by keeping these identities explicit:

• a pressure-index row is continuity, not a pressure conservation equation;

• \(\mathbf{A}\delta\mathbf{q}=\mathbf{B}\) already contains the Newton minus sign because \(\mathbf{B}=-\mathbf{R}\);

• the artificial-compressibility coefficient belongs to the pseudo-time mass matrix and must not be confused with the stored pressure state;

• a coarse-grid or machine-learning path must use the same row/column ordering and sign convention as the main residual if it reuses the CMPS matrix.