Turbulence closure

Steady manufactured solution of the compressible RANS equations coupled to the negative Spalart-Allmaras model on the unit box (M = 0.5, Re = 1000, nu_tilde / nu = 8 (1 + 0.4 sin), a smooth positive manufactured wall distance between 0.1 and 0.8), Dirichlet data from the exact state on all sides. The forcing is generated by automatic differentiation of the same flux and source kernels the solver uses (turbulence.manufactured_source). Errors: L2 norm of the mean-flow conserved variables, and of rho nu_tilde scaled by Re. The solve starts from the nodal values of the exact solution and removes the truncation-error residual with p-multigrid Newton-Krylov. The eddy viscosity of the manufactured state is 4.7 times the molecular one, so the mean flow sees an effective Reynolds number of about 175: between the diffusion- and the convection-dominated regimes of the laminar manufactured solutions, where the state converges at a rate between p and p + 1. The laminar manufactured solution with the same mean state is therefore solved at that Reynolds number on the same meshes, and the checks require the guaranteed rate p of both equations and that the model changes neither the order nor the size of the mean-flow error.

• p = 1: order of the mean-flow state (finest pair)1.914expected >= 0.75 (guaranteed rate p = 1, tolerance 0.25; design rate 2)pass
• p = 1: order of nu_tilde (finest pair)2.202expected >= 0.75 (guaranteed rate p = 1, tolerance 0.25; design rate 2)pass
• p = 1: mean-flow order with the model minus the laminar order at the effective Reynolds number0.006866expected > -0.15pass
• p = 1: mean-flow error with the model over the laminar error at the effective Reynolds number (finest mesh)1.024expected 0.5 .. 1.5pass
• p = 2: order of the mean-flow state (finest pair)2.477expected >= 1.75 (guaranteed rate p = 2, tolerance 0.25; design rate 3)pass
• p = 2: order of nu_tilde (finest pair)2.612expected >= 1.75 (guaranteed rate p = 2, tolerance 0.25; design rate 3)pass
• p = 2: mean-flow order with the model minus the laminar order at the effective Reynolds number-0.07206expected > -0.15pass
• p = 2: mean-flow error with the model over the laminar error at the effective Reynolds number (finest mesh)1.148expected 0.5 .. 1.5pass
• every model term is active at the domain centre (smallest |term| Re)0.00169expected > 1e-3pass

gyre 0.4.0 · 39979fcf90 · 2026-11-12 11:33 · cpu · quick suite · record (JSON)Same code as this build

Jacobian-vector products of the weak RANS-SA residual from automatic differentiation against central differences (best of three step sizes), the dot-product identity of tangent and adjoint products, the element-diagonal blocks from the split-argument residual against the dense Jacobian (2-D), and the derivative with respect to the geometry nodes, which includes the wall distance. States are random perturbations of the free stream with nu_tilde / nu between -2 and 12, so both branches of the model and its limiters are differentiated; every boundary-condition type appears.

• 2-D ring, adiabatic wall, far field: JVP vs central difference4.13e-10expected < 1e-6pass
• 2-D ring, adiabatic wall, far field: JVP vs central difference, turbulence row3.95e-09expected < 1e-6pass
• 2-D ring, adiabatic wall, far field: <v,Ju> = <J^T v,u>7.51e-15expected < 1e-12pass
• 2-D ring, adiabatic wall, far field: shape derivative of the residual vs central difference1.66e-10expected < 1e-6pass
• 2-D ring, adiabatic wall, far field: derivative of the wall distance w.r.t. the nodes1.45e-09expected < 1e-6pass
• 2-D ring, adiabatic wall, far field: share of nodes on the negative branch0.1759expected 0.05 .. 0.5pass
• 2-D ring, adiabatic wall, far field: block inverse times dense diagonal block9.06e-14expected < 1e-8pass
• 2-D channel: isothermal wall, inlet, outlet, slip wall: JVP vs central difference2.06e-10expected < 1e-6pass
• 2-D channel: isothermal wall, inlet, outlet, slip wall: JVP vs central difference, turbulence row1.38e-09expected < 1e-6pass
• 2-D channel: isothermal wall, inlet, outlet, slip wall: <v,Ju> = <J^T v,u>2.33e-14expected < 1e-12pass
• 2-D channel: isothermal wall, inlet, outlet, slip wall: shape derivative of the residual vs central difference7.91e-11expected < 1e-6pass
• 2-D channel: isothermal wall, inlet, outlet, slip wall: derivative of the wall distance w.r.t. the nodes3.73e-10expected < 1e-6pass
• 2-D channel: isothermal wall, inlet, outlet, slip wall: share of nodes on the negative branch0.2099expected 0.05 .. 0.5pass
• 2-D channel: isothermal wall, inlet, outlet, slip wall: block inverse times dense diagonal block8.23e-13expected < 1e-8pass
• 3-D box: adiabatic wall, far field, symmetry: JVP vs central difference2.66e-10expected < 1e-6pass
• 3-D box: adiabatic wall, far field, symmetry: JVP vs central difference, turbulence row6.48e-09expected < 1e-6pass
• 3-D box: adiabatic wall, far field, symmetry: <v,Ju> = <J^T v,u>1.77e-14expected < 1e-12pass
• 3-D box: adiabatic wall, far field, symmetry: shape derivative of the residual vs central difference6.76e-11expected < 1e-6pass
• 3-D box: adiabatic wall, far field, symmetry: derivative of the wall distance w.r.t. the nodes4.21e-10expected < 1e-6pass
• 3-D box: adiabatic wall, far field, symmetry: share of nodes on the negative branch0.1296expected 0.05 .. 0.5pass

