add upwind4

src/lib.rs

@@ -101,13 +101,24 @@ impl Universe {
 
 fn diffx(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
     for j in 0..prev.shape()[0] {
-        trad4(prev.slice(s!(j, ..)), fut.slice_mut(s!(j, ..)));
+        upwind4(prev.slice(s!(j, ..)), fut.slice_mut(s!(j, ..)));
     }
 }
 
 fn diffy(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
     for i in 0..prev.shape()[1] {
-        trad4(prev.slice(s!(.., i)), fut.slice_mut(s!(.., i)));
+        upwind4(prev.slice(s!(.., i)), fut.slice_mut(s!(.., i)));
+    }
+}
+
+fn dissx(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
+    for j in 0..prev.shape()[0] {
+        upwind4_diss(prev.slice(s!(j, ..)), fut.slice_mut(s!(j, ..)));
+    }
+}
+fn dissy(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
+    for i in 0..prev.shape()[1] {
+        upwind4_diss(prev.slice(s!(.., i)), fut.slice_mut(s!(.., i)));
     }
 }
 
@@ -152,3 +163,157 @@ fn trad4(prev: ArrayView1<f32>, mut fut: ArrayViewMut1<f32>) {
         + diag[4] * prev[(1)];
     fut[(nx - 1)] += diff / dx;
 }
+
+fn upwind4(prev: ArrayView1<f32>, mut fut: ArrayViewMut1<f32>) {
+    assert_eq!(prev.shape(), fut.shape());
+    let nx = prev.shape()[0];
+
+    let dx = 1.0 / (nx - 1) as f32;
+
+    let diag = [
+        -1.0 / 24.0,
+        1.0 / 4.0,
+        -7.0 / 8.0,
+        0.0,
+        7.0 / 8.0,
+        -1.0 / 4.0,
+        1.0 / 24.0,
+    ];
+
+    let diff = diag[0] * prev[(nx - 3)]
+        + diag[1] * prev[(nx - 2)]
+        + diag[2] * prev[(nx - 1)]
+        + diag[3] * prev[(0)]
+        + diag[4] * prev[(1)]
+        + diag[5] * prev[(2)]
+        + diag[6] * prev[(3)];
+    fut[0] += diff / dx;
+    let diff = diag[0] * prev[(nx - 2)]
+        + diag[1] * prev[(nx - 1)]
+        + diag[2] * prev[(0)]
+        + diag[3] * prev[(1)]
+        + diag[4] * prev[(2)]
+        + diag[5] * prev[(3)]
+        + diag[6] * prev[(4)];
+    fut[1] += diff / dx;
+    let diff = diag[0] * prev[(nx - 1)]
+        + diag[1] * prev[(0)]
+        + diag[2] * prev[(1)]
+        + diag[3] * prev[(2)]
+        + diag[4] * prev[(3)]
+        + diag[5] * prev[(4)]
+        + diag[6] * prev[(5)];
+    fut[2] += diff / dx;
+
+    for i in 3..nx - 3 {
+        let diff = diag[0] * prev[(i - 3)]
+            + diag[1] * prev[(i - 2)]
+            + diag[2] * prev[(i - 1)]
+            + diag[3] * prev[(i)]
+            + diag[4] * prev[(i + 1)]
+            + diag[5] * prev[(i + 2)]
+            + diag[6] * prev[(i + 3)];
+        fut[(i)] += diff / dx;
+    }
+    let diff = diag[0] * prev[(nx - 6)]
+        + diag[1] * prev[(nx - 5)]
+        + diag[2] * prev[(nx - 4)]
+        + diag[3] * prev[(nx - 3)]
+        + diag[4] * prev[(nx - 2)]
+        + diag[5] * prev[(nx - 1)]
+        + diag[6] * prev[(0)];
+    fut[(nx - 3)] += diff / dx;
+    let diff = diag[0] * prev[(nx - 5)]
+        + diag[1] * prev[(nx - 4)]
+        + diag[2] * prev[(nx - 3)]
+        + diag[3] * prev[(nx - 2)]
+        + diag[4] * prev[(nx - 1)]
+        + diag[5] * prev[(0)]
+        + diag[6] * prev[(1)];
+    fut[(nx - 2)] += diff / dx;
+    let diff = diag[0] * prev[(nx - 4)]
+        + diag[1] * prev[(nx - 3)]
+        + diag[2] * prev[(nx - 2)]
+        + diag[3] * prev[(nx - 1)]
+        + diag[4] * prev[(0)]
+        + diag[5] * prev[(1)]
+        + diag[6] * prev[(2)];
+    fut[(nx - 1)] += diff / dx;
+}
+
+fn upwind4_diss(prev: ArrayView1<f32>, mut fut: ArrayViewMut1<f32>) {
+    assert_eq!(prev.shape(), fut.shape());
+    let nx = prev.shape()[0];
+
+    let dx = 1.0 / (nx - 1) as f32;
+
+    let diag = [
+        1.0 / 24.0,
+        -1.0 / 4.0,
+        5.0 / 8.0,
+        -5.0 / 6.0,
+        5.0 / 8.0,
+        -1.0 / 4.0,
+        1.0 / 24.0,
+    ];
+
+    let diff = diag[0] * prev[(nx - 3)]
+        + diag[1] * prev[(nx - 2)]
+        + diag[2] * prev[(nx - 1)]
+        + diag[3] * prev[(0)]
+        + diag[4] * prev[(1)]
+        + diag[5] * prev[(2)]
+        + diag[6] * prev[(3)];
+    fut[0] += diff / dx;
+    let diff = diag[0] * prev[(nx - 2)]
+        + diag[1] * prev[(nx - 1)]
+        + diag[2] * prev[(0)]
+        + diag[3] * prev[(1)]
+        + diag[4] * prev[(2)]
+        + diag[5] * prev[(3)]
+        + diag[6] * prev[(4)];
+    fut[1] += diff / dx;
+    let diff = diag[0] * prev[(nx - 1)]
+        + diag[1] * prev[(0)]
+        + diag[2] * prev[(1)]
+        + diag[3] * prev[(2)]
+        + diag[4] * prev[(3)]
+        + diag[5] * prev[(4)]
+        + diag[6] * prev[(5)];
+    fut[2] += diff / dx;
+
+    for i in 3..nx - 3 {
+        let diff = diag[0] * prev[(i - 3)]
+            + diag[1] * prev[(i - 2)]
+            + diag[2] * prev[(i - 1)]
+            + diag[3] * prev[(i)]
+            + diag[4] * prev[(i + 1)]
+            + diag[5] * prev[(i + 2)]
+            + diag[6] * prev[(i + 3)];
+        fut[(i)] += diff / dx;
+    }
+    let diff = diag[0] * prev[(nx - 6)]
+        + diag[1] * prev[(nx - 5)]
+        + diag[2] * prev[(nx - 4)]
+        + diag[3] * prev[(nx - 3)]
+        + diag[4] * prev[(nx - 2)]
+        + diag[5] * prev[(nx - 1)]
+        + diag[6] * prev[(0)];
+    fut[(nx - 3)] += diff / dx;
+    let diff = diag[0] * prev[(nx - 5)]
+        + diag[1] * prev[(nx - 4)]
+        + diag[2] * prev[(nx - 3)]
+        + diag[3] * prev[(nx - 2)]
+        + diag[4] * prev[(nx - 1)]
+        + diag[5] * prev[(0)]
+        + diag[6] * prev[(1)];
+    fut[(nx - 2)] += diff / dx;
+    let diff = diag[0] * prev[(nx - 4)]
+        + diag[1] * prev[(nx - 3)]
+        + diag[2] * prev[(nx - 2)]
+        + diag[3] * prev[(nx - 1)]
+        + diag[4] * prev[(0)]
+        + diag[5] * prev[(1)]
+        + diag[6] * prev[(2)];
+    fut[(nx - 1)] += diff / dx;
+}