1#![allow(clippy::needless_range_loop)]
6#[allow(unused_imports)]
7use super::functions::*;
8#[cfg(test)]
9mod tests {
10 use super::*;
11 use crate::multiscale_methods::AdaptiveCriterion;
12 use crate::multiscale_methods::Atom;
13 use crate::multiscale_methods::CgBead;
14 use crate::multiscale_methods::CoarseGrid;
15 use crate::multiscale_methods::ConcurrentDomainDecomposition;
16 use crate::multiscale_methods::EquationFreeIntegrator;
17 use crate::multiscale_methods::HmmConfig;
18 use crate::multiscale_methods::MicroDataTable;
19 use crate::multiscale_methods::MultilevelErrorEstimator;
20 use crate::multiscale_methods::PhaseFieldParams;
21 use crate::multiscale_methods::ResolutionLevel;
22 use crate::multiscale_methods::RichardsonEstimator;
23 use crate::multiscale_methods::SubDomain;
24 use crate::multiscale_methods::UnitCell;
25 #[test]
26 fn test_unit_cell_uniform_voigt() {
27 let cell = UnitCell::uniform(4, 4, 4, 0.25, 2.5);
28 let v = cell.voigt_average();
29 assert!(
30 (v - 2.5).abs() < 1e-10,
31 "uniform Voigt should equal property"
32 );
33 }
34 #[test]
35 fn test_unit_cell_uniform_reuss() {
36 let cell = UnitCell::uniform(4, 4, 4, 0.25, 3.0);
37 let r = cell.reuss_average();
38 assert!((r - 3.0).abs() < 1e-10);
39 }
40 #[test]
41 fn test_unit_cell_geometric_average() {
42 let cell = UnitCell::uniform(2, 2, 2, 0.5, std::f64::consts::E);
43 let g = cell.geometric_average();
44 assert!((g - std::f64::consts::E).abs() < 1e-8);
45 }
46 #[test]
47 fn test_voigt_geq_reuss() {
48 let cell = UnitCell::checkerboard(4, 4, 4, 0.25, 10.0, 1.0);
49 let v = cell.voigt_average();
50 let r = cell.reuss_average();
51 assert!(v >= r - 1e-10, "Voigt >= Reuss");
52 }
53 #[test]
54 fn test_checkerboard_cell_idx() {
55 let cell = UnitCell::checkerboard(3, 3, 1, 1.0, 5.0, 2.0);
56 assert_eq!(cell.idx(0, 0, 0), 0);
57 assert_eq!(cell.idx(1, 0, 0), 1);
58 }
59 #[test]
60 fn test_homogenize_hs_bounds_ordered() {
61 let cell = UnitCell::checkerboard(4, 4, 4, 0.25, 10.0, 1.0);
62 let result = homogenize_two_phase(&cell, 10.0, 1.0);
63 assert!(
64 result.hs_lower <= result.hs_upper + 1e-10,
65 "HS lower <= HS upper"
66 );
67 }
68 #[test]
69 fn test_homogenize_bounds_contain_effective() {
70 let cell = UnitCell::checkerboard(4, 4, 4, 0.25, 10.0, 1.0);
71 let result = homogenize_two_phase(&cell, 10.0, 1.0);
72 assert!(result.effective_property >= result.reuss_bound - 1e-10);
73 assert!(result.effective_property <= result.voigt_bound + 1e-10);
74 }
75 #[test]
76 fn test_homogenize_uniform_trivial() {
77 let cell = UnitCell::uniform(4, 4, 4, 0.25, 5.0);
78 let result = homogenize_two_phase(&cell, 5.0, 5.0);
79 assert!((result.voigt_bound - 5.0).abs() < 1e-8);
80 }
81 #[test]
82 fn test_cell_problem_1d_mean_zero() {
83 let cell = UnitCell::uniform(8, 1, 1, 0.125, 1.0);
84 let chi = cell_problem_1d(&cell, 2000, 1e-12);
85 let mean: f64 = chi.iter().sum::<f64>() / chi.len() as f64;
86 assert!(mean.abs() < 0.5, "cell problem solution should be bounded");
87 }
88 #[test]
89 fn test_efi_step_scalar_decay() {
90 let efi = EquationFreeIntegrator::new(0.01, 5, 5, 0.1);
91 let coarse = vec![1.0f64];
92 let result = efi.step(
93 &coarse,
94 |c| c.to_vec(),
95 |fine| fine.iter().map(|&x| x * 0.99).collect(),
96 |fine| fine.to_vec(),
97 );
98 assert!(result[0] < 1.0, "decaying system should decrease");
99 assert!(result[0] > 0.0, "should remain positive");
100 }
101 #[test]
102 fn test_efi_step_preserves_length() {
103 let efi = EquationFreeIntegrator::new(0.01, 3, 3, 0.05);
104 let coarse = vec![1.0, 2.0, 3.0];
105 let result = efi.step(
106 &coarse,
107 |c| c.to_vec(),
108 |f| f.iter().map(|&x| x * 0.999).collect(),
109 |f| f.to_vec(),
110 );
111 assert_eq!(result.len(), 3);
112 }
113 #[test]
114 fn test_hmm_flux_positive_for_positive_strain() {
115 let config = HmmConfig {
116 n_macro_points: 4,
117 micro_steps: 10,
118 dt_micro: 0.01,
119 dt_macro: 0.1,
120 micro_domain_size: 1.0,
121 };
122 let flux = hmm_estimate_flux(&config, 1.0, |_x| 1.0);
123 assert!(flux.is_finite());
124 }
125 #[test]
126 fn test_hmm_flux_zero_for_zero_strain() {
127 let config = HmmConfig {
128 n_macro_points: 4,
129 micro_steps: 10,
130 dt_micro: 0.01,
131 dt_macro: 0.1,
132 micro_domain_size: 1.0,
133 };
134 let flux = hmm_estimate_flux(&config, 0.0, |_x| 1.0);
135 assert!(flux.abs() < 1e-6, "zero strain => near-zero flux");
136 }
137 #[test]
138 fn test_subdomain_contains() {
139 let dom = SubDomain::new(0, [0.0, 0.0, 0.0], [1.0, 1.0, 1.0], false);
140 assert!(dom.contains([0.5, 0.5, 0.5]));
141 assert!(!dom.contains([1.5, 0.5, 0.5]));
142 }
143 #[test]
144 fn test_subdomain_volume() {
145 let dom = SubDomain::new(0, [0.0, 0.0, 0.0], [2.0, 3.0, 4.0], false);
146 assert!((dom.volume() - 24.0).abs() < 1e-10);
147 }
148 #[test]
149 fn test_subdomain_overlap_detection() {
150 let d0 = SubDomain::new(0, [0.0, 0.0, 0.0], [2.0, 2.0, 2.0], false);
151 let d1 = SubDomain::new(1, [1.0, 1.0, 1.0], [3.0, 3.0, 3.0], true);
152 let d2 = SubDomain::new(2, [5.0, 5.0, 5.0], [6.0, 6.0, 6.0], false);
153 assert!(d0.overlaps(&d1), "d0 and d1 should overlap");
154 assert!(!d0.overlaps(&d2), "d0 and d2 should not overlap");
155 }
156 #[test]
157 fn test_domain_decomp_interpolation() {
158 let mut dd = ConcurrentDomainDecomposition::new();
159 dd.add_domain(SubDomain::new(0, [0.0, 0.0, 0.0], [1.0, 1.0, 1.0], false));
160 dd.add_domain(SubDomain::new(1, [1.0, 0.0, 0.0], [2.0, 1.0, 1.0], false));
161 let field = [10.0, 20.0];
162 let val = dd.interpolate_field([0.5, 0.5, 0.5], &field);
163 assert!(val.is_finite());
164 assert!(val > 0.0);
165 }
166 #[test]
167 fn test_coarse_grid_average_values() {
168 let fine = vec![1.0, 3.0, 5.0, 7.0];
169 let cg = CoarseGrid::from_fine(fine, 2);
170 assert!((cg.coarse_values[0] - 2.0).abs() < 1e-10);
171 assert!((cg.coarse_values[1] - 6.0).abs() < 1e-10);
172 }
173 #[test]
174 fn test_downscale_length() {
175 let fine = vec![1.0f64; 8];
176 let cg = CoarseGrid::from_fine(fine.clone(), 2);
177 let d = cg.downscale();
178 assert_eq!(d.len(), fine.len());
179 }
180 #[test]
181 fn test_downscale_linear_length() {
182 let fine = vec![0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0];
183 let cg = CoarseGrid::from_fine(fine.clone(), 2);
184 let d = cg.downscale_linear();
185 assert_eq!(d.len(), fine.len());
186 }
187 #[test]
188 fn test_lj_potential_at_min() {
189 let sigma = 1.0;
190 let epsilon = 1.0;
191 let r_min = 2.0_f64.powf(1.0 / 6.0) * sigma;
192 let u = lennard_jones_potential(r_min, epsilon, sigma);
193 assert!(
194 (u - (-epsilon)).abs() < 1e-8,
195 "LJ minimum should be -epsilon"
196 );
197 }
198 #[test]
199 fn test_lj_force_zero_at_min() {
200 let sigma = 1.0;
201 let epsilon = 1.0;
202 let r_min = 2.0_f64.powf(1.0 / 6.0) * sigma;
203 let f = lennard_jones_force(r_min, epsilon, sigma);
204 assert!(f.abs() < 1e-7, "LJ force at minimum should be ~0");
205 }
206 #[test]
207 fn test_atom_kinetic_energy() {
208 let atom = Atom::new([0.0; 3], [1.0, 0.0, 0.0], 2.0);
209 assert!((atom.kinetic_energy() - 1.0).abs() < 1e-10);
210 }
211 #[test]
212 fn test_atomistic_temperature_positive() {
213 let atoms = vec![
214 Atom::new([0.0; 3], [1.0, 0.0, 0.0], 1.0),
215 Atom::new([1.0; 3], [0.0, 1.0, 0.0], 1.0),
216 ];
217 let t = atomistic_temperature(&atoms, 1.0);
218 assert!(t > 0.0);
219 }
220 #[test]
221 fn test_virial_stress_tensor_shape() {
222 let atoms = vec![
223 Atom::new([0.0, 0.0, 0.0], [0.0; 3], 1.0),
224 Atom::new([1.5, 0.0, 0.0], [0.0; 3], 1.0),
225 ];
226 let pairs = vec![(0, 1)];
227 let sigma = virial_stress_tensor(&atoms, &pairs, 1.0, 1.0, 1.0);
228 assert_eq!(sigma.len(), 9);
229 }
230 #[test]
231 fn test_richardson_extrapolation() {
232 let re = RichardsonEstimator::new(2.0);
233 let u_c = 1.04;
234 let u_f = 1.01;
235 let extrap = re.extrapolate(u_c, u_f);
236 assert!((extrap - 1.0).abs() < 1e-8);
237 }
238 #[test]
239 fn test_richardson_error_estimate() {
240 let re = RichardsonEstimator::new(2.0);
241 let err = re.error_estimate(1.04, 1.01);
242 assert!((err - 0.01).abs() < 1e-10);
243 }
244 #[test]
245 fn test_multilevel_estimator_convergence() {
246 let mut est = MultilevelErrorEstimator::new(3, 1e-3);
247 est.record(0, 1e-5);
248 est.record(1, 1e-5);
249 est.record(2, 1e-5);
250 assert!(est.all_converged());
251 }
252 #[test]
253 fn test_multilevel_not_converged() {
254 let mut est = MultilevelErrorEstimator::new(2, 1e-6);
255 est.record(0, 1.0);
256 assert!(!est.all_converged());
257 }
258 #[test]
259 fn test_macro_diffusion_boundary_values() {
260 let n = 10;
261 let k_eff = vec![1.0f64; n];
262 let source = vec![0.0f64; n];
263 let result = solve_macro_diffusion_1d(n, 1.0, &k_eff, &source, 0.0, 1.0);
264 assert!((result[0] - 0.0).abs() < 1e-10);
265 assert!((result[n + 1] - 1.0).abs() < 1e-10);
266 }
267 #[test]
268 fn test_macro_diffusion_linear_solution() {
269 let n = 10;
270 let k_eff = vec![1.0f64; n];
271 let source = vec![0.0f64; n];
272 let result = solve_macro_diffusion_1d(n, 1.0, &k_eff, &source, 0.0, 1.0);
273 for i in 1..=n {
274 let x = i as f64 / (n + 1) as f64;
275 assert!(
276 (result[i] - x).abs() < 1e-6,
277 "linear solution expected at node {}",
278 i
279 );
280 }
281 }
282 #[test]
283 fn test_restriction_l2_sum_preserved() {
284 let fine = vec![1.0f64, 3.0, 5.0, 7.0];
285 let coarse = restriction_l2(&fine, 2);
286 assert_eq!(coarse.len(), 2);
287 assert!((coarse[0] - 2.0).abs() < 1e-10);
288 }
289 #[test]
290 fn test_prolongation_linear_length() {
291 let coarse = vec![0.0f64, 1.0];
292 let fine = prolongation_linear(&coarse, 4);
293 assert_eq!(fine.len(), 8);
294 }
295 #[test]
296 fn test_prolongation_preserves_values_at_centres() {
297 let coarse = vec![0.0f64, 2.0];
298 let fine = prolongation_linear(&coarse, 2);
299 assert!((fine[0] - 0.0).abs() < 1e-10);
300 assert!((fine[fine.len() - 1] - 2.0).abs() < 1e-10);
301 }
302 #[test]
303 fn test_allen_cahn_step_bounds() {
304 let params = PhaseFieldParams::new(0.01, 1.0, 1.0);
305 let mut phi = vec![0.2, 0.5, 0.8, 0.5, 0.2, 0.5, 0.8, 0.5];
306 allen_cahn_step(&mut phi, ¶ms, 0.1, 0.001);
307 for &p in &phi {
308 assert!((0.0..=1.0).contains(&p), "phi out of [0,1]: {}", p);
309 }
310 }
311 #[test]
312 fn test_phase_field_params_interface_width() {
313 let params = PhaseFieldParams::new(0.04, 1.0, 2.0);
314 let w = params.interface_width();
315 assert!(w > 0.0, "interface width should be positive");
316 }
317 #[test]
318 fn test_scale_separation_pass() {
319 assert!(check_scale_separation(0.001, 1.0, 100.0));
320 }
321 #[test]
322 fn test_scale_separation_fail() {
323 assert!(!check_scale_separation(0.1, 1.0, 100.0));
324 }
325 #[test]
326 fn test_table_interpolation_boundary() {
327 let table = MicroDataTable::new(vec![0.0, 1.0, 2.0], vec![0.0, 1.0, 4.0]);
328 assert!((table.query(0.0) - 0.0).abs() < 1e-10);
329 assert!((table.query(2.0) - 4.0).abs() < 1e-10);
330 }
331 #[test]
332 fn test_table_interpolation_midpoint() {
333 let table = MicroDataTable::new(vec![0.0, 2.0], vec![0.0, 4.0]);
334 let mid = table.query(1.0);
335 assert!(
336 (mid - 2.0).abs() < 1e-10,
337 "linear interpolation at midpoint"
338 );
339 }
340 #[test]
341 fn test_table_build() {
342 let table = MicroDataTable::build(0.0, 1.0, 5, |x| x * x);
343 assert_eq!(table.inputs.len(), 5);
344 assert!((table.outputs[0] - 0.0).abs() < 1e-10);
345 assert!((table.outputs[4] - 1.0).abs() < 1e-10);
346 }
347 #[test]
348 fn test_adaptive_refine_trigger() {
349 let crit = AdaptiveCriterion::new(0.01, 0.0001, 0.5);
350 let level = crit.decide(0.1, 0.0, ResolutionLevel::Coarse);
351 assert_eq!(level, ResolutionLevel::Meso);
352 }
353 #[test]
354 fn test_adaptive_coarsen_trigger() {
355 let crit = AdaptiveCriterion::new(0.01, 0.001, 0.5);
356 let level = crit.decide(1e-5, 1e-5, ResolutionLevel::Fine);
357 assert_eq!(level, ResolutionLevel::Meso);
358 }
359 #[test]
360 fn test_adaptive_no_change() {
361 let crit = AdaptiveCriterion::new(0.1, 0.001, 1.0);
362 let level = crit.decide(0.05, 0.0, ResolutionLevel::Meso);
363 assert_eq!(level, ResolutionLevel::Meso);
364 }
365 #[test]
366 fn test_cg_bond_force_magnitude() {
367 let pa = [0.0, 0.0, 0.0];
368 let pb = [2.0, 0.0, 0.0];
369 let f = cg_bond_force(pa, pb, 10.0, 1.5);
370 let mag = norm3(f);
371 assert!(
372 (mag - 10.0 * 0.5).abs() < 1e-8,
373 "bond force magnitude should be k*(r-r0)"
374 );
375 }
376 #[test]
377 fn test_cg_bead_kinetic_energy_zero() {
378 let bead = CgBead::new([0.0; 3], 1.0);
379 let ke = 0.5 * bead.mass * dot3(bead.vel, bead.vel);
380 assert!((ke - 0.0).abs() < 1e-14);
381 }
382 #[test]
383 fn test_multigrid_vcycle_convergence() {
384 let n = 17usize;
385 let h = 1.0 / (n - 1) as f64;
386 let f = vec![0.0f64; n];
387 let mut u = vec![0.0f64; n];
388 u[0] = 0.0;
389 u[n - 1] = 1.0;
390 multigrid_vcycle(&mut u, &f, h, 3, 3);
391 assert!(u[n / 2].abs() < 2.0, "multigrid solution should be bounded");
392 }
393 #[test]
394 fn test_arlequin_blend_extremes() {
395 assert!((arlequin_blend_energy(10.0, 2.0, 1.0) - 10.0).abs() < 1e-10);
396 assert!((arlequin_blend_energy(10.0, 2.0, 0.0) - 2.0).abs() < 1e-10);
397 }
398 #[test]
399 fn test_arlequin_blend_midpoint() {
400 let blend = arlequin_blend_energy(4.0, 2.0, 0.5);
401 assert!((blend - 3.0).abs() < 1e-10);
402 }
403 #[test]
404 fn test_hill_mandel_zero_for_consistent() {
405 let s = vec![1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
406 let e = vec![1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
407 let ms = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
408 let me = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
409 let err = hill_mandel_error(&s, &e, ms, me);
410 assert!(err.is_finite());
411 }
412 #[test]
413 fn test_msd_zero_for_same_positions() {
414 let atoms = vec![Atom::new([1.0, 2.0, 3.0], [0.0; 3], 1.0)];
415 let refs = vec![[1.0, 2.0, 3.0]];
416 let msd = mean_squared_displacement(&atoms, &refs);
417 assert!(msd.abs() < 1e-14);
418 }
419 #[test]
420 fn test_msd_unit_displacement() {
421 let atoms = vec![Atom::new([1.0, 0.0, 0.0], [0.0; 3], 1.0)];
422 let refs = vec![[0.0, 0.0, 0.0]];
423 let msd = mean_squared_displacement(&atoms, &refs);
424 assert!((msd - 1.0).abs() < 1e-10);
425 }
426 #[test]
427 fn test_cauchy_born_stress_zero_for_identity() {
428 let sigma = 1.0;
429 let r_min = 2.0_f64.powf(1.0 / 6.0) * sigma;
430 let s = cauchy_born_stress_1d(1.0, 1.0, sigma, r_min);
431 assert!(
432 s.abs() < 1e-5,
433 "stress at equilibrium should be ~0, got {}",
434 s
435 );
436 }
437 #[test]
438 fn test_cauchy_born_modulus_positive() {
439 let modulus = cauchy_born_modulus(&[10.0, 10.0], &[1.0, 1.0], 2.0);
440 assert!(modulus > 0.0);
441 }
442}