gyre 0.4.0 · 39979fcf90 · 2026-11-12 11:40 · cpu · quick suite · record (JSON)Same code as this build

Zero-pressure-gradient turbulent flat plate of the NASA Turbulence Modeling Resource (2D Zero Pressure Gradient Flat Plate Verification Case: M = 0.2, Re = 5e6 per unit length, plate from x = 0 to 2, symmetry upstream to x = -0.33, far field at y = 1, total-pressure inflow, pressure outflow, adiabatic wall, free-stream nu_tilde = 3 nu). Each run starts from the uniform free stream and is converged nine orders by the pseudo-transient Newton-Krylov solver. The skin friction, taken from the numerical wall flux of the residual, is compared with the TMR's SA results on its finest grid (545 x 385 points; file FlatPlate/SA/cf_plate.dat, mean of CFL3D and FUN3D, which differ by less than 0.1 % at the stations used; at x = 0.97008 their grid-converged value is cf = 2.705513e-03). The meshes here are far coarser: the coarse mesh has 15 x 10 elements with a first element height of 5e-05. The comparison is with the reference model implementation (standard SA with the ft2 term; this code solves SA-neg, identical for positive nu_tilde), not with experiment.

• coarse mesh, p = 1: converged from the free stream (residual drop)8.21e-10expected <= 1e-9pass
• coarse mesh, p = 1: Newton steps33expected <= 100 (cap 150)pass
• coarse mesh, p = 2: converged from the free stream (residual drop)2.31e-11expected <= 1e-9pass
• coarse mesh, p = 2: Newton steps34expected <= 100 (cap 150)pass
• p = 2: |cf / cf_TMR - 1| at x = 0.970.001656expected < 0.02pass
• p = 2: largest |cf / cf_TMR - 1| over the stations x >= 0.50.02153expected < 0.04pass
• p = 1: largest |cf / cf_TMR - 1| over the stations x >= 0.50.09989expected < 0.15pass
• raising the degree from 1 to 2 reduces the largest difference (ratio)0.2155expected < 0.5pass

Run notes The leading-edge region (x < 0.3) is the least resolved: the first element on the plate is 0.02 long. Negative nu_tilde occurs at the edge of the boundary layer, where the negative branch of the model takes over; it is reported, not clipped. Not covered: a mesh-converged result, the velocity profile in wall units, other Reynolds numbers.

gyre 0.4.0 · 39979fcf90 · 2026-11-12 11:36 · cpu · quick suite · record (JSON)Same code as this build

The turbulent flat plate of the preset flat_plate_sa (M = 0.2, Re = 5e6, degree 2, first cell 5e-5 high) is solved on its 2-D section and on a spanwise-periodic extrusion of that section into hexahedra (3 across a span of 3), both from the free stream with line Jacobi as the smoother of the p-multigrid levels 2-1-0. A correct 3-D implementation of the RANS equations (viscous and turbulent fluxes, the source terms of the model, the wall distance, the boundary conditions, the line smoother, force integration and the discrete adjoint) must give the drag of the section per unit span, its skin friction and the adjoint derivative of the drag. Section and extrusion have the same in-plane equations only at the same BR2 penalty, whose default is the number of faces of an element (4 for quadrilaterals, 6 for hexahedra): the section is therefore solved with 6, the default of the extrusion, and once more with its own default, which is reported without a check. Three kinds of statement are checked. Before the 3-D solve, the solution of the section is extruded and put into the 3-D residual: the equations the section has must be satisfied as well as the section satisfies them (below 1e-8 of the free-stream residual; the flow solves stop at a drop of 1e-9), which tests the 3-D operator without the 3-D solver. The spanwise momentum equation is not satisfied by the extruded solution, for the reason given in the notes, so that the 3-D solution has a small spanwise velocity that alternates within each element: it must leave the layers identical and mirror-symmetric and exert no net spanwise force (round-off only: below 1e-8). Being odd across the mid-plane of an element, it cannot change the mean over the span at first order: drag, mean skin friction and the derivative must agree with the section to 1e-6. The size of the spanwise velocity itself is reported, not checked.

