aeon-tk 0.3.3

Toolkit for running finite difference simulations with adaptive mesh refinement
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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150

use crate::kernel::{Border, Derivative, Kernel, SecondDerivative, Value};

// ************************************
// Convolution ************************
// ************************************

/// A N-dimensional tensor product of several seperable kernels.
pub trait Convolution<const N: usize> {
    fn border_width(&self, axis: usize) -> usize;
    fn interior(&self, axis: usize) -> &[f64];
    fn free(&self, border: Border, axis: usize) -> &[f64];
    fn scale(&self, spacing: [f64; N]) -> f64;
}

macro_rules! impl_convolution_for_tuples {
    ($($N:literal => $($T:ident $i:tt)*)*) => {
        $(
            impl<$($T: Kernel,)*> Convolution<$N> for ($($T,)*) {
                fn border_width(&self, axis: usize) -> usize {
                    match axis {
                        $($i => self.$i.border_width(),)*
                        _ => panic!("Invalid Axis")
                    }
                }

                fn interior(&self, axis: usize) -> &[f64] {
                    match axis {
                        $($i => self.$i.interior(),)*
                        _ => panic!("Invalid Axis")
                    }
                }

                fn free(&self, border: Border, axis: usize) -> &[f64] {
                    match axis {
                        $($i => self.$i.free(border),)*
                        _ => panic!("Invalid Axis")
                    }
                }

                fn scale(&self, spacing: [f64; $N]) -> f64 {
                    let mut result = 1.0;

                    $(
                        result *= self.$i.scale(spacing[$i]);
                    )*

                    result
                }
            }
        )*

    };
}

impl_convolution_for_tuples! {
   1 => K0 0
   2 => K0 0 K1 1
   3 => K0 0 K1 1 K2 2
   4 => K0 0 K1 1 K2 2 K3 3
}

/// Computes the gradient along the given axis.
#[derive(Clone)]
pub struct Gradient<const ORDER: usize>(pub usize);

impl<const N: usize, const ORDER: usize> Convolution<N> for Gradient<ORDER> {
    fn border_width(&self, axis: usize) -> usize {
        if axis == self.0 {
            Derivative::<ORDER>.border_width()
        } else {
            Value.border_width()
        }
    }

    fn interior(&self, axis: usize) -> &[f64] {
        if axis == self.0 {
            Derivative::<ORDER>.interior()
        } else {
            Value.interior()
        }
    }

    fn free(&self, border: Border, axis: usize) -> &[f64] {
        if axis == self.0 {
            Derivative::<ORDER>.free(border)
        } else {
            Value.free(border)
        }
    }

    fn scale(&self, spacing: [f64; N]) -> f64 {
        1.0 / spacing[self.0]
    }
}

/// Computes the mixed derivative of the given axes.
#[derive(Debug, Clone, Copy)]
pub struct Hessian<const ORDER: usize>(pub usize, pub usize);

impl<const ORDER: usize> Hessian<ORDER> {
    /// Constructs a convolution which computes the given entry of the hessian matrix.
    pub const fn new(i: usize, j: usize) -> Self {
        Self(i, j)
    }
}

impl<const ORDER: usize> Hessian<ORDER> {
    fn is_second(&self, axis: usize) -> bool {
        self.0 == self.1 && axis == self.0
    }

    fn is_first(&self, axis: usize) -> bool {
        self.0 == axis || self.1 == axis
    }
}

impl<const N: usize, const ORDER: usize> Convolution<N> for Hessian<ORDER> {
    fn border_width(&self, axis: usize) -> usize {
        if self.is_second(axis) {
            SecondDerivative::<ORDER>.border_width()
        } else if self.is_first(axis) {
            Derivative::<ORDER>.border_width()
        } else {
            Value.border_width()
        }
    }

    fn interior(&self, axis: usize) -> &[f64] {
        if self.is_second(axis) {
            SecondDerivative::<ORDER>.interior()
        } else if self.is_first(axis) {
            Derivative::<ORDER>.interior()
        } else {
            Value.interior()
        }
    }

    fn free(&self, border: Border, axis: usize) -> &[f64] {
        if self.is_second(axis) {
            SecondDerivative::<ORDER>.free(border)
        } else if self.is_first(axis) {
            Derivative::<ORDER>.free(border)
        } else {
            Value.free(border)
        }
    }
    fn scale(&self, spacing: [f64; N]) -> f64 {
        1.0 / (spacing[self.0] * spacing[self.1])
    }
}