noether-guard 0.1.0

Physics linter — verify conservation laws in simulations using Noether's theorem
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

use crate::symmetry::Symmetry;

/// A conservation law derived from a symmetry via Noether's theorem.
#[derive(Debug, Clone)]
pub struct ConservationLaw {
    pub symmetry: Symmetry,
    pub name: String,
    pub initial_value: f64,
    pub tolerance: f64,
}

impl ConservationLaw {
    /// Create a new conservation law.
    pub fn new(symmetry: Symmetry, initial_value: f64, tolerance: f64) -> Self {
        Self {
            name: symmetry.conserved_quantity().to_string(),
            symmetry,
            initial_value,
            tolerance,
        }
    }

    /// Check whether `value` is conserved within tolerance.
    pub fn is_conserved(&self, value: f64) -> bool {
        (value - self.initial_value).abs() <= self.tolerance
    }

    /// Compute the absolute drift from the initial value.
    pub fn drift(&self, value: f64) -> f64 {
        (value - self.initial_value).abs()
    }

    /// Compute the drift rate given a previous drift and elapsed time.
    pub fn drift_rate(&self, previous_drift: f64, dt: f64) -> f64 {
        if dt > 0.0 {
            (self.initial_value - previous_drift).abs() / dt
        } else {
            0.0
        }
    }
}

/// Tracks a single conserved quantity over time.
#[derive(Debug, Clone)]
pub struct ConservedQuantity {
    pub name: String,
    pub values: Vec<(f64, f64)>,
    pub drift_rate: f64,
}

impl ConservedQuantity {
    pub fn new(name: &str) -> Self {
        Self {
            name: name.to_string(),
            values: Vec::new(),
            drift_rate: 0.0,
        }
    }

    /// Record a (time, value) pair.
    pub fn record(&mut self, time: f64, value: f64) {
        self.values.push((time, value));
        if self.values.len() >= 2 {
            let (t0, v0) = self.values[self.values.len() - 2];
            let (t1, v1) = self.values[self.values.len() - 1];
            let dt = t1 - t0;
            if dt > 0.0 {
                self.drift_rate = (v1 - v0) / dt;
            }
        }
    }

    /// Current drift rate.
    pub fn current_drift_rate(&self) -> f64 {
        self.drift_rate
    }

    /// Total accumulated drift from first recorded value.
    pub fn total_drift(&self) -> f64 {
        if self.values.len() < 2 {
            return 0.0;
        }
        let first = self.values.first().unwrap().1;
        let last = self.values.last().unwrap().1;
        (last - first).abs()
    }

    /// Number of recorded data points.
    pub fn len(&self) -> usize {
        self.values.len()
    }

    /// Whether any data has been recorded.
    pub fn is_empty(&self) -> bool {
        self.values.is_empty()
    }
}