Apply SAT on boundaries

euler/src/lib.rs

@@ -1002,104 +1002,139 @@ fn vortexify(
 
 #[allow(non_snake_case)]
 /// Boundary conditions (SAT)
-fn SAT_characteristics(
+pub fn SAT_characteristics(
     op: &dyn SbpOperator2d,
     k: &mut Diff,
     y: &Field,
     metrics: &Metrics,
     boundaries: &BoundaryTerms,
+) {
+    SAT_north(op, k, y, metrics, boundaries.north);
+    SAT_south(op, k, y, metrics, boundaries.south);
+    SAT_east(op, k, y, metrics, boundaries.east);
+    SAT_west(op, k, y, metrics, boundaries.west);
+}
+
+#[allow(non_snake_case)]
+pub fn SAT_north(
+    op: &dyn SbpOperator2d,
+    k: &mut Diff,
+    y: &Field,
+    metrics: &Metrics,
+    boundary: ArrayView2<Float>,
 ) {
     let ny = y.ny();
+
+    let hi = if op.is_h2eta() {
+        (ny - 2) as Float / op.heta()[0]
+    } else {
+        (ny - 1) as Float / op.heta()[0]
+    };
+    let sign = -1.0;
+    let tau = 1.0;
+    let slice = s![y.ny() - 1, ..];
+    SAT_characteristic(
+        k.north_mut(),
+        y.north(),
+        boundary,
+        hi,
+        sign,
+        tau,
+        metrics.detj().slice(slice),
+        metrics.detj_deta_dx().slice(slice),
+        metrics.detj_deta_dy().slice(slice),
+    );
+}
+
+#[allow(non_snake_case)]
+pub fn SAT_south(
+    op: &dyn SbpOperator2d,
+    k: &mut Diff,
+    y: &Field,
+    metrics: &Metrics,
+    boundary: ArrayView2<Float>,
+) {
+    let ny = y.ny();
+    let hi = if op.is_h2eta() {
+        (ny - 2) as Float / op.heta()[0]
+    } else {
+        (ny - 1) as Float / op.heta()[0]
+    };
+    let sign = 1.0;
+    let tau = -1.0;
+    let slice = s![0, ..];
+    SAT_characteristic(
+        k.south_mut(),
+        y.south(),
+        boundary,
+        hi,
+        sign,
+        tau,
+        metrics.detj().slice(slice),
+        metrics.detj_deta_dx().slice(slice),
+        metrics.detj_deta_dy().slice(slice),
+    );
+}
+
+#[allow(non_snake_case)]
+pub fn SAT_west(
+    op: &dyn SbpOperator2d,
+    k: &mut Diff,
+    y: &Field,
+    metrics: &Metrics,
+    boundary: ArrayView2<Float>,
+) {
     let nx = y.nx();
 
-    // North boundary
+    let hi = if op.is_h2xi() {
-    {
+        (nx - 2) as Float / op.hxi()[0]
-        let hi = if op.is_h2eta() {
+    } else {
-            (ny - 2) as Float / op.heta()[0]
+        (nx - 1) as Float / op.hxi()[0]
-        } else {
+    };
-            (ny - 1) as Float / op.heta()[0]
+    let sign = 1.0;
-        };
+    let tau = -1.0;
-        let sign = -1.0;
+    let slice = s![.., 0];
-        let tau = 1.0;
+    SAT_characteristic(
-        let slice = s![y.ny() - 1, ..];
+        k.west_mut(),
-        SAT_characteristic(
+        y.west(),
-            k.north_mut(),
+        boundary,
-            y.north(),
+        hi,
-            boundaries.north,
+        sign,
-            hi,
+        tau,
-            sign,
+        metrics.detj().slice(slice),
-            tau,
+        metrics.detj_dxi_dx().slice(slice),
-            metrics.detj().slice(slice),
+        metrics.detj_dxi_dy().slice(slice),
-            metrics.detj_deta_dx().slice(slice),
+    );
-            metrics.detj_deta_dy().slice(slice),
+}
-        );
+
-    }
+#[allow(non_snake_case)]
-    // South boundary
+pub fn SAT_east(
-    {
+    op: &dyn SbpOperator2d,
-        let hi = if op.is_h2eta() {
+    k: &mut Diff,
-            (ny - 2) as Float / op.heta()[0]
+    y: &Field,
-        } else {
+    metrics: &Metrics,
-            (ny - 1) as Float / op.heta()[0]
+    boundary: ArrayView2<Float>,
-        };
+) {
-        let sign = 1.0;
+    let nx = y.nx();
-        let tau = -1.0;
+    let hi = if op.is_h2xi() {
-        let slice = s![0, ..];
+        (nx - 2) as Float / op.hxi()[0]
-        SAT_characteristic(
+    } else {
-            k.south_mut(),
+        (nx - 1) as Float / op.hxi()[0]
-            y.south(),
+    };
-            boundaries.south,
+    let sign = -1.0;
-            hi,
+    let tau = 1.0;
-            sign,
+    let slice = s![.., y.nx() - 1];
-            tau,
+    SAT_characteristic(
-            metrics.detj().slice(slice),
+        k.east_mut(),
-            metrics.detj_deta_dx().slice(slice),
+        y.east(),
-            metrics.detj_deta_dy().slice(slice),
+        boundary,
-        );
+        hi,
-    }
+        sign,
-    // West Boundary
+        tau,
-    {
+        metrics.detj().slice(slice),
-        let hi = if op.is_h2xi() {
+        metrics.detj_dxi_dx().slice(slice),
-            (nx - 2) as Float / op.hxi()[0]
+        metrics.detj_dxi_dy().slice(slice),
-        } else {
+    );
-            (nx - 1) as Float / op.hxi()[0]
-        };
-        let sign = 1.0;
-        let tau = -1.0;
-        let slice = s![.., 0];
-        SAT_characteristic(
-            k.west_mut(),
-            y.west(),
-            boundaries.west,
-            hi,
-            sign,
-            tau,
-            metrics.detj().slice(slice),
-            metrics.detj_dxi_dx().slice(slice),
-            metrics.detj_dxi_dy().slice(slice),
-        );
-    }
-    // East Boundary
-    {
-        let hi = if op.is_h2xi() {
-            (nx - 2) as Float / op.hxi()[0]
-        } else {
-            (nx - 1) as Float / op.hxi()[0]
-        };
-        let sign = -1.0;
-        let tau = 1.0;
-        let slice = s![.., y.nx() - 1];
-        SAT_characteristic(
-            k.east_mut(),
-            y.east(),
-            boundaries.east,
-            hi,
-            sign,
-            tau,
-            metrics.detj().slice(slice),
-            metrics.detj_dxi_dx().slice(slice),
-            metrics.detj_dxi_dy().slice(slice),
-        );
-    }
 }
 
 #[allow(non_snake_case)]

multigrid/src/main.rs

@@ -42,6 +42,9 @@ struct CliOptions {
     /// in json format
     #[argh(switch)]
     output_json: bool,
+    /// distribute the computation on multiple threads
+    #[argh(switch)]
+    distribute: bool,
 }
 
 #[derive(Default, serde::Serialize)]
@@ -71,7 +74,6 @@ fn main() {
             return;
         }
     };
-
     let parsing::RuntimeConfiguration {
         names,
         grids,
@@ -109,9 +111,19 @@ fn main() {
 
     let ntime = (integration_time / dt).round() as u64;
 
-    {
+    if opt.distribute {
-        let nthreads = opt.jobs.unwrap_or(1);
+        let sys = sys.distribute(ntime);
-        todo!("nthreads");
+        let timer = if opt.timings {
+            Some(std::time::Instant::now())
+        } else {
+            None
+        };
+        sys.run();
+        if let Some(timer) = timer {
+            let duration = timer.elapsed();
+            println!("Duration: {:?}", duration);
+        }
+        return;
     }
 
     let should_output = |itime| {

multigrid/src/parsing.rs

@@ -168,8 +168,7 @@ impl input::Configuration {
         let grid_connections = self
             .grids
             .iter()
-            .enumerate()
+            .map(|(name, g)| {
-            .map(|(i, (name, g))| {
                 let default_bc = default.boundary_conditions.clone().unwrap_or_default();
                 use input::BoundaryType;
                 g.boundary_conditions
@@ -190,7 +189,7 @@ impl input::Configuration {
                                 name
                             ),
                         },
-                        Some(BoundaryType::This) => euler::BoundaryCharacteristic::Grid(i),
+                        Some(BoundaryType::This) => euler::BoundaryCharacteristic::This,
                         Some(BoundaryType::Vortex) => euler::BoundaryCharacteristic::Vortex(
                             if let BoundaryConditions::Vortex(vortex) = &boundary_conditions {
                                 vortex.clone()

multigrid/src/system.rs

@@ -179,7 +179,7 @@ impl System {
 
     /// Spreads the computation over n threads, in a thread per grid way.
     /// This system can only be called once for ntime calls.
-    pub fn distribute(self, ntime: usize) -> DistributedSystem {
+    pub fn distribute(self, ntime: u64) -> DistributedSystem {
         let nthreads = self.grids.len();
         let time = 0.0;
         // alt: crossbeam::WaitGroup
@@ -199,28 +199,40 @@ impl System {
             | euler::BoundaryCharacteristic::Interpolate(i, _) = wb.north()
             {
                 let (s, r) = crossbeam_channel::bounded(1);
-                pull_channels[*i].south_mut().replace(r).unwrap();
+                pull_channels[*i]
+                    .south_mut()
+                    .replace(r)
+                    .and_then::<(), _>(|_| panic!("channel is already present"));
                 *local_push.north_mut() = Some(s);
             }
             if let euler::BoundaryCharacteristic::Grid(i)
             | euler::BoundaryCharacteristic::Interpolate(i, _) = wb.south()
             {
                 let (s, r) = crossbeam_channel::bounded(1);
-                pull_channels[*i].north_mut().replace(r).unwrap();
+                pull_channels[*i]
+                    .north_mut()
+                    .replace(r)
+                    .and_then::<(), _>(|_| panic!("channel is already present"));
                 *local_push.south_mut() = Some(s);
             }
             if let euler::BoundaryCharacteristic::Grid(i)
             | euler::BoundaryCharacteristic::Interpolate(i, _) = wb.east()
             {
                 let (s, r) = crossbeam_channel::bounded(1);
-                pull_channels[*i].west_mut().replace(r).unwrap();
+                pull_channels[*i]
+                    .west_mut()
+                    .replace(r)
+                    .and_then::<(), _>(|_| panic!("channel is already present"));
                 *local_push.east_mut() = Some(s);
             }
             if let euler::BoundaryCharacteristic::Grid(i)
             | euler::BoundaryCharacteristic::Interpolate(i, _) = wb.west()
             {
                 let (s, r) = crossbeam_channel::bounded(1);
-                pull_channels[*i].east_mut().replace(r).unwrap();
+                pull_channels[*i]
+                    .east_mut()
+                    .replace(r)
+                    .and_then::<(), _>(|_| panic!("channel is already present"));
                 *local_push.west_mut() = Some(s);
             }
 
@@ -279,7 +291,7 @@ impl System {
                             ],
                             boundary_conditions,
                             grid: (grid, metrics),
-                            output: (),
+                            _output: (),
                             push,
                             sbp,
                             t: time,
@@ -294,7 +306,7 @@ impl System {
         // Spawn a new communicator
 
         DistributedSystem {
-            ntime,
+            _ntime: ntime,
             start: b,
             sys: tids,
         }
@@ -335,13 +347,13 @@ pub struct DistributedSystem {
     /// collect, initialise, return something to main)
     /// This simply waits until all threads are ready, then starts the computation on all threads
     start: Arc<Barrier>,
-    ntime: usize,
+    _ntime: u64,
     /// These should be joined to mark the end of the computation
     sys: Vec<std::thread::JoinHandle<()>>,
 }
 
 impl DistributedSystem {
-    fn run(self) {
+    pub fn run(self) {
         // This should start as we have n thread, but barrier blocks on n+1
         self.start.wait();
         self.sys.into_iter().for_each(|tid| tid.join().unwrap());
@@ -383,16 +395,17 @@ struct DistributedSystemPart {
     k: [Diff; 4],
     wb: WorkBuffers,
     barrier: Arc<Barrier>,
-    output: (), // hdf5::Dataset eventually,
+    _output: (), // hdf5::Dataset eventually,
     t: Float,
     dt: Float,
-    ntime: usize,
+    ntime: u64,
 }
 
 impl DistributedSystemPart {
     fn advance(&mut self) {
         self.barrier.wait();
-        for i in 0..self.ntime {
+        for _i in 0..self.ntime {
+            println!("step: {}", _i);
             let metrics = &self.grid.1;
             let wb = &mut self.wb.0;
             let sbp = &self.sbp;
@@ -400,7 +413,7 @@ impl DistributedSystemPart {
             let boundary_conditions = &self.boundary_conditions;
             let grid = &self.grid.0;
 
-            let mut rhs = |k: &mut euler::Diff, y: &euler::Field, time: Float| {
+            let rhs = |k: &mut euler::Diff, y: &euler::Field, time: Float| {
                 // Send off the boundaries optimistically, in case some grid is ready
                 if let Some(s) = &push.north {
                     s.send(y.north().to_owned()).unwrap()
@@ -419,10 +432,6 @@ impl DistributedSystemPart {
                 // This computation does not depend on the boundaries
                 euler::RHS_no_SAT(sbp.deref(), k, y, metrics, wb);
 
-                fn north_sat() {
-                    todo!()
-                }
-
                 // Get boundaries, but be careful and maximise the amount of work which can be
                 // performed before we have all of them, whilst ensuring threads can sleep for as
                 // long as possible
@@ -430,9 +439,14 @@ impl DistributedSystemPart {
                 let mut selectable = 0;
                 let recv_north = match boundary_conditions.north() {
                     DistributedBoundaryConditions::Channel(r)
-                    | DistributedBoundaryConditions::Interpolate(r, _) => Some(r),
+                    | DistributedBoundaryConditions::Interpolate(r, _) => {
+                        selectable += 1;
+                        Some(select.recv(r))
+                    }
                     DistributedBoundaryConditions::This => {
-                        todo!()
+                        let data = y.south();
+                        euler::SAT_north(sbp.deref(), k, y, metrics, data.view());
+                        None
                     }
                     DistributedBoundaryConditions::Vortex(vp) => {
                         let mut data = y.north().to_owned();
@@ -443,37 +457,149 @@ impl DistributedSystemPart {
                             fiter.next().unwrap(),
                             fiter.next().unwrap(),
                         );
-                        let (x, y) = grid.north();
+                        let (gx, gy) = grid.north();
-                        vp.evaluate(time, x, y, rho, rhou, rhov, e);
+                        vp.evaluate(time, gx, gy, rho, rhou, rhov, e);
 
-                        north_sat();
+                        euler::SAT_north(sbp.deref(), k, y, metrics, data.view());
                         None
                     }
                     DistributedBoundaryConditions::Eval(eval) => {
-                        todo!()
+                        let mut data = y.north().to_owned();
+                        let mut fiter = data.outer_iter_mut();
+                        let (rho, rhou, rhov, e) = (
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                        );
+                        let (gx, gy) = grid.north();
+                        eval.evaluate(time, gx, gy, rho, rhou, rhov, e);
+                        euler::SAT_north(sbp.deref(), k, y, metrics, data.view());
+                        None
                     }
-                    _ => None,
                 };
-                let recv_south = if let Some(r) = boundary_conditions.south().channel() {
+                let recv_south = match boundary_conditions.south() {
-                    selectable += 1;
+                    DistributedBoundaryConditions::Channel(r)
-                    Some(select.recv(r))
+                    | DistributedBoundaryConditions::Interpolate(r, _) => {
-                } else {
+                        selectable += 1;
-                    // Do SAT boundary from other BC
+                        Some(select.recv(r))
-                    None
+                    }
+                    DistributedBoundaryConditions::This => {
+                        let data = y.north();
+                        euler::SAT_south(sbp.deref(), k, y, metrics, data.view());
+                        None
+                    }
+                    DistributedBoundaryConditions::Vortex(vp) => {
+                        let mut data = y.south().to_owned();
+                        let mut fiter = data.outer_iter_mut();
+                        let (rho, rhou, rhov, e) = (
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                        );
+                        let (gx, gy) = grid.south();
+                        vp.evaluate(time, gx, gy, rho, rhou, rhov, e);
+
+                        euler::SAT_south(sbp.deref(), k, y, metrics, data.view());
+                        None
+                    }
+                    DistributedBoundaryConditions::Eval(eval) => {
+                        let mut data = y.south().to_owned();
+                        let mut fiter = data.outer_iter_mut();
+                        let (rho, rhou, rhov, e) = (
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                        );
+                        let (gx, gy) = grid.south();
+                        eval.evaluate(time, gx, gy, rho, rhou, rhov, e);
+                        euler::SAT_south(sbp.deref(), k, y, metrics, data.view());
+                        None
+                    }
                 };
-                let recv_west = if let Some(r) = boundary_conditions.west().channel() {
+                let recv_east = match boundary_conditions.east() {
-                    selectable += 1;
+                    DistributedBoundaryConditions::Channel(r)
-                    Some(select.recv(r))
+                    | DistributedBoundaryConditions::Interpolate(r, _) => {
-                } else {
+                        selectable += 1;
-                    // Do SAT boundary from other BC
+                        Some(select.recv(r))
-                    None
+                    }
+                    DistributedBoundaryConditions::This => {
+                        let data = y.west();
+                        euler::SAT_east(sbp.deref(), k, y, metrics, data.view());
+                        None
+                    }
+                    DistributedBoundaryConditions::Vortex(vp) => {
+                        let mut data = y.east().to_owned();
+                        let mut fiter = data.outer_iter_mut();
+                        let (rho, rhou, rhov, e) = (
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                        );
+                        let (gx, gy) = grid.east();
+                        vp.evaluate(time, gx, gy, rho, rhou, rhov, e);
+
+                        euler::SAT_east(sbp.deref(), k, y, metrics, data.view());
+                        None
+                    }
+                    DistributedBoundaryConditions::Eval(eval) => {
+                        let mut data = y.east().to_owned();
+                        let mut fiter = data.outer_iter_mut();
+                        let (rho, rhou, rhov, e) = (
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                        );
+                        let (gx, gy) = grid.east();
+                        eval.evaluate(time, gx, gy, rho, rhou, rhov, e);
+                        euler::SAT_east(sbp.deref(), k, y, metrics, data.view());
+                        None
+                    }
                 };
-                let recv_east = if let Some(r) = boundary_conditions.east().channel() {
+                let recv_west = match boundary_conditions.west() {
-                    selectable += 1;
+                    DistributedBoundaryConditions::Channel(r)
-                    Some(select.recv(r))
+                    | DistributedBoundaryConditions::Interpolate(r, _) => {
-                } else {
+                        selectable += 1;
-                    // Do SAT boundary from other BC
+                        Some(select.recv(r))
-                    None
+                    }
+                    DistributedBoundaryConditions::This => {
+                        let data = y.east();
+                        euler::SAT_west(sbp.deref(), k, y, metrics, data.view());
+                        None
+                    }
+                    DistributedBoundaryConditions::Vortex(vp) => {
+                        let mut data = y.west().to_owned();
+                        let mut fiter = data.outer_iter_mut();
+                        let (rho, rhou, rhov, e) = (
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                        );
+                        let (gx, gy) = grid.west();
+                        vp.evaluate(time, gx, gy, rho, rhou, rhov, e);
+
+                        euler::SAT_west(sbp.deref(), k, y, metrics, data.view());
+                        None
+                    }
+                    DistributedBoundaryConditions::Eval(eval) => {
+                        let mut data = y.west().to_owned();
+                        let mut fiter = data.outer_iter_mut();
+                        let (rho, rhou, rhov, e) = (
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                            fiter.next().unwrap(),
+                        );
+                        let (gx, gy) = grid.west();
+                        eval.evaluate(time, gx, gy, rho, rhou, rhov, e);
+                        euler::SAT_west(sbp.deref(), k, y, metrics, data.view());
+                        None
+                    }
                 };
 
                 // Get an item off each channel, waiting minimally before processing that boundary.
@@ -484,21 +610,29 @@ impl DistributedSystemPart {
                     let s = select.select();
                     let sindex = s.index();
                     match Some(sindex) {
-                        recv_north => {
+                        x if x == recv_north => {
-                            let r = s.recv(boundary_conditions.north().channel().unwrap());
+                            let r = s
-                            // process into boundary SAT here
+                                .recv(boundary_conditions.north().channel().unwrap())
+                                .unwrap();
+                            euler::SAT_north(sbp.deref(), k, y, metrics, r.view());
                         }
-                        recv_south => {
+                        x if x == recv_south => {
-                            let r = s.recv(boundary_conditions.south().channel().unwrap());
+                            let r = s
-                            // process into boundary SAT here
+                                .recv(boundary_conditions.south().channel().unwrap())
+                                .unwrap();
+                            euler::SAT_south(sbp.deref(), k, y, metrics, r.view());
                         }
-                        recv_west => {
+                        x if x == recv_west => {
-                            let r = s.recv(boundary_conditions.west().channel().unwrap());
+                            let r = s
-                            // process into boundary SAT here
+                                .recv(boundary_conditions.west().channel().unwrap())
+                                .unwrap();
+                            euler::SAT_west(sbp.deref(), k, y, metrics, r.view());
                         }
-                        recv_east => {
+                        x if x == recv_east => {
-                            let r = s.recv(boundary_conditions.east().channel().unwrap());
+                            let r = s
-                            // process into boundary SAT here
+                                .recv(boundary_conditions.east().channel().unwrap())
+                                .unwrap();
+                            euler::SAT_east(sbp.deref(), k, y, metrics, r.view());
                         }
                         _ => unreachable!(),
                     }