qualia_core_db/solvers/calculus/
potential.rs

1//! Bounded finite-difference potential solver.
2
3#[derive(Debug, Clone, Copy, PartialEq)]
4pub enum PotentialError {
5    InvalidGrid,
6    DimensionMismatch,
7    NonFiniteInput,
8    ConvergenceFailed { residual: f64 },
9}
10
11#[repr(C)]
12#[derive(Debug, Clone, Copy, PartialEq)]
13pub struct PoissonGrid {
14    pub width: usize,
15    pub height: usize,
16    pub spacing: f64,
17}
18
19impl PoissonGrid {
20    pub fn point_count(self) -> Option<usize> {
21        self.width.checked_mul(self.height)
22    }
23}
24
25#[repr(C)]
26#[derive(Debug, Clone, Copy, PartialEq)]
27pub struct PoissonReport {
28    pub iterations: u32,
29    pub residual_inf: f64,
30    pub minimum: f64,
31    pub maximum: f64,
32}
33
34pub fn solve_poisson_dirichlet(
35    grid: PoissonGrid,
36    source: &[f64],
37    boundary_values: &[f64],
38    solution: &mut [f64],
39    tolerance: f64,
40    max_iterations: u32,
41) -> Result<PoissonReport, PotentialError> {
42    let count = grid.point_count().ok_or(PotentialError::InvalidGrid)?;
43    if grid.width < 3
44        || grid.height < 3
45        || !grid.spacing.is_finite()
46        || grid.spacing <= 0.0
47        || !tolerance.is_finite()
48        || tolerance <= 0.0
49        || max_iterations == 0
50    {
51        return Err(PotentialError::InvalidGrid);
52    }
53    if source.len() != count || boundary_values.len() != count || solution.len() != count {
54        return Err(PotentialError::DimensionMismatch);
55    }
56    if source
57        .iter()
58        .chain(boundary_values)
59        .any(|value| !value.is_finite())
60    {
61        return Err(PotentialError::NonFiniteInput);
62    }
63
64    for y in 0..grid.height {
65        for x in 0..grid.width {
66            let index = y * grid.width + x;
67            if x == 0 || y == 0 || x + 1 == grid.width || y + 1 == grid.height {
68                solution[index] = boundary_values[index];
69            } else if !solution[index].is_finite() {
70                solution[index] = 0.0;
71            }
72        }
73    }
74
75    let h2 = grid.spacing * grid.spacing;
76    let mut residual = f64::INFINITY;
77    for iteration in 1..=max_iterations {
78        for y in 1..grid.height - 1 {
79            for x in 1..grid.width - 1 {
80                let index = y * grid.width + x;
81                solution[index] = 0.25
82                    * (solution[index - 1]
83                        + solution[index + 1]
84                        + solution[index - grid.width]
85                        + solution[index + grid.width]
86                        + h2 * source[index]);
87            }
88        }
89
90        residual = 0.0;
91        for y in 1..grid.height - 1 {
92            for x in 1..grid.width - 1 {
93                let index = y * grid.width + x;
94                let discrete = (4.0 * solution[index]
95                    - solution[index - 1]
96                    - solution[index + 1]
97                    - solution[index - grid.width]
98                    - solution[index + grid.width])
99                    / h2;
100                residual = residual.max((discrete - source[index]).abs());
101            }
102        }
103        if residual <= tolerance {
104            let (minimum, maximum) = solution.iter().fold(
105                (f64::INFINITY, f64::NEG_INFINITY),
106                |(minimum, maximum), value| (minimum.min(*value), maximum.max(*value)),
107            );
108            return Ok(PoissonReport {
109                iterations: iteration,
110                residual_inf: residual,
111                minimum,
112                maximum,
113            });
114        }
115    }
116    Err(PotentialError::ConvergenceFailed { residual })
117}
118
119pub fn discrete_maximum_principle_holds(
120    grid: PoissonGrid,
121    source: &[f64],
122    boundary_values: &[f64],
123    solution: &[f64],
124    tolerance: f64,
125) -> Result<bool, PotentialError> {
126    let count = grid.point_count().ok_or(PotentialError::InvalidGrid)?;
127    if source.len() != count || boundary_values.len() != count || solution.len() != count {
128        return Err(PotentialError::DimensionMismatch);
129    }
130    if source.iter().any(|value| value.abs() > tolerance) {
131        return Ok(false);
132    }
133    let mut boundary_minimum = f64::INFINITY;
134    let mut boundary_maximum = f64::NEG_INFINITY;
135    for y in 0..grid.height {
136        for x in 0..grid.width {
137            if x == 0 || y == 0 || x + 1 == grid.width || y + 1 == grid.height {
138                let value = boundary_values[y * grid.width + x];
139                boundary_minimum = boundary_minimum.min(value);
140                boundary_maximum = boundary_maximum.max(value);
141            }
142        }
143    }
144    Ok(solution.iter().all(|value| {
145        *value >= boundary_minimum - tolerance && *value <= boundary_maximum + tolerance
146    }))
147}
148
149#[cfg(test)]
150mod tests {
151    use super::*;
152
153    #[test]
154    fn manufactured_poisson_solution_converges_with_residual() {
155        const N: usize = 25;
156        let grid = PoissonGrid {
157            width: N,
158            height: N,
159            spacing: 1.0 / (N - 1) as f64,
160        };
161        let mut source = vec![0.0; N * N];
162        let boundary = vec![0.0; N * N];
163        let mut solution = vec![0.0; N * N];
164        for y in 0..N {
165            let yy = y as f64 * grid.spacing;
166            for x in 0..N {
167                let xx = x as f64 * grid.spacing;
168                source[y * N + x] = 2.0 * (yy - yy * yy) + 2.0 * (xx - xx * xx);
169            }
170        }
171        let report =
172            solve_poisson_dirichlet(grid, &source, &boundary, &mut solution, 1e-8, 50_000).unwrap();
173        assert!(report.residual_inf <= 1e-8);
174        let mut maximum_error = 0.0_f64;
175        for y in 0..N {
176            let yy = y as f64 * grid.spacing;
177            for x in 0..N {
178                let xx = x as f64 * grid.spacing;
179                let exact = xx * (1.0 - xx) * yy * (1.0 - yy);
180                maximum_error = maximum_error.max((solution[y * N + x] - exact).abs());
181            }
182        }
183        assert!(maximum_error < 5e-4);
184    }
185
186    #[test]
187    fn harmonic_solution_obeys_discrete_maximum_principle() {
188        const N: usize = 9;
189        let grid = PoissonGrid {
190            width: N,
191            height: N,
192            spacing: 1.0 / (N - 1) as f64,
193        };
194        let source = [0.0; N * N];
195        let mut boundary = [0.0; N * N];
196        for y in 0..N {
197            boundary[y * N + N - 1] = 1.0;
198        }
199        let mut solution = [0.0; N * N];
200        solve_poisson_dirichlet(grid, &source, &boundary, &mut solution, 1e-9, 20_000).unwrap();
201        assert!(
202            discrete_maximum_principle_holds(grid, &source, &boundary, &solution, 1e-9).unwrap()
203        );
204    }
205}
qualia_core_db/solvers/calculus/potential.rs

qualia_core_db/solvers/calculus/
potential.rs
