sciforge_lib/maths/pde/
diffusion.rs1pub fn heat_equation_1d_explicit(u: &[f64], dx: f64, dt: f64, alpha: f64) -> Vec<f64> {
2 let n = u.len();
3 let r = alpha * dt / (dx * dx);
4 let mut u_new = u.to_vec();
5 for i in 1..n - 1 {
6 u_new[i] = u[i] + r * (u[i + 1] - 2.0 * u[i] + u[i - 1]);
7 }
8 u_new
9}
10
11pub fn heat_equation_1d_implicit(u: &[f64], dx: f64, dt: f64, alpha: f64) -> Vec<f64> {
12 let n = u.len();
13 let r = alpha * dt / (dx * dx);
14 let mut a = vec![0.0; n];
15 let mut b = vec![0.0; n];
16 let mut c = vec![0.0; n];
17 let mut d = u.to_vec();
18 b[0] = 1.0;
19 b[n - 1] = 1.0;
20 for i in 1..n - 1 {
21 a[i] = -r;
22 b[i] = 1.0 + 2.0 * r;
23 c[i] = -r;
24 }
25 thomas_solve(&a, &b, &c, &mut d);
26 d
27}
28
29pub fn heat_equation_1d_crank_nicolson(u: &[f64], dx: f64, dt: f64, alpha: f64) -> Vec<f64> {
30 let n = u.len();
31 let r = alpha * dt / (2.0 * dx * dx);
32 let mut a = vec![0.0; n];
33 let mut b = vec![0.0; n];
34 let mut c = vec![0.0; n];
35 let mut d = vec![0.0; n];
36 b[0] = 1.0;
37 d[0] = u[0];
38 b[n - 1] = 1.0;
39 d[n - 1] = u[n - 1];
40 for i in 1..n - 1 {
41 a[i] = -r;
42 b[i] = 1.0 + 2.0 * r;
43 c[i] = -r;
44 d[i] = r * u[i - 1] + (1.0 - 2.0 * r) * u[i] + r * u[i + 1];
45 }
46 thomas_solve(&a, &b, &c, &mut d);
47 d
48}
49
50pub fn diffusion_2d_explicit(
51 u: &[Vec<f64>],
52 dx: f64,
53 dy: f64,
54 dt: f64,
55 alpha: f64,
56) -> Vec<Vec<f64>> {
57 let ny = u.len();
58 let nx = u[0].len();
59 let rx = alpha * dt / (dx * dx);
60 let ry = alpha * dt / (dy * dy);
61 let mut u_new = u.to_vec();
62 for j in 1..ny - 1 {
63 for i in 1..nx - 1 {
64 u_new[j][i] = u[j][i]
65 + rx * (u[j][i + 1] - 2.0 * u[j][i] + u[j][i - 1])
66 + ry * (u[j + 1][i] - 2.0 * u[j][i] + u[j - 1][i]);
67 }
68 }
69 u_new
70}
71
72pub fn stability_criterion_explicit(dx: f64, alpha: f64) -> f64 {
73 dx * dx / (2.0 * alpha)
74}
75
76pub fn diffusion_green_function(x: f64, t: f64, alpha: f64) -> f64 {
77 if t <= 0.0 {
78 return 0.0;
79 }
80 (1.0 / (4.0 * std::f64::consts::PI * alpha * t).sqrt()) * (-x * x / (4.0 * alpha * t)).exp()
81}
82
83pub fn advection_diffusion_1d(u: &[f64], dx: f64, dt: f64, v: f64, alpha: f64) -> Vec<f64> {
84 let n = u.len();
85 let mut u_new = u.to_vec();
86 let r = alpha * dt / (dx * dx);
87 let c = v * dt / (2.0 * dx);
88 for i in 1..n - 1 {
89 u_new[i] = u[i] + r * (u[i + 1] - 2.0 * u[i] + u[i - 1]) - c * (u[i + 1] - u[i - 1]);
90 }
91 u_new
92}
93
94pub fn diffusion_reaction_1d(
95 u: &[f64],
96 dx: f64,
97 dt: f64,
98 alpha: f64,
99 reaction_rate: f64,
100) -> Vec<f64> {
101 let n = u.len();
102 let r = alpha * dt / (dx * dx);
103 let mut u_new = u.to_vec();
104 for i in 1..n - 1 {
105 u_new[i] = u[i]
106 + r * (u[i + 1] - 2.0 * u[i] + u[i - 1])
107 + dt * reaction_rate * u[i] * (1.0 - u[i]);
108 }
109 u_new
110}
111
112pub fn fisher_kpp_step(u: &[f64], dx: f64, dt: f64, d: f64, r: f64) -> Vec<f64> {
113 let n = u.len();
114 let ratio = d * dt / (dx * dx);
115 let mut u_new = u.to_vec();
116 for i in 1..n - 1 {
117 let diffusion = ratio * (u[i + 1] - 2.0 * u[i] + u[i - 1]);
118 let reaction = r * u[i] * (1.0 - u[i]);
119 u_new[i] = u[i] + diffusion + dt * reaction;
120 }
121 u_new
122}
123
124pub fn nonlinear_diffusion_1d(u: &[f64], dx: f64, dt: f64, m: f64) -> Vec<f64> {
125 let n = u.len();
126 let mut u_new = u.to_vec();
127 for i in 1..n - 1 {
128 let d_right = 0.5 * (u[i].abs().powf(m - 1.0) + u[i + 1].abs().powf(m - 1.0));
129 let d_left = 0.5 * (u[i].abs().powf(m - 1.0) + u[i - 1].abs().powf(m - 1.0));
130 u_new[i] =
131 u[i] + dt / (dx * dx) * (d_right * (u[i + 1] - u[i]) - d_left * (u[i] - u[i - 1]));
132 }
133 u_new
134}
135
136pub fn diffusion_2d_adi(u: &[Vec<f64>], dx: f64, dy: f64, dt: f64, alpha: f64) -> Vec<Vec<f64>> {
137 let ny = u.len();
138 let nx = u[0].len();
139 let rx = alpha * dt / (2.0 * dx * dx);
140 let ry = alpha * dt / (2.0 * dy * dy);
141 let mut half = u.to_vec();
142 for j in 1..ny - 1 {
143 let mut a = vec![0.0; nx];
144 let mut b = vec![0.0; nx];
145 let mut c = vec![0.0; nx];
146 let mut d = vec![0.0; nx];
147 b[0] = 1.0;
148 d[0] = u[j][0];
149 b[nx - 1] = 1.0;
150 d[nx - 1] = u[j][nx - 1];
151 for i in 1..nx - 1 {
152 a[i] = -rx;
153 b[i] = 1.0 + 2.0 * rx;
154 c[i] = -rx;
155 d[i] = u[j][i] + ry * (u[j + 1][i] - 2.0 * u[j][i] + u[j - 1][i]);
156 }
157 thomas_solve(&a, &b, &c, &mut d);
158 half[j] = d;
159 }
160 let mut result = half.clone();
161 for i in 1..nx - 1 {
162 let mut a = vec![0.0; ny];
163 let mut b = vec![0.0; ny];
164 let mut c = vec![0.0; ny];
165 let mut d = vec![0.0; ny];
166 b[0] = 1.0;
167 d[0] = half[0][i];
168 b[ny - 1] = 1.0;
169 d[ny - 1] = half[ny - 1][i];
170 for j in 1..ny - 1 {
171 a[j] = -ry;
172 b[j] = 1.0 + 2.0 * ry;
173 c[j] = -ry;
174 d[j] = half[j][i] + rx * (half[j][i + 1] - 2.0 * half[j][i] + half[j][i - 1]);
175 }
176 thomas_solve(&a, &b, &c, &mut d);
177 for j in 0..ny {
178 result[j][i] = d[j];
179 }
180 }
181 result
182}
183
184pub fn stability_criterion_2d(dx: f64, dy: f64, alpha: f64) -> f64 {
185 1.0 / (2.0 * alpha * (1.0 / (dx * dx) + 1.0 / (dy * dy)))
186}
187
188pub fn peclet_number(velocity: f64, length: f64, diffusivity: f64) -> f64 {
189 velocity * length / diffusivity
190}
191
192pub fn diffusion_steady_state_1d(n: usize, left_bc: f64, right_bc: f64) -> Vec<f64> {
193 (0..n)
194 .map(|i| left_bc + (right_bc - left_bc) * i as f64 / (n - 1) as f64)
195 .collect()
196}
197
198pub fn diffusion_analytical_finite(x: f64, t: f64, length: f64, alpha: f64, n_terms: usize) -> f64 {
199 let pi = std::f64::consts::PI;
200 let mut sum = 0.0;
201 for n in 1..=n_terms {
202 let nf = n as f64;
203 let coeff = if n % 2 == 0 { 0.0 } else { 4.0 / (nf * pi) };
204 sum +=
205 coeff * (nf * pi * x / length).sin() * (-(nf * pi / length).powi(2) * alpha * t).exp();
206 }
207 sum
208}
209
210pub fn mass_conservation_check(u: &[f64], dx: f64) -> f64 {
211 u.iter().sum::<f64>() * dx
212}
213
214pub fn diffusion_coefficient_from_msd(msd: f64, time: f64, dimensions: u32) -> f64 {
215 msd / (2.0 * dimensions as f64 * time)
216}
217
218fn thomas_solve(a: &[f64], b: &[f64], c: &[f64], d: &mut [f64]) {
219 let n = d.len();
220 let mut cp = vec![0.0; n];
221 let mut dp = d.to_vec();
222 cp[0] = c[0] / b[0];
223 dp[0] = d[0] / b[0];
224 for i in 1..n {
225 let m = b[i] - a[i] * cp[i - 1];
226 cp[i] = if i < n - 1 { c[i] / m } else { 0.0 };
227 dp[i] = (d[i] - a[i] * dp[i - 1]) / m;
228 }
229 d[n - 1] = dp[n - 1];
230 for i in (0..n - 1).rev() {
231 d[i] = dp[i] - cp[i] * d[i + 1];
232 }
233}