• the three flow solves converged from the free streamyesexpected truepass
• the solution of the section, extruded, in the 3-D discretization: residual of the equations of the section, relative to the residual of the free stream6.84e-10expected < 1e-8pass
• drag coefficient per unit span: |3-D - 2-D| / |2-D|5.23e-09expected < 1e-6pass
• skin friction at x = 0.97008, mean over the span: |3-D - 2-D| / |2-D|3.86e-09expected < 1e-6pass
• transposed solves reached their tolerance of 1e-10 (section and extrusion)yesexpected truepass
• adjoint derivative of the drag with respect to the Mach number: |3-D - 2-D| / |2-D|2.42e-11expected < 1e-6pass
• the spanwise layers repeat each other: largest difference of the state between a layer and the first7.69e-12expected < 1e-8pass
• every layer is symmetric about its mid-plane (spanwise momentum odd, the other variables even): largest defect1.08e-11expected < 1e-8pass
• no net spanwise force: |cs|3.08e-19expected < 1e-8pass
• lines of the extrusion: those of the section in each spanwise layer (elements in lines)393expected = 3 x 131, in 3 x 18 linespass

Run notes The 3-D solution is not exactly uniform along the span: its spanwise momentum reaches 6.9e-07 and the in-plane variables differ by up to 7.4e-07 between spanwise nodes (table 'span'). This is a property of the BR2 scheme, not of the extrusion: the viscous flux inside an element is evaluated with the gradient lifted from all its faces, the flux through a face with the local gradient and the lifting of that face alone. Where the in-plane velocity jumps across in-plane faces, the two differ in the normal stress on the spanwise faces, which leaves a residual in the spanwise momentum equation of every element (table 'span': the extruded solution of the section). It changes sign across the mid-plane of the element, which is why the layers repeat each other and there is no net spanwise force; how little the mean over the span is affected is what the comparison with the section shows. With the default penalty of quadrilaterals (4) the section gives cd = 5.753647e-03 and cf = 2.701031e-03 at the station; the extrusion, whose default is 6, differs from these by 0.26 % and 0.16 %. That is a difference between two discretizations on a mesh of 150 elements, not an error of the 3-D implementation, and it is why the comparison is made at equal penalty. The extruded mesh has 54 lines holding 393 of its 450 elements (the longest 12 elements, 3 colours); 72900 unknowns in blocks of 162. An extrusion has only faces whose two sides use the same local axes, and the flow has no spanwise component to speak of: the study does not exercise the face-orientation tables, the two further components of the vorticity in the model or a flow that varies along the span.

gyre 0.4.0 · 39979fcf90 · 2026-11-12 11:33 · gpu · NVIDIA GeForce RTX 3060 Ti · quick suite · record (JSON)Same code as this build

Total derivatives of the drag coefficient of the turbulent flat plate (NASA TMR conditions, M = 0.2, Re = 5e6, coarse mesh of 150 elements, p = 1, 3000 unknowns, cd = 6.284664e-03) from one discrete-adjoint solve of the coupled RANS-SA system, against central differences of re-converged solutions. Inputs: Mach number, Reynolds number, angle of the inflow, free-stream nu_tilde / nu, and a wall-normal bump of the plate surface (a shape derivative, which includes the sensitivity of the wall distance).

• dcd/dmach: adjoint vs central difference (rel.)5.41e-09expected < 1e-5pass
• dcd/dre: adjoint vs central difference (rel.)4.01e-09expected < 1e-5pass
• dcd/dalpha: adjoint vs central difference (rel.)6.55e-11expected < 1e-5pass
• dcd/dchi_inf: adjoint vs central difference (rel., best of 3 steps)3.18e-08expected < 1e-5pass
• dcd/dshape (wall bump, directional): adjoint vs central difference (rel.)1.79e-08expected < 1e-5pass
• adjoint linear solve reached its tolerance (residual reduction)9.21e-11expected < 1e-10pass

Run notes The mesh is wall resolved, with element aspect ratios of the order of 1000 at the wall. Element-block Jacobi does not see that anisotropy: as the only preconditioner it left the adjoint solve short of its tolerance (see the comment in the study); the piecewise-constant coarse level of the p-multigrid hierarchy is what makes it converge. There is no line-implicit or incomplete-factorisation preconditioner.

gyre 0.4.0 · 39979fcf90 · 2026-11-12 11:43 · cpu · quick suite · record (JSON)Same code as this build