numra_pde/
equations2d.rs

1//! Common 2D PDE equation builders.
2//!
3//! Convenience constructors for standard 2D PDEs using [`MOLSystem2D`].
4//!
5//! Author: Moussa Leblouba
6//! Date: 9 February 2026
7//! Modified: 2 May 2026
8
9use crate::boundary2d::BoundaryConditions2D;
10use crate::grid::Grid2D;
11use crate::mol2d::MOLSystem2D;
12use crate::sparse_assembly::{Operator2DCoefficients, SparseScalar};
13
14/// 2D Heat (diffusion) equation: u_t = alpha * (u_xx + u_yy).
15pub struct HeatEquation2D;
16
17impl HeatEquation2D {
18    /// Build a 2D heat equation MOL system.
19    pub fn build<S: SparseScalar>(
20        grid: Grid2D<S>,
21        alpha: S,
22        bc: &BoundaryConditions2D<S>,
23    ) -> MOLSystem2D<S> {
24        MOLSystem2D::heat(grid, alpha, bc)
25    }
26}
27
28/// 2D Advection-Diffusion: u_t = D*(u_xx + u_yy) - vx*u_x - vy*u_y.
29pub struct AdvectionDiffusion2D;
30
31impl AdvectionDiffusion2D {
32    /// Build a 2D advection-diffusion MOL system.
33    pub fn build<S: SparseScalar>(
34        grid: Grid2D<S>,
35        diffusion: S,
36        vx: S,
37        vy: S,
38        bc: &BoundaryConditions2D<S>,
39    ) -> MOLSystem2D<S> {
40        let coeffs = Operator2DCoefficients::advection_diffusion(diffusion, vx, vy);
41        MOLSystem2D::with_operator(grid, &coeffs, bc)
42    }
43}
44
45/// 2D Reaction-Diffusion: u_t = D*(u_xx + u_yy) + R(t, x, y, u).
46pub struct ReactionDiffusion2D;
47
48impl ReactionDiffusion2D {
49    /// Build a 2D reaction-diffusion MOL system.
50    pub fn build<S, R>(
51        grid: Grid2D<S>,
52        diffusion: S,
53        bc: &BoundaryConditions2D<S>,
54        reaction: R,
55    ) -> MOLSystem2D<S>
56    where
57        S: SparseScalar,
58        R: Fn(S, S, S, S) -> S + Send + Sync + 'static,
59    {
60        MOLSystem2D::heat(grid, diffusion, bc).with_reaction(reaction)
61    }
62
63    /// Build a 2D Fisher equation: u_t = D*(u_xx + u_yy) + r*u*(1-u).
64    pub fn fisher<S: SparseScalar>(
65        grid: Grid2D<S>,
66        diffusion: S,
67        growth_rate: S,
68        bc: &BoundaryConditions2D<S>,
69    ) -> MOLSystem2D<S> {
70        let r = growth_rate;
71        MOLSystem2D::heat(grid, diffusion, bc)
72            .with_reaction(move |_t, _x, _y, u| r * u * (S::ONE - u))
73    }
74}
75
76#[cfg(test)]
77mod tests {
78    use super::*;
79    use numra_ode::{DoPri5, OdeSystem, Solver, SolverOptions};
80
81    #[test]
82    fn test_heat_equation_2d() {
83        let grid = Grid2D::uniform(0.0, 1.0, 11, 0.0, 1.0, 11);
84        let bc = BoundaryConditions2D::all_zero_dirichlet();
85        let mol = HeatEquation2D::build(grid, 0.01_f64, &bc);
86        assert_eq!(mol.dim(), 81);
87    }
88
89    #[test]
90    fn test_advection_diffusion_2d() {
91        let grid = Grid2D::uniform(0.0, 1.0, 11, 0.0, 1.0, 11);
92        let bc = BoundaryConditions2D::all_zero_dirichlet();
93        let mol = AdvectionDiffusion2D::build(grid, 0.01, 1.0, 0.0, &bc);
94        assert_eq!(mol.dim(), 81);
95
96        // Should run without error
97        let u0 = vec![0.0; 81];
98        let options = SolverOptions::default().rtol(1e-4);
99        let result = DoPri5::solve(&mol, 0.0, 0.01, &u0, &options).unwrap();
100        assert!(result.success);
101    }
102
103    #[test]
104    fn test_fisher_2d() {
105        let n = 11;
106        let grid = Grid2D::uniform(0.0, 1.0, n, 0.0, 1.0, n);
107        let bc = BoundaryConditions2D::all_zero_dirichlet();
108        let mol = ReactionDiffusion2D::fisher(grid.clone(), 0.01, 1.0, &bc);
109
110        let nx_int = n - 2;
111        let ny_int = n - 2;
112        let n_int = nx_int * ny_int;
113
114        // IC: small perturbation in center
115        let mut u0 = vec![0.0; n_int];
116        for jj in 0..ny_int {
117            for ii in 0..nx_int {
118                let x = grid.x_grid.points()[ii + 1];
119                let y = grid.y_grid.points()[jj + 1];
120                let r2 = (x - 0.5) * (x - 0.5) + (y - 0.5) * (y - 0.5);
121                if r2 < 0.04 {
122                    u0[jj * nx_int + ii] = 0.5;
123                }
124            }
125        }
126
127        let options = SolverOptions::default().rtol(1e-4);
128        let result = DoPri5::solve(&mol, 0.0, 0.1, &u0, &options).unwrap();
129        assert!(result.success);
130    }
131}
numra_pde/equations2d.rs

numra_pde/
equations2d.rs
