From 95897777d662cb1381aed0e7459cc4aa7be95008 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Wed, 4 Aug 2021 18:30:15 +0000 Subject: [PATCH] Apply SAT on boundaries --- euler/src/lib.rs | 213 ++++++++++++++++++++-------------- multigrid/src/main.rs | 20 +++- multigrid/src/parsing.rs | 5 +- multigrid/src/system.rs | 242 ++++++++++++++++++++++++++++++--------- 4 files changed, 330 insertions(+), 150 deletions(-) diff --git a/euler/src/lib.rs b/euler/src/lib.rs index 8586bf6..350cb5f 100644 --- a/euler/src/lib.rs +++ b/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, ) { 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, +) { + 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, +) { let nx = y.nx(); - // North boundary - { - 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(), - boundaries.north, - hi, - sign, - tau, - metrics.detj().slice(slice), - metrics.detj_deta_dx().slice(slice), - metrics.detj_deta_dy().slice(slice), - ); - } - // South boundary - { - 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(), - boundaries.south, - hi, - sign, - tau, - metrics.detj().slice(slice), - metrics.detj_deta_dx().slice(slice), - metrics.detj_deta_dy().slice(slice), - ); - } - // West 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![.., 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), - ); - } + 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![.., 0]; + SAT_characteristic( + k.west_mut(), + y.west(), + boundary, + hi, + sign, + tau, + metrics.detj().slice(slice), + metrics.detj_dxi_dx().slice(slice), + metrics.detj_dxi_dy().slice(slice), + ); +} + +#[allow(non_snake_case)] +pub fn SAT_east( + op: &dyn SbpOperator2d, + k: &mut Diff, + y: &Field, + metrics: &Metrics, + boundary: ArrayView2, +) { + let nx = y.nx(); + 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(), + boundary, + hi, + sign, + tau, + metrics.detj().slice(slice), + metrics.detj_dxi_dx().slice(slice), + metrics.detj_dxi_dy().slice(slice), + ); } #[allow(non_snake_case)] diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index d2231a8..252a4fc 100644 --- a/multigrid/src/main.rs +++ b/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; - { - let nthreads = opt.jobs.unwrap_or(1); - todo!("nthreads"); + if opt.distribute { + let sys = sys.distribute(ntime); + 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| { diff --git a/multigrid/src/parsing.rs b/multigrid/src/parsing.rs index f0d5c14..5d87324 100644 --- a/multigrid/src/parsing.rs +++ b/multigrid/src/parsing.rs @@ -168,8 +168,7 @@ impl input::Configuration { let grid_connections = self .grids .iter() - .enumerate() - .map(|(i, (name, g))| { + .map(|(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() diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index d2a348f..27ca5b4 100644 --- a/multigrid/src/system.rs +++ b/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, - ntime: usize, + _ntime: u64, /// These should be joined to mark the end of the computation sys: Vec>, } 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, - 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(); - vp.evaluate(time, x, y, rho, rhou, rhov, e); + let (gx, gy) = grid.north(); + 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() { - selectable += 1; - Some(select.recv(r)) - } else { - // Do SAT boundary from other BC - None + let recv_south = match boundary_conditions.south() { + DistributedBoundaryConditions::Channel(r) + | DistributedBoundaryConditions::Interpolate(r, _) => { + selectable += 1; + Some(select.recv(r)) + } + 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() { - selectable += 1; - Some(select.recv(r)) - } else { - // Do SAT boundary from other BC - None + let recv_east = match boundary_conditions.east() { + DistributedBoundaryConditions::Channel(r) + | DistributedBoundaryConditions::Interpolate(r, _) => { + selectable += 1; + Some(select.recv(r)) + } + 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() { - selectable += 1; - Some(select.recv(r)) - } else { - // Do SAT boundary from other BC - None + let recv_west = match boundary_conditions.west() { + DistributedBoundaryConditions::Channel(r) + | DistributedBoundaryConditions::Interpolate(r, _) => { + selectable += 1; + Some(select.recv(r)) + } + 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 => { - let r = s.recv(boundary_conditions.north().channel().unwrap()); - // process into boundary SAT here + x if x == recv_north => { + let r = s + .recv(boundary_conditions.north().channel().unwrap()) + .unwrap(); + euler::SAT_north(sbp.deref(), k, y, metrics, r.view()); } - recv_south => { - let r = s.recv(boundary_conditions.south().channel().unwrap()); - // process into boundary SAT here + x if x == recv_south => { + let r = s + .recv(boundary_conditions.south().channel().unwrap()) + .unwrap(); + euler::SAT_south(sbp.deref(), k, y, metrics, r.view()); } - recv_west => { - let r = s.recv(boundary_conditions.west().channel().unwrap()); - // process into boundary SAT here + x if x == recv_west => { + let r = s + .recv(boundary_conditions.west().channel().unwrap()) + .unwrap(); + euler::SAT_west(sbp.deref(), k, y, metrics, r.view()); } - recv_east => { - let r = s.recv(boundary_conditions.east().channel().unwrap()); - // process into boundary SAT here + x if x == recv_east => { + let r = s + .recv(boundary_conditions.east().channel().unwrap()) + .unwrap(); + euler::SAT_east(sbp.deref(), k, y, metrics, r.view()); } _ => unreachable!(), }