1use serde::{Serialize, Deserialize};
4
5#[derive(Debug, Clone, Serialize, Deserialize)]
6pub struct ReactionDiffusionParams {
7 pub da: f64,
8 pub di: f64,
9 pub f: f64,
10 pub k: f64,
11 pub dt: f64,
12 pub dx: f64,
13}
14
15impl Default for ReactionDiffusionParams {
16 fn default() -> Self {
17 Self { da: 1.0, di: 0.5, f: 0.055, k: 0.062, dt: 1.0, dx: 1.0 }
18 }
19}
20
21pub struct GrayScott {
22 pub params: ReactionDiffusionParams,
23 pub n: usize,
24 pub u: Vec<f64>,
25 pub v: Vec<f64>,
26}
27
28impl GrayScott {
29 pub fn new(params: ReactionDiffusionParams, n: usize) -> Self {
30 let u = vec![1.0; n * n];
31 let v = vec![0.0; n * n];
32 Self { params, n, u, v }
33 }
34
35 pub fn seed_center(&mut self, radius: usize, u_val: f64, v_val: f64) {
36 let center = self.n / 2;
37 for i in (center - radius)..=(center + radius) {
38 for j in (center - radius)..=(center + radius) {
39 if i < self.n && j < self.n {
40 self.u[i * self.n + j] = u_val;
41 self.v[i * self.n + j] = v_val;
42 }
43 }
44 }
45 }
46
47 fn laplacian(&self, field: &[f64], i: usize, j: usize) -> f64 {
48 let n = self.n;
49 let dx2 = self.params.dx * self.params.dx;
50 let ip = if i + 1 < n { i + 1 } else { i };
51 let im = if i > 0 { i - 1 } else { 0 };
52 let jp = if j + 1 < n { j + 1 } else { j };
53 let jm = if j > 0 { j - 1 } else { 0 };
54 (field[ip * n + j] + field[im * n + j] + field[i * n + jp] + field[i * n + jm]
55 - 4.0 * field[i * n + j]) / dx2
56 }
57
58 pub fn step(&mut self) {
59 let dt = self.params.dt;
60 let f = self.params.f;
61 let k = self.params.k;
62 let da = self.params.da;
63 let di = self.params.di;
64 let mut u_new = self.u.clone();
65 let mut v_new = self.v.clone();
66
67 for i in 0..self.n {
68 for j in 0..self.n {
69 let idx = i * self.n + j;
70 let u = self.u[idx];
71 let v = self.v[idx];
72 let uvv = u * v * v;
73 u_new[idx] = u + dt * (da * self.laplacian(&self.u, i, j) - uvv + f * (1.0 - u));
74 v_new[idx] = v + dt * (di * self.laplacian(&self.v, i, j) + uvv - (f + k) * v);
75 }
76 }
77 self.u = u_new;
78 self.v = v_new;
79 }
80
81 pub fn evolve(&mut self, n_steps: usize) {
82 for _ in 0..n_steps { self.step(); }
83 }
84}
85
86pub struct FitzHughNagumo {
87 pub n: usize,
88 pub a: f64, pub b: f64, pub epsilon: f64,
89 pub du: f64, pub dv: f64, pub dt: f64, pub dx: f64,
90 pub u: Vec<f64>, pub v: Vec<f64>,
91}
92
93impl FitzHughNagumo {
94 pub fn new(n: usize, a: f64, b: f64, epsilon: f64, du: f64, dv: f64, dt: f64, dx: f64) -> Self {
95 Self { n, a, b, epsilon, du, dv, dt, dx, u: vec![0.0; n], v: vec![0.0; n] }
96 }
97
98 pub fn set_initial(&mut self, u: Vec<f64>, v: Vec<f64>) { self.u = u; self.v = v; }
99
100 fn laplacian_1d(&self, field: &[f64], i: usize) -> f64 {
101 let n = self.n;
102 let dx2 = self.dx * self.dx;
103 let ip = (i + 1).min(n - 1);
104 let im = if i > 0 { i - 1 } else { 0 };
105 (field[ip] + field[im] - 2.0 * field[i]) / dx2
106 }
107
108 pub fn step(&mut self) {
109 let mut u_new = self.u.clone();
110 let mut v_new = self.v.clone();
111 for i in 0..self.n {
112 let u = self.u[i]; let v = self.v[i];
113 u_new[i] = u + self.dt * (u - u.powi(3) / 3.0 - v + self.du * self.laplacian_1d(&self.u, i));
114 v_new[i] = v + self.dt * self.epsilon * (u + self.a - self.b * v + self.dv * self.laplacian_1d(&self.v, i));
115 }
116 self.u = u_new; self.v = v_new;
117 }
118
119 pub fn evolve(&mut self, n_steps: usize) { for _ in 0..n_steps { self.step(); } }
120}
121
122pub fn check_turing_instability(f_u: f64, f_v: f64, g_u: f64, g_v: f64, da: f64, di: f64) -> bool {
123 let trace = f_u + g_v;
124 let det = f_u * g_v - f_v * g_u;
125 let stable_no_diff = trace < 0.0 && det > 0.0;
126 let turing = da * g_v + di * f_u > 0.0 &&
127 (da * g_v + di * f_u).powi(2) - 4.0 * da * di * det > 0.0;
128 stable_no_diff && turing
129}
130
131pub fn turing_wavelength(f_u: f64, f_v: f64, g_u: f64, g_v: f64, da: f64, di: f64) -> f64 {
132 let k_c_sq = ((di * f_u + da * g_v) / (2.0 * da * di)).max(0.0);
133 if k_c_sq > 0.0 { 2.0 * std::f64::consts::PI / k_c_sq.sqrt() } else { f64::INFINITY }
134}
135
136#[cfg(test)]
137mod tests {
138 use super::*;
139 use approx::assert_relative_eq;
140
141 #[test]
142 fn test_gray_scott_init() {
143 let gs = GrayScott::new(ReactionDiffusionParams::default(), 10);
144 assert_eq!(gs.u.len(), 100);
145 assert!(gs.u.iter().all(|&u| u == 1.0));
146 assert!(gs.v.iter().all(|&v| v == 0.0));
147 }
148
149 #[test]
150 fn test_gray_scott_seed() {
151 let mut gs = GrayScott::new(ReactionDiffusionParams::default(), 20);
152 gs.seed_center(3, 0.5, 0.25);
153 assert!(gs.v[10 * 20 + 10] > 0.0);
154 }
155
156 #[test]
157 fn test_gray_scott_step() {
158 let mut gs = GrayScott::new(ReactionDiffusionParams::default(), 10);
159 gs.seed_center(2, 0.5, 0.25);
160 gs.step();
161 assert!(gs.u.iter().all(|&u| u.is_finite()));
162 assert!(gs.v.iter().all(|&v| v.is_finite()));
163 }
164
165 #[test]
166 fn test_gray_scott_evolve() {
167 let params = ReactionDiffusionParams { dt: 0.5, ..Default::default() };
168 let mut gs = GrayScott::new(params, 10);
169 gs.seed_center(2, 0.5, 0.25);
170 gs.evolve(10);
171 assert!(gs.u.iter().all(|&u| u.is_finite()));
172 }
173
174 #[test]
175 fn test_fitzhugh_nagumo() {
176 let mut fhn = FitzHughNagumo::new(50, 0.7, 0.8, 0.08, 1.0, 5.0, 0.01, 1.0);
177 let u: Vec<f64> = (0..50).map(|i| if i == 25 { 1.0 } else { 0.0 }).collect();
178 fhn.set_initial(u, vec![0.0; 50]);
179 fhn.evolve(100);
180 assert!(fhn.u.iter().all(|&u| u.is_finite()));
181 }
182
183 #[test]
184 fn test_turing_instability() {
185 assert!(check_turing_instability(0.5, -1.0, 1.0, -1.5, 1.0, 10.0));
190 }
191
192 #[test]
193 fn test_turing_stable() {
194 assert!(!check_turing_instability(-1.0, 1.0, -1.0, -1.0, 1.0, 1.0));
195 }
196
197 #[test]
198 fn test_turing_wavelength_finite() {
199 let lambda = turing_wavelength(1.0, -1.0, 2.0, -2.5, 1.0, 10.0);
200 assert!(lambda.is_finite() && lambda > 0.0);
201 }
202
203 #[test]
204 fn test_reaction_diffusion_params_default() {
205 let p = ReactionDiffusionParams::default();
206 assert_relative_eq!(p.f, 0.055);
207 assert_relative_eq!(p.k, 0.062);
208 }
209}