2dedc35turing-sphere: reaction-diffusion on the sphere, spectral solver on WebGPUJeremy Magland 1/**
2 * Reaction-diffusion model presets, ported from websph's
3 * SphericalReactionDiffusionDriver.m. Each species k solves
4 *
5 * d(u_k)/dt = D_k * lap_s(u_k) + f_k(t, x, y, z, u_1, ..., u_N)
6 *
7 * on the unit sphere. Reactions are vectorized over the grid.
8 */
10export type Params = Record<string, number>;
12export interface ParamSpec {
13 key: string;
14 label: string;
15 value: number;
16 min: number;
17 max: number;
18 step: number;
19}
21export interface ModelSpec {
22 key: string;
23 label: string;
24 blurb: string;
25 species: string[];
26 params: ParamSpec[];
27 /** Polynomial degree of the reaction in the fields (for dealiasing). */
28 pdeg: number;
29 /** Amplitude of the random perturbation seeded into the initial state. */
30 seedAmp: number;
31 diffusivities(p: Params): number[];
32 /** Fill out[k][i] with f_k evaluated at every grid point. */
33 reaction(
34 p: Params,
35 t: number,
36 x: Float64Array,
37 y: Float64Array,
38 z: Float64Array,
39 V: ArrayLike<number>[],
40 out: Float64Array[],
41 ): void;
42 /** Fill out[k][i] with the initial condition (noise added via randn). */
43 init(
44 p: Params,
45 x: Float64Array,
46 y: Float64Array,
47 z: Float64Array,
48 randn: () => number,
49 out: Float64Array[],
50 ): void;
51}
53const schnakenberg: ModelSpec = {
54 key: 'schnakenberg',
55 label: 'Schnakenberg',
56 blurb:
57 'Turing spots. The homogeneous state is stable to uniform perturbations ' +
58 'but unstable to degrees 14 ≤ l ≤ 40, most strongly at l = 24.',
59 species: ['u', 'v'],
60 params: [
61 { key: 'a', label: 'a', value: 0.1, min: 0.01, max: 0.5, step: 0.01 },
62 { key: 'b', label: 'b', value: 0.9, min: 0.1, max: 2, step: 0.05 },
63 { key: 'D1', label: 'D₁', value: 4e-4, min: 1e-5, max: 5e-3, step: 1e-5 },
64 { key: 'D2', label: 'D₂', value: 8e-3, min: 1e-4, max: 5e-2, step: 1e-4 },
65 { key: 'dt', label: 'dt', value: 0.05, min: 0.005, max: 0.5, step: 0.005 },
66 ],
67 pdeg: 3,
68 seedAmp: 1e-2,
69 diffusivities: (p) => [p.D1, p.D2],
70 reaction(p, _t, _x, _y, _z, V, out) {
71 const [u, v] = V;
72 const [fu, fv] = out;
73 const a = p.a, b = p.b;
74 const n = fu.length;
75 for (let i = 0; i < n; i++) {
76 const ui = u[i], vi = v[i];
77 const uuv = ui * ui * vi;
78 fu[i] = a - ui + uuv;
79 fv[i] = b - uuv;
80 }
81 },
82 init(p, x, _y, _z, randn, out) {
83 const [u, v] = out;
84 const a = p.a, b = p.b;
85 const us = a + b;
86 const vs = b / (us * us);
87 const n = x.length;
88 for (let i = 0; i < n; i++) {
89 u[i] = us + this.seedAmp * randn();
90 v[i] = vs;
91 }
92 },
93};
95const brusselator: ModelSpec = {
96 key: 'brusselator',
97 label: 'Brusselator',
98 blurb:
99 'Turing stripes and spots, from a smaller diffusivity contrast than ' +
100 'Schnakenberg but with a stiffer reaction.',
101 species: ['u', 'v'],
102 params: [
103 { key: 'A', label: 'A', value: 3, min: 0.5, max: 6, step: 0.1 },
104 { key: 'B', label: 'B', value: 9, min: 1, max: 15, step: 0.25 },
105 { key: 'D1', label: 'D₁', value: 3.33e-3, min: 1e-4, max: 2e-2, step: 1e-4 },
106 { key: 'D2', label: 'D₂', value: 1.67e-2, min: 1e-3, max: 1e-1, step: 1e-3 },
107 { key: 'dt', label: 'dt', value: 0.02, min: 0.002, max: 0.1, step: 0.002 },
108 ],
109 pdeg: 3,
110 seedAmp: 1e-2,
111 diffusivities: (p) => [p.D1, p.D2],
112 reaction(p, _t, _x, _y, _z, V, out) {
113 const [u, v] = V;
114 const [fu, fv] = out;
115 const A = p.A, B = p.B;
116 const n = fu.length;
117 for (let i = 0; i < n; i++) {
118 const ui = u[i], vi = v[i];
119 const uuv = ui * ui * vi;
120 fu[i] = A - (B + 1) * ui + uuv;
121 fv[i] = B * ui - uuv;
122 }
123 },
124 init(p, x, _y, _z, randn, out) {
125 const [u, v] = out;
126 const n = x.length;
127 for (let i = 0; i < n; i++) {
128 u[i] = p.A + this.seedAmp * randn();
129 v[i] = p.B / p.A;
130 }
131 },
132};
134const allenCahn: ModelSpec = {
135 key: 'allencahn',
136 label: 'Allen–Cahn',
137 blurb:
138 'A single species: interfaces form and then coarsen until one domain ' +
139 'swallows the sphere.',
140 species: ['u'],
141 params: [
142 { key: 'eps2', label: 'ε²', value: 1e-3, min: 1e-4, max: 1e-2, step: 1e-4 },
143 { key: 'dt', label: 'dt', value: 0.02, min: 0.002, max: 0.2, step: 0.002 },
144 ],
145 pdeg: 3,
146 seedAmp: 1e-2,
147 diffusivities: (p) => [p.eps2],
148 reaction(_p, _t, _x, _y, _z, V, out) {
149 const [u] = V;
150 const [fu] = out;
151 const n = fu.length;
152 for (let i = 0; i < n; i++) {
153 const ui = u[i];
154 fu[i] = ui - ui * ui * ui;
155 }
156 },
157 init(_p, x, _y, _z, randn, out) {
158 const [u] = out;
159 const n = x.length;
160 for (let i = 0; i < n; i++) {
161 u[i] = this.seedAmp * randn();
162 }
163 },
164};
166export const models: ModelSpec[] = [schnakenberg, brusselator, allenCahn];
168export const defaultParams = (m: ModelSpec): Params =>
169 Object.fromEntries(m.params.map((p) => [p.key, p.value]));
171/** Named parameter presets shown in the UI dropdown. The pattern length
172 * scale goes as 1/sqrt(D), so scaling both diffusivities moves the spot
173 * size without changing the dynamics. */
174export interface Preset {
175 key: string;
176 label: string;
177 modelKey: string;
178 /** Overrides applied on top of the model's default parameters. */
179 params?: Params;
180}
182export const presets: Preset[] = [
183 { key: 'schnak-spots', label: 'Schnakenberg — spots', modelKey: 'schnakenberg' },
184 {
185 key: 'schnak-coarse',
186 label: 'Schnakenberg — coarse spots',
187 modelKey: 'schnakenberg',
188 params: { D1: 1e-3, D2: 2e-2 },
189 },
190 {
191 key: 'schnak-fine',
192 label: 'Schnakenberg — fine spots',
193 modelKey: 'schnakenberg',
194 params: { D1: 1.6e-4, D2: 3.2e-3 },
195 },
196 { key: 'brussel', label: 'Brusselator — stripes & spots', modelKey: 'brusselator' },
197 { key: 'allencahn', label: 'Allen–Cahn — coarsening', modelKey: 'allencahn' },
198];