numra_pde/
equations3d.rs

1//! Common 3D PDE equation builders.
2//!
3//! Convenience constructors for standard 3D PDEs using [`MOLSystem3D`].
4//! Mirrors `equations2d.rs` with a third spatial axis added.
5//!
6//! Author: Moussa Leblouba
7//! Date: 7 May 2026
8
9use crate::boundary2d::BoundaryConditions3D;
10use crate::grid::Grid3D;
11use crate::mol3d::MOLSystem3D;
12use crate::sparse_assembly::{Operator3DCoefficients, SparseScalar};
13
14/// 3D Heat (diffusion) equation: u_t = alpha * (u_xx + u_yy + u_zz).
15pub struct HeatEquation3D;
16
17impl HeatEquation3D {
18    /// Build a 3D heat equation MOL system.
19    pub fn build<S: SparseScalar>(
20        grid: Grid3D<S>,
21        alpha: S,
22        bc: &BoundaryConditions3D<S>,
23    ) -> MOLSystem3D<S> {
24        MOLSystem3D::heat(grid, alpha, bc)
25    }
26}
27
28/// 3D Advection-Diffusion: u_t = D*(u_xx + u_yy + u_zz) - vx*u_x - vy*u_y - vz*u_z.
29pub struct AdvectionDiffusion3D;
30
31impl AdvectionDiffusion3D {
32    /// Build a 3D advection-diffusion MOL system.
33    #[allow(clippy::too_many_arguments)]
34    pub fn build<S: SparseScalar>(
35        grid: Grid3D<S>,
36        diffusion: S,
37        vx: S,
38        vy: S,
39        vz: S,
40        bc: &BoundaryConditions3D<S>,
41    ) -> MOLSystem3D<S> {
42        let coeffs = Operator3DCoefficients::advection_diffusion(diffusion, vx, vy, vz);
43        MOLSystem3D::with_operator(grid, &coeffs, bc)
44    }
45}
46
47/// 3D Reaction-Diffusion: u_t = D*(u_xx + u_yy + u_zz) + R(t, x, y, z, u).
48pub struct ReactionDiffusion3D;
49
50impl ReactionDiffusion3D {
51    /// Build a 3D reaction-diffusion MOL system.
52    pub fn build<S, R>(
53        grid: Grid3D<S>,
54        diffusion: S,
55        bc: &BoundaryConditions3D<S>,
56        reaction: R,
57    ) -> MOLSystem3D<S>
58    where
59        S: SparseScalar,
60        R: Fn(S, S, S, S, S) -> S + Send + Sync + 'static,
61    {
62        MOLSystem3D::heat(grid, diffusion, bc).with_reaction(reaction)
63    }
64
65    /// Build a 3D Fisher-KPP equation: u_t = D*(u_xx + u_yy + u_zz) + r*u*(1-u).
66    pub fn fisher<S: SparseScalar>(
67        grid: Grid3D<S>,
68        diffusion: S,
69        growth_rate: S,
70        bc: &BoundaryConditions3D<S>,
71    ) -> MOLSystem3D<S> {
72        let r = growth_rate;
73        MOLSystem3D::heat(grid, diffusion, bc)
74            .with_reaction(move |_t, _x, _y, _z, u| r * u * (S::ONE - u))
75    }
76}
77
78#[cfg(test)]
79mod tests {
80    use super::*;
81    use numra_ode::{DoPri5, OdeSystem, Solver, SolverOptions};
82
83    #[test]
84    fn test_heat_equation_3d() {
85        let grid = Grid3D::uniform(0.0, 1.0, 7, 0.0, 1.0, 7, 0.0, 1.0, 7);
86        let bc = BoundaryConditions3D::all_zero_dirichlet();
87        let mol = HeatEquation3D::build(grid, 0.01_f64, &bc);
88        assert_eq!(mol.dim(), 125); // 5*5*5
89    }
90
91    #[test]
92    fn test_advection_diffusion_3d() {
93        let grid = Grid3D::uniform(0.0, 1.0, 7, 0.0, 1.0, 7, 0.0, 1.0, 7);
94        let bc = BoundaryConditions3D::all_zero_dirichlet();
95        let mol = AdvectionDiffusion3D::build(grid, 0.01, 1.0, 0.0, 0.0, &bc);
96        assert_eq!(mol.dim(), 125);
97
98        // Should run without error
99        let u0 = vec![0.0; 125];
100        let options = SolverOptions::default().rtol(1e-4);
101        let result = DoPri5::solve(&mol, 0.0, 0.01, &u0, &options).unwrap();
102        assert!(result.success);
103    }
104
105    #[test]
106    fn test_reaction_diffusion_3d_custom() {
107        // u_t = D*Laplacian(u) - 0.5 * u (linear decay)
108        let n = 7;
109        let grid = Grid3D::uniform(0.0, 1.0, n, 0.0, 1.0, n, 0.0, 1.0, n);
110        let bc = BoundaryConditions3D::all_zero_dirichlet();
111        let mol = ReactionDiffusion3D::build(grid, 0.01, &bc, |_t, _x, _y, _z, u| -0.5 * u);
112
113        let nx_int = n - 2;
114        let n_int = nx_int * nx_int * nx_int;
115        let u0 = vec![0.1; n_int];
116
117        let options = SolverOptions::default().rtol(1e-4);
118        let result = DoPri5::solve(&mol, 0.0, 0.1, &u0, &options).unwrap();
119        assert!(result.success);
120
121        // With negative reaction + zero Dirichlet, solution should decay
122        let y_final = result.y_final().unwrap();
123        for &v in &y_final {
124            assert!(v < 0.1, "Expected decay, got {}", v);
125        }
126    }
127
128    #[test]
129    fn test_fisher_3d() {
130        let n = 9;
131        let grid = Grid3D::uniform(0.0, 1.0, n, 0.0, 1.0, n, 0.0, 1.0, n);
132        let bc = BoundaryConditions3D::all_zero_dirichlet();
133        let mol = ReactionDiffusion3D::fisher(grid.clone(), 0.01, 1.0, &bc);
134
135        let nx_int = n - 2;
136        let ny_int = n - 2;
137        let nz_int = n - 2;
138        let n_int = nx_int * ny_int * nz_int;
139
140        // IC: small spherical bump in centre
141        let mut u0 = vec![0.0; n_int];
142        for kk in 0..nz_int {
143            for jj in 0..ny_int {
144                for ii in 0..nx_int {
145                    let x = grid.x_grid.points()[ii + 1];
146                    let y = grid.y_grid.points()[jj + 1];
147                    let z = grid.z_grid.points()[kk + 1];
148                    let r2 = (x - 0.5) * (x - 0.5) + (y - 0.5) * (y - 0.5) + (z - 0.5) * (z - 0.5);
149                    if r2 < 0.05 {
150                        u0[kk * (nx_int * ny_int) + jj * nx_int + ii] = 0.5;
151                    }
152                }
153            }
154        }
155
156        let options = SolverOptions::default().rtol(1e-4);
157        let result = DoPri5::solve(&mol, 0.0, 0.1, &u0, &options).unwrap();
158        assert!(result.success);
159
160        // Solution stays in [-0.1, 1.1] for Fisher-KPP
161        let y_final = result.y_final().unwrap();
162        for &v in &y_final {
163            assert!(
164                (-0.1..=1.1).contains(&v),
165                "Solution out of expected range: {}",
166                v
167            );
168        }
169    }
170}
numra_pde/equations3d.rs

numra_pde/
equations3d.rs
