16  Laminar Cylinder Vortex Shedding — Time Integration and Precision Comparison

Temporal order of accuracy of implicit dual-time BDF1/BDF2 against explicit RK4, and FP32/FP64 arithmetic precision

16.1 Run Provenance and Artifacts

This case combines two independent campaigns that answer different questions and should not be cross-compared as if they were one controlled sweep.

Item BDF1/BDF2 Δt sweep SLAU2 FP64 / FP32
Wake CFD baseline version 0.22.2 0.22.2
Solver DBNS DBNS
Compute device CPU GPU (NVIDIA RTX A6000)
Parallel execution 4 MPI ranks 1 MPI rank
Arithmetic precision FP64 FP64 / FP32
Run date 2026-08-29 – 2026-08-31 2026-06-04
Note

Updated 2026-09-01. The time-integration comparison below replaces a single-point BDF1/BDF2 result from an earlier build (0.13.1, 2026-06-25) with a five-point Δt sweep on a current build. A single point cannot distinguish “close to the reference” from “converging at the right rate,” so the earlier numbers are superseded rather than compared directly. The FP32/FP64 precision comparison is a separate, unrelated campaign and is unchanged.

16.1.1 Downloadable Files

Not available: standalone input decks and solver result files have not yet been published for download.

16.2 Verification Objective

This case asks two questions about Wake CFD’s implicit dual-time integration (BDF1, first order; BDF2, second order), holding mesh, physics, and build fixed:

  • Does it converge to the same answer as an independent time-integration path? The explicit-RK4 case solves the same flow with explicit classical RK4 — no pseudo-time loop, so it carries no dual-time convergence question of its own — and serves as the time-axis reference here.
  • Does it converge at the theoretically expected rate? Refining the physical step Δt over a five-point sweep (T/32 through T/512, T = the shedding period) and fitting |q(Δt) − q_ref| against Δt on a log-log plot gives the observed order of accuracy directly. Theory predicts a slope of 1 for BDF1 and 2 for BDF2.

Comparing against the explicit-RK4 case rather than the literature values directly is intentional: literature comparison mixes spatial discretization error and physical-model error into the same number as temporal error, so it cannot isolate the effect being tested here. That case uses the same mesh and build with a time step (4.29×10^-6 s) fine enough that its own temporal error is negligible, making it a clean reference for the time axis alone.

A second, independent question — FP32 versus FP64 arithmetic precision — is examined in the same case document because it was run on the same physical setup, but the two should be read separately (see Run Provenance).

16.3 Reference Solution and Data

16.3.1 Numerical Reference

The explicit-RK4 case provides the independent time-integration reference for this Re_D=100 Kármán vortex street: St=0.1661, mean C_D=1.3473, and C_L amplitude 0.3327. The same mesh and physical conditions are used here for the BDF1/BDF2 time-step sweep.

16.4 Simulation Conditions

Item Setting
Re_D 100
Mach 0.05
Reference state ISA sea level, U=17.0387 m/s (as in the explicit-RK4 case)
Total simulated time 12 s per run
Readout window t=6→12 s (16 cycles), after discarding the transient

16.5 Mesh

Item Setting
Geometry / mesh Same as the explicit-RK4 case — 52,428 cells, 1.67% blockage

16.6 Boundary and Initial Conditions

The cylinder boundary conditions follow the explicit-RK4 case.

Item Setting
Initial condition Uniform free stream, cold start
Note

Why cold start rather than a restart from the explicit-RK4 case. Inheriting that case’s saturated solution would make agreement between the two partly a statement about shared initial conditions rather than independent convergence to the same limit cycle. Each BDF1/BDF2 run develops its own vortex street from a uniform free stream instead.

16.7 Numerical Methods

Item Setting
Convective scheme SLAU2, second order, weighted least-squares gradient, no limiter
Time integration Implicit dual-time BDF1 / BDF2; pseudo-time CFL=30 (constant, does not affect the physical step)
Δt sweep T/32, T/64, T/128, T/256, T/512 (T=0.353341 s, from the explicit-RK4 case’s measured St)

16.8 Results

16.8.1 Convergence to the reference solution

At the finest step (Δt=T/512), BDF2 reproduces the explicit-RK4 result to four significant figures:

Metric BDF2 (Δt=T/512) Explicit-RK4 reference Deviation
St 0.1661 0.1661 0.00%
Mean C_D 1.3473 1.3473 0.00%
C_L amplitude 0.3326 0.3327 −0.03%

BDF1, being only first order, has not converged at the same step size — even at Δt=T/512 it still differs from the reference:

Metric BDF1 (Δt=T/512) Deviation from reference
St 0.1645 −0.96%
Mean C_D 1.3359 −0.85%
C_L amplitude 0.3152 −5.26%

16.8.2 Order of accuracy

