legendre 0.1.0

Block-structured, deterministic, scheduler-driven PDE simulation framework
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

//! Second-order central finite differences on Cartesian grids.

use super::{
    operators::{Discretizes, Laplacian},
    stencil::Stencil,
};
use crate::{
    core::storage::Real,
    geometry::{
        cartesian::{CartesianGrid, for_each_interior},
        grid::{BlockId, Grid},
    },
};
/// Central finite-difference policy (second order) on Cartesian grids.
#[derive(Debug, Clone, Copy, Default)]
pub struct FiniteDifference;

impl<const D: usize> Discretizes<CartesianGrid<D>, Laplacian> for FiniteDifference {
    type Stencil = CentralLaplacian;

    fn build(&self, _grid: &CartesianGrid<D>, _op: Laplacian) -> CentralLaplacian {
        CentralLaplacian
    }
}

/// (2·D+1)-point central Laplacian. Holds no geometry: spacing is read from
/// `(grid, block)` at apply time, so the same stencil serves every
/// refinement level.
#[derive(Debug, Clone, Copy, Default)]
pub struct CentralLaplacian;

impl<const D: usize> Stencil<CartesianGrid<D>> for CentralLaplacian {
    fn ghost_width(&self) -> u32 {
        1
    }

    fn apply<T: Real>(
        &self,
        grid: &CartesianGrid<D>,
        block: BlockId,
        input: <CartesianGrid<D> as Grid>::View<'_, T>,
        output: &mut <CartesianGrid<D> as Grid>::ViewMut<'_, T>,
    ) {
        let h = grid.spacing(block);
        let inv_h2: [T; D] = std::array::from_fn(|d| T::from_f64(1.0 / (h[d] * h[d])));
        for_each_interior(input.interior(), |idx| {
            let center = input.get(idx);
            let mut acc = T::ZERO;
            for d in 0..D {
                let mut plus = idx;
                plus[d] += 1;
                let mut minus = idx;
                minus[d] -= 1;
                acc += (input.get(plus) - (center + center) + input.get(minus)) * inv_h2[d];
            }
            output.set(idx, acc);
        });
    }
}