Add D2 operator and heat equation

Cargo.toml

@@ -7,6 +7,7 @@ members = [
     "maxwell",
     "shallow_water",
     "gridgeneration",
+    "heat-equation",
 ]
 
 default-members = ["sbp", "euler", "maxwell", "shallow_water"]

heat-equation/Cargo.toml

@@ -0,0 +1,11 @@
+[package]
+name = "heat-equation"
+version = "0.1.0"
+authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
+edition = "2018"
+
+
+[dependencies]
+sbp = { path = "../sbp" }
+ndarray = "0.13.1"
+plotters = { version = "0.3.0", default-features = false, features = ["bitmap_gif", "bitmap_backend", "line_series"] }

heat-equation/src/main.rs

@@ -0,0 +1,68 @@
+use ndarray::Array1;
+use plotters::prelude::*;
+use sbp::{
+    integrate::{integrate, Rk4},
+    operators::{SbpOperator1d2, SBP4},
+    Float,
+};
+
+fn main() {
+    let drawing_area = BitMapBackend::gif("result.gif", (300, 300), 100)
+        .unwrap()
+        .into_drawing_area();
+    let mut chart = ChartBuilder::on(&drawing_area)
+        .x_label_area_size(40)
+        .y_label_area_size(40)
+        .build_cartesian_2d(0.0..1.0, -1.0..2.0)
+        .unwrap();
+
+    let nx: usize = 101;
+    let dt = 0.2 / nx.pow(2) as Float / 3.0;
+    let x = Array1::from_shape_fn((nx,), |i| i as Float / (nx - 1) as Float);
+    let v0 = x.map(|&x| 0.0 * x * x * x * x);
+
+    let op = SBP4;
+
+    let mut k = [v0.clone(), v0.clone(), v0.clone(), v0.clone()];
+    let rhs = move |fut: &mut Array1<Float>, prev: &Array1<Float>, _t: Float| {
+        fut.fill(0.0);
+        op.diff2(prev.view(), fut.view_mut());
+
+        let h = 1.0 / (nx - 1) as Float;
+        let tau = (1.0, -1.0);
+        let v0 = 1.0;
+        let vn = -1.0;
+
+        let d1 = op.d1();
+
+        for (d, fut) in d1.iter().zip(fut.iter_mut()) {
+            *fut += tau.0 / (h * h) * d * (prev[0] - v0);
+        }
+        for (d, fut) in d1.iter().zip(fut.iter_mut().rev().take(d1.len()).rev()) {
+            *fut += tau.1 / (h * h) * d * (prev[nx - 1] - vn);
+        }
+    };
+
+    let mut v1 = v0.clone();
+    let mut v2 = v0;
+    for i in 0..90 {
+        if i % 3 == 0 {
+            drawing_area.fill(&WHITE).unwrap();
+            chart
+                .configure_mesh()
+                .x_desc("x")
+                .y_desc("y")
+                .draw()
+                .unwrap();
+            chart
+                .draw_series(LineSeries::new(
+                    x.iter().zip(v1.iter()).map(|(&x, &y)| (x, y)),
+                    &BLACK,
+                ))
+                .unwrap();
+            drawing_area.present().unwrap();
+        }
+        integrate::<Rk4, _, _, _>(rhs, &v1, &mut v2, &mut 0.0, dt, &mut k);
+        std::mem::swap(&mut v1, &mut v2);
+    }
+}

sbp/src/operators.rs

@@ -17,6 +17,12 @@ pub trait SbpOperator1d: Send + Sync {
     fn h_matrix(&self, n: usize) -> sprs::CsMat<Float>;
 }
 
+pub trait SbpOperator1d2: SbpOperator1d {
+    fn diff2(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
+    /// Lacks a scaling of 1/h^2
+    fn d1(&self) -> &[Float];
+}
+
 pub trait SbpOperator2d: Send + Sync {
     fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
     fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);

sbp/src/operators/traditional4.rs

@@ -16,10 +16,22 @@ impl SBP4 {
     ];
     #[rustfmt::skip]
     const BLOCK: &'static [&'static [Float]] = &[
-        &[-1.41176470588235e+00, 1.73529411764706e+00, -2.35294117647059e-01, -8.82352941176471e-02],
-        &[-5.00000000000000e-01, 0.00000000000000e+00, 5.00000000000000e-01],
-        &[9.30232558139535e-02, -6.86046511627907e-01, 0.00000000000000e+00, 6.86046511627907e-01, -9.30232558139535e-02],
-        &[3.06122448979592e-02, 0.00000000000000e+00, -6.02040816326531e-01, 0.00000000000000e+00, 6.53061224489796e-01, -8.16326530612245e-02],
+        &[-24.0/17.0, 59.0/34.0, -4.0/17.0, -3.0/34.0],
+        &[-1.0/2.0, 0.0, 1.0/2.0],
+        &[4.0/43.0, -59.0/86.0, 0.0, 59.0/86.0, -4.0/43.0],
+        &[3.0/98.0, 0.0, -59.0/98.0, 0.0, 32.0/49.0, -4.0/49.0]
+    ];
+
+    #[rustfmt::skip]
+    const D2DIAG: &'static [Float] = &[
+        -1.0 / 12.0, 4.0 / 3.0, -5.0 / 2.0, 4.0 / 3.0, -1.0 / 12.0
+    ];
+    #[rustfmt::skip]
+    const D2BLOCK: &'static [&'static [Float]] = &[
+        &[2.0, -5.0, 4.0, -1.0],
+        &[1.0, -2.0, 1.0],
+        &[-4.0/43.0, 59.0/43.0, -110.0/43.0, 59.0/43.0, -4.0/43.0],
+        &[-1.0/49.0, 0.0, 59.0/49.0, -118.0/49.0, 64.0/49.0, -4.0/49.0]
     ];
 }
 
@@ -81,6 +93,24 @@ impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &SBP4) {
     }
 }
 
+impl super::SbpOperator1d2 for SBP4 {
+    fn diff2(&self, prev: ArrayView1<Float>, mut fut: ArrayViewMut1<Float>) {
+        super::diff_op_1d(
+            Self::D2BLOCK,
+            Self::D2DIAG,
+            super::Symmetry::Symmetric,
+            super::OperatorType::Normal,
+            prev,
+            fut.view_mut(),
+        );
+        let hi = (prev.len() - 1) as Float;
+        fut.map_inplace(|x| *x *= hi)
+    }
+    fn d1(&self) -> &[Float] {
+        &[-88.0 / 17.0, 144.0 / 17.0, -72.0 / 17.0, 16.0 / 17.0]
+    }
+}
+
 #[test]
 fn test_trad4() {
     use super::testing::*;