diff --git a/enzyme/test/Integration/ForwardMode/euler_homogeneity.c b/enzyme/test/Integration/ForwardMode/euler_homogeneity.c new file mode 100644 index 000000000000..5e74a6975123 --- /dev/null +++ b/enzyme/test/Integration/ForwardMode/euler_homogeneity.c @@ -0,0 +1,125 @@ +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - + +// Gate a full 5x5 flux Jacobian, built column by column with forward mode, +// against an exact algebraic identity of the primal -- no hardcoded derivative +// values and no finite differencing. +// +// For an ideal gas, the Euler normal flux is homogeneous by degree one in the +// conserved state, so Euler's homogeneous-function theorem gives +// +// A(U) * U == F(U), A = dF/dU, +// +// which validates all 25 entries at once. This is the identity underneath +// Steger-Warming flux-vector splitting, and it must hold to machine precision. +// +// The stiffened-gas branch is the matched counterexample. Writing p = p_ideal +// - gamma*p_inf splits the flux as F = F_ideal + dF, where dF is +// +// dF = (0, -gamma*p_inf*nx, -gamma*p_inf*ny, -gamma*p_inf*nz, -gamma*p_inf*un) +// +// and is homogeneous of degree ZERO (the momentum rows are constant in U, and +// un = (rho*u . n)/rho is a ratio of degree-one quantities). Euler's theorem +// applied degree-wise then predicts the failure exactly: +// +// A(U) * U - F(U) == -dF, +// +// so the identity breaks by a closed-form amount that Enzyme must reproduce to +// machine precision. A wrong Jacobian fails both halves. + +#include "../euler.h" +#include "../frechet.h" +#include "../test_utils.h" + +extern void __enzyme_fwddiff(void *, ...); +extern int enzyme_dup; +extern int enzyme_const; + +// A[row * NVARS + col] = dF[row]/dU[col], one seeded tangent per column. The +// blocks are 5x5, so tall-thin forward mode is the right shape here. +void flux_jacobian(const double *U, const double *n, const EulerEos *eos, + double *A) { + double dU[NVARS], F[NVARS], dF[NVARS]; + + for (int col = 0; col < NVARS; col++) { + for (int i = 0; i < NVARS; i++) + dU[i] = 0.0; + dU[col] = 1.0; + + __enzyme_fwddiff((void *)euler_flux, enzyme_dup, U, dU, enzyme_const, n, + enzyme_const, eos, enzyme_dup, F, dF); + + for (int row = 0; row < NVARS; row++) + A[row * NVARS + col] = dF[row]; + } +} + +void matvec(const double *A, const double *x, double *y) { + for (int row = 0; row < NVARS; row++) { + y[row] = 0.0; + for (int col = 0; col < NVARS; col++) + y[row] += A[row * NVARS + col] * x[col]; + } +} + +int main() { + for (int s = 0; s < EULER_NSTATES; s++) { + const double *q = euler_primitives[s]; + + for (int f = 0; f < EULER_NNORMALS; f++) { + const double *n = euler_normals[f]; + + // Ideal gas: A(U) * U == F(U). + { + EulerEos eos = {1.4, 0.0}; + double U[NVARS], F[NVARS], A[NVARS * NVARS], AU[NVARS]; + + euler_from_primitive(q[0], q[1], q[2], q[3], q[4], &eos, U); + euler_flux(U, n, &eos, F); + flux_jacobian(U, n, &eos, A); + matvec(A, U, AU); + + double err = frechet_rel_error(AU, F, NVARS); + printf("state %d normal %d ideal : rel err %g\n", s, f, err); + APPROX_EQ(err, 0.0, 1e-13); + } + + // Stiffened gas: A(U) * U - F(U) == -dF, exactly. p_inf is scaled to the + // state's own pressure so every state stays physical. + { + EulerEos eos = {1.4, 0.5 * q[4]}; + double U[NVARS], F[NVARS], A[NVARS * NVARS], AU[NVARS]; + double deviation[NVARS], expected[NVARS]; + + euler_from_primitive(q[0], q[1], q[2], q[3], q[4], &eos, U); + euler_flux(U, n, &eos, F); + flux_jacobian(U, n, &eos, A); + matvec(A, U, AU); + + double gp = eos.gamma * eos.p_inf; + double un = euler_normal_velocity(U, n); + expected[RHO] = 0.0; + expected[RHOU] = gp * n[0]; + expected[RHOV] = gp * n[1]; + expected[RHOW] = gp * n[2]; + expected[RHOE] = gp * un; + + for (int i = 0; i < NVARS; i++) + deviation[i] = AU[i] - F[i]; + + double err = frechet_rel_error(deviation, expected, NVARS); + printf("state %d normal %d stiffened : rel err %g\n", s, f, err); + APPROX_EQ(err, 0.0, 1e-13); + } + } + } + + printf("done\n"); + return 0; +} diff --git a/enzyme/test/Integration/ForwardMode/flux_consistency.c b/enzyme/test/Integration/ForwardMode/flux_consistency.c new file mode 100644 index 000000000000..fc2ad6367465 --- /dev/null +++ b/enzyme/test/Integration/ForwardMode/flux_consistency.c @@ -0,0 +1,194 @@ +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - + +// Differentiate the consistency property of an approximate Riemann solver. +// +// Every consistent numerical flux satisfies Fstar(U, U, n) == F(U, n) as an +// identity in U. Differentiating both sides gives a Frechet gate on the two +// face blocks a CFD Jacobian is assembled from: +// +// dFstar/dUL + dFstar/dUR == A(U).n at UL == UR, +// +// with A the physical flux Jacobian. It is exact, needs no finite difference, +// and holds independently of which branch the solver takes. The state space +// spans both subsonic and supersonic faces: +// +// * Subsonic faces land in the HLL star region, where the flux is a rational +// function of both states and both blocks are dense. +// * Supersonic faces land in the fully upwind branch, where Fstar == FL and +// the right block must come back EXACTLY zero -- a structural property, not +// an approximate one, so it is gated at zero tolerance. +// +// Rusanov additionally exercises a tie in fmax: at UL == UR the two wave-speed +// estimates are equal, so the dissipation coefficient sits exactly on the +// non-differentiable point. The identity survives anyway, because whatever +// one-sided value the branch picks multiplies (UR - UL) == 0; the surviving +// terms are -lambda/2 and +lambda/2, which cancel in the sum. Enzyme is free +// to choose either side of the tie and still has to pass. + +#include "../euler.h" +#include "../frechet.h" +#include "../test_utils.h" + +extern void __enzyme_fwddiff(void *, ...); +extern int enzyme_dup; +extern int enzyme_const; + +// Rusanov / local Lax-Friedrichs. +void rusanov_flux(const double *UL, const double *UR, const double *n, + const EulerEos *eos, double *F) { + double FL[NVARS], FR[NVARS]; + euler_flux(UL, n, eos, FL); + euler_flux(UR, n, eos, FR); + + double lam = fmax(fabs(euler_normal_velocity(UL, n)) + + euler_sound_speed(UL, eos), + fabs(euler_normal_velocity(UR, n)) + + euler_sound_speed(UR, eos)); + + for (int i = 0; i < NVARS; i++) + F[i] = 0.5 * (FL[i] + FR[i]) - 0.5 * lam * (UR[i] - UL[i]); +} + +// HLL with Davis wave-speed estimates. +void hll_flux(const double *UL, const double *UR, const double *n, + const EulerEos *eos, double *F) { + double FL[NVARS], FR[NVARS]; + euler_flux(UL, n, eos, FL); + euler_flux(UR, n, eos, FR); + + double unL = euler_normal_velocity(UL, n); + double unR = euler_normal_velocity(UR, n); + double cL = euler_sound_speed(UL, eos); + double cR = euler_sound_speed(UR, eos); + + double sl = fmin(unL - cL, unR - cR); + double sr = fmax(unL + cL, unR + cR); + + if (sl >= 0.0) { + for (int i = 0; i < NVARS; i++) + F[i] = FL[i]; + } else if (sr <= 0.0) { + for (int i = 0; i < NVARS; i++) + F[i] = FR[i]; + } else { + double inv = 1.0 / (sr - sl); + for (int i = 0; i < NVARS; i++) + F[i] = (sr * FL[i] - sl * FR[i] + sl * sr * (UR[i] - UL[i])) * inv; + } +} + +// The physical flux Jacobian A(U).n, one seeded tangent per column. +void flux_jacobian(const double *U, const double *n, const EulerEos *eos, + double *A) { + double dU[NVARS], F[NVARS], dF[NVARS]; + + for (int col = 0; col < NVARS; col++) { + for (int i = 0; i < NVARS; i++) + dU[i] = 0.0; + dU[col] = 1.0; + + __enzyme_fwddiff((void *)euler_flux, enzyme_dup, U, dU, enzyme_const, n, + enzyme_const, eos, enzyme_dup, F, dF); + + for (int row = 0; row < NVARS; row++) + A[row * NVARS + col] = dF[row]; + } +} + +// The two face blocks of a numerical flux. Seeding one side at a time keeps +// the left and right contributions separable, which the upwinding check needs; +// the consistency identity then uses their sum. +#define DEFINE_FACE_BLOCKS(NAME, FLUX) \ + void NAME(const double *UL, const double *UR, const double *n, \ + const EulerEos *eos, double *AL, double *AR) { \ + double dUL[NVARS], dUR[NVARS], F[NVARS], dF[NVARS]; \ + \ + for (int col = 0; col < NVARS; col++) { \ + for (int i = 0; i < NVARS; i++) { \ + dUL[i] = 0.0; \ + dUR[i] = 0.0; \ + } \ + dUL[col] = 1.0; \ + __enzyme_fwddiff((void *)FLUX, enzyme_dup, UL, dUL, enzyme_dup, UR, dUR, \ + enzyme_const, n, enzyme_const, eos, enzyme_dup, F, dF); \ + for (int row = 0; row < NVARS; row++) \ + AL[row * NVARS + col] = dF[row]; \ + } \ + \ + for (int col = 0; col < NVARS; col++) { \ + for (int i = 0; i < NVARS; i++) { \ + dUL[i] = 0.0; \ + dUR[i] = 0.0; \ + } \ + dUR[col] = 1.0; \ + __enzyme_fwddiff((void *)FLUX, enzyme_dup, UL, dUL, enzyme_dup, UR, dUR, \ + enzyme_const, n, enzyme_const, eos, enzyme_dup, F, dF); \ + for (int row = 0; row < NVARS; row++) \ + AR[row * NVARS + col] = dF[row]; \ + } \ + } + +DEFINE_FACE_BLOCKS(rusanov_face_blocks, rusanov_flux) +DEFINE_FACE_BLOCKS(hll_face_blocks, hll_flux) + +int main() { + EulerEos eos = {1.4, 0.0}; + + for (int s = 0; s < EULER_NSTATES; s++) { + const double *q = euler_primitives[s]; + + for (int f = 0; f < EULER_NNORMALS; f++) { + const double *n = euler_normals[f]; + + double U[NVARS], A[NVARS * NVARS]; + double AL[NVARS * NVARS], AR[NVARS * NVARS], sum[NVARS * NVARS]; + double Fstar[NVARS], F[NVARS]; + + euler_from_primitive(q[0], q[1], q[2], q[3], q[4], &eos, U); + euler_flux(U, n, &eos, F); + flux_jacobian(U, n, &eos, A); + + // The primal is consistent to begin with -- if this fails, the + // derivative identity below would be gating the wrong thing. + rusanov_flux(U, U, n, &eos, Fstar); + APPROX_EQ(frechet_rel_error(Fstar, F, NVARS), 0.0, 1e-14); + hll_flux(U, U, n, &eos, Fstar); + APPROX_EQ(frechet_rel_error(Fstar, F, NVARS), 0.0, 1e-14); + + rusanov_face_blocks(U, U, n, &eos, AL, AR); + for (int i = 0; i < NVARS * NVARS; i++) + sum[i] = AL[i] + AR[i]; + double err_rusanov = frechet_rel_error(sum, A, NVARS * NVARS); + + hll_face_blocks(U, U, n, &eos, AL, AR); + for (int i = 0; i < NVARS * NVARS; i++) + sum[i] = AL[i] + AR[i]; + double err_hll = frechet_rel_error(sum, A, NVARS * NVARS); + + printf("state %d normal %d: rusanov %g hll %g\n", s, f, err_rusanov, + err_hll); + APPROX_EQ(err_rusanov, 0.0, 1e-13); + APPROX_EQ(err_hll, 0.0, 1e-13); + + // On a supersonic face, HLL degenerates to pure upwinding, so the right + // block must be structurally absent rather than merely small. + double un = euler_normal_velocity(U, n); + if (un - euler_sound_speed(U, &eos) > 0.0) { + for (int i = 0; i < NVARS * NVARS; i++) + APPROX_EQ(AR[i], 0.0, 0.0); + printf("state %d normal %d: hll fully upwind, right block exact zero\n", + s, f); + } + } + } + + printf("done\n"); + return 0; +} diff --git a/enzyme/test/Integration/ForwardMode/guarded_normalize.c b/enzyme/test/Integration/ForwardMode/guarded_normalize.c new file mode 100644 index 000000000000..6eedfafe7c14 --- /dev/null +++ b/enzyme/test/Integration/ForwardMode/guarded_normalize.c @@ -0,0 +1,370 @@ +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - + +// Normalising a vector that can vanish, which is the geometric kernel under +// every rotated-hybrid Riemann solver and every surface-normal calculation. +// The primal hazard is well known and every code guards it; the derivative +// hazard is separate, survives the usual guards, and is what this file pins +// down. +// +// Away from the degeneracy the Jacobian is the closed-form projector +// +// d(v/|v|)/dv = (I - n n^T) / |v|, n = v/|v|, +// +// which is gated directly, along with three structural facts that need no +// reference values at all: the Jacobian is symmetric, it annihilates n on both +// sides, and -- for the full rotated frame below -- the orthonormality of the +// frame is preserved under differentiation. +// +// At the degeneracy the three ways of writing the guard stop agreeing, and the +// difference is invisible in the primal: +// +// * Guarding after the sqrt still EVALUATES sqrt(0), whose own derivative is +// 0/0. The branch nonetheless returns a constant, and Enzyme returns an +// exactly zero tangent rather than letting the NaN escape. +// * Guarding before the sqrt never reaches it. Same answer, and this is the +// version that is obviously correct by inspection. +// * Flooring the length branchlessly with fmax(|v|, eps) is the version a +// vectorising author reaches for, and it is equally NaN-free -- but its +// derivative at the origin is I/eps, which is 1e12 here. That is the right +// derivative of what was actually written, and it is a landmine: it is +// finite, so nothing traps, and it lands in an assembled Jacobian as a row +// twelve orders of magnitude out of scale. +// +// The last point is the reason to gate this rather than just document it. All +// three idioms look interchangeable, all three keep the primal finite, and only +// two of them keep the Jacobian usable. + +#include "../frechet.h" +#include "../test_utils.h" + +extern double __enzyme_fwddiff(void *, ...); +extern int enzyme_dup; +extern int enzyme_const; + +#define GUARD_EPS 1.0e-12 + +// Guard AFTER the sqrt: sqrt(0) is evaluated at the origin. +void normalize_guard_after(const double *v, double *n) { + double len = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + if (len <= GUARD_EPS) { + n[0] = 1.0; + n[1] = 0.0; + n[2] = 0.0; + return; + } + double inv = 1.0 / len; + n[0] = v[0] * inv; + n[1] = v[1] * inv; + n[2] = v[2] * inv; +} + +// Guard BEFORE the sqrt: the sqrt is never reached at the origin. +void normalize_guard_before(const double *v, double *n) { + double sq = v[0] * v[0] + v[1] * v[1] + v[2] * v[2]; + if (sq <= GUARD_EPS * GUARD_EPS) { + n[0] = 1.0; + n[1] = 0.0; + n[2] = 0.0; + return; + } + double inv = 1.0 / sqrt(sq); + n[0] = v[0] * inv; + n[1] = v[1] * inv; + n[2] = v[2] * inv; +} + +// Branchless floor: no branch to mispredict, no NaN, and a derivative that +// blows up as 1/eps at the origin instead of vanishing. +void normalize_guard_floor(const double *v, double *n) { + double len = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + double inv = 1.0 / fmax(len, GUARD_EPS); + n[0] = v[0] * inv; + n[1] = v[1] * inv; + n[2] = v[2] * inv; +} + +#define DEFINE_NORMALIZE_JACOBIAN(NAME, FN) \ + void NAME(const double *v, double *J) { \ + double dv[3], n[3], dn[3]; \ + for (int col = 0; col < 3; col++) { \ + for (int i = 0; i < 3; i++) \ + dv[i] = 0.0; \ + dv[col] = 1.0; \ + __enzyme_fwddiff((void *)FN, enzyme_dup, v, dv, enzyme_dup, n, dn); \ + for (int row = 0; row < 3; row++) \ + J[row * 3 + col] = dn[row]; \ + } \ + } + +DEFINE_NORMALIZE_JACOBIAN(jacobian_guard_after, normalize_guard_after) +DEFINE_NORMALIZE_JACOBIAN(jacobian_guard_before, normalize_guard_before) +DEFINE_NORMALIZE_JACOBIAN(jacobian_guard_floor, normalize_guard_floor) + +// (I - n n^T) / |v| +void projector(const double *v, double *J) { + double len = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + double n[3] = {v[0] / len, v[1] / len, v[2] / len}; + for (int row = 0; row < 3; row++) + for (int col = 0; col < 3; col++) + J[row * 3 + col] = ((row == col ? 1.0 : 0.0) - n[row] * n[col]) / len; +} + +int all_finite(const double *a, int n) { + for (int i = 0; i < n; i++) + if (!(a[i] == a[i]) || a[i] > 1e300 || a[i] < -1e300) + return 0; + return 1; +} + +// The rotated frame of Nishikawa & Kitamura: the principal direction follows +// the velocity jump, the transverse direction is what is left of the face +// normal after projecting that out. Two nested guarded normalisations and a +// sign flip, which is a fair sample of what real geometry code looks like. +// +// x = (uL, uR); out = (n1, n2, alpha1, alpha2). +void rotated_frame(const double *x, const double *nf, double *out) { + double dq[3] = {x[3] - x[0], x[4] - x[1], x[5] - x[2]}; + double len = sqrt(dq[0] * dq[0] + dq[1] * dq[1] + dq[2] * dq[2]); + + if (len <= GUARD_EPS) { + // Negligible velocity jump: fall back to the grid-aligned frame. + for (int i = 0; i < 3; i++) { + out[i] = nf[i]; + out[3 + i] = nf[i]; + } + out[6] = 0.0; + out[7] = 1.0; + return; + } + + double n1[3] = {dq[0] / len, dq[1] / len, dq[2] / len}; + double a1 = nf[0] * n1[0] + nf[1] * n1[1] + nf[2] * n1[2]; + + // Reorient so the principal weight is non-negative. + if (a1 < 0.0) { + for (int i = 0; i < 3; i++) + n1[i] = -n1[i]; + a1 = -a1; + } + + double t[3] = {nf[0] - a1 * n1[0], nf[1] - a1 * n1[1], nf[2] - a1 * n1[2]}; + double tlen = sqrt(t[0] * t[0] + t[1] * t[1] + t[2] * t[2]); + + if (tlen <= GUARD_EPS) { + // Velocity jump parallel to the face normal: the transverse direction is + // undefined and carries no weight. + for (int i = 0; i < 3; i++) { + out[i] = n1[i]; + out[3 + i] = nf[i]; + } + out[6] = 1.0; + out[7] = 0.0; + return; + } + + for (int i = 0; i < 3; i++) { + out[i] = n1[i]; + out[3 + i] = t[i] / tlen; + } + out[6] = a1; + out[7] = nf[0] * out[3] + nf[1] * out[4] + nf[2] * out[5]; +} + +// Three identities the frame satisfies for every input, hence three gradients +// that must vanish identically. +double frame_n1_normsq(const double *x, const double *nf) { + double o[8]; + rotated_frame(x, nf, o); + return o[0] * o[0] + o[1] * o[1] + o[2] * o[2]; +} + +double frame_n1_dot_n2(const double *x, const double *nf) { + double o[8]; + rotated_frame(x, nf, o); + return o[0] * o[3] + o[1] * o[4] + o[2] * o[5]; +} + +double frame_alpha_normsq(const double *x, const double *nf) { + double o[8]; + rotated_frame(x, nf, o); + return o[6] * o[6] + o[7] * o[7]; +} + +#define DEFINE_IDENTITY_GRADIENT(NAME, FN) \ + void NAME(const double *x, const double *nf, double *g) { \ + double dx[6]; \ + for (int col = 0; col < 6; col++) { \ + for (int i = 0; i < 6; i++) \ + dx[i] = 0.0; \ + dx[col] = 1.0; \ + g[col] = __enzyme_fwddiff((void *)FN, enzyme_dup, x, dx, enzyme_const, \ + nf); \ + } \ + } + +DEFINE_IDENTITY_GRADIENT(grad_n1_normsq, frame_n1_normsq) +DEFINE_IDENTITY_GRADIENT(grad_n1_dot_n2, frame_n1_dot_n2) +DEFINE_IDENTITY_GRADIENT(grad_alpha_normsq, frame_alpha_normsq) + +int main() { + double J[9], K[9], want[9]; + + // Away from the degeneracy every idiom gives the projector. + { + const double vectors[3][3] = { + {3.0, -4.0, 12.0}, {1.0, 0.0, 0.0}, {-0.02, 0.005, 0.031}}; + + for (int i = 0; i < 3; i++) { + const double *v = vectors[i]; + projector(v, want); + + jacobian_guard_after(v, J); + APPROX_EQ(frechet_rel_error(J, want, 9), 0.0, 1e-14); + jacobian_guard_before(v, J); + APPROX_EQ(frechet_rel_error(J, want, 9), 0.0, 1e-14); + jacobian_guard_floor(v, J); + APPROX_EQ(frechet_rel_error(J, want, 9), 0.0, 1e-14); + + // Structural: symmetric, and it annihilates n from both sides. + double len = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + double n[3] = {v[0] / len, v[1] / len, v[2] / len}; + double scale = 1.0 / len; + for (int r = 0; r < 3; r++) { + double row = 0.0, col = 0.0; + for (int c = 0; c < 3; c++) { + APPROX_EQ(J[r * 3 + c], J[c * 3 + r], 1e-15 * scale); + row += J[r * 3 + c] * n[c]; + col += J[c * 3 + r] * n[c]; + } + APPROX_EQ(row, 0.0, 1e-14 * scale); + APPROX_EQ(col, 0.0, 1e-14 * scale); + } + printf("normalize |v|=%.4g: projector, symmetry and null space all hold\n", + len); + } + } + + // At the origin, the two branch idioms return an exactly zero tangent -- the + // sqrt's 0/0 does not leak out of the untaken path. + { + double origin[3] = {0.0, 0.0, 0.0}; + + jacobian_guard_after(origin, J); + printf("origin, guard after sqrt : J[0][0] = %g\n", J[0]); + if (!all_finite(J, 9)) { + fprintf(stderr, "guard-after produced a non-finite tangent at the " + "origin\n"); + abort(); + } + for (int i = 0; i < 9; i++) + APPROX_EQ(J[i], 0.0, 0.0); + + jacobian_guard_before(origin, K); + printf("origin, guard before sqrt: J[0][0] = %g\n", K[0]); + if (!all_finite(K, 9)) { + fprintf(stderr, "guard-before produced a non-finite tangent at the " + "origin\n"); + abort(); + } + for (int i = 0; i < 9; i++) + APPROX_EQ(K[i], 0.0, 0.0); + + // The branchless floor is finite too, and that is exactly the problem: it + // is I/eps, a perfectly valid derivative of what was written and a row of + // an assembled Jacobian that is 1e12 out of scale. + jacobian_guard_floor(origin, K); + printf("origin, branchless floor : J[0][0] = %g (1/eps = %g)\n", K[0], + 1.0 / GUARD_EPS); + if (!all_finite(K, 9)) { + fprintf(stderr, "branchless floor produced a non-finite tangent\n"); + abort(); + } + for (int row = 0; row < 3; row++) + for (int col = 0; col < 3; col++) + APPROX_EQ(K[row * 3 + col], + (row == col ? 1.0 / GUARD_EPS : 0.0), 1e-3); + } + + // Just inside the guard the branch idioms are still zero; just outside, the + // stiffness the guard was hiding is fully present. The guard relocates the + // blow-up, it does not remove it. + { + double inside[3] = {0.5 * GUARD_EPS, 0.0, 0.0}; + double outside[3] = {2.0 * GUARD_EPS, 0.0, 0.0}; + + jacobian_guard_after(inside, J); + for (int i = 0; i < 9; i++) + APPROX_EQ(J[i], 0.0, 0.0); + + jacobian_guard_after(outside, J); + projector(outside, want); + APPROX_EQ(frechet_rel_error(J, want, 9), 0.0, 1e-14); + printf("just outside the guard: |J| ~ %.3g\n", fabs(J[4])); + if (!(fabs(J[4]) > 1e11)) { + fprintf(stderr, "expected the projector to be stiff just outside the " + "guard\n"); + abort(); + } + } + + // The rotated frame: orthonormality is an identity in the inputs, so its + // gradient has to vanish -- including on both degenerate branches, where the + // frame is assembled from entirely different expressions. + { + const double nf[3] = {0.6, -0.8, 0.0}; + const double cases[4][6] = { + {10.0, 3.0, -2.0, 40.0, -9.0, 5.0}, // generic + {10.0, 3.0, -2.0, 10.0, 3.0, -2.0}, // no velocity jump + {0.0, 0.0, 0.0, 0.6, -0.8, 0.0}, // jump parallel to the normal + {1.0, 1.0, 1.0, 1.0 - 0.6, 1.0 + 0.8, 1.0}, // antiparallel, flips sign + }; + const char *labels[4] = {"generic", "no jump", "parallel", "antiparallel"}; + + for (int c = 0; c < 4; c++) { + const double *x = cases[c]; + double o[8], g[6]; + + rotated_frame(x, nf, o); + APPROX_EQ(frame_n1_normsq(x, nf), 1.0, 1e-14); + APPROX_EQ(frame_alpha_normsq(x, nf), 1.0, 1e-14); + + grad_n1_normsq(x, nf, g); + if (!all_finite(g, 6)) { + fprintf(stderr, "non-finite d|n1|^2 on the %s case\n", labels[c]); + abort(); + } + for (int i = 0; i < 6; i++) + APPROX_EQ(g[i], 0.0, 1e-13); + + grad_n1_dot_n2(x, nf, g); + if (!all_finite(g, 6)) { + fprintf(stderr, "non-finite d(n1.n2) on the %s case\n", labels[c]); + abort(); + } + for (int i = 0; i < 6; i++) + APPROX_EQ(g[i], 0.0, 1e-13); + + grad_alpha_normsq(x, nf, g); + if (!all_finite(g, 6)) { + fprintf(stderr, "non-finite d(alpha1^2+alpha2^2) on the %s case\n", + labels[c]); + abort(); + } + for (int i = 0; i < 6; i++) + APPROX_EQ(g[i], 0.0, 1e-13); + + printf("rotated frame %-13s: orthonormality gradients all vanish\n", + labels[c]); + } + } + + printf("done\n"); + return 0; +} diff --git a/enzyme/test/Integration/ForwardMode/newton_inversion.c b/enzyme/test/Integration/ForwardMode/newton_inversion.c new file mode 100644 index 000000000000..633c25f87f90 --- /dev/null +++ b/enzyme/test/Integration/ForwardMode/newton_inversion.c @@ -0,0 +1,182 @@ +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - + +// Forward mode through a Newton solve whose trip count is data dependent. +// +// The loop runs up to 50 times, breaks early once the step is small, can break +// instead on a collapsed heat capacity, and clamps its iterate to a physical +// range every pass -- so the differentiated control flow is decided at runtime. +// The answer is fixed by the implicit function theorem and never mentions the +// iterates: +// +// dT/de = 1 / cv(T), +// dT/dY_s = -e_s(T) / cv(T), +// dT/dT_guess == 0. +// +// The last one is the sharpest gate in the file. A converged root cannot +// depend on where the iteration started, so the tangent with respect to the +// initial guess has to vanish identically -- not to the solver tolerance, to +// machine precision. It fails loudly for any AD that hands back the derivative +// of the final iterate in place of the derivative of the solution. +// +// The rest of the test measures what decides that. A Newton step applied at an +// already-converged iterate discards whatever tangent the loop accumulated (see +// thermo.h), so everything hinges on whether the last step lands there: +// +// 1. test-after-update -- quadratic convergence means the loop always runs +// one full step past the point a 1e-6 test could trip, so it performs the +// collapse itself and the tangent comes back at machine precision. +// 2. test-before-update -- that final step is skipped. Both the primal and +// the tangent are left at the step tolerance, the tangent roughly an order +// of magnitude worse, so the derivative is the least accurate thing a +// converged-looking solve returns. Gated as an inequality, so the file +// fails if it stops measuring that. +// 3. test-before-update + polish -- machine precision again, independent of +// where the test was placed. +// +// Finally, on the temperature clamps the derivative is exactly zero: fmax/fmin +// have selected a constant branch and no tangent flows through it. + +#include "../test_utils.h" +#include "../thermo.h" + +extern double __enzyme_fwddiff(void *, ...); +extern int enzyme_dup; +extern int enzyme_const; + +double dT_denergy(double e_target, const double *Y, double T_guess, + int test_before_update, int polish) { + return __enzyme_fwddiff((void *)thermo_temperature, e_target, 1.0, + enzyme_const, Y, enzyme_const, T_guess, enzyme_const, + test_before_update, enzyme_const, polish); +} + +double dT_dmassfraction(double e_target, const double *Y, double T_guess, + int test_before_update, int polish, int s) { + double dY[NSPECIES]; + for (int i = 0; i < NSPECIES; i++) + dY[i] = 0.0; + dY[s] = 1.0; + + return __enzyme_fwddiff((void *)thermo_temperature, e_target, 0.0, enzyme_dup, + Y, dY, enzyme_const, T_guess, enzyme_const, + test_before_update, enzyme_const, polish); +} + +double dT_dguess(double e_target, const double *Y, double T_guess, + int test_before_update, int polish) { + return __enzyme_fwddiff((void *)thermo_temperature, e_target, 0.0, + enzyme_const, Y, T_guess, 1.0, enzyme_const, + test_before_update, enzyme_const, polish); +} + +// Relative deviation of `got` from a nonzero reference. +double rel(double got, double want) { return fabs(got - want) / fabs(want); } + +int main() { + const double air[NSPECIES] = {0.767, 0.233}; + const double even[NSPECIES] = {0.5, 0.5}; + const double *mixtures[2] = {air, even}; + + const double temperatures[] = {200.0, 300.0, 1000.0, 3000.0, 8000.0, 15000.0}; + const int ntemp = sizeof(temperatures) / sizeof(temperatures[0]); + + // Deliberately poor initial guess so the loop actually iterates. + const double T_guess = 500.0; + + double worst_stale = 0.0; + + for (int m = 0; m < 2; m++) { + const double *Y = mixtures[m]; + + for (int t = 0; t < ntemp; t++) { + double T_ref = temperatures[t]; + double e_target = thermo_energy(T_ref, Y); + + double T = + thermo_temperature(e_target, Y, T_guess, THERMO_TEST_AFTER_UPDATE, 0); + APPROX_EQ(rel(T, T_ref), 0.0, 1e-12); + + double cv = thermo_cv(T, Y); + double want = 1.0 / cv; + + double after = + dT_denergy(e_target, Y, T_guess, THERMO_TEST_AFTER_UPDATE, 0); + double before = + dT_denergy(e_target, Y, T_guess, THERMO_TEST_BEFORE_UPDATE, 0); + double repaired = + dT_denergy(e_target, Y, T_guess, THERMO_TEST_BEFORE_UPDATE, 1); + + double err_after = rel(after, want); + double err_before = rel(before, want); + double err_repaired = rel(repaired, want); + + // The primal the stale tangent sits next to, for comparison. + double T_stale = thermo_temperature(e_target, Y, T_guess, + THERMO_TEST_BEFORE_UPDATE, 0); + double err_primal = rel(T_stale, T_ref); + + printf("mix %d T %7.1f: dT/de after %.2e before %.2e repaired %.2e | " + "primal %.2e\n", + m, T_ref, err_after, err_before, err_repaired, err_primal); + + APPROX_EQ(err_after, 0.0, 1e-13); + APPROX_EQ(err_repaired, 0.0, 1e-13); + APPROX_EQ(err_primal, 0.0, 1e-7); + + if (err_before > worst_stale) + worst_stale = err_before; + + // dT/dY_s against -e_s(T)/cv(T), through the same iterative kernel. + for (int s = 0; s < NSPECIES; s++) { + double got = dT_dmassfraction(e_target, Y, T_guess, + THERMO_TEST_AFTER_UPDATE, 0, s); + APPROX_EQ(rel(got, -thermo_energy_species(s, T) / cv), 0.0, 1e-11); + } + + // The converged root cannot depend on where the iteration started. + double dguess = + dT_dguess(e_target, Y, T_guess, THERMO_TEST_AFTER_UPDATE, 0); + APPROX_EQ(dguess, 0.0, 1e-13); + } + } + + // The stale tangent has to be visibly stale somewhere in the sweep, or this + // file has stopped testing what it claims to. It is worst where cv varies + // fastest, which is what the uncollapsed term is proportional to. + printf("worst stale tangent: %.3e\n", worst_stale); + if (worst_stale < 1.0e-9) { + fprintf(stderr, + "no measurably stale tangent anywhere in the sweep (worst %.3e) -- " + "the loop is no longer stopping on its step tolerance\n", + worst_stale); + abort(); + } + + // On the clamps the derivative is structurally zero, not merely small. + { + double T = + thermo_temperature(-1.0e5, air, T_guess, THERMO_TEST_AFTER_UPDATE, 1); + double dT = dT_denergy(-1.0e5, air, T_guess, THERMO_TEST_AFTER_UPDATE, 1); + printf("clamped low : T %g dT/de %g\n", T, dT); + APPROX_EQ(T, THERMO_T_MIN, 0.0); + APPROX_EQ(dT, 0.0, 0.0); + } + { + double T = + thermo_temperature(1.0e12, air, T_guess, THERMO_TEST_AFTER_UPDATE, 1); + double dT = dT_denergy(1.0e12, air, T_guess, THERMO_TEST_AFTER_UPDATE, 1); + printf("clamped high: T %g dT/de %g\n", T, dT); + APPROX_EQ(T, THERMO_T_MAX, 0.0); + APPROX_EQ(dT, 0.0, 0.0); + } + + printf("done\n"); + return 0; +} diff --git a/enzyme/test/Integration/ForwardMode/substencil_weights.c b/enzyme/test/Integration/ForwardMode/substencil_weights.c new file mode 100644 index 000000000000..6f21e0a10239 --- /dev/null +++ b/enzyme/test/Integration/ForwardMode/substencil_weights.c @@ -0,0 +1,410 @@ +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - + +// Shock-capturing reconstruction: five stencil values in, one face value out, +// through nonlinear weights that are deliberately non-smooth in the data. +// WENO5-JS, WENO5-Z and TENO5 differ only in how they weight the same three +// fifth-order candidates, and between them they cover three separate hazards +// that a plain "AD vs finite difference" check would either miss or misreport. +// +// 1. WENO5-JS is smooth: its 1/(beta+eps)^2 denominators are held off zero by +// eps = 1e-6, so AD and a central difference agree, and this is the control +// case that says the harness works. +// +// 2. WENO5-Z has a kink. Its tau5 = |beta0 - beta2| is exactly zero on any +// symmetric stencil -- beta0 and beta2 are equal there as an algebraic +// identity, not by luck -- so a symmetric stencil sits exactly on the +// non-differentiable point. The one-sided slopes really do differ (by +// around 20% of the value here), AD returns one of them, and a central +// difference returns their average, which is neither. This is the case +// where FD is wrong and AD is right, and the test pins down all three +// numbers rather than just asserting a disagreement. +// +// 3. TENO5 is discontinuous. Its sharp cutoff sets delta_k = 0 the moment a +// candidate's normalised strength falls below C_T = 1e-5, so the final +// weights delta_k * d_k are piecewise CONSTANT in the data. Two things +// follow, and both are gated: inside a branch the reconstruction is exactly +// linear in the stencil, so the gradient is exactly the frozen-weight +// stencil to machine precision; and across a branch the primal jumps, so a +// central difference straddling the cutoff is meaningless while AD's +// within-branch answer stays exact. The crossing is located by bisection +// rather than hardcoded, so the test keeps working if the constants move. +// +// The linear-profile check at the top is the one that ties the schemes +// together: when every candidate agrees, all three collapse to the optimal +// fifth-order stencil (1/30, -13/60, 47/60, 9/20, -1/20), and the |tau5| kink +// is annihilated because the weight perturbation multiplies (f_k - q) == 0. + +#include "../frechet.h" +#include "../test_utils.h" + +extern double __enzyme_fwddiff(void *, ...); +extern int enzyme_dup; +extern int enzyme_const; + +#define NSTENCIL 5 + +// Sharp cutoff threshold, sharpness exponent and denominator floors. +#define TENO_CT 1.0e-5 +#define EPS_TENO 1.0e-40 +#define EPS_WENOZ 1.0e-40 +#define EPS_WENOJS 1.0e-6 + +// Optimal sub-stencil combination weights. +static const double d_opt[3] = {0.1, 0.6, 0.3}; + +// The nine fifth-order Lagrange coefficients, candidate-major. +static const double sub_coeff[9] = { + 1.0 / 3.0, -7.0 / 6.0, 11.0 / 6.0, // f0: v[0], v[1], v[2] + -1.0 / 6.0, 5.0 / 6.0, 1.0 / 3.0, // f1: v[1], v[2], v[3] + 1.0 / 3.0, 5.0 / 6.0, -1.0 / 6.0, // f2: v[2], v[3], v[4] +}; + +// The optimal fifth-order stencil, which every scheme has to reduce to when +// its candidates agree: sum_k d_k * sub_coeff[k]. +static const double optimal[NSTENCIL] = {1.0 / 30.0, -13.0 / 60.0, 47.0 / 60.0, + 9.0 / 20.0, -1.0 / 20.0}; + +// Jiang-Shu smoothness indicators. +void smoothness(const double *v, double *beta) { + double a0 = v[0] - 2.0 * v[1] + v[2]; + double b0 = v[0] - 4.0 * v[1] + 3.0 * v[2]; + double a1 = v[1] - 2.0 * v[2] + v[3]; + double b1 = v[1] - v[3]; + double a2 = v[2] - 2.0 * v[3] + v[4]; + double b2 = 3.0 * v[2] - 4.0 * v[3] + v[4]; + + beta[0] = (13.0 / 12.0) * a0 * a0 + 0.25 * b0 * b0; + beta[1] = (13.0 / 12.0) * a1 * a1 + 0.25 * b1 * b1; + beta[2] = (13.0 / 12.0) * a2 * a2 + 0.25 * b2 * b2; +} + +// The three fifth-order candidates at the face i+1/2. +void candidates(const double *v, double *f) { + f[0] = sub_coeff[0] * v[0] + sub_coeff[1] * v[1] + sub_coeff[2] * v[2]; + f[1] = sub_coeff[3] * v[1] + sub_coeff[4] * v[2] + sub_coeff[5] * v[3]; + f[2] = sub_coeff[6] * v[2] + sub_coeff[7] * v[3] + sub_coeff[8] * v[4]; +} + +double normalise(const double *w, const double *f) { + return (w[0] * f[0] + w[1] * f[1] + w[2] * f[2]) / (w[0] + w[1] + w[2]); +} + +// WENO5-JS: alpha_k = d_k / (beta_k + eps)^2. +double weno5js(const double *v) { + double beta[3], f[3], w[3]; + smoothness(v, beta); + candidates(v, f); + for (int k = 0; k < 3; k++) { + double den = beta[k] + EPS_WENOJS; + w[k] = d_opt[k] / (den * den); + } + return normalise(w, f); +} + +// WENO5-Z with tau5 supplied by the caller, so the two smooth branches of +// |beta0 - beta2| can be differentiated separately. `sign` selects which one: +// on a symmetric stencil both give the same value and their derivatives are +// precisely the two one-sided limits of the real scheme. +double weno5z_branch(const double *v, double sign) { + double beta[3], f[3], w[3]; + smoothness(v, beta); + candidates(v, f); + double tau5 = sign * (beta[0] - beta[2]); + for (int k = 0; k < 3; k++) + w[k] = d_opt[k] * (1.0 + tau5 / (beta[k] + EPS_WENOZ)); + return normalise(w, f); +} + +// WENO5-Z proper: tau5 = |beta0 - beta2|. +double weno5z(const double *v) { + double beta[3], f[3], w[3]; + smoothness(v, beta); + candidates(v, f); + double tau5 = fabs(beta[0] - beta[2]); + for (int k = 0; k < 3; k++) + w[k] = d_opt[k] * (1.0 + tau5 / (beta[k] + EPS_WENOZ)); + return normalise(w, f); +} + +// TENO5's cutoff flags. Kept separate from the reconstruction because they are +// exactly what the frozen-weight identity freezes. +void teno5_cutoff(const double *v, double *delta) { + double beta[3], g[3]; + smoothness(v, beta); + double tau5 = fabs(beta[0] - beta[2]); + + double sum = 0.0; + for (int k = 0; k < 3; k++) { + double x = 1.0 + tau5 / (beta[k] + EPS_TENO); + double x3 = x * x * x; + g[k] = x3 * x3; + sum += g[k]; + } + for (int k = 0; k < 3; k++) + delta[k] = (g[k] / sum < TENO_CT) ? 0.0 : 1.0; +} + +double teno5(const double *v) { + double delta[3], f[3], w[3]; + teno5_cutoff(v, delta); + candidates(v, f); + for (int k = 0; k < 3; k++) + w[k] = delta[k] * d_opt[k]; + + double sum = w[0] + w[1] + w[2]; + if (sum < EPS_TENO) + return f[1]; // every candidate cut off: fall back to the central one + return (w[0] * f[0] + w[1] * f[1] + w[2] * f[2]) / sum; +} + +#define DEFINE_GRADIENT(NAME, SCHEME) \ + void NAME(const double *v, double *g) { \ + double dv[NSTENCIL]; \ + for (int col = 0; col < NSTENCIL; col++) { \ + for (int i = 0; i < NSTENCIL; i++) \ + dv[i] = 0.0; \ + dv[col] = 1.0; \ + g[col] = __enzyme_fwddiff((void *)SCHEME, enzyme_dup, v, dv); \ + } \ + } + +DEFINE_GRADIENT(weno5js_gradient, weno5js) +DEFINE_GRADIENT(weno5z_gradient, weno5z) +DEFINE_GRADIENT(teno5_gradient, teno5) + +// The gradient of one smooth branch of WENO5-Z. +void weno5z_branch_gradient(const double *v, double sign, double *g) { + double dv[NSTENCIL]; + for (int col = 0; col < NSTENCIL; col++) { + for (int i = 0; i < NSTENCIL; i++) + dv[i] = 0.0; + dv[col] = 1.0; + g[col] = __enzyme_fwddiff((void *)weno5z_branch, enzyme_dup, v, dv, + enzyme_const, sign); + } +} + +// Central-difference gradient, for the comparisons where FD is the thing under +// examination rather than the reference. +typedef double (*Scheme)(const double *); + +void fd_gradient(Scheme q, const double *v, double h, double *g) { + double vp[NSTENCIL], vm[NSTENCIL]; + for (int col = 0; col < NSTENCIL; col++) { + for (int i = 0; i < NSTENCIL; i++) { + vp[i] = v[i]; + vm[i] = v[i]; + } + vp[col] += h; + vm[col] -= h; + g[col] = (q(vp) - q(vm)) / (2.0 * h); + } +} + +// TENO5's exact gradient inside a branch: with delta frozen the reconstruction +// is linear in v, so the gradient is just the normalised candidate stencil. +void teno5_frozen_gradient(const double *v, double *g) { + double delta[3]; + teno5_cutoff(v, delta); + + double w[3], sum = 0.0; + for (int k = 0; k < 3; k++) { + w[k] = delta[k] * d_opt[k]; + sum += w[k]; + } + + for (int j = 0; j < NSTENCIL; j++) + g[j] = 0.0; + if (sum < EPS_TENO) { + g[1] = sub_coeff[3]; + g[2] = sub_coeff[4]; + g[3] = sub_coeff[5]; + return; + } + for (int k = 0; k < 3; k++) + for (int c = 0; c < 3; c++) + g[k + c] += (w[k] / sum) * sub_coeff[3 * k + c]; +} + +// A smooth asymmetric base with a step of amplitude s across the face, which +// drives the TENO cutoff through its flips as s grows. +void step_profile(double s, double *v) { + for (int j = 0; j < NSTENCIL; j++) + v[j] = sin(0.7 * (double)j) + (j >= 3 ? s : 0.0); +} + +int main() { + double g[NSTENCIL], g2[NSTENCIL], fd[NSTENCIL]; + + // Every scheme reduces to the optimal fifth-order stencil when the + // candidates agree, and the |tau5| kink cannot bite there because the weight + // perturbation multiplies (f_k - q) == 0. + { + double linear[NSTENCIL]; + for (int j = 0; j < NSTENCIL; j++) + linear[j] = 3.0 + 2.0 * (double)j; + + weno5js_gradient(linear, g); + APPROX_EQ(frechet_rel_error(g, optimal, NSTENCIL), 0.0, 1e-14); + weno5z_gradient(linear, g); + APPROX_EQ(frechet_rel_error(g, optimal, NSTENCIL), 0.0, 1e-14); + teno5_gradient(linear, g); + APPROX_EQ(frechet_rel_error(g, optimal, NSTENCIL), 0.0, 1e-14); + printf("linear profile: all three schemes give the optimal stencil\n"); + } + + // WENO5-JS is smooth, so AD and a central difference agree. Control case. + { + double v[NSTENCIL]; + for (int j = 0; j < NSTENCIL; j++) + v[j] = sin(0.7 * (double)j); + + weno5js_gradient(v, g); + fd_gradient(weno5js, v, 1e-6, fd); + double err = frechet_rel_error(fd, g, NSTENCIL); + printf("weno5-js smooth: AD vs FD %.3e\n", err); + APPROX_EQ(err, 0.0, 1e-8); + } + + // WENO5-Z on a symmetric stencil sits exactly on the |tau5| kink. + { + double v[NSTENCIL] = {1.0, 0.0, 1.0, 0.0, 1.0}; + double beta[3]; + smoothness(v, beta); + + // beta0 == beta2 is an algebraic identity for a symmetric stencil, so the + // kink is hit exactly rather than approached. + APPROX_EQ(beta[0] - beta[2], 0.0, 0.0); + + double plus[NSTENCIL], minus[NSTENCIL], mean[NSTENCIL]; + weno5z_gradient(v, g); + weno5z_branch_gradient(v, 1.0, plus); + weno5z_branch_gradient(v, -1.0, minus); + for (int j = 0; j < NSTENCIL; j++) + mean[j] = 0.5 * (plus[j] + minus[j]); + + double jump = 0.0; + for (int j = 0; j < NSTENCIL; j++) { + double diff = fabs(plus[j] - minus[j]); + if (diff > jump) + jump = diff; + } + printf("weno5-z kink: one-sided slopes differ by %.3e\n", jump); + + // The kink is real, not a rounding artefact. + if (jump < 1e-3) { + fprintf(stderr, "one-sided slopes agree -- no kink to test\n"); + abort(); + } + + // AD lands on one of the two one-sided limits, exactly. + double to_plus = frechet_rel_error(g, plus, NSTENCIL); + double to_minus = frechet_rel_error(g, minus, NSTENCIL); + printf("weno5-z kink: AD to (+) branch %.3e, to (-) branch %.3e\n", to_plus, + to_minus); + if (to_plus > 1e-14 && to_minus > 1e-14) { + fprintf(stderr, "AD matched neither one-sided limit\n"); + abort(); + } + + // A central difference returns their average, which is neither -- this is + // the case where finite differencing is simply the wrong instrument. It + // tracks the mean to O(h) rather than the usual O(h^2), since the second + // order term is exactly what the kink destroys, so the tolerance here is + // the step size and not its square. That is still four orders below the + // distance from AD, which is the comparison the test is making. + fd_gradient(weno5z, v, 1e-6, fd); + printf("weno5-z kink: FD to mean %.3e, FD to AD %.3e\n", + frechet_rel_error(fd, mean, NSTENCIL), + frechet_rel_error(fd, g, NSTENCIL)); + APPROX_EQ(frechet_rel_error(fd, mean, NSTENCIL), 0.0, 1e-5); + if (frechet_rel_error(fd, g, NSTENCIL) < 1e-3) { + fprintf(stderr, "FD agreed with AD at the kink -- expected it not to\n"); + abort(); + } + } + + // TENO5's weights are piecewise constant, so inside a branch its gradient is + // exactly the frozen-weight stencil, discontinuity or not. + { + const double amplitudes[] = {0.0, 0.2, 0.35, 0.8, 1.5, 3.0}; + const int n = sizeof(amplitudes) / sizeof(amplitudes[0]); + + for (int i = 0; i < n; i++) { + double v[NSTENCIL], delta[3]; + step_profile(amplitudes[i], v); + teno5_cutoff(v, delta); + teno5_gradient(v, g); + teno5_frozen_gradient(v, g2); + + double err = frechet_rel_error(g, g2, NSTENCIL); + printf("teno5 s=%.2f delta=(%.0f,%.0f,%.0f): AD vs frozen weights %.3e\n", + amplitudes[i], delta[0], delta[1], delta[2], err); + APPROX_EQ(err, 0.0, 1e-14); + } + } + + // Bisect to a cutoff crossing and show what each instrument reports there. + { + double lo = 0.20, hi = 1.00; // delta_2 flips 1 -> 0 somewhere inside + double v[NSTENCIL], delta[3]; + + step_profile(lo, v); + teno5_cutoff(v, delta); + double delta2_lo = delta[2]; + step_profile(hi, v); + teno5_cutoff(v, delta); + if (delta2_lo == delta[2]) { + fprintf(stderr, "no cutoff crossing bracketed\n"); + abort(); + } + + for (int it = 0; it < 200; it++) { + double mid = 0.5 * (lo + hi); + step_profile(mid, v); + teno5_cutoff(v, delta); + if (delta[2] == delta2_lo) + lo = mid; + else + hi = mid; + } + printf("teno5 cutoff crossing bracketed to [%.17g, %.17g]\n", lo, hi); + + // The primal really is discontinuous across it. + double vlo[NSTENCIL], vhi[NSTENCIL]; + step_profile(lo, vlo); + step_profile(hi, vhi); + double jump = fabs(teno5(vhi) - teno5(vlo)); + printf("teno5 primal jump across the cutoff: %.3e\n", jump); + if (jump < 1e-3) { + fprintf(stderr, "cutoff crossing is not discontinuous -- nothing to " + "distinguish AD from FD\n"); + abort(); + } + + // AD on the low side reports that side's exact linear stencil... + teno5_gradient(vlo, g); + teno5_frozen_gradient(vlo, g2); + APPROX_EQ(frechet_rel_error(g, g2, NSTENCIL), 0.0, 1e-14); + + // ...while a central difference wide enough to straddle the crossing is + // reporting the jump divided by the step, and agrees with nothing. + fd_gradient(teno5, vlo, 0.05, fd); + double err = frechet_rel_error(fd, g, NSTENCIL); + printf("teno5 at the cutoff: straddling FD vs AD %.3e\n", err); + if (err < 0.1) { + fprintf(stderr, "straddling FD tracked AD -- expected it not to\n"); + abort(); + } + } + + printf("done\n"); + return 0; +} diff --git a/enzyme/test/Integration/ReverseMode/entropy_hessian.c b/enzyme/test/Integration/ReverseMode/entropy_hessian.c new file mode 100644 index 000000000000..4a9bb6e9f8fa --- /dev/null +++ b/enzyme/test/Integration/ReverseMode/entropy_hessian.c @@ -0,0 +1,205 @@ +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - + +// Forward-over-reverse on the mathematical entropy of the Euler equations, +// which is the one place in a compressible solver where second derivatives are +// load bearing. Entropy-stable schemes exist because +// +// eta(U) = -rho * s / (gamma - 1), s = log(p) - gamma * log(rho) +// +// is convex, and every guarantee they offer follows from that -- so the Hessian +// this test builds is not an abstraction, it is the object whose definiteness +// makes the scheme provably stable. +// +// Nested AD is the fragile case in practice, and this file gates it four ways, +// none of which needs a hardcoded derivative: +// +// 1. The reverse-mode gradient must reproduce the entropy variables +// w = (( gamma - s)/(gamma-1) - beta|u|^2, 2 beta u, 2 beta v, 2 beta w, +// -2 beta), beta = rho/2p, in closed form. This is the vector an +// entropy-stable flux is actually written in terms of. +// 2. The Hessian must be symmetric. It is a Hessian, so this is free, exact, +// and impossible to satisfy by accident across 25 independently computed +// entries. +// 3. The Hessian must be positive definite, which is the convexity that the +// scheme's stability proof rests on. Cholesky either completes or it does +// not; no tolerance is involved. +// 4. Forward-over-reverse must agree with plain forward mode applied to the +// closed-form gradient. Two different AD compositions, one of them +// nested, over the same mathematics. +// +// Together these say the second-order tape is right without ever writing down +// a second derivative by hand. + +#include "../euler.h" +#include "../frechet.h" +#include "../test_utils.h" + +extern double __enzyme_autodiff(void *, ...); +extern void __enzyme_fwddiff(void *, ...); +extern int enzyme_dup; +extern int enzyme_const; + +// Harten's entropy for the compressible Euler equations. +double mathematical_entropy(const double *U, const EulerEos *eos) { + double p = euler_pressure(U, eos); + double s = log(p) - eos->gamma * log(U[RHO]); + return -U[RHO] * s / (eos->gamma - 1.0); +} + +// The entropy variables, in the form an entropy-stable flux consumes them. +void entropy_variables(const double *U, const EulerEos *eos, double *w) { + double rho = U[RHO]; + double inv = 1.0 / rho; + double u = U[RHOU] * inv; + double v = U[RHOV] * inv; + double z = U[RHOW] * inv; + + double p = euler_pressure(U, eos); + double beta = rho / (2.0 * p); + double s = log(p) - eos->gamma * log(rho); + + w[RHO] = (eos->gamma - s) / (eos->gamma - 1.0) - beta * (u * u + v * v + z * z); + w[RHOU] = 2.0 * beta * u; + w[RHOV] = 2.0 * beta * v; + w[RHOW] = 2.0 * beta * z; + w[RHOE] = -2.0 * beta; +} + +void entropy_gradient(const double *U, const EulerEos *eos, double *g) { + for (int i = 0; i < NVARS; i++) + g[i] = 0.0; + + __enzyme_autodiff((void *)mathematical_entropy, enzyme_dup, U, g, + enzyme_const, eos); +} + +// Forward over reverse: seed one input direction, differentiate the whole +// reverse sweep, and read off a column of the Hessian. +void entropy_hessian(const double *U, const EulerEos *eos, double *H) { + double dU[NVARS], g[NVARS], dg[NVARS]; + + for (int col = 0; col < NVARS; col++) { + for (int i = 0; i < NVARS; i++) { + dU[i] = 0.0; + dg[i] = 0.0; + } + dU[col] = 1.0; + + __enzyme_fwddiff((void *)entropy_gradient, enzyme_dup, U, dU, enzyme_const, + eos, enzyme_dup, g, dg); + + for (int row = 0; row < NVARS; row++) + H[row * NVARS + col] = dg[row]; + } +} + +// Plain forward mode over the closed-form gradient, for comparison. +void entropy_variables_jacobian(const double *U, const EulerEos *eos, + double *J) { + double dU[NVARS], w[NVARS], dw[NVARS]; + + for (int col = 0; col < NVARS; col++) { + for (int i = 0; i < NVARS; i++) + dU[i] = 0.0; + dU[col] = 1.0; + + __enzyme_fwddiff((void *)entropy_variables, enzyme_dup, U, dU, enzyme_const, + eos, enzyme_dup, w, dw); + + for (int row = 0; row < NVARS; row++) + J[row * NVARS + col] = dw[row]; + } +} + +// Cholesky without pivoting: succeeds exactly when A is positive definite. +// Returns the smallest pivot seen, or a negative value if the factorisation +// broke down. +double cholesky_min_pivot(const double *A, int n) { + double L[NVARS * NVARS]; + double min_pivot = 0.0; + + for (int i = 0; i < n * n; i++) + L[i] = 0.0; + + for (int i = 0; i < n; i++) { + for (int j = 0; j <= i; j++) { + double sum = A[i * n + j]; + for (int k = 0; k < j; k++) + sum -= L[i * n + k] * L[j * n + k]; + + if (i == j) { + if (sum <= 0.0) + return -1.0; + if (min_pivot == 0.0 || sum < min_pivot) + min_pivot = sum; + L[i * n + j] = sqrt(sum); + } else { + L[i * n + j] = sum / L[j * n + j]; + } + } + } + return min_pivot; +} + +int main() { + EulerEos eos = {1.4, 0.0}; + + for (int s = 0; s < EULER_NSTATES; s++) { + const double *q = euler_primitives[s]; + double U[NVARS], w[NVARS], g[NVARS]; + double H[NVARS * NVARS], J[NVARS * NVARS]; + + euler_from_primitive(q[0], q[1], q[2], q[3], q[4], &eos, U); + + // 1. The reverse gradient is the entropy-variable vector. + entropy_variables(U, &eos, w); + entropy_gradient(U, &eos, g); + double grad_err = frechet_rel_error(g, w, NVARS); + APPROX_EQ(grad_err, 0.0, 1e-13); + + entropy_hessian(U, &eos, H); + + // 2. Symmetry, across 25 separately computed entries. + double scale = 0.0; + for (int i = 0; i < NVARS * NVARS; i++) + if (fabs(H[i]) > scale) + scale = fabs(H[i]); + + double asym = 0.0; + for (int row = 0; row < NVARS; row++) + for (int col = 0; col < NVARS; col++) { + double d = fabs(H[row * NVARS + col] - H[col * NVARS + row]); + if (d > asym) + asym = d; + } + APPROX_EQ(asym / scale, 0.0, 1e-13); + + // 3. Convexity: Cholesky completes, so the Hessian is positive definite. + double pivot = cholesky_min_pivot(H, NVARS); + if (pivot < 0.0) { + fprintf(stderr, "entropy Hessian is not positive definite at state %d -- " + "the convexity the scheme relies on is broken\n", + s); + abort(); + } + + // 4. Nested AD against single-level AD of the closed-form gradient. + entropy_variables_jacobian(U, &eos, J); + double nest_err = frechet_rel_error(H, J, NVARS * NVARS); + APPROX_EQ(nest_err, 0.0, 1e-12); + + printf("state %d: grad %.2e asym %.2e min pivot %.3e nested vs flat " + "%.2e\n", + s, grad_err, asym / scale, pivot, nest_err); + } + + printf("done\n"); + return 0; +} diff --git a/enzyme/test/Integration/ReverseMode/newton_inversion.c b/enzyme/test/Integration/ReverseMode/newton_inversion.c new file mode 100644 index 000000000000..df1395510f71 --- /dev/null +++ b/enzyme/test/Integration/ReverseMode/newton_inversion.c @@ -0,0 +1,151 @@ +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -S | %lli - +// RUN: %clang -std=c11 -O0 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O1 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O2 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - +// RUN: %clang -std=c11 -O3 %s -S -emit-llvm -o - | %opt - %OPloadEnzyme %enzyme -enzyme-inline=1 -S | %lli - + +// Reverse mode through the same data-dependent Newton solve that +// ForwardMode/newton_inversion.c covers, which means taping a loop whose trip +// count is not known until it runs: it breaks on a converged step, can break +// instead on a collapsed heat capacity, and clamps its iterate every pass. +// +// The whole thermodynamic state goes in as one vector so a single reverse sweep +// returns every sensitivity at once, +// +// x = (e, Y_0, Y_1, T_guess), +// dT/dx = (1/cv(T), -e_0(T)/cv(T), -e_1(T)/cv(T), 0), +// +// each entry known in closed form from the implicit function theorem. The +// trailing zero is the interesting one: it says the taped adjoint of a +// converged iteration carries no memory of where the iteration started. +// +// The file also cross-checks reverse against forward on the identical kernel. +// That comparison needs no reference values and no finite differencing at all +// -- the two modes traverse the loop in opposite directions, over a tape in one +// case and not the other, so agreement to machine precision is a much tighter +// statement than either mode agreeing with a closed form. + +#include "../frechet.h" +#include "../test_utils.h" +#include "../thermo.h" + +extern double __enzyme_autodiff(void *, ...); +extern double __enzyme_fwddiff(void *, ...); +extern int enzyme_dup; +extern int enzyme_const; + +#define NINPUT 4 + +// x = (e_target, Y_0, Y_1, T_guess). +double temperature_of(const double *x, int test_before_update, int polish) { + return thermo_temperature(x[0], x + 1, x[3], test_before_update, polish); +} + +void gradient_reverse(const double *x, int test_before_update, int polish, + double *g) { + for (int i = 0; i < NINPUT; i++) + g[i] = 0.0; + + __enzyme_autodiff((void *)temperature_of, enzyme_dup, x, g, enzyme_const, + test_before_update, enzyme_const, polish); +} + +void gradient_forward(const double *x, int test_before_update, int polish, + double *g) { + double dx[NINPUT]; + + for (int col = 0; col < NINPUT; col++) { + for (int i = 0; i < NINPUT; i++) + dx[i] = 0.0; + dx[col] = 1.0; + + g[col] = __enzyme_fwddiff((void *)temperature_of, enzyme_dup, x, dx, + enzyme_const, test_before_update, enzyme_const, + polish); + } +} + +double rel(double got, double want) { return fabs(got - want) / fabs(want); } + +int main() { + const double mixtures[2][NSPECIES] = {{0.767, 0.233}, {0.5, 0.5}}; + const double temperatures[] = {200.0, 300.0, 1000.0, 3000.0, 8000.0, 15000.0}; + const int ntemp = sizeof(temperatures) / sizeof(temperatures[0]); + + for (int m = 0; m < 2; m++) { + for (int t = 0; t < ntemp; t++) { + double T_ref = temperatures[t]; + double x[NINPUT]; + + x[0] = thermo_energy(T_ref, mixtures[m]); + x[1] = mixtures[m][0]; + x[2] = mixtures[m][1]; + x[3] = 500.0; // deliberately poor guess, so the loop iterates + + double T = temperature_of(x, THERMO_TEST_AFTER_UPDATE, 0); + APPROX_EQ(rel(T, T_ref), 0.0, 1e-12); + + double cv = thermo_cv(T, x + 1); + double want[NINPUT]; + want[0] = 1.0 / cv; + want[1] = -thermo_energy_species(0, T) / cv; + want[2] = -thermo_energy_species(1, T) / cv; + want[3] = 0.0; + + double rev[NINPUT], fwd[NINPUT]; + gradient_reverse(x, THERMO_TEST_AFTER_UPDATE, 0, rev); + gradient_forward(x, THERMO_TEST_AFTER_UPDATE, 0, fwd); + + printf("mix %d T %7.1f: dT/de %.6e dT/dY0 %.6e dT/dY1 %.6e " + "dT/dTguess %.3e\n", + m, T_ref, rev[0], rev[1], rev[2], rev[3]); + + // Against the implicit-function answer. + for (int i = 0; i < 3; i++) + APPROX_EQ(rel(rev[i], want[i]), 0.0, 1e-11); + + // The adjoint keeps no memory of the initial guess. + APPROX_EQ(rev[3], 0.0, 1e-13); + + // Reverse and forward have to agree with each other, not just with the + // closed form. Scaled against the gradient as a whole, since one of its + // entries is legitimately zero. + APPROX_EQ(frechet_rel_error(rev, fwd, NINPUT), 0.0, 1e-13); + + // The polish leaves a converged adjoint alone rather than perturbing it. + double polished[NINPUT]; + gradient_reverse(x, THERMO_TEST_AFTER_UPDATE, 1, polished); + for (int i = 0; i < 3; i++) + APPROX_EQ(rel(polished[i], want[i]), 0.0, 1e-11); + } + } + + // The clamped branches tape correctly too: no adjoint flows out of a + // saturated fmax/fmin, on any input. + { + double x[NINPUT] = {-1.0e5, 0.767, 0.233, 500.0}; + double g[NINPUT]; + double T = temperature_of(x, THERMO_TEST_AFTER_UPDATE, 1); + gradient_reverse(x, THERMO_TEST_AFTER_UPDATE, 1, g); + printf("clamped low : T %g grad %g %g %g %g\n", T, g[0], g[1], g[2], g[3]); + APPROX_EQ(T, THERMO_T_MIN, 0.0); + for (int i = 0; i < NINPUT; i++) + APPROX_EQ(g[i], 0.0, 0.0); + } + { + double x[NINPUT] = {1.0e12, 0.767, 0.233, 500.0}; + double g[NINPUT]; + double T = temperature_of(x, THERMO_TEST_AFTER_UPDATE, 1); + gradient_reverse(x, THERMO_TEST_AFTER_UPDATE, 1, g); + printf("clamped high: T %g grad %g %g %g %g\n", T, g[0], g[1], g[2], g[3]); + APPROX_EQ(T, THERMO_T_MAX, 0.0); + for (int i = 0; i < NINPUT; i++) + APPROX_EQ(g[i], 0.0, 0.0); + } + + printf("done\n"); + return 0; +} diff --git a/enzyme/test/Integration/euler.h b/enzyme/test/Integration/euler.h new file mode 100644 index 000000000000..f19380238c16 --- /dev/null +++ b/enzyme/test/Integration/euler.h @@ -0,0 +1,84 @@ +// Compressible Euler physics shared by the CFD-derived integration tests. +// +// The conserved state is U = (rho, rho*u, rho*v, rho*w, rho*E) and the closure +// is the stiffened-gas (Tammann) equation of state +// +// p = (gamma - 1) * (rho*E - 0.5*rho*|u|^2) - gamma*p_inf, +// +// which reduces to the ideal gas when p_inf == 0. Both branches live in one +// closure deliberately: the ideal gas makes the normal flux homogeneous of +// degree one in U, and p_inf breaks that homogeneity by an exactly known +// amount. ForwardMode/euler_homogeneity.c gates an Enzyme-built 5x5 flux +// Jacobian against both facts. + +#include + +enum { RHO = 0, RHOU = 1, RHOV = 2, RHOW = 3, RHOE = 4, NVARS = 5 }; + +typedef struct { + double gamma; + double p_inf; +} EulerEos; + +// Static pressure from the conserved state. +double euler_pressure(const double *U, const EulerEos *eos) { + double inv = 1.0 / U[RHO]; + double ke = + 0.5 * inv * (U[RHOU] * U[RHOU] + U[RHOV] * U[RHOV] + U[RHOW] * U[RHOW]); + return (eos->gamma - 1.0) * (U[RHOE] - ke) - eos->gamma * eos->p_inf; +} + +// Frozen speed of sound. +double euler_sound_speed(const double *U, const EulerEos *eos) { + double p = euler_pressure(U, eos); + return sqrt(eos->gamma * (p + eos->p_inf) / U[RHO]); +} + +// Velocity projected onto the face normal. +double euler_normal_velocity(const double *U, const double *n) { + return (U[RHOU] * n[0] + U[RHOV] * n[1] + U[RHOW] * n[2]) / U[RHO]; +} + +// Normal flux F(U).n through a face with unit normal n. +void euler_flux(const double *U, const double *n, const EulerEos *eos, + double *F) { + double un = euler_normal_velocity(U, n); + double p = euler_pressure(U, eos); + F[RHO] = U[RHO] * un; + F[RHOU] = U[RHOU] * un + p * n[0]; + F[RHOV] = U[RHOV] * un + p * n[1]; + F[RHOW] = U[RHOW] * un + p * n[2]; + F[RHOE] = (U[RHOE] + p) * un; +} + +// Assemble a conserved state from the primitives (rho, u, v, w, p). +void euler_from_primitive(double rho, double u, double v, double w, double p, + const EulerEos *eos, double *U) { + U[RHO] = rho; + U[RHOU] = rho * u; + U[RHOV] = rho * v; + U[RHOW] = rho * w; + U[RHOE] = (p + eos->gamma * eos->p_inf) / (eos->gamma - 1.0) + + 0.5 * rho * (u * u + v * v + w * w); +} + +// Representative states, in primitive form (rho, u, v, w, p). These are the +// real operating range of a hypersonic solver rather than round numbers: the +// Mach-8 freestream carries rho ~ 1e-2 against rho*E ~ 1e5, so a probe +// direction that ignores the per-variable scale leaves the linear regime in +// one component while barely moving another. +#define EULER_NSTATES 4 +const double euler_primitives[EULER_NSTATES][5] = { + {1.225, 50.0, 10.0, -5.0, 101325.0}, // sea-level subsonic + {0.0184, 2400.0, 60.0, 0.0, 1197.0}, // Mach 8 at 30 km + {0.1400, 300.0, -120.0, 45.0, 1.6e5}, // post-shock, hot and slow + {1.0e-3, 1500.0, 0.0, -200.0, 25.0}, // strong expansion, near vacuum +}; + +// Face normals: axis-aligned, oblique, and fully three-dimensional. +#define EULER_NNORMALS 3 +const double euler_normals[EULER_NNORMALS][3] = { + {1.0, 0.0, 0.0}, + {0.6, -0.8, 0.0}, + {0.4242640687119285, 0.5656854249492380, -0.7071067811865476}, +}; diff --git a/enzyme/test/Integration/frechet.h b/enzyme/test/Integration/frechet.h new file mode 100644 index 000000000000..ff7dc8350d77 --- /dev/null +++ b/enzyme/test/Integration/frechet.h @@ -0,0 +1,86 @@ +// Frechet finite-difference gate for Enzyme-built Jacobians. +// +// Rather than compare a derivative against a hardcoded constant, compare the +// Jacobian's action against a central difference of the primal itself, +// +// J*v == (R(U + eps*v) - R(U - eps*v)) / (2*eps) + O(eps^2), +// +// so a kernel too involved to differentiate by hand can still be gated. The +// probe direction v is the part that needs care: conserved states span many +// decades (rho ~ 1e-2 against rho*E ~ 1e5 for a hypersonic freestream), so an +// unscaled direction leaves the linear regime in the small components while +// barely perturbing the large ones. frechet_direction gives every variable its +// own scale, taken from the state, and is deterministic in `seed` so a failure +// reproduces exactly. +// +// State vectors are interlaced: U[cell * nvars + var]. + +#include +#include + +// Residual R(U) evaluated into `out`; `ctx` carries whatever the kernel needs. +typedef void (*FrechetResidual)(const double *U, double *out, void *ctx); + +// Deterministic per-variable-scaled probe direction. A variable that is +// identically zero across the state falls back to unit scale so its columns are +// still probed. +void frechet_direction(int nvars, int ncells, const double *U, unsigned seed, + double *v) { + for (int var = 0; var < nvars; var++) { + double scale = 0.0; + for (int cell = 0; cell < ncells; cell++) { + double mag = fabs(U[cell * nvars + var]); + if (mag > scale) + scale = mag; + } + if (scale == 0.0) + scale = 1.0; + for (int cell = 0; cell < ncells; cell++) { + int i = cell * nvars + var; + v[i] = scale * sin(12.9898 * (double)(i + 1) + 78.233 * (double)seed); + } + } +} + +// Central-difference directional derivative of R at U along v. +void frechet_apply(FrechetResidual R, void *ctx, const double *U, int n, int m, + const double *v, double eps, double *out) { + double *Up = (double *)malloc(sizeof(double) * n); + double *Um = (double *)malloc(sizeof(double) * n); + double *Rp = (double *)malloc(sizeof(double) * m); + double *Rm = (double *)malloc(sizeof(double) * m); + + for (int i = 0; i < n; i++) { + Up[i] = U[i] + eps * v[i]; + Um[i] = U[i] - eps * v[i]; + } + R(Up, Rp, ctx); + R(Um, Rm, ctx); + for (int i = 0; i < m; i++) + out[i] = (Rp[i] - Rm[i]) / (2.0 * eps); + + free(Up); + free(Um); + free(Rp); + free(Rm); +} + +// Largest componentwise deviation of `a` from `b`, relative to the scale of b. +// Comparing against the vector scale rather than each component keeps a +// near-cancelling component from dominating the report. +double frechet_rel_error(const double *a, const double *b, int m) { + double scale = 0.0, err = 0.0; + for (int i = 0; i < m; i++) { + double mag = fabs(b[i]); + if (mag > scale) + scale = mag; + } + if (scale == 0.0) + scale = 1.0; + for (int i = 0; i < m; i++) { + double rel = fabs(a[i] - b[i]) / scale; + if (rel > err) + err = rel; + } + return err; +} diff --git a/enzyme/test/Integration/thermo.h b/enzyme/test/Integration/thermo.h new file mode 100644 index 000000000000..31d12e1af3f9 --- /dev/null +++ b/enzyme/test/Integration/thermo.h @@ -0,0 +1,132 @@ +// Thermally perfect thermodynamics and the temperature inversion that goes with +// it, shared by the forward- and reverse-mode Newton-inversion tests. +// +// Each species carries a rigid-rotor translational-rotational energy plus one +// simple-harmonic-oscillator vibrational mode, +// +// e_s(T) = R_s * (5/2 * T + theta_s / (exp(theta_s/T) - 1)), +// +// so the mixture energy is a genuinely nonlinear, non-invertible-in-closed-form +// function of temperature -- exactly the situation that forces a solver into a +// Newton loop, and the reason a CFD residual has an iterative kernel sitting in +// the middle of its differentiable path. +// +// The useful property for AD is that the loop has an exact answer that never +// mentions the iterates. Implicit differentiation of e_mix(T, Y) == e gives +// +// dT/de = 1 / cv_mix(T, Y), +// dT/dY_s = -e_s(T) / cv_mix(T, Y), +// +// so the derivative can be gated with no finite differencing and no dependence +// on how the loop happened to terminate. + +#include + +#define NSPECIES 2 + +#define THERMO_T_MIN 50.0 +#define THERMO_T_MAX 50000.0 +#define THERMO_MAX_ITER 50 +#define THERMO_REL_TOL 1.0e-6 +#define THERMO_CV_FLOOR 1.0e-8 + +// N2 and O2: gas constant [J/(kg K)] and characteristic vibrational +// temperature [K]. +const double thermo_R[NSPECIES] = {296.8035, 259.8367}; +const double thermo_theta[NSPECIES] = {3395.0, 2239.0}; + +// Species internal energy [J/kg]. +double thermo_energy_species(int s, double T) { + double x = thermo_theta[s] / T; + return thermo_R[s] * (2.5 * T + thermo_theta[s] / (exp(x) - 1.0)); +} + +// Species specific heat at constant volume [J/(kg K)] -- d/dT of the above. +double thermo_cv_species(int s, double T) { + double x = thermo_theta[s] / T; + double ex = exp(x); + double den = ex - 1.0; + return thermo_R[s] * (2.5 + x * x * ex / (den * den)); +} + +double thermo_energy(double T, const double *Y) { + double e = 0.0; + for (int s = 0; s < NSPECIES; s++) + e += Y[s] * thermo_energy_species(s, T); + return e; +} + +double thermo_cv(double T, const double *Y) { + double cv = 0.0; + for (int s = 0; s < NSPECIES; s++) + cv += Y[s] * thermo_cv_species(s, T); + return cv; +} + +// Where the convergence test sits relative to the Newton update. The two +// placements are equally common in solver code and indistinguishable for the +// primal, but they are NOT equivalent under differentiation -- see below. +#define THERMO_TEST_AFTER_UPDATE 0 +#define THERMO_TEST_BEFORE_UPDATE 1 + +// Newton solve for T given the mixture internal energy, clamped to a physical +// range. The trip count is data dependent: the loop breaks early on a +// converged step, and can break on a collapsed heat capacity instead. +// +// Differentiating the update +// +// T' = T - (e_mix(T, Y) - e) / cv_mix(T, Y) +// +// gives dT' = 1/cv + (e_mix(T,Y) - e) * cv'/cv^2 * dT, so the incoming tangent +// survives only multiplied by the residual. One Newton step applied at an +// already-converged iterate therefore DISCARDS whatever the loop accumulated +// and lands on the implicit-function answer. Everything about this kernel's AD +// accuracy follows from where that last step falls: +// +// * test-after-update runs one full iteration past convergence -- quadratic +// convergence means the step that trips a 1e-6 test is itself around 1e-13 +// -- so the loop performs the collapse on its own and the tangent comes out +// at machine precision. `polish` is then redundant. +// * test-before-update skips that final step and returns the previous +// iterate. Both the primal and the tangent are left at the step tolerance, +// and the tangent is the worse of the two by roughly an order of magnitude, +// because it carries the residual's error on top of the state's. +// * `polish` restores both to machine precision, making the result +// independent of where the test was placed. +// +// So the derivative is the least accurate thing a converged-looking solve +// returns, which is why an iterative kernel wants its tangent collapsed +// explicitly rather than by luck of loop structure. +double thermo_temperature(double e_target, const double *Y, double T_guess, + int test_before_update, int polish) { + double T = fmax(T_guess, 100.0); + + for (int iter = 0; iter < THERMO_MAX_ITER; iter++) { + double res = thermo_energy(T, Y) - e_target; + double cv = thermo_cv(T, Y); + if (cv < THERMO_CV_FLOOR) + break; + + double dT = res / cv; + if (test_before_update && fabs(dT) < THERMO_REL_TOL * T) + break; + + T -= dT; + T = fmax(T, THERMO_T_MIN); + T = fmin(T, THERMO_T_MAX); + + if (!test_before_update && fabs(dT) < THERMO_REL_TOL * T) + break; + } + + if (polish) { + double cv = thermo_cv(T, Y); + if (cv >= THERMO_CV_FLOOR) { + T -= (thermo_energy(T, Y) - e_target) / cv; + T = fmax(T, THERMO_T_MIN); + T = fmin(T, THERMO_T_MAX); + } + } + + return T; +}