Temporal order of accuracy — |error| vs Δt on a log-log plot, fitted slope is the observed order
Metric BDF1 observed order BDF2 observed order Theoretical
St 0.94 2.08 1 / 2
Mean C_D 0.92 1.38 1 / 2
C_L amplitude 0.91 1.71 1 / 2

BDF1 matches its theoretical first order on every metric. BDF2 matches its theoretical second order on St (2.08) but comes out lower for mean C_D and C_L amplitude (1.38, 1.71).

Note

The lower BDF2 orders on C_D/C_L are not a sign of reduced accuracy. By Δt=T/128, BDF2’s error on these two metrics has already fallen to the level of the explicit-RK4 case’s own cycle-to-cycle variability (±0.0067 on mean C_D), so the two finer points in the sweep are fitting the reference’s statistical noise floor rather than remaining truncation error — the log-log slope flattens once the signal being fit is noise. St does not show this problem because its cycle-to-cycle variability is two orders of magnitude smaller.

Warning

One data point does not fit the pattern. BDF1 at the coarsest step (Δt=T/32, 32 steps per shedding period) fails a physical self-consistency check: the drag signal should oscillate at exactly twice the lift-signal frequency, and every other configuration in this sweep reproduces that to within 0.3% — but this one point gives a ratio of 0.026, and its cycle-to-cycle variability is an order of magnitude higher than any other point in the sweep. Excluding it changes the fitted BDF1 order by less than 0.01 on every metric, so it does not affect the order-of-accuracy conclusion. It is, however, a practical warning: on this mesh, 32 steps per period is too coarse for BDF1 to produce a trustworthy force signal, while BDF2 at the same step count is unaffected.

16.8.3 FP64 versus FP32

In the integrated force coefficients the two are indistinguishable: mean C_D is 1.3604 (FP64) versus 1.3603 (FP32), a difference below 0.01%, and both give a C_L amplitude of 0.3354 and the same St.

Important

“Indistinguishable in C_D” is not “indistinguishable”. The same pair of runs differs clearly in two other places:

Metric FP64 FP32
Final flow residual (log₁₀) −8.98 −7.09
Drag spectrum noise floor decays cleanly flat floor near 1e-5
Wall-clock time 5.16 h 2.75 h

FP32 stalls about 1.9 orders higher in residual — it reaches its precision floor and cannot converge further. The reported 1.9× speed-up therefore includes doing fewer useful iterations, not only faster arithmetic.

Practical reading: FP32 is adequate when the quantity of interest is an integrated force averaged over a saturated limit cycle. It is not adequate when you need deep residual convergence, or when you need weak spectral content such as higher harmonics, which the FP32 noise floor buries.

The SLAU2 FP64 result agrees with the explicit-RK4 case’s reference within 1% for every metric (St +0.66%, mean C_D +0.97%, C_L amplitude +0.81%). Note that this pair also differs in solver version and CPU/GPU execution (see Run Provenance), so this is a useful consistency check rather than a controlled comparison.

BDF1 lift-coefficient time history — amplitude saturates below the literature value

BDF2 lift-coefficient time history — amplitude agrees closely with the literature value

16.9 Conclusions and Applicable Scope

  • Both BDF1 and BDF2 implicit dual-time integration converge to the explicit-RK4 reference as Δt is refined, at close to their theoretical order of accuracy — first order for BDF1, second order for BDF2 (cleanly confirmed via St; the C_D/C_L-amplitude orders are limited by the reference’s own statistical noise floor rather than by the scheme).
  • BDF2 at Δt=T/512 already matches the RK4 reference to four significant figures. BDF1 needs a substantially finer step for comparable accuracy — even at Δt=T/512 it is still 5.3% off in C_L amplitude.
  • For BDF1, keep at least on the order of 64 steps per period of the unsteady phenomenon of interest; 32 steps per period was shown empirically to corrupt the force signal on this mesh. BDF2 remained clean at 32 steps per period. These step counts are specific to this Reynolds number and mesh and should be re-verified under other conditions.
  • FP32 matches FP64 in integrated force coefficients (below 0.01%) but converges about 1.9 orders less deeply and raises the spectral noise floor. Choose FP32 for integrated loads on a saturated limit cycle, not for deep convergence or weak spectral content.
  • These results were obtained for one Reynolds number (Re_D=100), one mesh, and one convective scheme (SLAU2). They cannot be generalized to other Reynolds numbers, mesh resolutions, or flux schemes.

16.10 References

  • Williamson, C. H. K., “Defining a Universal and Continuous Strouhal–Reynolds Number Relationship for the Laminar Vortex Shedding of a Circular Cylinder,” Physics of Fluids, 31(10), 1988, pp. 2742–2744.
  • Henderson, R. D., “Details of the Drag Curve near the Onset of Vortex Shedding,” Physics of Fluids, 7(9), 1995, pp. 2102–2104.