SAT in x direction

src/maxwell.rs

@@ -98,32 +98,32 @@ impl System {
             diffy_periodic(y.0.view(), k[i].1.view_mut());
 
             // Boundary conditions (SAT)
-            let ny = y.1.shape()[0];
+            let ny = y.0.shape()[0];
-            let nx = y.1.shape()[1];
+            let nx = y.0.shape()[1];
-            let h = 0.3; // Get from schema
+            let h = 49.0 / 144.0 / (nx - 1) as f32; // TODO: Get from schema
-            let hinv = 1.0 / h / (nx - 1) as f32;
+            let hinv = 1.0 / h;
 
             // East boundary
             for j in 0..ny {
-                let tau = (0.0, -1.0, -1.0);
+                let tau = -1.0;
                 let g = (y.0[(j, 0)], y.1[(j, 0)], y.2[(j, 0)]);
                 let v = (y.0[(j, nx - 1)], y.1[(j, nx - 1)], y.2[(j, nx - 1)]);
 
                 // A+ = (0, 0, 0; 0, 1/2, -1/2; 0, -1/2, 1/2);
                 k[i].0[(j, nx - 1)] += 0.0;
-                k[i].1[(j, nx - 1)] += tau.1 * hinv * (0.5 * (v.1 - g.1) - 0.5 * (v.2 - g.2));
+                k[i].1[(j, nx - 1)] += tau * hinv * (0.5 * (v.1 - g.1) - 0.5 * (v.2 - g.2));
-                k[i].2[(j, nx - 1)] += tau.2 * hinv * (-0.5 * (v.1 - g.1) + 0.5 * (v.2 - g.2));
+                k[i].2[(j, nx - 1)] += tau * hinv * (-0.5 * (v.1 - g.1) + 0.5 * (v.2 - g.2));
             }
             // West boundary
             for j in 0..ny {
-                let tau = (0.0, 1.0, -1.0);
+                let tau = 1.0;
                 let g = (y.0[(j, nx - 1)], y.1[(j, nx - 1)], y.2[(j, nx - 1)]);
                 let v = (y.0[(j, 0)], y.1[(j, 0)], y.2[(j, 0)]);
 
                 // A- = (0, 0, 0; 0, -1/2, -1/2; 0, -1/2, -1/2);
                 k[i].0[(j, 0)] += 0.0;
-                k[i].1[(j, 0)] += tau.1 * hinv * (-0.5 * (v.1 - g.1) - 0.5 * (v.2 - g.2));
+                k[i].1[(j, 0)] += tau * hinv * (-0.5 * (v.1 - g.1) - 0.5 * (v.2 - g.2));
-                k[i].2[(j, 0)] += tau.2 * hinv * (-0.5 * (v.1 - g.1) - 0.5 * (v.2 - g.2));
+                k[i].2[(j, 0)] += tau * hinv * (-0.5 * (v.1 - g.1) - 0.5 * (v.2 - g.2));
             }
         }
 

src/operators.rs

@@ -132,25 +132,24 @@ fn upwind4(prev: ArrayView1<f32>, mut fut: ArrayViewMut1<f32>) {
         ],
     ]);
 
-    let h_block = [49.0 / 144.0, 61.0 / 48.0, 41.0 / 48.0, 149.0 / 144.0];
+    // let h_block = [49.0 / 144.0, 61.0 / 48.0, 41.0 / 48.0, 149.0 / 144.0];
 
     let first_elems = prev.slice(s!(..7));
     for i in 0..4 {
         let diff = first_elems.dot(&block.slice(s!(i, ..)));
-        fut[i] += diff / dx * h_block[i];
+        fut[i] += diff / dx;
     }
 
     for i in 4..nx - 4 {
         let diff = diag.dot(&prev.slice(s!(i - 3..i + 3 + 1)));
         fut[(i)] += diff / dx;
     }
-
     let last_elems = prev.slice(s!(nx - 7..));
     for i in 0..4 {
         let ii = nx - 4 + i;
         let block = block.slice(s!(3 - i, ..;-1));
         let diff = last_elems.dot(&block);
-        fut[ii] += diff / dx * h_block[3 - i];
+        fut[ii] += -diff / dx;
     }
 }
 
@@ -231,6 +230,41 @@ fn upwind4_periodic(prev: ArrayView1<f32>, mut fut: ArrayViewMut1<f32>) {
     fut[(nx - 1)] += diff / dx;
 }
 
+#[test]
+fn upwind4_test() {
+    let nx = 20;
+    let dx = 1.0 / (nx - 1) as f32;
+    let mut source: ndarray::Array1<f32> = ndarray::Array1::zeros((nx));
+    let mut res = ndarray::Array1::zeros((nx));
+    let mut target = ndarray::Array1::zeros((nx));
+
+    for i in 0..nx {
+        source[i] = i as f32 * dx;
+        target[i] = 1.0;
+    }
+    res.fill(0.0);
+    upwind4(source.view(), res.view_mut());
+    assert!(res.all_close(&target, 1e-4));
+
+    for i in 0..nx {
+        let x = i as f32 * dx;
+        source[i] = x * x;
+        target[i] = 2.0 * x;
+    }
+    res.fill(0.0);
+    upwind4(source.view(), res.view_mut());
+    assert!(res.all_close(&target, 1e-4));
+
+    for i in 0..nx {
+        let x = i as f32 * dx;
+        source[i] = x * x * x;
+        target[i] = 3.0 * x * x;
+    }
+    res.fill(0.0);
+    upwind4(source.view(), res.view_mut());
+    assert!(res.all_close(&target, 1e-2));
+}
+
 fn upwind4_diss(prev: ArrayView1<f32>, mut fut: ArrayViewMut1<f32>) {
     assert_eq!(prev.shape(), fut.shape());
     let nx = prev.shape()[0];