1/**
2 * Solver correctness tests against the f64 CPU transform backend.
3 *
4 * A. Linear reaction + diffusion, single mode: every (l,m) mode of
5 * f = c*u with implicit diffusion follows the exact scalar recurrence
6 * g = (1 + dt*c) / (1 + dt*D*l(l+1)).
7 * B. Uniform state, nonlinear reaction: the l=0 mode follows the explicit
8 * Euler map of the reaction ODE exactly.
9 * C. Turing linear stability: a small single-mode perturbation of the
10 * Schnakenberg fixed point follows the 2x2 linearized IMEX recurrence,
11 * and the (24, 7) mode lies in the unstable band.
12 *
13 * Run: node scripts/test-node.ts
14 */
15import { CpuBackend } from '../src/solver/backend.ts';
16import { Simulation, gridForLmax } from '../src/solver/simulation.ts';
17import { models, defaultParams } from '../src/solver/models.ts';
18import type { ModelSpec } from '../src/solver/models.ts';
19import { lmIndex } from '../src/sht/layout.ts';
21let failures = 0;
22function check(name: string, ok: boolean, detail: string): void {
23 console.log(`${ok ? 'PASS' : 'FAIL'} ${name} ${detail}`);
24 if (!ok) failures++;
25}
27// ---------------------------------------------------------------- test A
28{
29 const lmax = 15;
30 const { nlat, nphi } = gridForLmax(lmax, 1);
31 const backend = new CpuBackend({ lmax, mmax: lmax, nlat, nphi });
32 const c = -0.3;
33 const D = 0.01;
34 const model: ModelSpec = {
35 key: 'linear', label: 'linear', blurb: '', species: ['u'],
36 params: [], pdeg: 1, seedAmp: 0,
37 diffusivities: () => [D],
38 reaction(_p, _t, _x, _y, _z, V, out) {
39 for (let i = 0; i < out[0].length; i++) out[0][i] = c * V[0][i];
40 },
41 init() {},
42 };
43 const sim = new Simulation(backend, model, { dt: 0.1 });
44 const l = 5, m = 2;
45 const idx = lmIndex(lmax, l, m);
46 sim.U[0][2 * idx] = 0.8;
47 sim.U[0][2 * idx + 1] = -0.35;
49 const nsteps = 20;
50 for (let s = 0; s < nsteps; s++) await sim.step();
52 const g = (1 + 0.1 * c) / (1 + 0.1 * D * l * (l + 1));
53 const gn = Math.pow(g, nsteps);
54 const errRe = Math.abs(sim.U[0][2 * idx] - 0.8 * gn);
55 const errIm = Math.abs(sim.U[0][2 * idx + 1] - -0.35 * gn);
56 let leak = 0;
57 for (let i = 0; i < backend.nlm; i++) {
58 if (i === idx) continue;
59 leak = Math.max(leak, Math.abs(sim.U[0][2 * i]), Math.abs(sim.U[0][2 * i + 1]));
60 }
61 check('A: single-mode linear recurrence', errRe < 1e-12 && errIm < 1e-12,
62 `err=(${errRe.toExponential(2)}, ${errIm.toExponential(2)})`);
63 check('A: no leakage into other modes', leak < 1e-12, `leak=${leak.toExponential(2)}`);
64}
66// ---------------------------------------------------------------- test B
67{
68 const schnak = models[0];
69 const p = defaultParams(schnak);
70 const lmax = 15;
71 const { nlat, nphi } = gridForLmax(lmax, schnak.pdeg);
72 const backend = new CpuBackend({ lmax, mmax: lmax, nlat, nphi });
73 const uniform: ModelSpec = {
74 ...schnak,
75 seedAmp: 0,
76 init(pp, x, _y, _z, _randn, out) {
77 out[0].fill(1.2);
78 out[1].fill(0.8);
79 void pp; void x;
80 },
81 };
82 const sim = new Simulation(backend, uniform, p);
83 await sim.init(1);
84 const nsteps = 50;
85 for (let s = 0; s < nsteps; s++) await sim.step();
87 // reference: explicit Euler on the 2-species ODE (l=0 is untouched by diffusion)
88 let u = 1.2, v = 0.8;
89 for (let s = 0; s < nsteps; s++) {
90 const fu = p.a - u + u * u * v;
91 const fv = p.b - u * u * v;
92 u += p.dt * fu;
93 v += p.dt * fv;
94 }
95 // the area mean is the l=0 coefficient of U (V lags U by one step)
96 const sqrt4pi = Math.sqrt(4 * Math.PI);
97 const i00 = 2 * lmIndex(lmax, 0, 0);
98 const errU = Math.abs(sim.U[0][i00] / sqrt4pi - u);
99 const errV = Math.abs(sim.U[1][i00] / sqrt4pi - v);
100 check('B: uniform nonlinear reaction ODE', errU < 1e-10 && errV < 1e-10,
101 `err=(${errU.toExponential(2)}, ${errV.toExponential(2)}) u=${u.toFixed(6)} v=${v.toFixed(6)}`);
102}
104// ---------------------------------------------------------------- test C
105{
106 const schnak = models[0];
107 const p = defaultParams(schnak);
108 const lmax = 31;
109 const { nlat, nphi } = gridForLmax(lmax, schnak.pdeg);
110 const backend = new CpuBackend({ lmax, mmax: lmax, nlat, nphi });
111 const sim = new Simulation(backend, schnak, p);
113 const us = p.a + p.b; // 1.0
114 const vs = p.b / (us * us); // 0.9
115 const sqrt4pi = Math.sqrt(4 * Math.PI);
116 const l = 24, m = 7;
117 const idx = lmIndex(lmax, l, m);
118 const eps = 1e-6;
119 const c0 = [eps, 0.5 * eps];
120 // fixed point + single-mode perturbation, set directly in spectral space
121 sim.U[0][2 * lmIndex(lmax, 0, 0)] = us * sqrt4pi;
122 sim.U[1][2 * lmIndex(lmax, 0, 0)] = vs * sqrt4pi;
123 sim.U[0][2 * idx] = c0[0];
124 sim.U[1][2 * idx] = c0[1];
126 const nsteps = 20;
127 for (let s = 0; s < nsteps; s++) await sim.step();
129 // linearized IMEX recurrence: c' = diag(1/(1+dt*Dk*lam)) * (I + dt*J) * c
130 const lam = l * (l + 1);
131 const J = [
132 [-1 + 2 * us * vs, us * us],
133 [-2 * us * vs, -us * us],
134 ];
135 let c = [...c0];
136 for (let s = 0; s < nsteps; s++) {
137 const r0 = c[0] + p.dt * (J[0][0] * c[0] + J[0][1] * c[1]);
138 const r1 = c[1] + p.dt * (J[1][0] * c[0] + J[1][1] * c[1]);
139 c = [r0 / (1 + p.dt * p.D1 * lam), r1 / (1 + p.dt * p.D2 * lam)];
140 }
141 const got = [sim.U[0][2 * idx], sim.U[1][2 * idx]];
142 const errU = Math.abs(got[0] - c[0]) / Math.abs(c[0]);
143 const errV = Math.abs(got[1] - c[1]) / Math.abs(c[1]);
144 check('C: linearized Turing-mode recurrence', errU < 1e-4 && errV < 1e-4,
145 `rel err=(${errU.toExponential(2)}, ${errV.toExponential(2)})`);
146 check('C: (l=24, m=7) is growing', Math.abs(got[0]) > Math.abs(c0[0]),
147 `|c|: ${Math.abs(c0[0]).toExponential(2)} -> ${Math.abs(got[0]).toExponential(2)}`);
148}
150console.log(failures === 0 ? '\nAll tests passed.' : `\n${failures} test(s) FAILED.`);
151process.exit(failures === 0 ? 0 : 1);