Ensim5
Ensim5 explores compiler friendly SIMD expressions and cache locality for isentropic flow and kinematic equations.
One can model a piston flow lane from plenum intake to exhaust as nine chambers
(atmosphere, plenum, throttle, runner, piston, runner, exhaust, exhaust, atmosphere) and with eight nozzles,
mapping parallel computation perfectly to floating point AVX256 or doubling point AVX512.
The latter is required for numerical stability, but sharing at least the nozzle velocity equation with nozzle count N = 8:
/*
* ______________
* /
* / y * Pt / pt
* / -------------
* / y - 1 2
* u = M _ / 1 + ----- M
* \/ 2
*
* Pt
* pt = -------
* Rs * Tt
*
*/
fn void calc_nozzle_velocities()
{
for(size_t i = 0; i < N; i++)
{
const real Rs = g_specific_gas_constant_j_per_kg_k;
const real Tt = chamber_total_temperature_k[i];
const real M = nozzle_mach[i];
const real X = g_gamma * Rs * Tt;
const real Y = 0.5_r * (g_gamma - 1.0_r) * M * M;
const real u = M * sqrt(X / (1.0_r + Y));
const real A = chamber_nozzle_real_flow_area_m2[i];
const real mute = A == 0.0_r ? 0.0_r : 1.0_r;
nozzle_velocity_m_per_s[i] = u * mute;
}
}
The floating point (four bytes per, AVX256) output yields:
<ensim::flow<9ul, 4ul>::calc_nozzle_velocities()>:
vmovups 0x1d4(%rdi),%ymm0
vbroadcastss 0x0(%rip),%ymm1
vmulps 0x168(%rdi),%ymm1,%ymm1
vmulps %ymm0,%ymm0,%ymm2
vbroadcastss 0x0(%rip),%ymm3
vbroadcastss 0x0(%rip),%ymm4
vfmadd231ps %ymm3,%ymm2,%ymm4
vrcpps %ymm4,%ymm2
vmulps %ymm2,%ymm1,%ymm3
vfmsub213ps %ymm1,%ymm3,%ymm4
vfnmadd213ps %ymm3,%ymm2,%ymm4
vsqrtps %ymm4,%ymm1
vmulps %ymm0,%ymm1,%ymm0
vxorps %xmm1,%xmm1,%xmm1
vcmpneqps 0x6c(%rdi),%ymm1,%ymm1
vandps %ymm0,%ymm1,%ymm0
vmovups %ymm0,0x1f8(%rdi)
vzeroupper
ret
This allows for a real time engine to compute an entire piston’s chamber gas flow speed from intake to exhaust with a couple cache line reads and ~fifteen instructions. Multiply this instruction count by eight, and a straight eight engine can compute its flow field in roughly one hundred instructions.
Piston inertia torque can also be calculated with some clever sinusoid identity swaps and precomputed sin and cos thetas with this straight eight configuration:
/*
* Hailemariam Nigus. Kinematics and Load Formulation of Engine Crank Mechanism. Mechanics, Materials/
* Science & Engineering Journal, 2015, ⟨10.13140/RG.2.1.3257.1928⟩. ⟨hal-01305936⟩
*
* 2 r 1 3r
* Ti = I * w * [ ---- sin(t) - --- sin(2t) - ---- * sin(3t) ]
* 4l 2 4l
*
* These identities free up the SIMD lanes:
*
* sin(2t) = 2 sin(t) * 1 cos(t)
*
* 3
* sin(3t) = 3 sin(t) − 4 sin (t)
*
*/
fn void calc_inertia_torques()
{
for(size_t i = 0; i < W; i++)
{
const real r = crank_throw_length_m[i];
const real l = connecting_rod_length_m[i];
const real I = moment_of_inertia_kg_m2[i];
const real w = crankshaft_angular_velocity_r_per_s;
const real rl = r / l;
const real s = sint[i];
const real c = cost[i];
const real X = 0.25_r * rl * s;
const real Y = s * c;
const real Z = 0.75_r * rl * (3.0_r * s - 4.0_r * s * s * s);
inertia_torque_n_m[i] = I * w * w * (X - Y - Z);
}
}
See similar SIMD performance and aesthetics:
<ensim::inline_pistons<4ul>::calc_inertia_torques()>:
vmovss 0x180(%rdi),%xmm0
vmulss %xmm0,%xmm0,%xmm0
vmovups 0x10(%rdi),%xmm1
vmovups 0x20(%rdi),%xmm2
vmovups 0x90(%rdi),%xmm3
vrcpps %xmm2,%xmm4
vmulps %xmm4,%xmm1,%xmm5
vfmsub213ps %xmm1,%xmm5,%xmm2
vfnmadd213ps %xmm5,%xmm4,%xmm2
vbroadcastss 0x0(%rip),%xmm1
vmulps %xmm3,%xmm3,%xmm4
vbroadcastss 0x0(%rip),%xmm5
vbroadcastss 0x0(%rip),%xmm6
vfmadd231ps %xmm5,%xmm4,%xmm6
vbroadcastss 0x0(%rip),%xmm4
vmulps %xmm3,%xmm4,%xmm4
vmulps %xmm6,%xmm4,%xmm4
vmulps %xmm4,%xmm2,%xmm4
vfnmadd213ps 0xa0(%rdi),%xmm2,%xmm1
vfnmadd213ps %xmm4,%xmm3,%xmm1
vbroadcastss %xmm0,%xmm0
vmulps 0x110(%rdi),%xmm0,%xmm0
vmulps %xmm1,%xmm0,%xmm0
vmovups %xmm0,0x130(%rdi)
ret
We are not limited to the above equations, however, so see the implementation of ensim.cc. Compressible ideal gas flow is modelled in full, as are piston kinematics.
perf-stat
Powering our flow velocity model with our piston kinematic model above, we output to a branchless and cache efficient one dimensional Lax-Friedrichs-Reimann CFD pipe for audio generation, generating 48000 audio samples in 0.2 seconds:
Intel(R) Core(TM) i7-8665U CPU @ 1.90GHz
0 context-switches:u # 0.0 cs_per_second
0 cpu-migrations:u # 0.0 migrations_per_second
2,376,758 L1-dcache-load-misses:u # 0.7 % l1d_miss_rate
83,108 branch-misses:u # 0.1 % branch_miss_rate
1,627,262,496 instructions:u # 2.2 insn_per_cycle
362,993,860 dTLB-loads:u # 0.0 % dtlb_miss_rate
While single threaded performance is excellent even on a T490, pipe audio fidelity is reduced due to incompatible one dimensional CFD boundary conditions when connected to the zero dimensional isentropic flow model. I suppose combining two interpretations of the laws-of-physics will do that. I suspect a full 1D LFR pipe model in-lieu of the hybrid zero dimensional and one dimensional CFD hybrid is required to generate the most realistic sounding exhaust pulse waves.