From 35b8af8b2dd364b0b77ea2d97d0e6f4a7a79bf32 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Tue, 3 Aug 2021 16:18:20 +0200 Subject: [PATCH 01/25] Partial implementation --- euler/src/lib.rs | 53 ++++ multigrid/Cargo.toml | 1 + multigrid/src/main.rs | 213 +-------------- multigrid/src/system.rs | 579 ++++++++++++++++++++++++++++++++++++++++ 4 files changed, 636 insertions(+), 210 deletions(-) create mode 100644 multigrid/src/system.rs diff --git a/euler/src/lib.rs b/euler/src/lib.rs index 32e299c..8586bf6 100644 --- a/euler/src/lib.rs +++ b/euler/src/lib.rs @@ -517,6 +517,59 @@ pub fn RHS_trad( SAT_characteristics(op, k, y, metrics, boundaries); } +#[allow(non_snake_case)] +pub fn RHS_no_SAT( + op: &dyn SbpOperator2d, + k: &mut Diff, + y: &Field, + metrics: &Metrics, + tmp: &mut (Field, Field, Field, Field, Field, Field), +) { + let ehat = &mut tmp.0; + let fhat = &mut tmp.1; + fluxes((ehat, fhat), y, metrics, &mut tmp.2); + let dE = &mut tmp.2; + let dF = &mut tmp.3; + + op.diffxi(ehat.rho(), dE.rho_mut()); + op.diffxi(ehat.rhou(), dE.rhou_mut()); + op.diffxi(ehat.rhov(), dE.rhov_mut()); + op.diffxi(ehat.e(), dE.e_mut()); + + op.diffeta(fhat.rho(), dF.rho_mut()); + op.diffeta(fhat.rhou(), dF.rhou_mut()); + op.diffeta(fhat.rhov(), dF.rhov_mut()); + op.diffeta(fhat.e(), dF.e_mut()); + + if let Some(diss_op) = op.upwind() { + let ad_xi = &mut tmp.4; + let ad_eta = &mut tmp.5; + upwind_dissipation( + &*diss_op, + (ad_xi, ad_eta), + y, + metrics, + (&mut tmp.0, &mut tmp.1), + ); + + azip!((out in &mut k.0, + eflux in &dE.0, + fflux in &dF.0, + ad_xi in &ad_xi.0, + ad_eta in &ad_eta.0, + detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) { + *out = (-eflux - fflux + ad_xi + ad_eta)/detj + }); + } else { + azip!((out in &mut k.0, + eflux in &dE.0, + fflux in &dF.0, + detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) { + *out = (-eflux - fflux )/detj + }); + } +} + #[allow(non_snake_case)] pub fn RHS_upwind( op: &dyn SbpOperator2d, diff --git a/multigrid/Cargo.toml b/multigrid/Cargo.toml index 91490c4..52ae7c9 100644 --- a/multigrid/Cargo.toml +++ b/multigrid/Cargo.toml @@ -18,3 +18,4 @@ json5 = "0.3.0" indexmap = { version = "1.5.2", features = ["serde-1"] } argh = "0.1.4" evalexpr = "6.3.0" +crossbeam-channel = "0.5.0" diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index c3b3894..d2231a8 100644 --- a/multigrid/src/main.rs +++ b/multigrid/src/main.rs @@ -1,8 +1,6 @@ use argh::FromArgs; -use rayon::prelude::*; use euler::eval::Evaluator; -use sbp::operators::SbpOperator2d; use sbp::*; mod file; @@ -10,207 +8,8 @@ mod input; mod parsing; use file::*; mod eval; - -struct System { - fnow: Vec, - fnext: Vec, - wb: Vec, - k: [Vec; 4], - grids: Vec, - metrics: Vec, - bt: Vec, - eb: Vec, - time: Float, - operators: Vec>, -} - -use std::sync::atomic::{AtomicBool, Ordering}; -pub(crate) static MULTITHREAD: AtomicBool = AtomicBool::new(false); - -impl integrate::Integrable for System { - type State = Vec; - type Diff = Vec; - - fn scaled_add(s: &mut Self::State, o: &Self::Diff, scale: Float) { - if MULTITHREAD.load(Ordering::Acquire) { - s.par_iter_mut() - .zip(o.par_iter()) - .for_each(|(s, o)| euler::Field::scaled_add(s, o, scale)) - } else { - s.iter_mut() - .zip(o.iter()) - .for_each(|(s, o)| euler::Field::scaled_add(s, o, scale)) - } - } -} - -impl System { - fn new( - grids: Vec, - bt: Vec, - operators: Vec>, - ) -> Self { - let fnow = grids - .iter() - .map(|g| euler::Field::new(g.ny(), g.nx())) - .collect::>(); - let fnext = fnow.clone(); - let wb = grids - .iter() - .map(|g| euler::WorkBuffers::new(g.ny(), g.nx())) - .collect(); - let k = grids - .iter() - .map(|g| euler::Diff::zeros((g.ny(), g.nx()))) - .collect::>(); - let k = [k.clone(), k.clone(), k.clone(), k]; - let metrics = grids - .iter() - .zip(&operators) - .map(|(g, op)| g.metrics(&**op).unwrap()) - .collect::>(); - - let eb = bt - .iter() - .zip(&grids) - .map(|(bt, grid)| euler::BoundaryStorage::new(bt, grid)) - .collect(); - - Self { - fnow, - fnext, - k, - wb, - grids, - metrics, - bt, - eb, - time: 0.0, - operators, - } - } - - fn vortex(&mut self, t: Float, vortex_params: &euler::VortexParameters) { - for (f, g) in self.fnow.iter_mut().zip(&self.grids) { - f.vortex(g.x(), g.y(), t, &vortex_params); - } - } - - fn advance(&mut self, dt: Float) { - let metrics = &self.metrics; - let grids = &self.grids; - let bt = &self.bt; - let wb = &mut self.wb; - let eb = &mut self.eb; - let operators = &self.operators; - - let rhs = move |fut: &mut Vec, prev: &Vec, time: Float| { - let prev_all = &prev; - if MULTITHREAD.load(Ordering::Acquire) { - rayon::scope(|s| { - for (((((((fut, prev), wb), grid), metrics), op), bt), eb) in fut - .iter_mut() - .zip(prev.iter()) - .zip(wb.iter_mut()) - .zip(grids) - .zip(metrics.iter()) - .zip(operators.iter()) - .zip(bt.iter()) - .zip(eb.iter_mut()) - { - s.spawn(move |_| { - let bc = euler::boundary_extracts(prev_all, bt, prev, grid, eb, time); - if op.upwind().is_some() { - euler::RHS_upwind(&**op, fut, prev, metrics, &bc, &mut wb.0); - } else { - euler::RHS_trad(&**op, fut, prev, metrics, &bc, &mut wb.0); - } - }) - } - }); - } else { - for (((((((fut, prev), wb), grid), metrics), op), bt), eb) in fut - .iter_mut() - .zip(prev.iter()) - .zip(wb.iter_mut()) - .zip(grids) - .zip(metrics.iter()) - .zip(operators.iter()) - .zip(bt.iter()) - .zip(eb.iter_mut()) - { - let bc = euler::boundary_extracts(prev_all, bt, prev, grid, eb, time); - if op.upwind().is_some() { - euler::RHS_upwind(&**op, fut, prev, metrics, &bc, &mut wb.0); - } else { - euler::RHS_trad(&**op, fut, prev, metrics, &bc, &mut wb.0); - } - } - } - }; - - integrate::integrate::( - rhs, - &self.fnow, - &mut self.fnext, - &mut self.time, - dt, - &mut self.k, - ); - - std::mem::swap(&mut self.fnow, &mut self.fnext); - } - - /// Suggested maximum dt for this problem - fn max_dt(&self) -> Float { - let is_h2 = self - .operators - .iter() - .any(|op| op.is_h2xi() || op.is_h2eta()); - let c_max = if is_h2 { 0.5 } else { 1.0 }; - let mut max_dt: Float = Float::INFINITY; - - for (field, metrics) in self.fnow.iter().zip(self.metrics.iter()) { - let nx = field.nx(); - let ny = field.ny(); - - let rho = field.rho(); - let rhou = field.rhou(); - let rhov = field.rhov(); - - let mut max_u: Float = 0.0; - let mut max_v: Float = 0.0; - - for ((((((rho, rhou), rhov), detj_dxi_dx), detj_dxi_dy), detj_deta_dx), detj_deta_dy) in - rho.iter() - .zip(rhou.iter()) - .zip(rhov.iter()) - .zip(metrics.detj_dxi_dx()) - .zip(metrics.detj_dxi_dy()) - .zip(metrics.detj_deta_dx()) - .zip(metrics.detj_deta_dy()) - { - let u = rhou / rho; - let v = rhov / rho; - - let uhat: Float = detj_dxi_dx * u + detj_dxi_dy * v; - let vhat: Float = detj_deta_dx * u + detj_deta_dy * v; - - max_u = max_u.max(uhat.abs()); - max_v = max_v.max(vhat.abs()); - } - - let dx = 1.0 / nx as Float; - let dy = 1.0 / ny as Float; - - let c_dt = Float::max(max_u / dx, max_v / dy); - - max_dt = Float::min(max_dt, c_max / c_dt); - } - - max_dt - } -} +mod system; +use system::*; #[derive(Debug, FromArgs)] /// Options for configuring and running the solver @@ -312,13 +111,7 @@ fn main() { { let nthreads = opt.jobs.unwrap_or(1); - if nthreads > 1 { - MULTITHREAD.store(true, Ordering::Release); - rayon::ThreadPoolBuilder::new() - .num_threads(nthreads) - .build_global() - .unwrap(); - } + todo!("nthreads"); } let should_output = |itime| { diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs new file mode 100644 index 0000000..380517c --- /dev/null +++ b/multigrid/src/system.rs @@ -0,0 +1,579 @@ +use crate::utils::Direction; +use crossbeam_channel::{Receiver, Select, Sender}; +use euler::{ + eval::{self, Evaluator}, + Diff, Field, VortexParameters, WorkBuffers, +}; +use ndarray::Array2; +use sbp::grid::{Grid, Metrics}; +use sbp::operators::{InterpolationOperator, SbpOperator2d}; +use sbp::*; +use std::sync::{Arc, Barrier}; + +pub struct System { + pub fnow: Vec, + pub fnext: Vec, + pub wb: Vec, + pub k: [Vec; 4], + pub grids: Vec, + pub metrics: Vec, + pub bt: Vec, + pub eb: Vec, + pub time: Float, + pub operators: Vec>, +} + +impl integrate::Integrable for System { + type State = Vec; + type Diff = Vec; + + fn scaled_add(s: &mut Self::State, o: &Self::Diff, scale: Float) { + s.iter_mut() + .zip(o.iter()) + .for_each(|(s, o)| euler::Field::scaled_add(s, o, scale)) + } +} + +impl System { + pub fn new( + grids: Vec, + bt: Vec, + operators: Vec>, + ) -> Self { + let fnow = grids + .iter() + .map(|g| euler::Field::new(g.ny(), g.nx())) + .collect::>(); + let fnext = fnow.clone(); + let wb = grids + .iter() + .map(|g| euler::WorkBuffers::new(g.ny(), g.nx())) + .collect(); + let k = grids + .iter() + .map(|g| euler::Diff::zeros((g.ny(), g.nx()))) + .collect::>(); + let k = [k.clone(), k.clone(), k.clone(), k]; + let metrics = grids + .iter() + .zip(&operators) + .map(|(g, op)| g.metrics(&**op).unwrap()) + .collect::>(); + + let eb = bt + .iter() + .zip(&grids) + .map(|(bt, grid)| euler::BoundaryStorage::new(bt, grid)) + .collect(); + + Self { + fnow, + fnext, + k, + wb, + grids, + metrics, + bt, + eb, + time: 0.0, + operators, + } + } + + pub fn vortex(&mut self, t: Float, vortex_params: &euler::VortexParameters) { + for (f, g) in self.fnow.iter_mut().zip(&self.grids) { + f.vortex(g.x(), g.y(), t, &vortex_params); + } + } + + pub fn advance(&mut self, dt: Float) { + let metrics = &self.metrics; + let grids = &self.grids; + let bt = &self.bt; + let wb = &mut self.wb; + let eb = &mut self.eb; + let operators = &self.operators; + + let rhs = move |fut: &mut Vec, prev: &Vec, time: Float| { + let prev_all = &prev; + for (((((((fut, prev), wb), grid), metrics), op), bt), eb) in fut + .iter_mut() + .zip(prev.iter()) + .zip(wb.iter_mut()) + .zip(grids) + .zip(metrics.iter()) + .zip(operators.iter()) + .zip(bt.iter()) + .zip(eb.iter_mut()) + { + let bc = euler::boundary_extracts(prev_all, bt, prev, grid, eb, time); + if op.upwind().is_some() { + euler::RHS_upwind(&**op, fut, prev, metrics, &bc, &mut wb.0); + } else { + euler::RHS_trad(&**op, fut, prev, metrics, &bc, &mut wb.0); + } + } + }; + + integrate::integrate::( + rhs, + &self.fnow, + &mut self.fnext, + &mut self.time, + dt, + &mut self.k, + ); + + std::mem::swap(&mut self.fnow, &mut self.fnext); + } + + /// Suggested maximum dt for this problem + pub fn max_dt(&self) -> Float { + let is_h2 = self + .operators + .iter() + .any(|op| op.is_h2xi() || op.is_h2eta()); + let c_max = if is_h2 { 0.5 } else { 1.0 }; + let mut max_dt: Float = Float::INFINITY; + + for (field, metrics) in self.fnow.iter().zip(self.metrics.iter()) { + let nx = field.nx(); + let ny = field.ny(); + + let rho = field.rho(); + let rhou = field.rhou(); + let rhov = field.rhov(); + + let mut max_u: Float = 0.0; + let mut max_v: Float = 0.0; + + for ((((((rho, rhou), rhov), detj_dxi_dx), detj_dxi_dy), detj_deta_dx), detj_deta_dy) in + rho.iter() + .zip(rhou.iter()) + .zip(rhov.iter()) + .zip(metrics.detj_dxi_dx()) + .zip(metrics.detj_dxi_dy()) + .zip(metrics.detj_deta_dx()) + .zip(metrics.detj_deta_dy()) + { + let u = rhou / rho; + let v = rhov / rho; + + let uhat: Float = detj_dxi_dx * u + detj_dxi_dy * v; + let vhat: Float = detj_deta_dx * u + detj_deta_dy * v; + + max_u = max_u.max(uhat.abs()); + max_v = max_v.max(vhat.abs()); + } + + let dx = 1.0 / nx as Float; + let dy = 1.0 / ny as Float; + + let c_dt = Float::max(max_u / dx, max_v / dy); + + max_dt = Float::min(max_dt, c_max / c_dt); + } + + max_dt + } + + /// 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 { + let nthreads = self.grids.len(); + let time = 0.0; + // alt: crossbeam::WaitGroup + let b = Arc::new(Barrier::new(nthreads + 1)); + let dt = self.max_dt(); + + // Build up the boundary conditions + // Assume all boundaries are push/pull through channels + let channels = (0..nthreads) + .map(|_| { + use crossbeam_channel::unbounded; + Direction { + north: unbounded(), + south: unbounded(), + west: unbounded(), + east: unbounded(), + } + }) + .collect::>(); + + // TODO: Iterate through all grids and see if they need ourself to push + let mut requested_channels = (0..nthreads) + .map(|_| Direction::default()) + .collect::>>(); + + let mut tids = Vec::new(); + for (((((((current, fut), grid), metrics), sbp), wb), bt), req_channel) in self + .fnow + .into_iter() + .zip(self.fnext.into_iter()) + .zip(self.grids.into_iter()) + .zip(self.metrics.into_iter()) + .zip(self.operators.into_iter()) + .zip(self.wb.into_iter()) + .zip(self.bt) + .zip(requested_channels) + { + let builder = std::thread::Builder::new().name(format!("eulersolver: {}", "smth")); + let barrier = b.clone(); + + let Direction { + north: bt_north, + south: bt_south, + west: bt_west, + east: bt_east, + } = bt; + let boundary_conditions = Direction { + north: match bt_north { + euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, + euler::BoundaryCharacteristic::Grid(i) => { + *requested_channels[i].south_mut() = true; + DistributedBoundaryConditions::Channel(channels[i].south().1.clone()) + } + euler::BoundaryCharacteristic::Interpolate(i, int_op) => { + *requested_channels[i].south_mut() = true; + DistributedBoundaryConditions::Interpolate( + channels[i].south().1.clone(), + int_op, + ) + } + euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), + euler::BoundaryCharacteristic::Vortex(vp) => { + DistributedBoundaryConditions::Vortex(vp) + } + euler::BoundaryCharacteristic::Eval(eval) => { + DistributedBoundaryConditions::Eval(eval) + } + }, + south: match bt_south { + euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, + euler::BoundaryCharacteristic::Grid(i) => { + *requested_channels[i].north_mut() = true; + DistributedBoundaryConditions::Channel(channels[i].north().1.clone()) + } + euler::BoundaryCharacteristic::Interpolate(i, int_op) => { + *requested_channels[i].north_mut() = true; + DistributedBoundaryConditions::Interpolate( + channels[i].north().1.clone(), + int_op, + ) + } + euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), + euler::BoundaryCharacteristic::Vortex(vp) => { + DistributedBoundaryConditions::Vortex(vp) + } + euler::BoundaryCharacteristic::Eval(eval) => { + DistributedBoundaryConditions::Eval(eval) + } + }, + east: match bt_east { + euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, + euler::BoundaryCharacteristic::Grid(i) => { + *requested_channels[i].west_mut() = true; + DistributedBoundaryConditions::Channel(channels[i].west().1.clone()) + } + euler::BoundaryCharacteristic::Interpolate(i, int_op) => { + *requested_channels[i].west_mut() = true; + DistributedBoundaryConditions::Interpolate( + channels[i].west().1.clone(), + int_op, + ) + } + euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), + euler::BoundaryCharacteristic::Vortex(vp) => { + DistributedBoundaryConditions::Vortex(vp) + } + euler::BoundaryCharacteristic::Eval(eval) => { + DistributedBoundaryConditions::Eval(eval) + } + }, + west: match bt_west { + euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, + euler::BoundaryCharacteristic::Grid(i) => { + *requested_channels[i].east_mut() = true; + DistributedBoundaryConditions::Channel(channels[i].east().1.clone()) + } + euler::BoundaryCharacteristic::Interpolate(i, int_op) => { + *requested_channels[i].east_mut() = true; + DistributedBoundaryConditions::Interpolate( + channels[i].east().1.clone(), + int_op, + ) + } + euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), + euler::BoundaryCharacteristic::Vortex(vp) => { + DistributedBoundaryConditions::Vortex(vp) + } + euler::BoundaryCharacteristic::Eval(eval) => { + DistributedBoundaryConditions::Eval(eval) + } + }, + }; + + tids.push( + builder + .spawn(move || { + let mut sys = DistributedSystemPart { + barrier, + ntime, + dt, + current, + fut, + k: [todo!(); 4], + boundary_conditions, + grid: (grid, metrics), + output: (), + push: todo!(), + sbp, + t: time, + wb, + }; + sys.advance(); + }) + .unwrap(), + ); + } + // Set up communicators + // Spawn a new communicator + + DistributedSystem { + ntime, + start: b, + sys: tids, + } + } +} + +// single-threaded items: clone_from +// sync points: scaled_add, rhs +// +// Could instead make every thread (of a threadpool?) carry one grid themselves, +// and instead only wait on obtaining bc, which would make synchronisation +// be somewhat minimal +// +// Difficulties/implementation notes: +// * How to get output from each thread? Push to the multi-producer, single consumer queue +// with capacity of 2*N (which ensures internal synchronisation stays active), or many +// channels which are select!'ed. Choose instead a (name, itime, data) tuple, which is +// communicated. Also initialise the output file all at the start, since we know the +// number of time steps, and which timesteps will be used (should_output). +// * How to balance the number of threads/jobs? +// Each grid will be pushed to a thread, and a thread pool created for use by all threads +// combined. Set affinity to limit the number of cores available (could have side effects +// on heat etc, although this would also be the case for multi-threaded code). A +// thread-per-core architecture + the output thread if we have enough cores, else let os +// control the scheduling for us. +// Use cgroups to artificially limit the number of cores. +// Use taskset to artificially limit the number of cores. +// * Mechanism to wait for bc, is select available? Use a channel with 0-4 capacity, each +// other grid additonally pushes their direction to always use the available +// Use select or similar from crossbeam_channel +// * Mechanism to push bc: try_push and then move to next bc if not available. Spawn a +// thread into a threadpool to handle this. +// * Each grid is spawned onto an async threadpool with M threads allocated to itself, this should +// allow async waiting on the channel (or in the future from process/internet/bc processor) +// * BC: A local thread can have a handle +pub struct DistributedSystem { + /// Simple messaging system to be replaced by a more sophisticated system (e.g. run 5 steps, + /// 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, + /// These should be joined to mark the end of the computation + sys: Vec>, +} + +impl DistributedSystem { + 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()); + } +} + +// #[derive(Debug)] +pub enum DistributedBoundaryConditions { + This, + + Vortex(VortexParameters), + Eval(std::sync::Arc>), + + Interpolate(Receiver>, Box), + Channel(Receiver>), +} + +impl DistributedBoundaryConditions { + fn channel(&self) -> Option<&Receiver>> { + match self { + Self::Interpolate(r, _) | Self::Channel(r) => Some(r), + _ => None, + } + } +} + +#[derive(Debug, Clone)] +enum PushCommunicator { + Channel(Sender>), + None, +} + +impl Default for PushCommunicator { + fn default() -> Self { + Self::None + } +} + +struct DistributedSystemPart { + grid: (Grid, Metrics), + sbp: Box, + + boundary_conditions: Direction, + /// Subscribers to the boundaries of self + push: Direction, + + current: Field, + fut: Field, + k: [Diff; 4], + wb: WorkBuffers, + barrier: Arc, + output: (), // hdf5::Dataset eventually, + t: Float, + dt: Float, + ntime: usize, +} + +impl DistributedSystemPart { + fn advance(&mut self) { + self.barrier.wait(); + for i in 0..self.ntime { + let metrics = &self.grid.1; + let wb = &mut self.wb.0; + let sbp = &self.sbp; + let push = &self.push; + let boundary_conditions = &self.boundary_conditions; + let grid = &self.grid.0; + + let mut rhs = |k: &mut euler::Diff, y: &euler::Field, time: Float| { + // Send off the boundaries optimistically, in case some grid is ready + match &push.north { + PushCommunicator::None => (), + PushCommunicator::Channel(s) => s.send(y.north().to_owned()).unwrap(), + } + match &push.south { + PushCommunicator::None => (), + PushCommunicator::Channel(s) => s.send(y.south().to_owned()).unwrap(), + } + match &push.east { + PushCommunicator::None => (), + PushCommunicator::Channel(s) => s.send(y.east().to_owned()).unwrap(), + } + match &push.west { + PushCommunicator::None => (), + PushCommunicator::Channel(s) => s.send(y.west().to_owned()).unwrap(), + } + + use std::ops::Deref; + // 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 + let mut select = Select::new(); + let mut selectable = 0; + let recv_north = match boundary_conditions.north() { + DistributedBoundaryConditions::Channel(r) + | DistributedBoundaryConditions::Interpolate(r, _) => Some(r), + DistributedBoundaryConditions::This => { + todo!() + } + DistributedBoundaryConditions::Vortex(vp) => { + 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 (x, y) = grid.north(); + vp.evaluate(time, x, y, rho, rhou, rhov, e); + + north_sat(); + None + } + DistributedBoundaryConditions::Eval(eval) => { + todo!() + } + _ => 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_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 = if let Some(r) = boundary_conditions.east().channel() { + selectable += 1; + Some(select.recv(r)) + } else { + // Do SAT boundary from other BC + None + }; + + // Get an item off each channel, waiting minimally before processing that boundary. + // The waiting ensures other grids can be processed by the core in case of + // oversubscription (in case of a more grids than core scenario) + // This minimises the amount of time waiting on boundary conditions + while selectable != 0 { + 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 + } + recv_south => { + let r = s.recv(boundary_conditions.south().channel().unwrap()); + // process into boundary SAT here + } + recv_west => { + let r = s.recv(boundary_conditions.west().channel().unwrap()); + // process into boundary SAT here + } + recv_east => { + let r = s.recv(boundary_conditions.east().channel().unwrap()); + // process into boundary SAT here + } + _ => unreachable!(), + } + select.remove(sindex); + selectable -= 1; + } + }; + integrate::integrate::( + rhs, + &self.current, + &mut self.fut, + &mut self.t, + self.dt, + &mut self.k, + ) + } + } +} From b11f3c9abb2b884942e136f7bc0374cf9b3d0469 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Tue, 3 Aug 2021 17:17:39 +0000 Subject: [PATCH 02/25] Improve channel distribution --- multigrid/src/system.rs | 210 ++++++++++++++-------------------------- 1 file changed, 75 insertions(+), 135 deletions(-) diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 380517c..d2a348f 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -187,26 +187,48 @@ impl System { let dt = self.max_dt(); // Build up the boundary conditions - // Assume all boundaries are push/pull through channels - let channels = (0..nthreads) - .map(|_| { - use crossbeam_channel::unbounded; - Direction { - north: unbounded(), - south: unbounded(), - west: unbounded(), - east: unbounded(), - } - }) - .collect::>(); + let mut push_channels: Vec>>>> = + Vec::with_capacity(nthreads); + let mut pull_channels: Vec>>>> = + vec![Direction::default(); nthreads]; - // TODO: Iterate through all grids and see if they need ourself to push - let mut requested_channels = (0..nthreads) - .map(|_| Direction::default()) - .collect::>>(); + // Build the set of communicators between boundaries + for wb in &self.bt { + let mut local_push = Direction::default(); + if let euler::BoundaryCharacteristic::Grid(i) + | euler::BoundaryCharacteristic::Interpolate(i, _) = wb.north() + { + let (s, r) = crossbeam_channel::bounded(1); + pull_channels[*i].south_mut().replace(r).unwrap(); + *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(); + *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(); + *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(); + *local_push.west_mut() = Some(s); + } + + push_channels.push(local_push); + } let mut tids = Vec::new(); - for (((((((current, fut), grid), metrics), sbp), wb), bt), req_channel) in self + for ((((((((current, fut), grid), metrics), sbp), wb), bt), chan), push) in self .fnow .into_iter() .zip(self.fnext.into_iter()) @@ -215,103 +237,30 @@ impl System { .zip(self.operators.into_iter()) .zip(self.wb.into_iter()) .zip(self.bt) - .zip(requested_channels) + .zip(pull_channels) + .zip(push_channels) { let builder = std::thread::Builder::new().name(format!("eulersolver: {}", "smth")); let barrier = b.clone(); - let Direction { - north: bt_north, - south: bt_south, - west: bt_west, - east: bt_east, - } = bt; - let boundary_conditions = Direction { - north: match bt_north { - euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, - euler::BoundaryCharacteristic::Grid(i) => { - *requested_channels[i].south_mut() = true; - DistributedBoundaryConditions::Channel(channels[i].south().1.clone()) - } - euler::BoundaryCharacteristic::Interpolate(i, int_op) => { - *requested_channels[i].south_mut() = true; - DistributedBoundaryConditions::Interpolate( - channels[i].south().1.clone(), - int_op, - ) - } - euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), - euler::BoundaryCharacteristic::Vortex(vp) => { - DistributedBoundaryConditions::Vortex(vp) - } - euler::BoundaryCharacteristic::Eval(eval) => { - DistributedBoundaryConditions::Eval(eval) - } - }, - south: match bt_south { - euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, - euler::BoundaryCharacteristic::Grid(i) => { - *requested_channels[i].north_mut() = true; - DistributedBoundaryConditions::Channel(channels[i].north().1.clone()) - } - euler::BoundaryCharacteristic::Interpolate(i, int_op) => { - *requested_channels[i].north_mut() = true; - DistributedBoundaryConditions::Interpolate( - channels[i].north().1.clone(), - int_op, - ) - } - euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), - euler::BoundaryCharacteristic::Vortex(vp) => { - DistributedBoundaryConditions::Vortex(vp) - } - euler::BoundaryCharacteristic::Eval(eval) => { - DistributedBoundaryConditions::Eval(eval) - } - }, - east: match bt_east { - euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, - euler::BoundaryCharacteristic::Grid(i) => { - *requested_channels[i].west_mut() = true; - DistributedBoundaryConditions::Channel(channels[i].west().1.clone()) - } - euler::BoundaryCharacteristic::Interpolate(i, int_op) => { - *requested_channels[i].west_mut() = true; - DistributedBoundaryConditions::Interpolate( - channels[i].west().1.clone(), - int_op, - ) - } - euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), - euler::BoundaryCharacteristic::Vortex(vp) => { - DistributedBoundaryConditions::Vortex(vp) - } - euler::BoundaryCharacteristic::Eval(eval) => { - DistributedBoundaryConditions::Eval(eval) - } - }, - west: match bt_west { - euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, - euler::BoundaryCharacteristic::Grid(i) => { - *requested_channels[i].east_mut() = true; - DistributedBoundaryConditions::Channel(channels[i].east().1.clone()) - } - euler::BoundaryCharacteristic::Interpolate(i, int_op) => { - *requested_channels[i].east_mut() = true; - DistributedBoundaryConditions::Interpolate( - channels[i].east().1.clone(), - int_op, - ) - } - euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), - euler::BoundaryCharacteristic::Vortex(vp) => { - DistributedBoundaryConditions::Vortex(vp) - } - euler::BoundaryCharacteristic::Eval(eval) => { - DistributedBoundaryConditions::Eval(eval) - } - }, - }; + let boundary_conditions = bt.zip(chan).map(|(bt, chan)| match bt { + euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, + euler::BoundaryCharacteristic::Grid(_) => { + DistributedBoundaryConditions::Channel(chan.unwrap()) + } + euler::BoundaryCharacteristic::Interpolate(_, int_op) => { + DistributedBoundaryConditions::Interpolate(chan.unwrap(), int_op) + } + euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), + euler::BoundaryCharacteristic::Vortex(vp) => { + DistributedBoundaryConditions::Vortex(vp) + } + euler::BoundaryCharacteristic::Eval(eval) => { + DistributedBoundaryConditions::Eval(eval) + } + }); + + let (ny, nx) = (grid.nx(), grid.ny()); tids.push( builder @@ -322,11 +271,16 @@ impl System { dt, current, fut, - k: [todo!(); 4], + k: [ + Diff::zeros((ny, nx)), + Diff::zeros((ny, nx)), + Diff::zeros((ny, nx)), + Diff::zeros((ny, nx)), + ], boundary_conditions, grid: (grid, metrics), output: (), - push: todo!(), + push, sbp, t: time, wb, @@ -414,17 +368,7 @@ impl DistributedBoundaryConditions { } } -#[derive(Debug, Clone)] -enum PushCommunicator { - Channel(Sender>), - None, -} - -impl Default for PushCommunicator { - fn default() -> Self { - Self::None - } -} +type PushCommunicator = Option>>; struct DistributedSystemPart { grid: (Grid, Metrics), @@ -458,21 +402,17 @@ impl DistributedSystemPart { let mut rhs = |k: &mut euler::Diff, y: &euler::Field, time: Float| { // Send off the boundaries optimistically, in case some grid is ready - match &push.north { - PushCommunicator::None => (), - PushCommunicator::Channel(s) => s.send(y.north().to_owned()).unwrap(), + if let Some(s) = &push.north { + s.send(y.north().to_owned()).unwrap() } - match &push.south { - PushCommunicator::None => (), - PushCommunicator::Channel(s) => s.send(y.south().to_owned()).unwrap(), + if let Some(s) = &push.south { + s.send(y.south().to_owned()).unwrap() } - match &push.east { - PushCommunicator::None => (), - PushCommunicator::Channel(s) => s.send(y.east().to_owned()).unwrap(), + if let Some(s) = &push.east { + s.send(y.east().to_owned()).unwrap() } - match &push.west { - PushCommunicator::None => (), - PushCommunicator::Channel(s) => s.send(y.west().to_owned()).unwrap(), + if let Some(s) = &push.west { + s.send(y.west().to_owned()).unwrap() } use std::ops::Deref; From 95897777d662cb1381aed0e7459cc4aa7be95008 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Wed, 4 Aug 2021 18:30:15 +0000 Subject: [PATCH 03/25] 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!(), } From 26159d5ffb9966bbff3461321af49f7a81d7e0bf Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Thu, 5 Aug 2021 19:28:14 +0000 Subject: [PATCH 04/25] Introduce precursor system --- multigrid/src/main.rs | 36 ++++++++------------------ multigrid/src/parsing.rs | 10 +++----- multigrid/src/system.rs | 55 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 70 insertions(+), 31 deletions(-) diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index 252a4fc..6507f99 100644 --- a/multigrid/src/main.rs +++ b/multigrid/src/main.rs @@ -9,7 +9,6 @@ mod parsing; use file::*; mod eval; mod system; -use system::*; #[derive(Debug, FromArgs)] /// Options for configuring and running the solver @@ -77,35 +76,22 @@ fn main() { let parsing::RuntimeConfiguration { names, grids, - grid_connections, + boundary_conditions, op: operators, integration_time, initial_conditions, - boundary_conditions: _, } = config.into_runtime(); - let mut sys = System::new(grids, grid_connections, operators); - match &initial_conditions { - /* - parsing::InitialConditions::File(f) => { - for grid in &sys.grids { - // Copy initial conditions from file, requires name of field - todo!() - } - } - */ - parsing::InitialConditions::Vortex(vortexparams) => sys.vortex(0.0, &vortexparams), - parsing::InitialConditions::Expressions(expr) => { - let t = 0.0; - for (grid, field) in sys.grids.iter().zip(sys.fnow.iter_mut()) { - // Evaluate the expressions on all variables - let x = grid.x(); - let y = grid.y(); - let (rho, rhou, rhov, e) = field.components_mut(); - (*expr).evaluate(t, x, y, rho, rhou, rhov, e); - } - } - } + let basesystem = system::BaseSystem::new( + names.clone(), + grids, + 0.0, + operators, + boundary_conditions, + initial_conditions.clone(), + ); + let mut sys = basesystem.create(); + // System::new(grids, grid_connections, operators); let dt = sys.max_dt(); diff --git a/multigrid/src/parsing.rs b/multigrid/src/parsing.rs index 5d87324..035ff70 100644 --- a/multigrid/src/parsing.rs +++ b/multigrid/src/parsing.rs @@ -27,7 +27,7 @@ pub enum InitialConditions { } #[derive(Clone, Debug)] -pub enum BoundaryConditions { +enum BoundaryConditions { Vortex(euler::VortexParameters), Expressions(std::sync::Arc), NotNeeded, @@ -36,11 +36,10 @@ pub enum BoundaryConditions { pub struct RuntimeConfiguration { pub names: Vec, pub grids: Vec, - pub grid_connections: Vec, + pub boundary_conditions: Vec, pub op: Vec>, pub integration_time: Float, pub initial_conditions: InitialConditions, - pub boundary_conditions: BoundaryConditions, } impl input::Configuration { @@ -165,7 +164,7 @@ impl input::Configuration { Box::new((matcher(eta), matcher(xi))) as Box }) .collect(); - let grid_connections = self + let boundary_conditions = self .grids .iter() .map(|(name, g)| { @@ -226,11 +225,10 @@ impl input::Configuration { RuntimeConfiguration { names, grids, - grid_connections, + boundary_conditions, op, integration_time: self.integration_time, initial_conditions, - boundary_conditions, } } } diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 27ca5b4..034a7d8 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -1,3 +1,4 @@ +use crate::parsing; use crate::utils::Direction; use crossbeam_channel::{Receiver, Select, Sender}; use euler::{ @@ -10,6 +11,60 @@ use sbp::operators::{InterpolationOperator, SbpOperator2d}; use sbp::*; use std::sync::{Arc, Barrier}; +pub struct BaseSystem { + pub names: Vec, + pub grids: Vec, + pub time: Float, + pub boundary_conditions: Vec, + pub initial_conditions: crate::parsing::InitialConditions, + pub operators: Vec>, +} + +impl BaseSystem { + pub fn new( + names: Vec, + grids: Vec, + time: Float, + operators: Vec>, + boundary_conditions: Vec, + initial_conditions: crate::parsing::InitialConditions, + ) -> Self { + Self { + names, + grids, + time, + boundary_conditions, + initial_conditions, + operators, + } + } + pub fn create(self) -> System { + let mut sys = System::new(self.grids, self.boundary_conditions, self.operators); + match &self.initial_conditions { + /* + parsing::InitialConditions::File(f) => { + for grid in &sys.grids { + // Copy initial conditions from file, requires name of field + todo!() + } + } + */ + parsing::InitialConditions::Vortex(vortexparams) => sys.vortex(0.0, &vortexparams), + parsing::InitialConditions::Expressions(expr) => { + let t = 0.0; + for (grid, field) in sys.grids.iter().zip(sys.fnow.iter_mut()) { + // Evaluate the expressions on all variables + let x = grid.x(); + let y = grid.y(); + let (rho, rhou, rhov, e) = field.components_mut(); + (*expr).evaluate(t, x, y, rho, rhou, rhov, e); + } + } + } + sys + } +} + pub struct System { pub fnow: Vec, pub fnext: Vec, From b142bb63e499679ea5070f96887273c126154614 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Fri, 6 Aug 2021 15:01:48 +0000 Subject: [PATCH 05/25] Add workbuffer for boundaries --- multigrid/src/system.rs | 66 ++++++++++++++++++++--------------------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 034a7d8..8e1c2b7 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -350,7 +350,10 @@ impl System { push, sbp, t: time, + wb, + wb_ns: Array2::zeros((4, nx)), + wb_ew: Array2::zeros((4, ny)), }; sys.advance(); }) @@ -454,19 +457,24 @@ struct DistributedSystemPart { t: Float, dt: Float, ntime: u64, + + /// Work buffer for boundaries + wb_ns: Array2, + wb_ew: Array2, } impl DistributedSystemPart { fn advance(&mut self) { self.barrier.wait(); for _i in 0..self.ntime { - println!("step: {}", _i); let metrics = &self.grid.1; let wb = &mut self.wb.0; let sbp = &self.sbp; let push = &self.push; let boundary_conditions = &self.boundary_conditions; let grid = &self.grid.0; + let wb_ns = &mut self.wb_ns; + let wb_ew = &mut self.wb_ew; let rhs = |k: &mut euler::Diff, y: &euler::Field, time: Float| { // Send off the boundaries optimistically, in case some grid is ready @@ -499,13 +507,11 @@ impl DistributedSystemPart { Some(select.recv(r)) } DistributedBoundaryConditions::This => { - let data = y.south(); - euler::SAT_north(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_north(sbp.deref(), k, y, metrics, y.south()); None } DistributedBoundaryConditions::Vortex(vp) => { - let mut data = y.north().to_owned(); - let mut fiter = data.outer_iter_mut(); + let mut fiter = wb_ns.outer_iter_mut(); let (rho, rhou, rhov, e) = ( fiter.next().unwrap(), fiter.next().unwrap(), @@ -515,12 +521,11 @@ impl DistributedSystemPart { let (gx, gy) = grid.north(); vp.evaluate(time, gx, gy, rho, rhou, rhov, e); - euler::SAT_north(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_north(sbp.deref(), k, y, metrics, wb_ns.view()); None } DistributedBoundaryConditions::Eval(eval) => { - let mut data = y.north().to_owned(); - let mut fiter = data.outer_iter_mut(); + let mut fiter = wb_ns.outer_iter_mut(); let (rho, rhou, rhov, e) = ( fiter.next().unwrap(), fiter.next().unwrap(), @@ -529,7 +534,7 @@ impl DistributedSystemPart { ); let (gx, gy) = grid.north(); eval.evaluate(time, gx, gy, rho, rhou, rhov, e); - euler::SAT_north(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_north(sbp.deref(), k, y, metrics, wb_ns.view()); None } }; @@ -540,13 +545,11 @@ impl DistributedSystemPart { Some(select.recv(r)) } DistributedBoundaryConditions::This => { - let data = y.north(); - euler::SAT_south(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_south(sbp.deref(), k, y, metrics, y.north()); None } DistributedBoundaryConditions::Vortex(vp) => { - let mut data = y.south().to_owned(); - let mut fiter = data.outer_iter_mut(); + let mut fiter = wb_ns.outer_iter_mut(); let (rho, rhou, rhov, e) = ( fiter.next().unwrap(), fiter.next().unwrap(), @@ -556,12 +559,11 @@ impl DistributedSystemPart { let (gx, gy) = grid.south(); vp.evaluate(time, gx, gy, rho, rhou, rhov, e); - euler::SAT_south(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_south(sbp.deref(), k, y, metrics, wb_ns.view()); None } DistributedBoundaryConditions::Eval(eval) => { - let mut data = y.south().to_owned(); - let mut fiter = data.outer_iter_mut(); + let mut fiter = wb_ns.outer_iter_mut(); let (rho, rhou, rhov, e) = ( fiter.next().unwrap(), fiter.next().unwrap(), @@ -570,7 +572,7 @@ impl DistributedSystemPart { ); let (gx, gy) = grid.south(); eval.evaluate(time, gx, gy, rho, rhou, rhov, e); - euler::SAT_south(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_south(sbp.deref(), k, y, metrics, wb_ns.view()); None } }; @@ -581,13 +583,11 @@ impl DistributedSystemPart { Some(select.recv(r)) } DistributedBoundaryConditions::This => { - let data = y.west(); - euler::SAT_east(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_east(sbp.deref(), k, y, metrics, y.west()); None } DistributedBoundaryConditions::Vortex(vp) => { - let mut data = y.east().to_owned(); - let mut fiter = data.outer_iter_mut(); + let mut fiter = wb_ew.outer_iter_mut(); let (rho, rhou, rhov, e) = ( fiter.next().unwrap(), fiter.next().unwrap(), @@ -597,12 +597,11 @@ impl DistributedSystemPart { let (gx, gy) = grid.east(); vp.evaluate(time, gx, gy, rho, rhou, rhov, e); - euler::SAT_east(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_east(sbp.deref(), k, y, metrics, wb_ew.view()); None } DistributedBoundaryConditions::Eval(eval) => { - let mut data = y.east().to_owned(); - let mut fiter = data.outer_iter_mut(); + let mut fiter = wb_ew.outer_iter_mut(); let (rho, rhou, rhov, e) = ( fiter.next().unwrap(), fiter.next().unwrap(), @@ -611,7 +610,7 @@ impl DistributedSystemPart { ); let (gx, gy) = grid.east(); eval.evaluate(time, gx, gy, rho, rhou, rhov, e); - euler::SAT_east(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_east(sbp.deref(), k, y, metrics, wb_ew.view()); None } }; @@ -622,13 +621,11 @@ impl DistributedSystemPart { Some(select.recv(r)) } DistributedBoundaryConditions::This => { - let data = y.east(); - euler::SAT_west(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_west(sbp.deref(), k, y, metrics, y.east()); None } DistributedBoundaryConditions::Vortex(vp) => { - let mut data = y.west().to_owned(); - let mut fiter = data.outer_iter_mut(); + let mut fiter = wb_ew.outer_iter_mut(); let (rho, rhou, rhov, e) = ( fiter.next().unwrap(), fiter.next().unwrap(), @@ -638,12 +635,11 @@ impl DistributedSystemPart { let (gx, gy) = grid.west(); vp.evaluate(time, gx, gy, rho, rhou, rhov, e); - euler::SAT_west(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_west(sbp.deref(), k, y, metrics, wb_ew.view()); None } DistributedBoundaryConditions::Eval(eval) => { - let mut data = y.west().to_owned(); - let mut fiter = data.outer_iter_mut(); + let mut fiter = wb_ew.outer_iter_mut(); let (rho, rhou, rhov, e) = ( fiter.next().unwrap(), fiter.next().unwrap(), @@ -652,7 +648,7 @@ impl DistributedSystemPart { ); let (gx, gy) = grid.west(); eval.evaluate(time, gx, gy, rho, rhou, rhov, e); - euler::SAT_west(sbp.deref(), k, y, metrics, data.view()); + euler::SAT_west(sbp.deref(), k, y, metrics, wb_ew.view()); None } }; @@ -669,24 +665,28 @@ impl DistributedSystemPart { let r = s .recv(boundary_conditions.north().channel().unwrap()) .unwrap(); + // TODO: Interpolation euler::SAT_north(sbp.deref(), k, y, metrics, r.view()); } x if x == recv_south => { let r = s .recv(boundary_conditions.south().channel().unwrap()) .unwrap(); + // TODO: Interpolation euler::SAT_south(sbp.deref(), k, y, metrics, r.view()); } x if x == recv_west => { let r = s .recv(boundary_conditions.west().channel().unwrap()) .unwrap(); + // TODO: Interpolation euler::SAT_west(sbp.deref(), k, y, metrics, r.view()); } x if x == recv_east => { let r = s .recv(boundary_conditions.east().channel().unwrap()) .unwrap(); + // TODO: Interpolation euler::SAT_east(sbp.deref(), k, y, metrics, r.view()); } _ => unreachable!(), From 1e363c05087b5bdd193ab59c406c63b47e519345 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Fri, 6 Aug 2021 17:24:15 +0000 Subject: [PATCH 06/25] Reimplement interpolation for distributed system --- multigrid/src/system.rs | 115 +++++++++++++++++++++++++++------------- 1 file changed, 77 insertions(+), 38 deletions(-) diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 8e1c2b7..1b778c6 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -327,7 +327,7 @@ impl System { } }); - let (ny, nx) = (grid.nx(), grid.ny()); + let (ny, nx) = (grid.ny(), grid.nx()); tids.push( builder @@ -429,15 +429,6 @@ pub enum DistributedBoundaryConditions { Channel(Receiver>), } -impl DistributedBoundaryConditions { - fn channel(&self) -> Option<&Receiver>> { - match self { - Self::Interpolate(r, _) | Self::Channel(r) => Some(r), - _ => None, - } - } -} - type PushCommunicator = Option>>; struct DistributedSystemPart { @@ -661,34 +652,82 @@ impl DistributedSystemPart { let s = select.select(); let sindex = s.index(); match Some(sindex) { - x if x == recv_north => { - let r = s - .recv(boundary_conditions.north().channel().unwrap()) - .unwrap(); - // TODO: Interpolation - euler::SAT_north(sbp.deref(), k, y, metrics, r.view()); - } - x if x == recv_south => { - let r = s - .recv(boundary_conditions.south().channel().unwrap()) - .unwrap(); - // TODO: Interpolation - euler::SAT_south(sbp.deref(), k, y, metrics, r.view()); - } - x if x == recv_west => { - let r = s - .recv(boundary_conditions.west().channel().unwrap()) - .unwrap(); - // TODO: Interpolation - euler::SAT_west(sbp.deref(), k, y, metrics, r.view()); - } - x if x == recv_east => { - let r = s - .recv(boundary_conditions.east().channel().unwrap()) - .unwrap(); - // TODO: Interpolation - euler::SAT_east(sbp.deref(), k, y, metrics, r.view()); - } + x if x == recv_north => match boundary_conditions.north() { + DistributedBoundaryConditions::Channel(r) => { + let r = s.recv(r).unwrap(); + euler::SAT_north(sbp.deref(), k, y, metrics, r.view()); + } + DistributedBoundaryConditions::Interpolate(r, int_op) => { + let r = s.recv(r).unwrap(); + let is_fine2coarse = r.shape()[1] > wb_ns.shape()[1]; + for (mut to, from) in wb_ns.outer_iter_mut().zip(r.outer_iter()) { + if is_fine2coarse { + int_op.fine2coarse(from.view(), to.view_mut()); + } else { + int_op.coarse2fine(from.view(), to.view_mut()); + } + } + euler::SAT_north(sbp.deref(), k, y, metrics, wb_ns.view()); + } + _ => unreachable!(), + }, + x if x == recv_south => match boundary_conditions.south() { + DistributedBoundaryConditions::Channel(r) => { + let r = s.recv(r).unwrap(); + euler::SAT_south(sbp.deref(), k, y, metrics, r.view()); + } + DistributedBoundaryConditions::Interpolate(r, int_op) => { + let r = s.recv(r).unwrap(); + let is_fine2coarse = r.shape()[1] > wb_ns.shape()[1]; + for (mut to, from) in wb_ns.outer_iter_mut().zip(r.outer_iter()) { + if is_fine2coarse { + int_op.fine2coarse(from.view(), to.view_mut()); + } else { + int_op.coarse2fine(from.view(), to.view_mut()); + } + } + euler::SAT_south(sbp.deref(), k, y, metrics, wb_ns.view()); + } + _ => unreachable!(), + }, + x if x == recv_west => match boundary_conditions.west() { + DistributedBoundaryConditions::Channel(r) => { + let r = s.recv(r).unwrap(); + euler::SAT_west(sbp.deref(), k, y, metrics, r.view()); + } + DistributedBoundaryConditions::Interpolate(r, int_op) => { + let r = s.recv(r).unwrap(); + let is_fine2coarse = r.shape()[1] > wb_ew.shape()[1]; + for (mut to, from) in wb_ew.outer_iter_mut().zip(r.outer_iter()) { + if is_fine2coarse { + int_op.fine2coarse(from.view(), to.view_mut()); + } else { + int_op.coarse2fine(from.view(), to.view_mut()); + } + } + euler::SAT_west(sbp.deref(), k, y, metrics, wb_ew.view()); + } + _ => unreachable!(), + }, + x if x == recv_east => match boundary_conditions.east() { + DistributedBoundaryConditions::Channel(r) => { + let r = s.recv(r).unwrap(); + euler::SAT_east(sbp.deref(), k, y, metrics, r.view()); + } + DistributedBoundaryConditions::Interpolate(r, int_op) => { + let r = s.recv(r).unwrap(); + let is_fine2coarse = r.shape()[1] > wb_ew.shape()[1]; + for (mut to, from) in wb_ew.outer_iter_mut().zip(r.outer_iter()) { + if is_fine2coarse { + int_op.fine2coarse(from.view(), to.view_mut()); + } else { + int_op.coarse2fine(from.view(), to.view_mut()); + } + } + euler::SAT_east(sbp.deref(), k, y, metrics, wb_ew.view()); + } + _ => unreachable!(), + }, _ => unreachable!(), } select.remove(sindex); From 4319e403a5f2c8579ad06850c5fb910e35d84255 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Fri, 6 Aug 2021 17:24:59 +0000 Subject: [PATCH 07/25] Remove option for number of threads --- multigrid/src/main.rs | 3 --- 1 file changed, 3 deletions(-) diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index 6507f99..01e37e4 100644 --- a/multigrid/src/main.rs +++ b/multigrid/src/main.rs @@ -15,9 +15,6 @@ mod system; struct CliOptions { #[argh(positional)] json: std::path::PathBuf, - /// number of simultaneous threads - #[argh(option, short = 'j')] - jobs: Option, /// name of output file #[argh( option, From 0a5647df3a9cd36f5a272526dd9bbe23518eb2ad Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Sat, 7 Aug 2021 17:30:00 +0000 Subject: [PATCH 08/25] Compilable state --- multigrid/src/main.rs | 63 +++-- multigrid/src/system.rs | 528 ++++++++++++++++++++++------------------ 2 files changed, 320 insertions(+), 271 deletions(-) diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index 01e37e4..ed16796 100644 --- a/multigrid/src/main.rs +++ b/multigrid/src/main.rs @@ -86,42 +86,28 @@ fn main() { operators, boundary_conditions, initial_conditions.clone(), + opt.output.clone(), ); - let mut sys = basesystem.create(); // System::new(grids, grid_connections, operators); - let dt = sys.max_dt(); - - let ntime = (integration_time / dt).round() as u64; - - 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| { - opt.number_of_outputs.map_or(false, |num_out| { - if num_out == 0 { - false - } else { - itime % (std::cmp::max(ntime / (num_out - 1), 1)) == 0 - } - }) + let sys = if opt.distribute { + basesystem.create_distributed() + } else { + basesystem.create() }; - let output = File::create(&opt.output, sys.grids.as_slice(), names).unwrap(); - let mut output = OutputThread::new(output); - output.add_timestep(0, &sys.fnow); + let dt = sys.max_dt(); + let ntime = (integration_time / dt).round() as u64; + let steps_between_outputs = if let Some(n) = opt.number_of_outputs { + std::cmp::max(n / ntime, 1) + } else { + ntime + }; + + sys.output(0); + //let output = File::create(&opt.output, sys.grids.as_slice(), names).unwrap(); + //let mut output = OutputThread::new(output); + //output.add_timestep(0, &sys.fnow); let progressbar = progressbar(opt.no_progressbar, ntime); @@ -131,6 +117,16 @@ fn main() { None }; + let mut itime = 0; + while itime < ntime { + sys.advance(steps_between_outputs); + progressbar.inc(1); + + itime += steps_between_outputs; + sys.output(itime); + } + + /* for itime in 0..ntime { if should_output(itime) { output.add_timestep(itime, &sys.fnow); @@ -138,6 +134,7 @@ fn main() { progressbar.inc(1); sys.advance(dt); } + */ progressbar.finish_and_clear(); let mut outinfo = OutputInformation { @@ -150,8 +147,9 @@ fn main() { outinfo.time_elapsed = Some(duration); } - output.add_timestep(ntime, &sys.fnow); + //output.add_timestep(ntime, &sys.fnow); + /* if opt.error { let time = ntime as Float * dt; let mut e = 0.0; @@ -171,6 +169,7 @@ fn main() { } outinfo.error = Some(e); } + */ if opt.output_json { println!("{}", json5::to_string(&outinfo).unwrap()); diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 1b778c6..8058b4b 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -1,5 +1,6 @@ use crate::parsing; use crate::utils::Direction; +use core::ops::Deref; use crossbeam_channel::{Receiver, Select, Sender}; use euler::{ eval::{self, Evaluator}, @@ -9,7 +10,6 @@ use ndarray::Array2; use sbp::grid::{Grid, Metrics}; use sbp::operators::{InterpolationOperator, SbpOperator2d}; use sbp::*; -use std::sync::{Arc, Barrier}; pub struct BaseSystem { pub names: Vec, @@ -18,6 +18,7 @@ pub struct BaseSystem { pub boundary_conditions: Vec, pub initial_conditions: crate::parsing::InitialConditions, pub operators: Vec>, + pub output: hdf5::File, } impl BaseSystem { @@ -28,7 +29,15 @@ impl BaseSystem { operators: Vec>, boundary_conditions: Vec, initial_conditions: crate::parsing::InitialConditions, + output: std::path::PathBuf, ) -> Self { + let output = hdf5::File::create(output).unwrap(); + output + .new_dataset::() + .resizable(true) + .chunk((1,)) + .create("t", (0,)) + .unwrap(); Self { names, grids, @@ -36,10 +45,54 @@ impl BaseSystem { boundary_conditions, initial_conditions, operators, + output, } } pub fn create(self) -> System { - let mut sys = System::new(self.grids, self.boundary_conditions, self.operators); + let fnow = self + .grids + .iter() + .map(|g| euler::Field::new(g.ny(), g.nx())) + .collect::>(); + let fnext = fnow.clone(); + let wb = self + .grids + .iter() + .map(|g| euler::WorkBuffers::new(g.ny(), g.nx())) + .collect(); + let k = self + .grids + .iter() + .map(|g| euler::Diff::zeros((g.ny(), g.nx()))) + .collect::>(); + let k = [k.clone(), k.clone(), k.clone(), k]; + let metrics = self + .grids + .iter() + .zip(&self.operators) + .map(|(g, op)| g.metrics(&**op).unwrap()) + .collect::>(); + + let eb = self + .boundary_conditions + .iter() + .zip(&self.grids) + .map(|(bt, grid)| euler::BoundaryStorage::new(bt, grid)) + .collect(); + + let sys = SingleThreadedSystem { + fnow, + fnext, + k, + wb, + grids: self.grids, + metrics, + bt: self.boundary_conditions, + eb, + time: self.time, + operators: self.operators, + output: todo!(), + }; match &self.initial_conditions { /* parsing::InitialConditions::File(f) => { @@ -61,11 +114,180 @@ impl BaseSystem { } } } - sys + System::SingleThreaded(sys) + } + + /// Spreads the computation over n threads, in a thread per grid way. + pub fn create_distributed(self) -> System { + let nthreads = self.grids.len(); + // let dt = self.max_dt(); + + // Build up the boundary conditions + let mut push_channels: Vec>>>> = + Vec::with_capacity(nthreads); + let mut pull_channels: Vec>>>> = + vec![Direction::default(); nthreads]; + + // Build the set of communicators between boundaries + for wb in &self.boundary_conditions { + let mut local_push = Direction::default(); + if let euler::BoundaryCharacteristic::Grid(i) + | euler::BoundaryCharacteristic::Interpolate(i, _) = wb.north() + { + let (s, r) = crossbeam_channel::bounded(1); + 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) + .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) + .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) + .and_then::<(), _>(|_| panic!("channel is already present")); + *local_push.west_mut() = Some(s); + } + + push_channels.push(local_push); + } + + let (master_send, master_recv) = crossbeam_channel::unbounded(); + + let mut tids = Vec::with_capacity(nthreads); + let mut communicators = Vec::with_capacity(nthreads); + + for (id, (((((name, grid), sbp), bt), chan), push)) in self + .names + .into_iter() + .zip(self.grids.into_iter()) + .zip(self.operators.into_iter()) + .zip(self.boundary_conditions) + .zip(pull_channels) + .zip(push_channels) + .enumerate() + { + let builder = std::thread::Builder::new().name(format!("eulersolver: {}", name)); + + let boundary_conditions = bt.zip(chan).map(|(bt, chan)| match bt { + euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, + euler::BoundaryCharacteristic::Grid(_) => { + DistributedBoundaryConditions::Channel(chan.unwrap()) + } + euler::BoundaryCharacteristic::Interpolate(_, int_op) => { + DistributedBoundaryConditions::Interpolate(chan.unwrap(), int_op) + } + euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), + euler::BoundaryCharacteristic::Vortex(vp) => { + DistributedBoundaryConditions::Vortex(vp) + } + euler::BoundaryCharacteristic::Eval(eval) => { + DistributedBoundaryConditions::Eval(eval) + } + }); + + let master_send = master_send.clone(); + + let (t_send, t_recv) = crossbeam_channel::unbounded(); + communicators.push(t_send); + let time = self.time; + + tids.push( + builder + .spawn(move || { + let (ny, nx) = (grid.ny(), grid.nx()); + let current = Field::new(ny, nx); + let fut = current.clone(); + let k = [ + Diff::zeros((ny, nx)), + Diff::zeros((ny, nx)), + Diff::zeros((ny, nx)), + Diff::zeros((ny, nx)), + ]; + let metrics = grid.metrics(sbp.deref()).unwrap(); + + let wb = WorkBuffers::new(ny, nx); + + let mut sys = DistributedSystemPart { + current, + fut, + k, + boundary_conditions, + grid: (grid, metrics), + output: todo!(), + push, + sbp, + t: time, + dt: Float::NAN, + + name, + id, + + recv: t_recv, + send: master_send, + + wb, + wb_ns: Array2::zeros((4, nx)), + wb_ew: Array2::zeros((4, ny)), + }; + + // init and send maxdt + // receive maxdt + sys.run(); + }) + .unwrap(), + ); + } + + System::MultiThreaded(DistributedSystem { + sys: tids, + recv: master_recv, + send: communicators, + }) } } -pub struct System { +pub enum System { + SingleThreaded(SingleThreadedSystem), + MultiThreaded(DistributedSystem), +} + +impl System { + pub fn max_dt(&self) -> Float { + todo!() + } + pub fn advance(&self, nsteps: u64) { + todo!() + } + pub fn output(&self, ntime: u64) { + todo!() + } +} + +pub struct SingleThreadedSystem { pub fnow: Vec, pub fnext: Vec, pub wb: Vec, @@ -76,9 +298,10 @@ pub struct System { pub eb: Vec, pub time: Float, pub operators: Vec>, + pub output: Vec, } -impl integrate::Integrable for System { +impl integrate::Integrable for SingleThreadedSystem { type State = Vec; type Diff = Vec; @@ -89,52 +312,7 @@ impl integrate::Integrable for System { } } -impl System { - pub fn new( - grids: Vec, - bt: Vec, - operators: Vec>, - ) -> Self { - let fnow = grids - .iter() - .map(|g| euler::Field::new(g.ny(), g.nx())) - .collect::>(); - let fnext = fnow.clone(); - let wb = grids - .iter() - .map(|g| euler::WorkBuffers::new(g.ny(), g.nx())) - .collect(); - let k = grids - .iter() - .map(|g| euler::Diff::zeros((g.ny(), g.nx()))) - .collect::>(); - let k = [k.clone(), k.clone(), k.clone(), k]; - let metrics = grids - .iter() - .zip(&operators) - .map(|(g, op)| g.metrics(&**op).unwrap()) - .collect::>(); - - let eb = bt - .iter() - .zip(&grids) - .map(|(bt, grid)| euler::BoundaryStorage::new(bt, grid)) - .collect(); - - Self { - fnow, - fnext, - k, - wb, - grids, - metrics, - bt, - eb, - time: 0.0, - operators, - } - } - +impl SingleThreadedSystem { pub fn vortex(&mut self, t: Float, vortex_params: &euler::VortexParameters) { for (f, g) in self.fnow.iter_mut().zip(&self.grids) { f.vortex(g.x(), g.y(), t, &vortex_params); @@ -170,7 +348,7 @@ impl System { } }; - integrate::integrate::( + integrate::integrate::( rhs, &self.fnow, &mut self.fnext, @@ -231,191 +409,48 @@ impl System { max_dt } - - /// 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: u64) -> DistributedSystem { - let nthreads = self.grids.len(); - let time = 0.0; - // alt: crossbeam::WaitGroup - let b = Arc::new(Barrier::new(nthreads + 1)); - let dt = self.max_dt(); - - // Build up the boundary conditions - let mut push_channels: Vec>>>> = - Vec::with_capacity(nthreads); - let mut pull_channels: Vec>>>> = - vec![Direction::default(); nthreads]; - - // Build the set of communicators between boundaries - for wb in &self.bt { - let mut local_push = Direction::default(); - if let euler::BoundaryCharacteristic::Grid(i) - | euler::BoundaryCharacteristic::Interpolate(i, _) = wb.north() - { - let (s, r) = crossbeam_channel::bounded(1); - 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) - .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) - .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) - .and_then::<(), _>(|_| panic!("channel is already present")); - *local_push.west_mut() = Some(s); - } - - push_channels.push(local_push); - } - - let mut tids = Vec::new(); - for ((((((((current, fut), grid), metrics), sbp), wb), bt), chan), push) in self - .fnow - .into_iter() - .zip(self.fnext.into_iter()) - .zip(self.grids.into_iter()) - .zip(self.metrics.into_iter()) - .zip(self.operators.into_iter()) - .zip(self.wb.into_iter()) - .zip(self.bt) - .zip(pull_channels) - .zip(push_channels) - { - let builder = std::thread::Builder::new().name(format!("eulersolver: {}", "smth")); - let barrier = b.clone(); - - let boundary_conditions = bt.zip(chan).map(|(bt, chan)| match bt { - euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, - euler::BoundaryCharacteristic::Grid(_) => { - DistributedBoundaryConditions::Channel(chan.unwrap()) - } - euler::BoundaryCharacteristic::Interpolate(_, int_op) => { - DistributedBoundaryConditions::Interpolate(chan.unwrap(), int_op) - } - euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), - euler::BoundaryCharacteristic::Vortex(vp) => { - DistributedBoundaryConditions::Vortex(vp) - } - euler::BoundaryCharacteristic::Eval(eval) => { - DistributedBoundaryConditions::Eval(eval) - } - }); - - let (ny, nx) = (grid.ny(), grid.nx()); - - tids.push( - builder - .spawn(move || { - let mut sys = DistributedSystemPart { - barrier, - ntime, - dt, - current, - fut, - k: [ - Diff::zeros((ny, nx)), - Diff::zeros((ny, nx)), - Diff::zeros((ny, nx)), - Diff::zeros((ny, nx)), - ], - boundary_conditions, - grid: (grid, metrics), - _output: (), - push, - sbp, - t: time, - - wb, - wb_ns: Array2::zeros((4, nx)), - wb_ew: Array2::zeros((4, ny)), - }; - sys.advance(); - }) - .unwrap(), - ); - } - // Set up communicators - // Spawn a new communicator - - DistributedSystem { - _ntime: ntime, - start: b, - sys: tids, - } - } } -// single-threaded items: clone_from -// sync points: scaled_add, rhs -// -// Could instead make every thread (of a threadpool?) carry one grid themselves, -// and instead only wait on obtaining bc, which would make synchronisation -// be somewhat minimal -// -// Difficulties/implementation notes: -// * How to get output from each thread? Push to the multi-producer, single consumer queue -// with capacity of 2*N (which ensures internal synchronisation stays active), or many -// channels which are select!'ed. Choose instead a (name, itime, data) tuple, which is -// communicated. Also initialise the output file all at the start, since we know the -// number of time steps, and which timesteps will be used (should_output). -// * How to balance the number of threads/jobs? -// Each grid will be pushed to a thread, and a thread pool created for use by all threads -// combined. Set affinity to limit the number of cores available (could have side effects -// on heat etc, although this would also be the case for multi-threaded code). A -// thread-per-core architecture + the output thread if we have enough cores, else let os -// control the scheduling for us. -// Use cgroups to artificially limit the number of cores. -// Use taskset to artificially limit the number of cores. -// * Mechanism to wait for bc, is select available? Use a channel with 0-4 capacity, each -// other grid additonally pushes their direction to always use the available -// Use select or similar from crossbeam_channel -// * Mechanism to push bc: try_push and then move to next bc if not available. Spawn a -// thread into a threadpool to handle this. -// * Each grid is spawned onto an async threadpool with M threads allocated to itself, this should -// allow async waiting on the channel (or in the future from process/internet/bc processor) -// * BC: A local thread can have a handle pub struct DistributedSystem { - /// Simple messaging system to be replaced by a more sophisticated system (e.g. run 5 steps, - /// collect, initialise, return something to main) - /// This simply waits until all threads are ready, then starts the computation on all threads - start: Arc, - _ntime: u64, - /// These should be joined to mark the end of the computation + recv: Receiver<(usize, MsgToHost)>, + send: Vec>, + /// All threads should be joined to mark the end of the computation sys: Vec>, } impl DistributedSystem { - 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()); + pub fn advance(&mut self, ntime: u64) { + for tid in &self.send { + tid.send(MsgFromHost::Advance(ntime)); + } } + pub fn output(&self) { + todo!() + } +} + +impl Drop for DistributedSystem { + fn drop(&mut self) { + for tid in &self.send { + tid.send(MsgFromHost::Stop); + } + let handles = std::mem::take(&mut self.sys); + for tid in handles { + tid.join().unwrap() + } + } +} + +enum MsgFromHost { + Advance(u64), + MaxDt(Float), + Output(u64), + Stop, +} + +enum MsgToHost { + MaxDt(Float), + CurrentTimestep(u64), } // #[derive(Debug)] @@ -441,23 +476,39 @@ struct DistributedSystemPart { current: Field, fut: Field, - k: [Diff; 4], - wb: WorkBuffers, - barrier: Arc, - _output: (), // hdf5::Dataset eventually, + t: Float, dt: Float, - ntime: u64, - /// Work buffer for boundaries + name: String, + id: usize, + recv: Receiver, + send: Sender<(usize, MsgToHost)>, + + output: hdf5::Dataset, + + k: [Diff; 4], + wb: WorkBuffers, + /// Work buffer for north/south boundary wb_ns: Array2, + /// Work buffer for east/west boundary wb_ew: Array2, } impl DistributedSystemPart { - fn advance(&mut self) { - self.barrier.wait(); - for _i in 0..self.ntime { + fn run(&mut self) { + loop { + match self.recv.recv().unwrap() { + MsgFromHost::MaxDt(dt) => self.dt = dt, + MsgFromHost::Advance(ntime) => self.advance(ntime), + MsgFromHost::Output(ntime) => todo!(), + MsgFromHost::Stop => return, + } + } + } + fn advance(&mut self, ntime: u64) { + for ntime in 0..ntime { + self.send.send((self.id, MsgToHost::CurrentTimestep(ntime))); let metrics = &self.grid.1; let wb = &mut self.wb.0; let sbp = &self.sbp; @@ -482,7 +533,6 @@ impl DistributedSystemPart { s.send(y.west().to_owned()).unwrap() } - use std::ops::Deref; // This computation does not depend on the boundaries euler::RHS_no_SAT(sbp.deref(), k, y, metrics, wb); From 4d44b4a74a2a034e27eb7963d89025b8baeab135 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Sun, 8 Aug 2021 20:28:37 +0000 Subject: [PATCH 09/25] Fix ouptut for single/multi backend --- multigrid/src/file.rs | 155 ----------------------------- multigrid/src/main.rs | 8 +- multigrid/src/system.rs | 209 ++++++++++++++++++++++++++++++++++++---- 3 files changed, 191 insertions(+), 181 deletions(-) delete mode 100644 multigrid/src/file.rs diff --git a/multigrid/src/file.rs b/multigrid/src/file.rs deleted file mode 100644 index 6be9569..0000000 --- a/multigrid/src/file.rs +++ /dev/null @@ -1,155 +0,0 @@ -use super::*; - -pub struct OutputThread { - rx: Option>>, - tx: Option)>>, - thread: Option>, -} - -impl OutputThread { - pub fn new(file: File) -> Self { - // Pingpong back and forth a number of Vec to be used for the - // output. The sync_channel applies some backpressure - let (tx_thread, rx) = std::sync::mpsc::channel::>(); - let (tx, rx_thread) = std::sync::mpsc::sync_channel::<(u64, Vec)>(3); - let thread = std::thread::Builder::new() - .name("multigrid_output".to_owned()) - .spawn(move || { - let mut times = Vec::::new(); - - for (ntime, fields) in rx_thread.iter() { - if !times.contains(&ntime) { - file.add_timestep(ntime, fields.as_slice()).unwrap(); - times.push(ntime); - } - tx_thread.send(fields).unwrap(); - } - }) - .unwrap(); - - Self { - tx: Some(tx), - rx: Some(rx), - thread: Some(thread), - } - } - - pub fn add_timestep(&mut self, ntime: u64, fields: &[euler::Field]) { - match self.rx.as_ref().unwrap().try_recv() { - Ok(mut copy_fields) => { - for (from, to) in fields.iter().zip(copy_fields.iter_mut()) { - to.clone_from(&from); - } - self.tx - .as_ref() - .unwrap() - .send((ntime, copy_fields)) - .unwrap(); - } - Err(std::sync::mpsc::TryRecvError::Empty) => { - let fields = fields.to_vec(); - self.tx.as_ref().unwrap().send((ntime, fields)).unwrap(); - } - Err(e) => panic!("{:?}", e), - }; - } -} - -impl Drop for OutputThread { - fn drop(&mut self) { - let tx = self.tx.take(); - std::mem::drop(tx); - let thread = self.thread.take().unwrap(); - thread.join().unwrap(); - } -} - -#[derive(Debug, Clone)] -pub struct File(hdf5::File, Vec); - -impl File { - pub fn create>( - path: P, - grids: &[sbp::grid::Grid], - names: Vec, - ) -> Result> { - assert_eq!(grids.len(), names.len()); - let file = hdf5::File::create(path.as_ref())?; - let _tds = file - .new_dataset::() - .resizable(true) - .chunk((1,)) - .create("t", (0,))?; - - for (name, grid) in names.iter().zip(grids.iter()) { - let g = file.create_group(name)?; - g.link_soft("/t", "t").unwrap(); - - let add_dim = |name| { - g.new_dataset::() - .chunk((grid.ny(), grid.nx())) - .gzip(9) - .create(name, (grid.ny(), grid.nx())) - }; - let xds = add_dim("x")?; - xds.write(grid.x())?; - let yds = add_dim("y")?; - yds.write(grid.y())?; - - let add_var = |name| { - g.new_dataset::() - .gzip(3) - .shuffle(true) - .chunk((1, grid.ny(), grid.nx())) - .resizable(true) - .create(name, (0, grid.ny(), grid.nx())) - }; - add_var("rho")?; - add_var("rhou")?; - add_var("rhov")?; - add_var("e")?; - } - - Ok(Self(file, names)) - } - - pub fn add_timestep( - &self, - t: u64, - fields: &[euler::Field], - ) -> Result<(), Box> { - let file = &self.0; - let tds = file.dataset("t")?; - let tpos = tds.size(); - tds.resize((tpos + 1,))?; - tds.write_slice(&[t], ndarray::s![tpos..tpos + 1])?; - - for (groupname, fnow) in self.1.iter().zip(fields.iter()) { - let g = file.group(groupname)?; - let (tpos, ny, nx) = { - let ds = g.dataset("rho")?; - let shape = ds.shape(); - (shape[0], shape[1], shape[2]) - }; - - let rhods = g.dataset("rho")?; - let rhouds = g.dataset("rhou")?; - let rhovds = g.dataset("rhov")?; - let eds = g.dataset("e")?; - - let (rho, rhou, rhov, e) = fnow.components(); - rhods.resize((tpos + 1, ny, nx))?; - rhods.write_slice(rho, ndarray::s![tpos, .., ..])?; - - rhouds.resize((tpos + 1, ny, nx))?; - rhouds.write_slice(rhou, ndarray::s![tpos, .., ..])?; - - rhovds.resize((tpos + 1, ny, nx))?; - rhovds.write_slice(rhov, ndarray::s![tpos, .., ..])?; - - eds.resize((tpos + 1, ny, nx))?; - eds.write_slice(e, ndarray::s![tpos, .., ..])?; - } - Ok(()) - } -} diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index ed16796..2f7a718 100644 --- a/multigrid/src/main.rs +++ b/multigrid/src/main.rs @@ -1,13 +1,10 @@ use argh::FromArgs; -use euler::eval::Evaluator; use sbp::*; -mod file; +mod eval; mod input; mod parsing; -use file::*; -mod eval; mod system; #[derive(Debug, FromArgs)] @@ -90,13 +87,14 @@ fn main() { ); // System::new(grids, grid_connections, operators); - let sys = if opt.distribute { + let mut sys = if opt.distribute { basesystem.create_distributed() } else { basesystem.create() }; let dt = sys.max_dt(); + sys.set_dt(dt); let ntime = (integration_time / dt).round() as u64; let steps_between_outputs = if let Some(n) = opt.number_of_outputs { std::cmp::max(n / ntime, 1) diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 8058b4b..cf92604 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -80,7 +80,38 @@ impl BaseSystem { .map(|(bt, grid)| euler::BoundaryStorage::new(bt, grid)) .collect(); - let sys = SingleThreadedSystem { + let mut outputs = Vec::with_capacity(self.grids.len()); + for (name, grid) in self.names.iter().zip(&self.grids) { + let g = self.output.create_group(name).unwrap(); + g.link_soft("/t", "t").unwrap(); + + let add_dim = |name| { + g.new_dataset::() + .chunk((grid.ny(), grid.nx())) + .gzip(9) + .create(name, (grid.ny(), grid.nx())) + }; + let xds = add_dim("x").unwrap(); + xds.write(grid.x()).unwrap(); + let yds = add_dim("y").unwrap(); + yds.write(grid.y()).unwrap(); + + let add_var = |name| { + g.new_dataset::() + .gzip(3) + .shuffle(true) + .chunk((1, grid.ny(), grid.nx())) + .resizable(true) + .create(name, (0, grid.ny(), grid.nx())) + }; + add_var("rho").unwrap(); + add_var("rhou").unwrap(); + add_var("rhov").unwrap(); + add_var("e").unwrap(); + outputs.push(g); + } + + let mut sys = SingleThreadedSystem { fnow, fnext, k, @@ -90,8 +121,9 @@ impl BaseSystem { bt: self.boundary_conditions, eb, time: self.time, + dt: Float::NAN, operators: self.operators, - output: todo!(), + output: (self.output, outputs), }; match &self.initial_conditions { /* @@ -120,7 +152,6 @@ impl BaseSystem { /// Spreads the computation over n threads, in a thread per grid way. pub fn create_distributed(self) -> System { let nthreads = self.grids.len(); - // let dt = self.max_dt(); // Build up the boundary conditions let mut push_channels: Vec>>>> = @@ -215,6 +246,8 @@ impl BaseSystem { communicators.push(t_send); let time = self.time; + let output = self.output.clone(); + tids.push( builder .spawn(move || { @@ -231,19 +264,48 @@ impl BaseSystem { let wb = WorkBuffers::new(ny, nx); + let g = { + let g = output.create_group(&name).unwrap(); + g.link_soft("/t", "t").unwrap(); + let add_dim = |name| { + g.new_dataset::() + .chunk((grid.ny(), grid.nx())) + .gzip(9) + .create(name, (grid.ny(), grid.nx())) + }; + let xds = add_dim("x").unwrap(); + xds.write(grid.x()).unwrap(); + let yds = add_dim("y").unwrap(); + yds.write(grid.y()).unwrap(); + + let add_var = |name| { + g.new_dataset::() + .gzip(3) + .shuffle(true) + .chunk((1, grid.ny(), grid.nx())) + .resizable(true) + .create(name, (0, grid.ny(), grid.nx())) + }; + add_var("rho").unwrap(); + add_var("rhou").unwrap(); + add_var("rhov").unwrap(); + add_var("e").unwrap(); + g + }; + let mut sys = DistributedSystemPart { current, fut, k, boundary_conditions, grid: (grid, metrics), - output: todo!(), + output: g, push, sbp, t: time, dt: Float::NAN, - name, + _name: name, id, recv: t_recv, @@ -266,6 +328,7 @@ impl BaseSystem { sys: tids, recv: master_recv, send: communicators, + output: self.output, }) } } @@ -277,13 +340,47 @@ pub enum System { impl System { pub fn max_dt(&self) -> Float { - todo!() + match self { + Self::SingleThreaded(sys) => sys.max_dt(), + Self::MultiThreaded(sys) => { + for send in &sys.send { + send.send(MsgFromHost::DtRequest).unwrap(); + } + let mut max_dt = Float::MAX; + let mut to_receive = sys.sys.len(); + while to_receive != 0 { + let dt = match sys.recv.recv().unwrap() { + (_, MsgToHost::MaxDt(dt)) => dt, + _ => unreachable!(), + }; + max_dt = max_dt.min(dt); + to_receive -= 1; + } + max_dt + } + } } - pub fn advance(&self, nsteps: u64) { - todo!() + pub fn set_dt(&mut self, dt: Float) { + match self { + Self::SingleThreaded(sys) => sys.dt = dt, + Self::MultiThreaded(sys) => { + for tid in &sys.send { + tid.send(MsgFromHost::DtSet(dt)).unwrap() + } + } + } + } + pub fn advance(&mut self, nsteps: u64) { + match self { + Self::SingleThreaded(sys) => sys.advance(nsteps), + Self::MultiThreaded(sys) => sys.advance(nsteps), + } } pub fn output(&self, ntime: u64) { - todo!() + match self { + Self::SingleThreaded(sys) => sys.output(ntime), + Self::MultiThreaded(sys) => sys.output(ntime), + } } } @@ -297,8 +394,9 @@ pub struct SingleThreadedSystem { pub bt: Vec, pub eb: Vec, pub time: Float, + pub dt: Float, pub operators: Vec>, - pub output: Vec, + pub output: (hdf5::File, Vec), } impl integrate::Integrable for SingleThreadedSystem { @@ -319,7 +417,13 @@ impl SingleThreadedSystem { } } - pub fn advance(&mut self, dt: Float) { + pub fn advance(&mut self, nsteps: u64) { + for _ in 0..nsteps { + self.advance_single_step(self.dt) + } + } + + pub fn advance_single_step(&mut self, dt: Float) { let metrics = &self.metrics; let grids = &self.grids; let bt = &self.bt; @@ -409,6 +513,34 @@ impl SingleThreadedSystem { max_dt } + + pub fn output(&self, ntime: u64) { + let tds = self.output.0.dataset("t").unwrap(); + let tpos = tds.size(); + tds.resize((tpos + 1,)).unwrap(); + tds.write_slice(&[ntime], ndarray::s![tpos..tpos + 1]) + .unwrap(); + for (group, fnow) in self.output.1.iter().zip(&self.fnow) { + let (ny, nx) = (fnow.ny(), fnow.nx()); + let rhods = group.dataset("rho").unwrap(); + let rhouds = group.dataset("rhou").unwrap(); + let rhovds = group.dataset("rhov").unwrap(); + let eds = group.dataset("e").unwrap(); + + let (rho, rhou, rhov, e) = fnow.components(); + rhods.resize((tpos + 1, ny, nx)).unwrap(); + rhods.write_slice(rho, ndarray::s![tpos, .., ..]).unwrap(); + + rhouds.resize((tpos + 1, ny, nx)).unwrap(); + rhouds.write_slice(rhou, ndarray::s![tpos, .., ..]).unwrap(); + + rhovds.resize((tpos + 1, ny, nx)).unwrap(); + rhovds.write_slice(rhov, ndarray::s![tpos, .., ..]).unwrap(); + + eds.resize((tpos + 1, ny, nx)).unwrap(); + eds.write_slice(e, ndarray::s![tpos, .., ..]).unwrap(); + } + } } pub struct DistributedSystem { @@ -416,23 +548,31 @@ pub struct DistributedSystem { send: Vec>, /// All threads should be joined to mark the end of the computation sys: Vec>, + output: hdf5::File, } impl DistributedSystem { pub fn advance(&mut self, ntime: u64) { for tid in &self.send { - tid.send(MsgFromHost::Advance(ntime)); + tid.send(MsgFromHost::Advance(ntime)).unwrap(); } } - pub fn output(&self) { - todo!() + pub fn output(&self, ntime: u64) { + for tid in &self.send { + tid.send(MsgFromHost::Output(ntime)).unwrap(); + } + let tds = self.output.dataset("t").unwrap(); + let tpos = tds.size(); + tds.resize((tpos + 1,)).unwrap(); + tds.write_slice(&[ntime], ndarray::s![tpos..tpos + 1]) + .unwrap(); } } impl Drop for DistributedSystem { fn drop(&mut self) { for tid in &self.send { - tid.send(MsgFromHost::Stop); + tid.send(MsgFromHost::Stop).unwrap(); } let handles = std::mem::take(&mut self.sys); for tid in handles { @@ -443,7 +583,8 @@ impl Drop for DistributedSystem { enum MsgFromHost { Advance(u64), - MaxDt(Float), + DtRequest, + DtSet(Float), Output(u64), Stop, } @@ -480,12 +621,12 @@ struct DistributedSystemPart { t: Float, dt: Float, - name: String, + _name: String, id: usize, recv: Receiver, send: Sender<(usize, MsgToHost)>, - output: hdf5::Dataset, + output: hdf5::Group, k: [Diff; 4], wb: WorkBuffers, @@ -499,16 +640,42 @@ impl DistributedSystemPart { fn run(&mut self) { loop { match self.recv.recv().unwrap() { - MsgFromHost::MaxDt(dt) => self.dt = dt, + MsgFromHost::DtSet(dt) => self.dt = dt, + MsgFromHost::DtRequest => todo!(), MsgFromHost::Advance(ntime) => self.advance(ntime), - MsgFromHost::Output(ntime) => todo!(), + MsgFromHost::Output(ntime) => self.output(ntime), MsgFromHost::Stop => return, } } } + + fn output(&mut self, _ntime: u64) { + let (ny, nx) = (self.current.ny(), self.current.nx()); + let rhods = self.output.dataset("rho").unwrap(); + let rhouds = self.output.dataset("rhou").unwrap(); + let rhovds = self.output.dataset("rhov").unwrap(); + let eds = self.output.dataset("e").unwrap(); + + let (rho, rhou, rhov, e) = self.current.components(); + let tpos = rhods.size() / (ny * nx); + rhods.resize((tpos + 1, ny, nx)).unwrap(); + rhods.write_slice(rho, ndarray::s![tpos, .., ..]).unwrap(); + + rhouds.resize((tpos + 1, ny, nx)).unwrap(); + rhouds.write_slice(rhou, ndarray::s![tpos, .., ..]).unwrap(); + + rhovds.resize((tpos + 1, ny, nx)).unwrap(); + rhovds.write_slice(rhov, ndarray::s![tpos, .., ..]).unwrap(); + + eds.resize((tpos + 1, ny, nx)).unwrap(); + eds.write_slice(e, ndarray::s![tpos, .., ..]).unwrap(); + } + fn advance(&mut self, ntime: u64) { for ntime in 0..ntime { - self.send.send((self.id, MsgToHost::CurrentTimestep(ntime))); + self.send + .send((self.id, MsgToHost::CurrentTimestep(ntime))) + .unwrap(); let metrics = &self.grid.1; let wb = &mut self.wb.0; let sbp = &self.sbp; From 2d473b825558d4fe71da8aea7a8854c4684acdca Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Mon, 16 Aug 2021 20:33:57 +0000 Subject: [PATCH 10/25] Clippy lints --- euler/src/eval.rs | 3 +++ euler/src/lib.rs | 1 + multigrid/src/eval.rs | 1 + multigrid/src/eval/evalexpr.rs | 2 ++ multigrid/src/main.rs | 4 ++-- multigrid/src/system.rs | 5 +++-- sbp/src/operators/algos.rs | 1 + utils/integrate/src/lib.rs | 2 +- 8 files changed, 14 insertions(+), 5 deletions(-) diff --git a/euler/src/eval.rs b/euler/src/eval.rs index 5e4c7cc..02344f6 100644 --- a/euler/src/eval.rs +++ b/euler/src/eval.rs @@ -5,6 +5,7 @@ use super::GAMMA; use ndarray::{azip, ArrayView, ArrayViewMut, Dimension}; pub trait Evaluator: Send + Sync { + #[allow(clippy::too_many_arguments)] fn evaluate( &self, t: Float, @@ -55,6 +56,7 @@ pub trait EvaluatorPressure: Send + Sync { rho: ArrayView, out: ArrayViewMut, ); + #[allow(clippy::too_many_arguments)] fn p( &self, t: Float, @@ -70,6 +72,7 @@ pub trait EvaluatorPressure: Send + Sync { impl<'a, D: Dimension, BP: EvaluatorPressure> Evaluator for EvaluatorPressureWrapper<'a, D, BP> { + #[allow(clippy::many_single_char_names)] fn evaluate( &self, t: Float, diff --git a/euler/src/lib.rs b/euler/src/lib.rs index 350cb5f..f712fb1 100644 --- a/euler/src/lib.rs +++ b/euler/src/lib.rs @@ -330,6 +330,7 @@ impl Field { let (rho, rhou, rhov, e) = self.components_mut(); vortex_param.evaluate(time, x, y, rho, rhou, rhov, e) } + #[allow(clippy::erasing_op, clippy::identity_op)] fn iter(&self) -> impl ExactSizeIterator + '_ { let n = self.nx() * self.ny(); let slice = self.0.as_slice().unwrap(); diff --git a/multigrid/src/eval.rs b/multigrid/src/eval.rs index 8b4b903..52c0857 100644 --- a/multigrid/src/eval.rs +++ b/multigrid/src/eval.rs @@ -10,6 +10,7 @@ pub enum Evaluator { } impl euler::eval::Evaluator for Evaluator { + #[allow(clippy::many_single_char_names)] fn evaluate( &self, t: Float, diff --git a/multigrid/src/eval/evalexpr.rs b/multigrid/src/eval/evalexpr.rs index e999252..ea50b71 100644 --- a/multigrid/src/eval/evalexpr.rs +++ b/multigrid/src/eval/evalexpr.rs @@ -79,6 +79,7 @@ pub struct EvaluatorConservation { } impl euler::eval::Evaluator for Evaluator { + #[allow(clippy::many_single_char_names)] fn evaluate( &self, t: Float, @@ -267,6 +268,7 @@ impl euler::eval::EvaluatorPressure for EvaluatorPressure { }) } + #[allow(clippy::many_single_char_names)] fn p( &self, t: Float, diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index 2f7a718..ad205fa 100644 --- a/multigrid/src/main.rs +++ b/multigrid/src/main.rs @@ -77,12 +77,12 @@ fn main() { } = config.into_runtime(); let basesystem = system::BaseSystem::new( - names.clone(), + names, grids, 0.0, operators, boundary_conditions, - initial_conditions.clone(), + initial_conditions, opt.output.clone(), ); // System::new(grids, grid_connections, operators); diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index cf92604..c698203 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -48,6 +48,7 @@ impl BaseSystem { output, } } + #[allow(clippy::many_single_char_names)] pub fn create(self) -> System { let fnow = self .grids @@ -134,7 +135,7 @@ impl BaseSystem { } } */ - parsing::InitialConditions::Vortex(vortexparams) => sys.vortex(0.0, &vortexparams), + parsing::InitialConditions::Vortex(vortexparams) => sys.vortex(0.0, vortexparams), parsing::InitialConditions::Expressions(expr) => { let t = 0.0; for (grid, field) in sys.grids.iter().zip(sys.fnow.iter_mut()) { @@ -413,7 +414,7 @@ impl integrate::Integrable for SingleThreadedSystem { impl SingleThreadedSystem { pub fn vortex(&mut self, t: Float, vortex_params: &euler::VortexParameters) { for (f, g) in self.fnow.iter_mut().zip(&self.grids) { - f.vortex(g.x(), g.y(), t, &vortex_params); + f.vortex(g.x(), g.y(), t, vortex_params); } } diff --git a/sbp/src/operators/algos.rs b/sbp/src/operators/algos.rs index 20ae1d5..74c351f 100644 --- a/sbp/src/operators/algos.rs +++ b/sbp/src/operators/algos.rs @@ -281,6 +281,7 @@ pub(crate) fn diff_op_2d_fallback( matrix: &BlockMatrix, optype: OperatorType, diff --git a/utils/integrate/src/lib.rs b/utils/integrate/src/lib.rs index 3785a77..d91afce 100644 --- a/utils/integrate/src/lib.rs +++ b/utils/integrate/src/lib.rs @@ -215,7 +215,7 @@ pub fn integrate( } }; - rhs(&mut k[i], &fut, simtime); + rhs(&mut k[i], fut, simtime); } } From 296dc98e017d4d9e335a6ef421c17c9a03cd3b8f Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Tue, 17 Aug 2021 11:48:52 +0000 Subject: [PATCH 11/25] Add progressbar to single/multi system --- multigrid/src/main.rs | 42 ++++++++++++++---------------------- multigrid/src/system.rs | 48 ++++++++++++++++++++++++++++++++++++++++- 2 files changed, 63 insertions(+), 27 deletions(-) diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index ad205fa..99d33eb 100644 --- a/multigrid/src/main.rs +++ b/multigrid/src/main.rs @@ -85,7 +85,6 @@ fn main() { initial_conditions, opt.output.clone(), ); - // System::new(grids, grid_connections, operators); let mut sys = if opt.distribute { basesystem.create_distributed() @@ -103,11 +102,10 @@ fn main() { }; sys.output(0); - //let output = File::create(&opt.output, sys.grids.as_slice(), names).unwrap(); - //let mut output = OutputThread::new(output); - //output.add_timestep(0, &sys.fnow); - let progressbar = progressbar(opt.no_progressbar, ntime); + if !opt.no_progressbar { + sys.add_progressbar(ntime) + } let timer = if opt.timings { Some(std::time::Instant::now()) @@ -118,22 +116,18 @@ fn main() { let mut itime = 0; while itime < ntime { sys.advance(steps_between_outputs); - progressbar.inc(1); itime += steps_between_outputs; sys.output(itime); } - - /* - for itime in 0..ntime { - if should_output(itime) { - output.add_timestep(itime, &sys.fnow); - } - progressbar.inc(1); - sys.advance(dt); + if itime != ntime { + sys.advance(ntime - itime); + sys.output(ntime); + } + + if !opt.no_progressbar { + sys.finish_progressbar(); } - */ - progressbar.finish_and_clear(); let mut outinfo = OutputInformation { filename: opt.output, @@ -181,14 +175,10 @@ fn main() { } } -fn progressbar(dummy: bool, ntime: u64) -> indicatif::ProgressBar { - if dummy { - indicatif::ProgressBar::hidden() - } else { - let progressbar = indicatif::ProgressBar::new(ntime); - progressbar.with_style( - indicatif::ProgressStyle::default_bar() - .template("{wide_bar:.cyan/blue} {pos}/{len} ({eta})"), - ) - } +fn progressbar(ntime: u64) -> indicatif::ProgressBar { + let progressbar = indicatif::ProgressBar::new(ntime); + progressbar.with_style( + indicatif::ProgressStyle::default_bar() + .template("{wide_bar:.cyan/blue} {pos}/{len} ({eta})"), + ) } diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index c698203..bcc14dd 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -125,6 +125,7 @@ impl BaseSystem { dt: Float::NAN, operators: self.operators, output: (self.output, outputs), + progressbar: None, }; match &self.initial_conditions { /* @@ -330,6 +331,7 @@ impl BaseSystem { recv: master_recv, send: communicators, output: self.output, + progressbar: None, }) } } @@ -383,6 +385,24 @@ impl System { Self::MultiThreaded(sys) => sys.output(ntime), } } + pub fn add_progressbar(&mut self, ntime: u64) { + match self { + Self::SingleThreaded(sys) => sys.progressbar = Some(super::progressbar(ntime)), + Self::MultiThreaded(sys) => sys.attach_progressbar(ntime), + } + } + pub fn finish_progressbar(&mut self) { + match self { + Self::SingleThreaded(sys) => sys.progressbar.take().unwrap().finish_and_clear(), + Self::MultiThreaded(sys) => { + let (target, pbs) = sys.progressbar.take().unwrap(); + for pb in pbs.into_iter() { + pb.finish_and_clear() + } + target.join_and_clear().unwrap(); + } + } + } } pub struct SingleThreadedSystem { @@ -398,6 +418,7 @@ pub struct SingleThreadedSystem { pub dt: Float, pub operators: Vec>, pub output: (hdf5::File, Vec), + pub progressbar: Option, } impl integrate::Integrable for SingleThreadedSystem { @@ -420,7 +441,10 @@ impl SingleThreadedSystem { pub fn advance(&mut self, nsteps: u64) { for _ in 0..nsteps { - self.advance_single_step(self.dt) + self.advance_single_step(self.dt); + if let Some(pbar) = &self.progressbar { + pbar.inc(1) + } } } @@ -550,6 +574,7 @@ pub struct DistributedSystem { /// All threads should be joined to mark the end of the computation sys: Vec>, output: hdf5::File, + progressbar: Option<(indicatif::MultiProgress, Vec)>, } impl DistributedSystem { @@ -557,6 +582,17 @@ impl DistributedSystem { for tid in &self.send { tid.send(MsgFromHost::Advance(ntime)).unwrap(); } + if let Some(pbar) = &self.progressbar { + let expected_messages = ntime * self.sys.len() as u64; + for _i in 0..expected_messages { + match self.recv.recv().unwrap() { + (i, MsgToHost::CurrentTimestep(_)) => { + pbar.1[i].inc(1); + } + _ => unreachable!(), + } + } + } } pub fn output(&self, ntime: u64) { for tid in &self.send { @@ -568,6 +604,16 @@ impl DistributedSystem { tds.write_slice(&[ntime], ndarray::s![tpos..tpos + 1]) .unwrap(); } + pub fn attach_progressbar(&mut self, ntime: u64) { + let target = indicatif::MultiProgress::new(); + let mut progressbars = Vec::with_capacity(self.sys.len()); + for _ in 0..self.sys.len() { + let pb = super::progressbar(ntime); + progressbars.push(target.add(pb)); + } + target.set_move_cursor(true); + self.progressbar = Some((target, progressbars)); + } } impl Drop for DistributedSystem { From 1bfd37b1644c12ec5a98df5b420b8c015fe795cb Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Tue, 17 Aug 2021 12:35:05 +0000 Subject: [PATCH 12/25] Fixup max dt --- multigrid/src/system.rs | 49 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 46 insertions(+), 3 deletions(-) diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index bcc14dd..dd8bcb7 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -318,8 +318,6 @@ impl BaseSystem { wb_ew: Array2::zeros((4, ny)), }; - // init and send maxdt - // receive maxdt sys.run(); }) .unwrap(), @@ -688,7 +686,10 @@ impl DistributedSystemPart { loop { match self.recv.recv().unwrap() { MsgFromHost::DtSet(dt) => self.dt = dt, - MsgFromHost::DtRequest => todo!(), + MsgFromHost::DtRequest => { + let dt = self.max_dt(); + self.send.send((self.id, MsgToHost::MaxDt(dt))).unwrap(); + } MsgFromHost::Advance(ntime) => self.advance(ntime), MsgFromHost::Output(ntime) => self.output(ntime), MsgFromHost::Stop => return, @@ -696,6 +697,48 @@ impl DistributedSystemPart { } } + fn max_dt(&self) -> Float { + let nx = self.current.nx(); + let ny = self.current.ny(); + + let (rho, rhou, rhov, _e) = self.current.components(); + + let mut max_u: Float = 0.0; + let mut max_v: Float = 0.0; + + for ((((((rho, rhou), rhov), detj_dxi_dx), detj_dxi_dy), detj_deta_dx), detj_deta_dy) in rho + .iter() + .zip(rhou.iter()) + .zip(rhov.iter()) + .zip(self.grid.1.detj_dxi_dx()) + .zip(self.grid.1.detj_dxi_dy()) + .zip(self.grid.1.detj_deta_dx()) + .zip(self.grid.1.detj_deta_dy()) + { + let u = rhou / rho; + let v = rhov / rho; + + let uhat: Float = detj_dxi_dx * u + detj_dxi_dy * v; + let vhat: Float = detj_deta_dx * u + detj_deta_dy * v; + + max_u = max_u.max(uhat.abs()); + max_v = max_v.max(vhat.abs()); + } + + let dx = 1.0 / nx as Float; + let dy = 1.0 / ny as Float; + + let c_dt = Float::max(max_u / dx, max_v / dy); + + let c_max = if self.sbp.is_h2xi() || self.sbp.is_h2eta() { + 0.5 + } else { + 1.0 + }; + + c_max / c_dt + } + fn output(&mut self, _ntime: u64) { let (ny, nx) = (self.current.ny(), self.current.nx()); let rhods = self.output.dataset("rho").unwrap(); From 05cb4551085445168c047d59a3f4b280ae8e819b Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Tue, 17 Aug 2021 12:42:24 +0000 Subject: [PATCH 13/25] Apply initial conditions to distributed --- multigrid/src/system.rs | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index dd8bcb7..016aab2 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -250,11 +250,26 @@ impl BaseSystem { let output = self.output.clone(); + let initial = self.initial_conditions.clone(); + tids.push( builder .spawn(move || { let (ny, nx) = (grid.ny(), grid.nx()); - let current = Field::new(ny, nx); + let mut current = Field::new(ny, nx); + + match &initial { + parsing::InitialConditions::Vortex(vortexparams) => { + current.vortex(grid.x(), grid.y(), time, vortexparams) + } + parsing::InitialConditions::Expressions(expr) => { + // Evaluate the expressions on all variables + let x = grid.x(); + let y = grid.y(); + let (rho, rhou, rhov, e) = current.components_mut(); + (*expr).evaluate(time, x, y, rho, rhou, rhov, e); + } + } let fut = current.clone(); let k = [ Diff::zeros((ny, nx)), From ed81ba995f5180e239d48d2649ae7f9c5e6fbce7 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Tue, 17 Aug 2021 13:15:36 +0000 Subject: [PATCH 14/25] Update indicatif for multi-progressbar --- multigrid/Cargo.toml | 2 +- multigrid/src/system.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/multigrid/Cargo.toml b/multigrid/Cargo.toml index 52ae7c9..067c0bc 100644 --- a/multigrid/Cargo.toml +++ b/multigrid/Cargo.toml @@ -11,7 +11,7 @@ euler = { path = "../euler", features = ["serde1"] } hdf5 = "0.7.0" integrate = { path = "../utils/integrate" } rayon = "1.3.0" -indicatif = "0.15.0" +indicatif = "0.17.0-beta.1" ndarray = { version = "0.14.0", features = ["serde"] } serde = { version = "1.0.115", features = ["derive"] } json5 = "0.3.0" diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 016aab2..2f25453 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -412,7 +412,7 @@ impl System { for pb in pbs.into_iter() { pb.finish_and_clear() } - target.join_and_clear().unwrap(); + target.clear().unwrap(); } } } From 67884c38c639acab8e8f6649f3d4c340de8e2c2b Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Tue, 17 Aug 2021 13:40:23 +0000 Subject: [PATCH 15/25] Simplify time tracking loop --- multigrid/src/main.rs | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index 99d33eb..fe748aa 100644 --- a/multigrid/src/main.rs +++ b/multigrid/src/main.rs @@ -115,15 +115,12 @@ fn main() { let mut itime = 0; while itime < ntime { - sys.advance(steps_between_outputs); + let nexttime = (itime + steps_between_outputs).max(ntime); + sys.advance(nexttime - itime); - itime += steps_between_outputs; + itime = nexttime; sys.output(itime); } - if itime != ntime { - sys.advance(ntime - itime); - sys.output(ntime); - } if !opt.no_progressbar { sys.finish_progressbar(); From a843ad9974ddb56ef150a214843fd267f2cf9eec Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Wed, 18 Aug 2021 10:21:38 +0000 Subject: [PATCH 16/25] Add error output for single/multi --- multigrid/src/main.rs | 20 +---------- multigrid/src/system.rs | 79 +++++++++++++++++++++++++++++++++++++++-- 2 files changed, 78 insertions(+), 21 deletions(-) diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index fe748aa..e57be84 100644 --- a/multigrid/src/main.rs +++ b/multigrid/src/main.rs @@ -138,27 +138,9 @@ fn main() { //output.add_timestep(ntime, &sys.fnow); - /* if opt.error { - let time = ntime as Float * dt; - let mut e = 0.0; - for ((fmod, grid), op) in sys.fnow.iter().zip(&sys.grids).zip(&sys.operators) { - let mut fvort = fmod.clone(); - match &initial_conditions { - parsing::InitialConditions::Vortex(vortexparams) => { - fvort.vortex(grid.x(), grid.y(), time, &vortexparams); - } - parsing::InitialConditions::Expressions(expr) => { - let (rho, rhou, rhov, e) = fvort.components_mut(); - expr.as_ref() - .evaluate(time, grid.x(), grid.y(), rho, rhou, rhov, e) - } - } - e += fmod.h2_err(&fvort, &**op); - } - outinfo.error = Some(e); + outinfo.error = Some(sys.error()) } - */ if opt.output_json { println!("{}", json5::to_string(&outinfo).unwrap()); diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 2f25453..3f65b7b 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -126,6 +126,7 @@ impl BaseSystem { operators: self.operators, output: (self.output, outputs), progressbar: None, + initial_conditions: self.initial_conditions.clone(), }; match &self.initial_conditions { /* @@ -250,7 +251,7 @@ impl BaseSystem { let output = self.output.clone(); - let initial = self.initial_conditions.clone(); + let initial_conditions = self.initial_conditions.clone(); tids.push( builder @@ -258,7 +259,7 @@ impl BaseSystem { let (ny, nx) = (grid.ny(), grid.nx()); let mut current = Field::new(ny, nx); - match &initial { + match &initial_conditions { parsing::InitialConditions::Vortex(vortexparams) => { current.vortex(grid.x(), grid.y(), time, vortexparams) } @@ -321,6 +322,7 @@ impl BaseSystem { sbp, t: time, dt: Float::NAN, + initial_conditions, _name: name, id, @@ -416,6 +418,24 @@ impl System { } } } + pub fn error(&self) -> Float { + match self { + Self::SingleThreaded(sys) => sys.error(), + Self::MultiThreaded(sys) => { + for sender in &sys.send { + sender.send(MsgFromHost::Error).unwrap(); + } + let mut e = 0.0; + for _ in 0..sys.sys.len() { + e += match sys.recv.recv().unwrap() { + (_, MsgToHost::Error(e)) => e, + (_, m) => panic!("Unexpected message: {:?}", m), + } + } + e + } + } + } } pub struct SingleThreadedSystem { @@ -432,6 +452,7 @@ pub struct SingleThreadedSystem { pub operators: Vec>, pub output: (hdf5::File, Vec), pub progressbar: Option, + pub initial_conditions: parsing::InitialConditions, } impl integrate::Integrable for SingleThreadedSystem { @@ -579,6 +600,25 @@ impl SingleThreadedSystem { eds.write_slice(e, ndarray::s![tpos, .., ..]).unwrap(); } } + + pub fn error(&self) -> Float { + let mut e = 0.0; + for ((fmod, grid), op) in self.fnow.iter().zip(&self.grids).zip(&self.operators) { + let mut fvort = fmod.clone(); + match &self.initial_conditions { + parsing::InitialConditions::Vortex(vortexparams) => { + fvort.vortex(grid.x(), grid.y(), self.time, &vortexparams); + } + parsing::InitialConditions::Expressions(expr) => { + let (rho, rhou, rhov, e) = fvort.components_mut(); + expr.as_ref() + .evaluate(self.time, grid.x(), grid.y(), rho, rhou, rhov, e) + } + } + e += fmod.h2_err(&fvort, &**op); + } + e + } } pub struct DistributedSystem { @@ -641,17 +681,32 @@ impl Drop for DistributedSystem { } } +/// Messages sent from the host to each compute thread +#[derive(Debug)] enum MsgFromHost { + /// Advance n steps Advance(u64), + /// Compute the maximum dt allowed by this grid DtRequest, + /// Set dt DtSet(Float), + /// Output the current time to file Output(u64), + /// Stop all computing Stop, + /// Request the current error + Error, } +/// Messages sent back to the host +#[derive(Debug)] enum MsgToHost { + /// Maximum dt allowed by the current grid MaxDt(Float), + /// Timestep which we have currently computed CurrentTimestep(u64), + /// Error from the current grid + Error(Float), } // #[derive(Debug)] @@ -687,6 +742,7 @@ struct DistributedSystemPart { send: Sender<(usize, MsgToHost)>, output: hdf5::Group, + initial_conditions: crate::parsing::InitialConditions, k: [Diff; 4], wb: WorkBuffers, @@ -708,6 +764,10 @@ impl DistributedSystemPart { MsgFromHost::Advance(ntime) => self.advance(ntime), MsgFromHost::Output(ntime) => self.output(ntime), MsgFromHost::Stop => return, + MsgFromHost::Error => self + .send + .send((self.id, MsgToHost::Error(self.error()))) + .unwrap(), } } } @@ -1066,4 +1126,19 @@ impl DistributedSystemPart { ) } } + + fn error(&self) -> Float { + let mut fvort = self.current.clone(); + match &self.initial_conditions { + parsing::InitialConditions::Vortex(vortexparams) => { + fvort.vortex(self.grid.0.x(), self.grid.0.y(), self.t, &vortexparams); + } + parsing::InitialConditions::Expressions(expr) => { + let (rho, rhou, rhov, e) = fvort.components_mut(); + expr.as_ref() + .evaluate(self.t, self.grid.0.x(), self.grid.0.y(), rho, rhou, rhov, e) + } + } + self.current.h2_err(&fvort, &*self.sbp) + } } From d6356da39380ab36db699cf5cc4ac0b487280a9b Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Wed, 18 Aug 2021 10:32:29 +0000 Subject: [PATCH 17/25] Swap solutions to propagate the solution --- multigrid/src/system.rs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 3f65b7b..1ba6bd1 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -1123,7 +1123,8 @@ impl DistributedSystemPart { &mut self.t, self.dt, &mut self.k, - ) + ); + std::mem::swap(&mut self.current, &mut self.fut); } } From 4b6ee18491032427ecc89510ddbb0a1e0512403f Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Wed, 18 Aug 2021 11:06:29 +0000 Subject: [PATCH 18/25] Add progressbar inside multi-sys --- multigrid/src/system.rs | 61 +++++++++++++++++++++-------------------- 1 file changed, 31 insertions(+), 30 deletions(-) diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 1ba6bd1..add187f 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -333,6 +333,8 @@ impl BaseSystem { wb, wb_ns: Array2::zeros((4, nx)), wb_ew: Array2::zeros((4, ny)), + + progressbar: None, }; sys.run(); @@ -410,10 +412,10 @@ impl System { match self { Self::SingleThreaded(sys) => sys.progressbar.take().unwrap().finish_and_clear(), Self::MultiThreaded(sys) => { - let (target, pbs) = sys.progressbar.take().unwrap(); - for pb in pbs.into_iter() { - pb.finish_and_clear() + for tid in &sys.send { + tid.send(MsgFromHost::ProgressbarDrop).unwrap(); } + let target = sys.progressbar.take().unwrap(); target.clear().unwrap(); } } @@ -627,7 +629,7 @@ pub struct DistributedSystem { /// All threads should be joined to mark the end of the computation sys: Vec>, output: hdf5::File, - progressbar: Option<(indicatif::MultiProgress, Vec)>, + progressbar: Option, } impl DistributedSystem { @@ -635,17 +637,6 @@ impl DistributedSystem { for tid in &self.send { tid.send(MsgFromHost::Advance(ntime)).unwrap(); } - if let Some(pbar) = &self.progressbar { - let expected_messages = ntime * self.sys.len() as u64; - for _i in 0..expected_messages { - match self.recv.recv().unwrap() { - (i, MsgToHost::CurrentTimestep(_)) => { - pbar.1[i].inc(1); - } - _ => unreachable!(), - } - } - } } pub fn output(&self, ntime: u64) { for tid in &self.send { @@ -659,13 +650,13 @@ impl DistributedSystem { } pub fn attach_progressbar(&mut self, ntime: u64) { let target = indicatif::MultiProgress::new(); - let mut progressbars = Vec::with_capacity(self.sys.len()); - for _ in 0..self.sys.len() { + for tid in &self.send { let pb = super::progressbar(ntime); - progressbars.push(target.add(pb)); + let pb = target.add(pb); + tid.send(MsgFromHost::Progressbar(pb)).unwrap(); } target.set_move_cursor(true); - self.progressbar = Some((target, progressbars)); + self.progressbar = Some(target); } } @@ -696,6 +687,10 @@ enum MsgFromHost { Stop, /// Request the current error Error, + /// Progressbar to report progress + Progressbar(indicatif::ProgressBar), + /// Clear and remove the progressbar + ProgressbarDrop, } /// Messages sent back to the host @@ -703,8 +698,6 @@ enum MsgFromHost { enum MsgToHost { /// Maximum dt allowed by the current grid MaxDt(Float), - /// Timestep which we have currently computed - CurrentTimestep(u64), /// Error from the current grid Error(Float), } @@ -750,6 +743,8 @@ struct DistributedSystemPart { wb_ns: Array2, /// Work buffer for east/west boundary wb_ew: Array2, + + progressbar: Option, } impl DistributedSystemPart { @@ -759,15 +754,17 @@ impl DistributedSystemPart { MsgFromHost::DtSet(dt) => self.dt = dt, MsgFromHost::DtRequest => { let dt = self.max_dt(); - self.send.send((self.id, MsgToHost::MaxDt(dt))).unwrap(); + self.send(MsgToHost::MaxDt(dt)).unwrap(); } MsgFromHost::Advance(ntime) => self.advance(ntime), MsgFromHost::Output(ntime) => self.output(ntime), MsgFromHost::Stop => return, - MsgFromHost::Error => self - .send - .send((self.id, MsgToHost::Error(self.error()))) - .unwrap(), + MsgFromHost::Error => self.send(MsgToHost::Error(self.error())).unwrap(), + MsgFromHost::Progressbar(pbar) => self.progressbar = Some(pbar), + MsgFromHost::ProgressbarDrop => { + let pb = self.progressbar.take().unwrap(); + pb.finish_and_clear() + } } } } @@ -837,10 +834,10 @@ impl DistributedSystemPart { } fn advance(&mut self, ntime: u64) { - for ntime in 0..ntime { - self.send - .send((self.id, MsgToHost::CurrentTimestep(ntime))) - .unwrap(); + for _itime in 0..ntime { + if let Some(pbar) = &self.progressbar { + pbar.inc(1) + } let metrics = &self.grid.1; let wb = &mut self.wb.0; let sbp = &self.sbp; @@ -1128,6 +1125,10 @@ impl DistributedSystemPart { } } + fn send(&self, msg: MsgToHost) -> Result<(), crossbeam_channel::SendError<(usize, MsgToHost)>> { + self.send.send((self.id, msg)) + } + fn error(&self) -> Float { let mut fvort = self.current.clone(); match &self.initial_conditions { From 4f0af1f6c1faa8faec6e59af499689ec8cf4797b Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Wed, 18 Aug 2021 11:45:55 +0000 Subject: [PATCH 19/25] Add synchronisation of multi-threaded system --- multigrid/Cargo.toml | 1 + multigrid/src/main.rs | 34 +++++++++++++++++++++------------- multigrid/src/system.rs | 15 +++++++++++++++ 3 files changed, 37 insertions(+), 13 deletions(-) diff --git a/multigrid/Cargo.toml b/multigrid/Cargo.toml index 067c0bc..112eaa6 100644 --- a/multigrid/Cargo.toml +++ b/multigrid/Cargo.toml @@ -19,3 +19,4 @@ indexmap = { version = "1.5.2", features = ["serde-1"] } argh = "0.1.4" evalexpr = "6.3.0" crossbeam-channel = "0.5.0" +crossbeam-utils = "0.8.5" diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs index e57be84..9b83ecd 100644 --- a/multigrid/src/main.rs +++ b/multigrid/src/main.rs @@ -108,6 +108,9 @@ fn main() { } let timer = if opt.timings { + if let system::System::MultiThreaded(sys) = &sys { + sys.synchronise() + } Some(std::time::Instant::now()) } else { None @@ -119,25 +122,30 @@ fn main() { sys.advance(nexttime - itime); itime = nexttime; - sys.output(itime); + if itime != ntime { + sys.output(itime); + } } + let timer = timer.map(|timer| { + if let system::System::MultiThreaded(sys) = &sys { + sys.synchronise(); + } + + timer.elapsed() + }); + sys.output(ntime); + + let mut outinfo = OutputInformation { + filename: opt.output, + time_elapsed: timer, + ..Default::default() + }; + if !opt.no_progressbar { sys.finish_progressbar(); } - let mut outinfo = OutputInformation { - filename: opt.output, - ..Default::default() - }; - - if let Some(timer) = timer { - let duration = timer.elapsed(); - outinfo.time_elapsed = Some(duration); - } - - //output.add_timestep(ntime, &sys.fnow); - if opt.error { outinfo.error = Some(sys.error()) } diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index add187f..f756c3c 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -415,6 +415,7 @@ impl System { for tid in &sys.send { tid.send(MsgFromHost::ProgressbarDrop).unwrap(); } + sys.synchronise(); let target = sys.progressbar.take().unwrap(); target.clear().unwrap(); } @@ -658,6 +659,17 @@ impl DistributedSystem { target.set_move_cursor(true); self.progressbar = Some(target); } + fn send_barrier(&self, barrier: &crossbeam_utils::sync::WaitGroup) { + for tid in &self.send { + tid.send(MsgFromHost::Barrier(barrier.clone())).unwrap() + } + } + pub fn synchronise(&self) { + // Syncronise before starting the timer + let barrier = crossbeam_utils::sync::WaitGroup::new(); + self.send_barrier(&barrier); + barrier.wait(); + } } impl Drop for DistributedSystem { @@ -687,6 +699,8 @@ enum MsgFromHost { Stop, /// Request the current error Error, + /// A barrier that must be waited on + Barrier(crossbeam_utils::sync::WaitGroup), /// Progressbar to report progress Progressbar(indicatif::ProgressBar), /// Clear and remove the progressbar @@ -760,6 +774,7 @@ impl DistributedSystemPart { MsgFromHost::Output(ntime) => self.output(ntime), MsgFromHost::Stop => return, MsgFromHost::Error => self.send(MsgToHost::Error(self.error())).unwrap(), + MsgFromHost::Barrier(barrier) => barrier.wait(), MsgFromHost::Progressbar(pbar) => self.progressbar = Some(pbar), MsgFromHost::ProgressbarDrop => { let pb = self.progressbar.take().unwrap(); From e7222a99b54363ade4cf14bc50f75ebf8bacd3e6 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Wed, 18 Aug 2021 12:39:08 +0000 Subject: [PATCH 20/25] Shorten thread name --- multigrid/src/system.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index f756c3c..5bb2d8b 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -224,7 +224,7 @@ impl BaseSystem { .zip(push_channels) .enumerate() { - let builder = std::thread::Builder::new().name(format!("eulersolver: {}", name)); + let builder = std::thread::Builder::new().name(format!("mg: {}", name)); let boundary_conditions = bt.zip(chan).map(|(bt, chan)| match bt { euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, From d2c811d3af009bb907dd20229278419c28cedf4d Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Fri, 20 Aug 2021 15:57:12 +0000 Subject: [PATCH 21/25] Rework wait primitive to condvar --- euler/src/lib.rs | 53 ++-- multigrid/Cargo.toml | 2 + multigrid/src/system.rs | 685 ++++++++++++++++++++++------------------ sbp/src/utils.rs | 19 ++ 4 files changed, 421 insertions(+), 338 deletions(-) diff --git a/euler/src/lib.rs b/euler/src/lib.rs index f712fb1..e4cf538 100644 --- a/euler/src/lib.rs +++ b/euler/src/lib.rs @@ -301,22 +301,18 @@ impl Field { pub fn west(&self) -> ArrayView2 { self.slice(s![.., .., 0]) } - #[allow(unused)] - fn north_mut(&mut self) -> ArrayViewMut2 { + pub fn north_mut(&mut self) -> ArrayViewMut2 { let ny = self.ny(); self.slice_mut(s![.., ny - 1, ..]) } - #[allow(unused)] - fn south_mut(&mut self) -> ArrayViewMut2 { + pub fn south_mut(&mut self) -> ArrayViewMut2 { self.slice_mut(s![.., 0, ..]) } - #[allow(unused)] - fn east_mut(&mut self) -> ArrayViewMut2 { + pub fn east_mut(&mut self) -> ArrayViewMut2 { let nx = self.nx(); self.slice_mut(s![.., .., nx - 1]) } - #[allow(unused)] - fn west_mut(&mut self) -> ArrayViewMut2 { + pub fn west_mut(&mut self) -> ArrayViewMut2 { self.slice_mut(s![.., .., 0]) } @@ -463,18 +459,18 @@ impl Diff { pub fn zeros((ny, nx): (usize, usize)) -> Self { Self(Array3::zeros((4, ny, nx))) } - fn north_mut(&mut self) -> ArrayViewMut2 { + pub fn north_mut(&mut self) -> ArrayViewMut2 { let ny = self.shape()[1]; self.0.slice_mut(s![.., ny - 1, ..]) } - fn south_mut(&mut self) -> ArrayViewMut2 { + pub fn south_mut(&mut self) -> ArrayViewMut2 { self.0.slice_mut(s![.., 0, ..]) } - fn east_mut(&mut self) -> ArrayViewMut2 { + pub fn east_mut(&mut self) -> ArrayViewMut2 { let nx = self.shape()[2]; self.0.slice_mut(s![.., .., nx - 1]) } - fn west_mut(&mut self) -> ArrayViewMut2 { + pub fn west_mut(&mut self) -> ArrayViewMut2 { self.0.slice_mut(s![.., .., 0]) } } @@ -1010,16 +1006,25 @@ pub fn SAT_characteristics( 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); + SAT_north(op, k.north_mut(), y, metrics, boundaries.north); + SAT_south(op, k.south_mut(), y, metrics, boundaries.south); + SAT_east(op, k.east_mut(), y, metrics, boundaries.east); + SAT_west(op, k.west_mut(), y, metrics, boundaries.west); } +pub const SAT_FUNCTIONS: Direction< + fn(&dyn SbpOperator2d, ArrayViewMut2, &Field, &Metrics, ArrayView2), +> = Direction { + north: SAT_north, + south: SAT_south, + west: SAT_west, + east: SAT_east, +}; + #[allow(non_snake_case)] pub fn SAT_north( op: &dyn SbpOperator2d, - k: &mut Diff, + k: ArrayViewMut2, y: &Field, metrics: &Metrics, boundary: ArrayView2, @@ -1035,7 +1040,7 @@ pub fn SAT_north( let tau = 1.0; let slice = s![y.ny() - 1, ..]; SAT_characteristic( - k.north_mut(), + k, y.north(), boundary, hi, @@ -1050,7 +1055,7 @@ pub fn SAT_north( #[allow(non_snake_case)] pub fn SAT_south( op: &dyn SbpOperator2d, - k: &mut Diff, + k: ArrayViewMut2, y: &Field, metrics: &Metrics, boundary: ArrayView2, @@ -1065,7 +1070,7 @@ pub fn SAT_south( let tau = -1.0; let slice = s![0, ..]; SAT_characteristic( - k.south_mut(), + k, y.south(), boundary, hi, @@ -1080,7 +1085,7 @@ pub fn SAT_south( #[allow(non_snake_case)] pub fn SAT_west( op: &dyn SbpOperator2d, - k: &mut Diff, + k: ArrayViewMut2, y: &Field, metrics: &Metrics, boundary: ArrayView2, @@ -1096,7 +1101,7 @@ pub fn SAT_west( let tau = -1.0; let slice = s![.., 0]; SAT_characteristic( - k.west_mut(), + k, y.west(), boundary, hi, @@ -1111,7 +1116,7 @@ pub fn SAT_west( #[allow(non_snake_case)] pub fn SAT_east( op: &dyn SbpOperator2d, - k: &mut Diff, + k: ArrayViewMut2, y: &Field, metrics: &Metrics, boundary: ArrayView2, @@ -1126,7 +1131,7 @@ pub fn SAT_east( let tau = 1.0; let slice = s![.., y.nx() - 1]; SAT_characteristic( - k.east_mut(), + k, y.east(), boundary, hi, diff --git a/multigrid/Cargo.toml b/multigrid/Cargo.toml index 112eaa6..66b7177 100644 --- a/multigrid/Cargo.toml +++ b/multigrid/Cargo.toml @@ -20,3 +20,5 @@ argh = "0.1.4" evalexpr = "6.3.0" crossbeam-channel = "0.5.0" crossbeam-utils = "0.8.5" +parking_lot = "0.11.1" +lock_api = "0.4.4" diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 5bb2d8b..fdb9e81 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -1,15 +1,17 @@ use crate::parsing; use crate::utils::Direction; use core::ops::Deref; -use crossbeam_channel::{Receiver, Select, Sender}; +use crossbeam_channel::{Receiver, Sender}; use euler::{ eval::{self, Evaluator}, Diff, Field, VortexParameters, WorkBuffers, }; use ndarray::Array2; +use parking_lot::{Condvar, Mutex}; use sbp::grid::{Grid, Metrics}; use sbp::operators::{InterpolationOperator, SbpOperator2d}; use sbp::*; +use std::sync::Arc; pub struct BaseSystem { pub names: Vec, @@ -157,56 +159,37 @@ impl BaseSystem { let nthreads = self.grids.len(); // Build up the boundary conditions - let mut push_channels: Vec>>>> = - Vec::with_capacity(nthreads); - let mut pull_channels: Vec>>>> = - vec![Direction::default(); nthreads]; + let mut pull = Vec::>::with_capacity(nthreads); + for wb in &self.boundary_conditions { + let data = wb.as_ref().map(|bc| { + if let euler::BoundaryCharacteristic::Grid(_) + | euler::BoundaryCharacteristic::Interpolate(_, _) = bc + { + CommunicatorData::NotAvailable + } else { + CommunicatorData::NotApplicable + } + }); + pull.push(Arc::new(Communicator { + cvar: Condvar::new(), + data: Mutex::new(data), + })); + } // Build the set of communicators between boundaries + let mut push = Vec::>>>::with_capacity(nthreads); for wb in &self.boundary_conditions { - let mut local_push = Direction::default(); - if let euler::BoundaryCharacteristic::Grid(i) - | euler::BoundaryCharacteristic::Interpolate(i, _) = wb.north() - { - let (s, r) = crossbeam_channel::bounded(1); - 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) - .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) - .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) - .and_then::<(), _>(|_| panic!("channel is already present")); - *local_push.west_mut() = Some(s); - } + let local_push = wb.as_ref().map(|bc| { + if let euler::BoundaryCharacteristic::Grid(i) + | euler::BoundaryCharacteristic::Interpolate(i, _) = bc + { + Some(pull[*i].clone()) + } else { + None + } + }); - push_channels.push(local_push); + push.push(local_push); } let (master_send, master_recv) = crossbeam_channel::unbounded(); @@ -214,25 +197,23 @@ impl BaseSystem { let mut tids = Vec::with_capacity(nthreads); let mut communicators = Vec::with_capacity(nthreads); - for (id, (((((name, grid), sbp), bt), chan), push)) in self + for (id, (((((name, grid), sbp), bt), pull), push)) in self .names .into_iter() .zip(self.grids.into_iter()) .zip(self.operators.into_iter()) .zip(self.boundary_conditions) - .zip(pull_channels) - .zip(push_channels) + .zip(pull) + .zip(push) .enumerate() { let builder = std::thread::Builder::new().name(format!("mg: {}", name)); - let boundary_conditions = bt.zip(chan).map(|(bt, chan)| match bt { + let boundary_conditions = bt.map(|bt| match bt { euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This, - euler::BoundaryCharacteristic::Grid(_) => { - DistributedBoundaryConditions::Channel(chan.unwrap()) - } + euler::BoundaryCharacteristic::Grid(_) => DistributedBoundaryConditions::Channel, euler::BoundaryCharacteristic::Interpolate(_, int_op) => { - DistributedBoundaryConditions::Interpolate(chan.unwrap(), int_op) + DistributedBoundaryConditions::Interpolate(int_op) } euler::BoundaryCharacteristic::MultiGrid(_) => unimplemented!(), euler::BoundaryCharacteristic::Vortex(vp) => { @@ -319,6 +300,7 @@ impl BaseSystem { grid: (grid, metrics), output: g, push, + pull, sbp, t: time, dt: Float::NAN, @@ -331,8 +313,20 @@ impl BaseSystem { send: master_send, wb, - wb_ns: Array2::zeros((4, nx)), - wb_ew: Array2::zeros((4, ny)), + workbuffer_edges: ( + Direction { + north: Array2::zeros((4, nx)), + south: Array2::zeros((4, nx)), + east: Array2::zeros((4, ny)), + west: Array2::zeros((4, ny)), + }, + Direction { + north: Some(Array2::zeros((4, nx))), + south: Some(Array2::zeros((4, nx))), + east: Some(Array2::zeros((4, ny))), + west: Some(Array2::zeros((4, ny))), + }, + ), progressbar: None, }; @@ -723,19 +717,34 @@ pub enum DistributedBoundaryConditions { Vortex(VortexParameters), Eval(std::sync::Arc>), - Interpolate(Receiver>, Box), - Channel(Receiver>), + Interpolate(Box), + Channel, } -type PushCommunicator = Option>>; +enum CommunicatorData { + NotApplicable, + NotAvailable, + Some(Array2), +} + +struct Communicator { + /// Waker for this grid, neighbours should have a reference + /// and notify when a boundary has been put + cvar: Condvar, + /// Internal data exchange, is None on missing data, inner type + /// can be set to None when grabbing the boundary + data: Mutex>, +} struct DistributedSystemPart { grid: (Grid, Metrics), sbp: Box, boundary_conditions: Direction, + /// Channel pullers + pull: Arc, /// Subscribers to the boundaries of self - push: Direction, + push: Direction>>, current: Field, fut: Field, @@ -753,10 +762,13 @@ struct DistributedSystemPart { k: [Diff; 4], wb: WorkBuffers, - /// Work buffer for north/south boundary - wb_ns: Array2, - /// Work buffer for east/west boundary - wb_ew: Array2, + /// Work buffer for boundaries + /// + // Option: This can be sent from the current thread to another, + // Will be replenished by arriving boundary conditions for + // zero-allocation in loop (no global locks). + // Should never be None on entry to loop + workbuffer_edges: (Direction>, Direction>>), progressbar: Option, } @@ -857,25 +869,42 @@ impl DistributedSystemPart { let wb = &mut self.wb.0; let sbp = &self.sbp; let push = &self.push; + let pull = &self.pull; let boundary_conditions = &self.boundary_conditions; let grid = &self.grid.0; - let wb_ns = &mut self.wb_ns; - let wb_ew = &mut self.wb_ew; + let workbuffer_edges = &mut self.workbuffer_edges; 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() - } - if let Some(s) = &push.south { - s.send(y.south().to_owned()).unwrap() - } - if let Some(s) = &push.east { - s.send(y.east().to_owned()).unwrap() - } - if let Some(s) = &push.west { - s.send(y.west().to_owned()).unwrap() - } + // Send off the boundaries eagerly, in case neighbouring grid is ready + push.as_ref() + .zip(workbuffer_edges.1.as_mut()) + .zip( + Direction::) -> &mut CommunicatorData> { + north: |x| x.south_mut(), + south: |x| x.north_mut(), + east: |x| x.west_mut(), + west: |x| x.east_mut(), + }, + ) + .zip(Direction { + north: y.north(), + south: y.south(), + west: y.west(), + east: y.east(), + }) + .map(|(((push, wb), sel), this)| { + if let Some(s) = push { + let mut wb = wb.take().unwrap(); + wb.assign(&this); + { + let mut s = s.data.lock(); + let none = + std::mem::replace(sel(&mut s), CommunicatorData::Some(wb)); + assert!(matches!(none, CommunicatorData::NotAvailable)); + } + s.cvar.notify_one(); + } + }); // This computation does not depend on the boundaries euler::RHS_no_SAT(sbp.deref(), k, y, metrics, wb); @@ -883,249 +912,277 @@ impl DistributedSystemPart { // 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 - let mut select = Select::new(); - let mut selectable = 0; - let recv_north = match boundary_conditions.north() { - DistributedBoundaryConditions::Channel(r) - | DistributedBoundaryConditions::Interpolate(r, _) => { - selectable += 1; - Some(select.recv(r)) - } - DistributedBoundaryConditions::This => { - euler::SAT_north(sbp.deref(), k, y, metrics, y.south()); - None - } - DistributedBoundaryConditions::Vortex(vp) => { - let mut fiter = wb_ns.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(); - vp.evaluate(time, gx, gy, rho, rhou, rhov, e); - - euler::SAT_north(sbp.deref(), k, y, metrics, wb_ns.view()); - None - } - DistributedBoundaryConditions::Eval(eval) => { - let mut fiter = wb_ns.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, wb_ns.view()); - None - } - }; - let recv_south = match boundary_conditions.south() { - DistributedBoundaryConditions::Channel(r) - | DistributedBoundaryConditions::Interpolate(r, _) => { - selectable += 1; - Some(select.recv(r)) - } - DistributedBoundaryConditions::This => { - euler::SAT_south(sbp.deref(), k, y, metrics, y.north()); - None - } - DistributedBoundaryConditions::Vortex(vp) => { - let mut fiter = wb_ns.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, wb_ns.view()); - None - } - DistributedBoundaryConditions::Eval(eval) => { - let mut fiter = wb_ns.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, wb_ns.view()); - None - } - }; - let recv_east = match boundary_conditions.east() { - DistributedBoundaryConditions::Channel(r) - | DistributedBoundaryConditions::Interpolate(r, _) => { - selectable += 1; - Some(select.recv(r)) - } - DistributedBoundaryConditions::This => { - euler::SAT_east(sbp.deref(), k, y, metrics, y.west()); - None - } - DistributedBoundaryConditions::Vortex(vp) => { - let mut fiter = wb_ew.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, wb_ew.view()); - None - } - DistributedBoundaryConditions::Eval(eval) => { - let mut fiter = wb_ew.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, wb_ew.view()); - None - } - }; - let recv_west = match boundary_conditions.west() { - DistributedBoundaryConditions::Channel(r) - | DistributedBoundaryConditions::Interpolate(r, _) => { - selectable += 1; - Some(select.recv(r)) - } - DistributedBoundaryConditions::This => { - euler::SAT_west(sbp.deref(), k, y, metrics, y.east()); - None - } - DistributedBoundaryConditions::Vortex(vp) => { - let mut fiter = wb_ew.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, wb_ew.view()); - None - } - DistributedBoundaryConditions::Eval(eval) => { - let mut fiter = wb_ew.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, wb_ew.view()); - None - } - }; - - // Get an item off each channel, waiting minimally before processing that boundary. - // The waiting ensures other grids can be processed by the core in case of - // oversubscription (in case of a more grids than core scenario) - // This minimises the amount of time waiting on boundary conditions - while selectable != 0 { - let s = select.select(); - let sindex = s.index(); - match Some(sindex) { - x if x == recv_north => match boundary_conditions.north() { - DistributedBoundaryConditions::Channel(r) => { - let r = s.recv(r).unwrap(); - euler::SAT_north(sbp.deref(), k, y, metrics, r.view()); + let computed = boundary_conditions + .as_ref() + .zip(euler::SAT_FUNCTIONS) + .zip(workbuffer_edges.0.as_mut()) + .zip(workbuffer_edges.1.as_mut()) + .zip(Direction { + north: y.south(), + south: y.north(), + east: y.west(), + west: y.east(), + }) + .zip(Direction { + north: grid.north(), + south: grid.south(), + east: grid.east(), + west: grid.west(), + }) + .map(|(((((bc, sat), wb0), wb1), self_edge), grid)| { + wb0.fill(0.0); + match bc { + DistributedBoundaryConditions::Channel + | DistributedBoundaryConditions::Interpolate(_) => false, + DistributedBoundaryConditions::This => { + sat(sbp.deref(), wb0.view_mut(), y, metrics, self_edge); + true } - DistributedBoundaryConditions::Interpolate(r, int_op) => { - let r = s.recv(r).unwrap(); - let is_fine2coarse = r.shape()[1] > wb_ns.shape()[1]; - for (mut to, from) in wb_ns.outer_iter_mut().zip(r.outer_iter()) { - if is_fine2coarse { - int_op.fine2coarse(from.view(), to.view_mut()); - } else { - int_op.coarse2fine(from.view(), to.view_mut()); + DistributedBoundaryConditions::Vortex(vp) => { + let wb1 = wb1.as_mut().unwrap(); + let mut fiter = wb1.outer_iter_mut(); + let (rho, rhou, rhov, e) = ( + fiter.next().unwrap(), + fiter.next().unwrap(), + fiter.next().unwrap(), + fiter.next().unwrap(), + ); + let (gx, gy) = grid; + vp.evaluate(time, gx, gy, rho, rhou, rhov, e); + + sat(sbp.deref(), wb0.view_mut(), y, metrics, wb1.view()); + true + } + DistributedBoundaryConditions::Eval(eval) => { + let wb1 = wb1.as_mut().unwrap(); + let mut fiter = wb1.outer_iter_mut(); + let (rho, rhou, rhov, e) = ( + fiter.next().unwrap(), + fiter.next().unwrap(), + fiter.next().unwrap(), + fiter.next().unwrap(), + ); + let (gx, gy) = grid; + eval.evaluate(time, gx, gy, rho, rhou, rhov, e); + sat(sbp.deref(), wb0.view_mut(), y, metrics, wb1.view()); + true + } + } + }); + + if computed.north { + k.north_mut().scaled_add(1.0, &workbuffer_edges.0.north()); + } + if computed.south { + k.south_mut().scaled_add(1.0, &workbuffer_edges.0.south()); + } + if computed.east { + k.east_mut().scaled_add(1.0, &workbuffer_edges.0.east()); + } + if computed.west { + k.west_mut().scaled_add(1.0, &workbuffer_edges.0.west()); + } + + let mut boundaries_remaining = computed.map(|b| !b); + + { + let mut data = pull.data.lock(); + 'check_boundaries: loop { + if boundaries_remaining.north { + if let CommunicatorData::Some(boundary) = std::mem::replace( + &mut (*data).north, + CommunicatorData::NotAvailable, + ) { + lock_api::MutexGuard::unlocked(&mut data, || { + let wb0 = workbuffer_edges.0.north_mut(); + let wb1 = workbuffer_edges.1.north.insert(boundary); + match boundary_conditions.north() { + DistributedBoundaryConditions::Channel => { + std::mem::swap(wb0, wb1); + } + DistributedBoundaryConditions::Interpolate(int_op) => { + let is_fine2coarse = wb1.shape()[1] > wb0.shape()[2]; + for (to, from) in + wb0.outer_iter_mut().zip(wb1.outer_iter()) + { + if is_fine2coarse { + int_op.fine2coarse(from, to); + } else { + int_op.coarse2fine(from, to); + } + } + // Reshape edge buffer to correct size + let wb = workbuffer_edges.1.north.take().unwrap(); + let mut vec = wb.into_raw_vec(); + vec.resize(wb0.len(), 0.0); + let wb = + Array2::from_shape_vec(wb0.raw_dim(), vec).unwrap(); + workbuffer_edges.1.north = Some(wb); + } + _ => unreachable!(), } + euler::SAT_north( + sbp.deref(), + k.north_mut(), + y, + metrics, + wb0.view(), + ); + boundaries_remaining.north = false; + }); + continue 'check_boundaries; + } + + if boundaries_remaining.south { + if let CommunicatorData::Some(boundary) = std::mem::replace( + &mut (*data).south, + CommunicatorData::NotAvailable, + ) { + lock_api::MutexGuard::unlocked(&mut data, || { + let wb0 = workbuffer_edges.0.south_mut(); + let wb1 = workbuffer_edges.1.south.insert(boundary); + match boundary_conditions.south() { + DistributedBoundaryConditions::Channel => { + std::mem::swap(wb0, wb1); + } + DistributedBoundaryConditions::Interpolate(int_op) => { + let is_fine2coarse = + wb1.shape()[1] > wb0.shape()[2]; + for (to, from) in + wb0.outer_iter_mut().zip(wb1.outer_iter()) + { + if is_fine2coarse { + int_op.fine2coarse(from, to); + } else { + int_op.coarse2fine(from, to); + } + } + // Reshape edge buffer to correct size + let wb = workbuffer_edges.1.south.take().unwrap(); + let mut vec = wb.into_raw_vec(); + vec.resize(wb0.len(), 0.0); + let wb = Array2::from_shape_vec(wb0.raw_dim(), vec) + .unwrap(); + workbuffer_edges.1.south = Some(wb); + } + _ => unreachable!(), + } + euler::SAT_south( + sbp.deref(), + k.south_mut(), + y, + metrics, + wb0.view(), + ); + boundaries_remaining.south = false; + }); + continue 'check_boundaries; } - euler::SAT_north(sbp.deref(), k, y, metrics, wb_ns.view()); } - _ => unreachable!(), - }, - x if x == recv_south => match boundary_conditions.south() { - DistributedBoundaryConditions::Channel(r) => { - let r = s.recv(r).unwrap(); - euler::SAT_south(sbp.deref(), k, y, metrics, r.view()); - } - DistributedBoundaryConditions::Interpolate(r, int_op) => { - let r = s.recv(r).unwrap(); - let is_fine2coarse = r.shape()[1] > wb_ns.shape()[1]; - for (mut to, from) in wb_ns.outer_iter_mut().zip(r.outer_iter()) { - if is_fine2coarse { - int_op.fine2coarse(from.view(), to.view_mut()); - } else { - int_op.coarse2fine(from.view(), to.view_mut()); - } + + if boundaries_remaining.east { + if let CommunicatorData::Some(boundary) = std::mem::replace( + &mut (*data).east, + CommunicatorData::NotAvailable, + ) { + lock_api::MutexGuard::unlocked(&mut data, || { + let wb0 = workbuffer_edges.0.east_mut(); + let wb1 = workbuffer_edges.1.east.insert(boundary); + match boundary_conditions.east() { + DistributedBoundaryConditions::Channel => { + std::mem::swap(wb0, wb1); + } + DistributedBoundaryConditions::Interpolate(int_op) => { + let is_fine2coarse = + wb1.shape()[1] > wb0.shape()[2]; + for (to, from) in + wb0.outer_iter_mut().zip(wb1.outer_iter()) + { + if is_fine2coarse { + int_op.fine2coarse(from, to); + } else { + int_op.coarse2fine(from, to); + } + } + // Reshape edge buffer to correct size + let wb = workbuffer_edges.1.east.take().unwrap(); + let mut vec = wb.into_raw_vec(); + vec.resize(wb0.len(), 0.0); + let wb = Array2::from_shape_vec(wb0.raw_dim(), vec) + .unwrap(); + workbuffer_edges.1.east = Some(wb); + } + _ => unreachable!(), + } + euler::SAT_east( + sbp.deref(), + k.east_mut(), + y, + metrics, + wb0.view(), + ); + boundaries_remaining.east = false; + }); + continue 'check_boundaries; } - euler::SAT_south(sbp.deref(), k, y, metrics, wb_ns.view()); } - _ => unreachable!(), - }, - x if x == recv_west => match boundary_conditions.west() { - DistributedBoundaryConditions::Channel(r) => { - let r = s.recv(r).unwrap(); - euler::SAT_west(sbp.deref(), k, y, metrics, r.view()); - } - DistributedBoundaryConditions::Interpolate(r, int_op) => { - let r = s.recv(r).unwrap(); - let is_fine2coarse = r.shape()[1] > wb_ew.shape()[1]; - for (mut to, from) in wb_ew.outer_iter_mut().zip(r.outer_iter()) { - if is_fine2coarse { - int_op.fine2coarse(from.view(), to.view_mut()); - } else { - int_op.coarse2fine(from.view(), to.view_mut()); - } + + if boundaries_remaining.west { + if let CommunicatorData::Some(boundary) = std::mem::replace( + &mut (*data).west, + CommunicatorData::NotAvailable, + ) { + lock_api::MutexGuard::unlocked(&mut data, || { + let wb0 = workbuffer_edges.0.west_mut(); + let wb1 = workbuffer_edges.1.west.insert(boundary); + match boundary_conditions.west() { + DistributedBoundaryConditions::Channel => { + std::mem::swap(wb0, wb1); + } + DistributedBoundaryConditions::Interpolate(int_op) => { + let is_fine2coarse = + wb1.shape()[1] > wb0.shape()[2]; + for (to, from) in + wb0.outer_iter_mut().zip(wb1.outer_iter()) + { + if is_fine2coarse { + int_op.fine2coarse(from, to); + } else { + int_op.coarse2fine(from, to); + } + } + // Reshape edge buffer to correct size + let wb = workbuffer_edges.1.west.take().unwrap(); + let mut vec = wb.into_raw_vec(); + vec.resize(wb0.len(), 0.0); + let wb = Array2::from_shape_vec(wb0.raw_dim(), vec) + .unwrap(); + workbuffer_edges.1.west = Some(wb); + } + _ => unreachable!(), + } + euler::SAT_west( + sbp.deref(), + k.west_mut(), + y, + metrics, + wb0.view(), + ); + boundaries_remaining.west = false; + }); + continue 'check_boundaries; } - euler::SAT_west(sbp.deref(), k, y, metrics, wb_ew.view()); } - _ => unreachable!(), - }, - x if x == recv_east => match boundary_conditions.east() { - DistributedBoundaryConditions::Channel(r) => { - let r = s.recv(r).unwrap(); - euler::SAT_east(sbp.deref(), k, y, metrics, r.view()); - } - DistributedBoundaryConditions::Interpolate(r, int_op) => { - let r = s.recv(r).unwrap(); - let is_fine2coarse = r.shape()[1] > wb_ew.shape()[1]; - for (mut to, from) in wb_ew.outer_iter_mut().zip(r.outer_iter()) { - if is_fine2coarse { - int_op.fine2coarse(from.view(), to.view_mut()); - } else { - int_op.coarse2fine(from.view(), to.view_mut()); - } - } - euler::SAT_east(sbp.deref(), k, y, metrics, wb_ew.view()); - } - _ => unreachable!(), - }, - _ => unreachable!(), + } + if !boundaries_remaining.any() { + break 'check_boundaries; + } + println!("{:?}", boundaries_remaining); + // We have no available boundaries yet, can wait until + // notified by other threads. Early continues + // ensures the boundary has not been notified earlier + pull.cvar.wait(&mut data); } - select.remove(sindex); - selectable -= 1; } }; integrate::integrate::( diff --git a/sbp/src/utils.rs b/sbp/src/utils.rs index 8d1f025..168a9ba 100644 --- a/sbp/src/utils.rs +++ b/sbp/src/utils.rs @@ -53,6 +53,16 @@ impl Direction { east: (self.east, other.east), } } + + /// Flips all direction through origo + pub fn opposite(self) -> Self { + Self { + north: self.south, + south: self.north, + east: self.west, + west: self.east, + } + } } impl Direction> { @@ -96,6 +106,15 @@ impl Direction { } } +impl Direction { + pub fn all(&self) -> bool { + self.north && self.south && self.east && self.west + } + pub fn any(&self) -> bool { + self.north || self.south || self.east || self.west + } +} + /// Linearly spaced parameters, apart from the boundaries which /// only have a distance of `h/2` from the boundary pub fn h2linspace(start: Float, end: Float, n: usize) -> ndarray::Array1 { From d0901f5755ed6e4f4f3c09c869d76b3da1fccfe8 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Sat, 21 Aug 2021 09:29:45 +0000 Subject: [PATCH 22/25] SAT boundaries for multi-thread fixing --- multigrid/Cargo.toml | 1 + multigrid/src/system.rs | 305 +++++++++++++++++----------------------- sbp/src/utils.rs | 12 ++ 3 files changed, 140 insertions(+), 178 deletions(-) diff --git a/multigrid/Cargo.toml b/multigrid/Cargo.toml index 66b7177..01dc0df 100644 --- a/multigrid/Cargo.toml +++ b/multigrid/Cargo.toml @@ -22,3 +22,4 @@ crossbeam-channel = "0.5.0" crossbeam-utils = "0.8.5" parking_lot = "0.11.1" lock_api = "0.4.4" +arrayvec = "0.7.1" diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index fdb9e81..c4c96cc 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -160,19 +160,10 @@ impl BaseSystem { // Build up the boundary conditions let mut pull = Vec::>::with_capacity(nthreads); - for wb in &self.boundary_conditions { - let data = wb.as_ref().map(|bc| { - if let euler::BoundaryCharacteristic::Grid(_) - | euler::BoundaryCharacteristic::Interpolate(_, _) = bc - { - CommunicatorData::NotAvailable - } else { - CommunicatorData::NotApplicable - } - }); + for _ in 0..nthreads { pull.push(Arc::new(Communicator { cvar: Condvar::new(), - data: Mutex::new(data), + data: Mutex::new(Direction::splat(()).map(|_| arrayvec::ArrayVec::new())), })); } @@ -721,11 +712,7 @@ pub enum DistributedBoundaryConditions { Channel, } -enum CommunicatorData { - NotApplicable, - NotAvailable, - Some(Array2), -} +type CommunicatorData = arrayvec::ArrayVec, 2>; struct Communicator { /// Waker for this grid, neighbours should have a reference @@ -898,9 +885,7 @@ impl DistributedSystemPart { wb.assign(&this); { let mut s = s.data.lock(); - let none = - std::mem::replace(sel(&mut s), CommunicatorData::Some(wb)); - assert!(matches!(none, CommunicatorData::NotAvailable)); + sel(&mut s).push(wb); } s.cvar.notify_one(); } @@ -987,16 +972,69 @@ impl DistributedSystemPart { { let mut data = pull.data.lock(); - 'check_boundaries: loop { - if boundaries_remaining.north { - if let CommunicatorData::Some(boundary) = std::mem::replace( - &mut (*data).north, - CommunicatorData::NotAvailable, - ) { - lock_api::MutexGuard::unlocked(&mut data, || { - let wb0 = workbuffer_edges.0.north_mut(); - let wb1 = workbuffer_edges.1.north.insert(boundary); - match boundary_conditions.north() { + 'check_boundaries: while boundaries_remaining.any() { + let boundaries = + boundaries_remaining + .zip(data.as_mut()) + .map( + |(remains, data)| { + if remains { + data.pop_at(0) + } else { + None + } + }, + ); + if boundaries.as_ref().map(Option::is_none).all() { + // Park thread while waiting for boundaries + pull.cvar.wait(&mut data); + continue 'check_boundaries; + } + // While we are waiting we can unlock mutex + lock_api::MutexGuard::unlocked(&mut data, || { + if let Some(boundary) = boundaries.north { + boundaries_remaining.north = false; + let wb0 = workbuffer_edges.0.north_mut(); + let wb1 = workbuffer_edges.1.north.insert(boundary); + match boundary_conditions.north() { + DistributedBoundaryConditions::Channel => { + std::mem::swap(wb0, wb1); + } + DistributedBoundaryConditions::Interpolate(int_op) => { + let is_fine2coarse = wb1.shape()[1] > wb0.shape()[2]; + for (to, from) in wb0.outer_iter_mut().zip(wb1.outer_iter()) + { + if is_fine2coarse { + int_op.fine2coarse(from, to); + } else { + int_op.coarse2fine(from, to); + } + } + // Reshape edge buffer to correct size + let wb = workbuffer_edges.1.north.take().unwrap(); + let mut vec = wb.into_raw_vec(); + vec.resize(wb0.len(), 0.0); + let wb = + Array2::from_shape_vec(wb0.raw_dim(), vec).unwrap(); + workbuffer_edges.1.north = Some(wb); + } + _ => unreachable!(), + } + euler::SAT_north( + sbp.deref(), + k.north_mut(), + y, + metrics, + wb0.view(), + ); + }; + + if boundaries_remaining.south { + boundaries_remaining.south = false; + if let Some(boundary) = boundaries.south { + let wb0 = workbuffer_edges.0.south_mut(); + let wb1 = workbuffer_edges.1.south.insert(boundary); + match boundary_conditions.south() { DistributedBoundaryConditions::Channel => { std::mem::swap(wb0, wb1); } @@ -1012,176 +1050,87 @@ impl DistributedSystemPart { } } // Reshape edge buffer to correct size - let wb = workbuffer_edges.1.north.take().unwrap(); + let wb = workbuffer_edges.1.south.take().unwrap(); let mut vec = wb.into_raw_vec(); vec.resize(wb0.len(), 0.0); let wb = Array2::from_shape_vec(wb0.raw_dim(), vec).unwrap(); - workbuffer_edges.1.north = Some(wb); + workbuffer_edges.1.south = Some(wb); } _ => unreachable!(), } - euler::SAT_north( + euler::SAT_south( sbp.deref(), - k.north_mut(), + k.south_mut(), y, metrics, wb0.view(), ); - boundaries_remaining.north = false; - }); - continue 'check_boundaries; + }; } - if boundaries_remaining.south { - if let CommunicatorData::Some(boundary) = std::mem::replace( - &mut (*data).south, - CommunicatorData::NotAvailable, - ) { - lock_api::MutexGuard::unlocked(&mut data, || { - let wb0 = workbuffer_edges.0.south_mut(); - let wb1 = workbuffer_edges.1.south.insert(boundary); - match boundary_conditions.south() { - DistributedBoundaryConditions::Channel => { - std::mem::swap(wb0, wb1); + if let Some(boundary) = boundaries.east { + boundaries_remaining.east = false; + let wb0 = workbuffer_edges.0.east_mut(); + let wb1 = workbuffer_edges.1.east.insert(boundary); + match boundary_conditions.east() { + DistributedBoundaryConditions::Channel => { + std::mem::swap(wb0, wb1); + } + DistributedBoundaryConditions::Interpolate(int_op) => { + let is_fine2coarse = wb1.shape()[1] > wb0.shape()[2]; + for (to, from) in wb0.outer_iter_mut().zip(wb1.outer_iter()) + { + if is_fine2coarse { + int_op.fine2coarse(from, to); + } else { + int_op.coarse2fine(from, to); } - DistributedBoundaryConditions::Interpolate(int_op) => { - let is_fine2coarse = - wb1.shape()[1] > wb0.shape()[2]; - for (to, from) in - wb0.outer_iter_mut().zip(wb1.outer_iter()) - { - if is_fine2coarse { - int_op.fine2coarse(from, to); - } else { - int_op.coarse2fine(from, to); - } - } - // Reshape edge buffer to correct size - let wb = workbuffer_edges.1.south.take().unwrap(); - let mut vec = wb.into_raw_vec(); - vec.resize(wb0.len(), 0.0); - let wb = Array2::from_shape_vec(wb0.raw_dim(), vec) - .unwrap(); - workbuffer_edges.1.south = Some(wb); - } - _ => unreachable!(), } - euler::SAT_south( - sbp.deref(), - k.south_mut(), - y, - metrics, - wb0.view(), - ); - boundaries_remaining.south = false; - }); - continue 'check_boundaries; + // Reshape edge buffer to correct size + let wb = workbuffer_edges.1.east.take().unwrap(); + let mut vec = wb.into_raw_vec(); + vec.resize(wb0.len(), 0.0); + let wb = + Array2::from_shape_vec(wb0.raw_dim(), vec).unwrap(); + workbuffer_edges.1.east = Some(wb); + } + _ => unreachable!(), } - } + euler::SAT_east(sbp.deref(), k.east_mut(), y, metrics, wb0.view()); + }; - if boundaries_remaining.east { - if let CommunicatorData::Some(boundary) = std::mem::replace( - &mut (*data).east, - CommunicatorData::NotAvailable, - ) { - lock_api::MutexGuard::unlocked(&mut data, || { - let wb0 = workbuffer_edges.0.east_mut(); - let wb1 = workbuffer_edges.1.east.insert(boundary); - match boundary_conditions.east() { - DistributedBoundaryConditions::Channel => { - std::mem::swap(wb0, wb1); + if let Some(boundary) = boundaries.west { + boundaries_remaining.west = false; + let wb0 = workbuffer_edges.0.west_mut(); + let wb1 = workbuffer_edges.1.west.insert(boundary); + match boundary_conditions.west() { + DistributedBoundaryConditions::Channel => { + std::mem::swap(wb0, wb1); + } + DistributedBoundaryConditions::Interpolate(int_op) => { + let is_fine2coarse = wb1.shape()[1] > wb0.shape()[2]; + for (to, from) in wb0.outer_iter_mut().zip(wb1.outer_iter()) + { + if is_fine2coarse { + int_op.fine2coarse(from, to); + } else { + int_op.coarse2fine(from, to); } - DistributedBoundaryConditions::Interpolate(int_op) => { - let is_fine2coarse = - wb1.shape()[1] > wb0.shape()[2]; - for (to, from) in - wb0.outer_iter_mut().zip(wb1.outer_iter()) - { - if is_fine2coarse { - int_op.fine2coarse(from, to); - } else { - int_op.coarse2fine(from, to); - } - } - // Reshape edge buffer to correct size - let wb = workbuffer_edges.1.east.take().unwrap(); - let mut vec = wb.into_raw_vec(); - vec.resize(wb0.len(), 0.0); - let wb = Array2::from_shape_vec(wb0.raw_dim(), vec) - .unwrap(); - workbuffer_edges.1.east = Some(wb); - } - _ => unreachable!(), } - euler::SAT_east( - sbp.deref(), - k.east_mut(), - y, - metrics, - wb0.view(), - ); - boundaries_remaining.east = false; - }); - continue 'check_boundaries; + // Reshape edge buffer to correct size + let wb = workbuffer_edges.1.west.take().unwrap(); + let mut vec = wb.into_raw_vec(); + vec.resize(wb0.len(), 0.0); + let wb = + Array2::from_shape_vec(wb0.raw_dim(), vec).unwrap(); + workbuffer_edges.1.west = Some(wb); + } + _ => unreachable!(), } - } - - if boundaries_remaining.west { - if let CommunicatorData::Some(boundary) = std::mem::replace( - &mut (*data).west, - CommunicatorData::NotAvailable, - ) { - lock_api::MutexGuard::unlocked(&mut data, || { - let wb0 = workbuffer_edges.0.west_mut(); - let wb1 = workbuffer_edges.1.west.insert(boundary); - match boundary_conditions.west() { - DistributedBoundaryConditions::Channel => { - std::mem::swap(wb0, wb1); - } - DistributedBoundaryConditions::Interpolate(int_op) => { - let is_fine2coarse = - wb1.shape()[1] > wb0.shape()[2]; - for (to, from) in - wb0.outer_iter_mut().zip(wb1.outer_iter()) - { - if is_fine2coarse { - int_op.fine2coarse(from, to); - } else { - int_op.coarse2fine(from, to); - } - } - // Reshape edge buffer to correct size - let wb = workbuffer_edges.1.west.take().unwrap(); - let mut vec = wb.into_raw_vec(); - vec.resize(wb0.len(), 0.0); - let wb = Array2::from_shape_vec(wb0.raw_dim(), vec) - .unwrap(); - workbuffer_edges.1.west = Some(wb); - } - _ => unreachable!(), - } - euler::SAT_west( - sbp.deref(), - k.west_mut(), - y, - metrics, - wb0.view(), - ); - boundaries_remaining.west = false; - }); - continue 'check_boundaries; - } - } - } - if !boundaries_remaining.any() { - break 'check_boundaries; - } - println!("{:?}", boundaries_remaining); - // We have no available boundaries yet, can wait until - // notified by other threads. Early continues - // ensures the boundary has not been notified earlier - pull.cvar.wait(&mut data); + euler::SAT_west(sbp.deref(), k.west_mut(), y, metrics, wb0.view()); + }; + }); } } }; diff --git a/sbp/src/utils.rs b/sbp/src/utils.rs index 168a9ba..bd87a7b 100644 --- a/sbp/src/utils.rs +++ b/sbp/src/utils.rs @@ -63,6 +63,18 @@ impl Direction { west: self.east, } } + + pub fn splat(t: T) -> Self + where + T: Copy, + { + Self { + north: t, + south: t, + east: t, + west: t, + } + } } impl Direction> { From 44e0eb98f36dbc22f513480a03bbd917cffddde5 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Fri, 24 Sep 2021 17:02:48 +0000 Subject: [PATCH 23/25] checkpoint --- multigrid/src/system.rs | 155 ++++++++++++++++++++++++---------------- 1 file changed, 92 insertions(+), 63 deletions(-) diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index c4c96cc..4ac66da 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -1,5 +1,6 @@ use crate::parsing; use crate::utils::Direction; +use arrayvec::ArrayVec; use core::ops::Deref; use crossbeam_channel::{Receiver, Sender}; use euler::{ @@ -163,7 +164,7 @@ impl BaseSystem { for _ in 0..nthreads { pull.push(Arc::new(Communicator { cvar: Condvar::new(), - data: Mutex::new(Direction::splat(()).map(|_| arrayvec::ArrayVec::new())), + data: Mutex::new(Direction::splat(()).map(|_| ArrayVec::new())), })); } @@ -304,20 +305,44 @@ impl BaseSystem { send: master_send, wb, - workbuffer_edges: ( + workbuffer_edges: { Direction { - north: Array2::zeros((4, nx)), - south: Array2::zeros((4, nx)), - east: Array2::zeros((4, ny)), - west: Array2::zeros((4, ny)), + north: (Array2::zeros((4, nx)), Array2::zeros((4, nx))), + south: (Array2::zeros((4, nx)), Array2::zeros((4, nx))), + east: (Array2::zeros((4, ny)), Array2::zeros((4, ny))), + west: (Array2::zeros((4, ny)), Array2::zeros((4, ny))), + } + }, + workbuffer_free: Direction { + north: { + let mut arr = ArrayVec::new(); + for _ in 0..2 { + arr.push(Array2::zeros((4, nx))) + } + arr }, - Direction { - north: Some(Array2::zeros((4, nx))), - south: Some(Array2::zeros((4, nx))), - east: Some(Array2::zeros((4, ny))), - west: Some(Array2::zeros((4, ny))), + south: { + let mut arr = ArrayVec::new(); + for _ in 0..2 { + arr.push(Array2::zeros((4, nx))) + } + arr }, - ), + east: { + let mut arr = ArrayVec::new(); + for _ in 0..2 { + arr.push(Array2::zeros((4, ny))) + } + arr + }, + west: { + let mut arr = ArrayVec::new(); + for _ in 0..2 { + arr.push(Array2::zeros((4, ny))) + } + arr + }, + }, progressbar: None, }; @@ -712,7 +737,7 @@ pub enum DistributedBoundaryConditions { Channel, } -type CommunicatorData = arrayvec::ArrayVec, 2>; +type CommunicatorData = ArrayVec, 2>; struct Communicator { /// Waker for this grid, neighbours should have a reference @@ -750,12 +775,9 @@ struct DistributedSystemPart { k: [Diff; 4], wb: WorkBuffers, /// Work buffer for boundaries - /// - // Option: This can be sent from the current thread to another, - // Will be replenished by arriving boundary conditions for - // zero-allocation in loop (no global locks). - // Should never be None on entry to loop - workbuffer_edges: (Direction>, Direction>>), + workbuffer_edges: Direction<(Array2, Array2)>, + /// These can be popped and pushed as we communicate data + workbuffer_free: Direction, progressbar: Option, } @@ -860,11 +882,12 @@ impl DistributedSystemPart { let boundary_conditions = &self.boundary_conditions; let grid = &self.grid.0; let workbuffer_edges = &mut self.workbuffer_edges; + let workbuffer_free = &mut self.workbuffer_free; let rhs = |k: &mut euler::Diff, y: &euler::Field, time: Float| { // Send off the boundaries eagerly, in case neighbouring grid is ready push.as_ref() - .zip(workbuffer_edges.1.as_mut()) + .zip(workbuffer_free.as_mut()) .zip( Direction::) -> &mut CommunicatorData> { north: |x| x.south_mut(), @@ -881,7 +904,7 @@ impl DistributedSystemPart { }) .map(|(((push, wb), sel), this)| { if let Some(s) = push { - let mut wb = wb.take().unwrap(); + let mut wb = wb.pop().unwrap(); wb.assign(&this); { let mut s = s.data.lock(); @@ -900,8 +923,7 @@ impl DistributedSystemPart { let computed = boundary_conditions .as_ref() .zip(euler::SAT_FUNCTIONS) - .zip(workbuffer_edges.0.as_mut()) - .zip(workbuffer_edges.1.as_mut()) + .zip(workbuffer_edges.as_mut()) .zip(Direction { north: y.south(), south: y.north(), @@ -914,18 +936,17 @@ impl DistributedSystemPart { east: grid.east(), west: grid.west(), }) - .map(|(((((bc, sat), wb0), wb1), self_edge), grid)| { - wb0.fill(0.0); + .map(|((((bc, sat), wb), self_edge), grid)| { + wb.0.fill(0.0); match bc { DistributedBoundaryConditions::Channel | DistributedBoundaryConditions::Interpolate(_) => false, DistributedBoundaryConditions::This => { - sat(sbp.deref(), wb0.view_mut(), y, metrics, self_edge); + sat(sbp.deref(), wb.0.view_mut(), y, metrics, self_edge); true } DistributedBoundaryConditions::Vortex(vp) => { - let wb1 = wb1.as_mut().unwrap(); - let mut fiter = wb1.outer_iter_mut(); + let mut fiter = wb.1.outer_iter_mut(); let (rho, rhou, rhov, e) = ( fiter.next().unwrap(), fiter.next().unwrap(), @@ -935,12 +956,11 @@ impl DistributedSystemPart { let (gx, gy) = grid; vp.evaluate(time, gx, gy, rho, rhou, rhov, e); - sat(sbp.deref(), wb0.view_mut(), y, metrics, wb1.view()); + sat(sbp.deref(), wb.0.view_mut(), y, metrics, wb.1.view()); true } DistributedBoundaryConditions::Eval(eval) => { - let wb1 = wb1.as_mut().unwrap(); - let mut fiter = wb1.outer_iter_mut(); + let mut fiter = wb.1.outer_iter_mut(); let (rho, rhou, rhov, e) = ( fiter.next().unwrap(), fiter.next().unwrap(), @@ -949,23 +969,27 @@ impl DistributedSystemPart { ); let (gx, gy) = grid; eval.evaluate(time, gx, gy, rho, rhou, rhov, e); - sat(sbp.deref(), wb0.view_mut(), y, metrics, wb1.view()); + sat(sbp.deref(), wb.0.view_mut(), y, metrics, wb.1.view()); true } } }); if computed.north { - k.north_mut().scaled_add(1.0, &workbuffer_edges.0.north()); + k.north_mut() + .scaled_add(1.0, &workbuffer_edges.north().0.view()); } if computed.south { - k.south_mut().scaled_add(1.0, &workbuffer_edges.0.south()); + k.south_mut() + .scaled_add(1.0, &workbuffer_edges.south().0.view()); } if computed.east { - k.east_mut().scaled_add(1.0, &workbuffer_edges.0.east()); + k.east_mut() + .scaled_add(1.0, &workbuffer_edges.east().0.view()); } if computed.west { - k.west_mut().scaled_add(1.0, &workbuffer_edges.0.west()); + k.west_mut() + .scaled_add(1.0, &workbuffer_edges.west().0.view()); } let mut boundaries_remaining = computed.map(|b| !b); @@ -994,15 +1018,17 @@ impl DistributedSystemPart { lock_api::MutexGuard::unlocked(&mut data, || { if let Some(boundary) = boundaries.north { boundaries_remaining.north = false; - let wb0 = workbuffer_edges.0.north_mut(); - let wb1 = workbuffer_edges.1.north.insert(boundary); + let wb = workbuffer_edges.north_mut(); + let wb_push = workbuffer_free.north_mut(); match boundary_conditions.north() { DistributedBoundaryConditions::Channel => { - std::mem::swap(wb0, wb1); + std::mem::swap(&mut wb.0, &mut boundary); + wb_push.push(boundary); } DistributedBoundaryConditions::Interpolate(int_op) => { - let is_fine2coarse = wb1.shape()[1] > wb0.shape()[2]; - for (to, from) in wb0.outer_iter_mut().zip(wb1.outer_iter()) + let is_fine2coarse = boundary.shape()[1] > wb.0.shape()[2]; + for (to, from) in + boundary.outer_iter_mut().zip(boundary.outer_iter()) { if is_fine2coarse { int_op.fine2coarse(from, to); @@ -1011,12 +1037,11 @@ impl DistributedSystemPart { } } // Reshape edge buffer to correct size - let wb = workbuffer_edges.1.north.take().unwrap(); - let mut vec = wb.into_raw_vec(); - vec.resize(wb0.len(), 0.0); - let wb = - Array2::from_shape_vec(wb0.raw_dim(), vec).unwrap(); - workbuffer_edges.1.north = Some(wb); + let mut vec = boundary.into_raw_vec(); + vec.resize(wb.0.len(), 0.0); + let boundary = + Array2::from_shape_vec(wb.0.raw_dim(), vec).unwrap(); + wb_push.push(boundary) } _ => unreachable!(), } @@ -1025,23 +1050,25 @@ impl DistributedSystemPart { k.north_mut(), y, metrics, - wb0.view(), + wb.0.view(), ); }; if boundaries_remaining.south { boundaries_remaining.south = false; if let Some(boundary) = boundaries.south { - let wb0 = workbuffer_edges.0.south_mut(); - let wb1 = workbuffer_edges.1.south.insert(boundary); + let wb = workbuffer_edges.north_mut(); + let wb_push = workbuffer_free.south_mut(); match boundary_conditions.south() { DistributedBoundaryConditions::Channel => { - std::mem::swap(wb0, wb1); + std::mem::swap(&mut wb.0, &mut boundary); + wb_push.push(boundary); } DistributedBoundaryConditions::Interpolate(int_op) => { - let is_fine2coarse = wb1.shape()[1] > wb0.shape()[2]; + let is_fine2coarse = + boundary.shape()[1] > wb.0.shape()[2]; for (to, from) in - wb0.outer_iter_mut().zip(wb1.outer_iter()) + wb.0.outer_iter_mut().zip(boundary.outer_iter()) { if is_fine2coarse { int_op.fine2coarse(from, to); @@ -1050,12 +1077,12 @@ impl DistributedSystemPart { } } // Reshape edge buffer to correct size - let wb = workbuffer_edges.1.south.take().unwrap(); - let mut vec = wb.into_raw_vec(); - vec.resize(wb0.len(), 0.0); - let wb = - Array2::from_shape_vec(wb0.raw_dim(), vec).unwrap(); - workbuffer_edges.1.south = Some(wb); + let mut vec = boundary.into_raw_vec(); + vec.resize(wb.0.len(), 0.0); + let boundary = + Array2::from_shape_vec(wb.0.raw_dim(), vec) + .unwrap(); + wb_push.push(boundary); } _ => unreachable!(), } @@ -1064,19 +1091,21 @@ impl DistributedSystemPart { k.south_mut(), y, metrics, - wb0.view(), + wb.0.view(), ); }; } if let Some(boundary) = boundaries.east { boundaries_remaining.east = false; - let wb0 = workbuffer_edges.0.east_mut(); - let wb1 = workbuffer_edges.1.east.insert(boundary); + let wb = workbuffer_edges.east_mut(); + let wb_push = workbuffer_free.east_mut(); match boundary_conditions.east() { DistributedBoundaryConditions::Channel => { - std::mem::swap(wb0, wb1); + std::mem::swap(&mut wb.0, &mut boundary); + wb_push.push(boundary); } + // TODO: From this point down DistributedBoundaryConditions::Interpolate(int_op) => { let is_fine2coarse = wb1.shape()[1] > wb0.shape()[2]; for (to, from) in wb0.outer_iter_mut().zip(wb1.outer_iter()) From 42bfdd1ca191d7c65c686b89e4a8d945c0398433 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Mon, 27 Sep 2021 20:28:07 +0000 Subject: [PATCH 24/25] checkpoint --- multigrid/src/system.rs | 128 +++++++++++++++++++--------------------- 1 file changed, 62 insertions(+), 66 deletions(-) diff --git a/multigrid/src/system.rs b/multigrid/src/system.rs index 4ac66da..d81fa06 100644 --- a/multigrid/src/system.rs +++ b/multigrid/src/system.rs @@ -737,7 +737,7 @@ pub enum DistributedBoundaryConditions { Channel, } -type CommunicatorData = ArrayVec, 2>; +type CommunicatorData = ArrayVec, 4>; struct Communicator { /// Waker for this grid, neighbours should have a reference @@ -1016,7 +1016,7 @@ impl DistributedSystemPart { } // While we are waiting we can unlock mutex lock_api::MutexGuard::unlocked(&mut data, || { - if let Some(boundary) = boundaries.north { + if let Some(mut boundary) = boundaries.north { boundaries_remaining.north = false; let wb = workbuffer_edges.north_mut(); let wb_push = workbuffer_free.north_mut(); @@ -1028,7 +1028,7 @@ impl DistributedSystemPart { DistributedBoundaryConditions::Interpolate(int_op) => { let is_fine2coarse = boundary.shape()[1] > wb.0.shape()[2]; for (to, from) in - boundary.outer_iter_mut().zip(boundary.outer_iter()) + wb.0.outer_iter_mut().zip(boundary.outer_iter()) { if is_fine2coarse { int_op.fine2coarse(from, to); @@ -1054,49 +1054,45 @@ impl DistributedSystemPart { ); }; - if boundaries_remaining.south { + if let Some(mut boundary) = boundaries.south { boundaries_remaining.south = false; - if let Some(boundary) = boundaries.south { - let wb = workbuffer_edges.north_mut(); - let wb_push = workbuffer_free.south_mut(); - match boundary_conditions.south() { - DistributedBoundaryConditions::Channel => { - std::mem::swap(&mut wb.0, &mut boundary); - wb_push.push(boundary); - } - DistributedBoundaryConditions::Interpolate(int_op) => { - let is_fine2coarse = - boundary.shape()[1] > wb.0.shape()[2]; - for (to, from) in - wb.0.outer_iter_mut().zip(boundary.outer_iter()) - { - if is_fine2coarse { - int_op.fine2coarse(from, to); - } else { - int_op.coarse2fine(from, to); - } - } - // Reshape edge buffer to correct size - let mut vec = boundary.into_raw_vec(); - vec.resize(wb.0.len(), 0.0); - let boundary = - Array2::from_shape_vec(wb.0.raw_dim(), vec) - .unwrap(); - wb_push.push(boundary); - } - _ => unreachable!(), + let wb = workbuffer_edges.north_mut(); + let wb_push = workbuffer_free.south_mut(); + match boundary_conditions.south() { + DistributedBoundaryConditions::Channel => { + std::mem::swap(&mut wb.0, &mut boundary); + wb_push.push(boundary); } - euler::SAT_south( - sbp.deref(), - k.south_mut(), - y, - metrics, - wb.0.view(), - ); - }; - } + DistributedBoundaryConditions::Interpolate(int_op) => { + let is_fine2coarse = boundary.shape()[1] > wb.0.shape()[2]; + for (to, from) in + wb.0.outer_iter_mut().zip(boundary.outer_iter()) + { + if is_fine2coarse { + int_op.fine2coarse(from, to); + } else { + int_op.coarse2fine(from, to); + } + } + // Reshape edge buffer to correct size + let mut vec = boundary.into_raw_vec(); + vec.resize(wb.0.len(), 0.0); + let boundary = + Array2::from_shape_vec(wb.0.raw_dim(), vec).unwrap(); + wb_push.push(boundary); + } + _ => unreachable!(), + } + euler::SAT_south( + sbp.deref(), + k.south_mut(), + y, + metrics, + wb.0.view(), + ); + }; - if let Some(boundary) = boundaries.east { + if let Some(mut boundary) = boundaries.east { boundaries_remaining.east = false; let wb = workbuffer_edges.east_mut(); let wb_push = workbuffer_free.east_mut(); @@ -1105,10 +1101,10 @@ impl DistributedSystemPart { std::mem::swap(&mut wb.0, &mut boundary); wb_push.push(boundary); } - // TODO: From this point down DistributedBoundaryConditions::Interpolate(int_op) => { - let is_fine2coarse = wb1.shape()[1] > wb0.shape()[2]; - for (to, from) in wb0.outer_iter_mut().zip(wb1.outer_iter()) + let is_fine2coarse = boundary.shape()[1] > wb.0.shape()[2]; + for (to, from) in + wb.0.outer_iter_mut().zip(boundary.outer_iter()) { if is_fine2coarse { int_op.fine2coarse(from, to); @@ -1117,29 +1113,30 @@ impl DistributedSystemPart { } } // Reshape edge buffer to correct size - let wb = workbuffer_edges.1.east.take().unwrap(); - let mut vec = wb.into_raw_vec(); - vec.resize(wb0.len(), 0.0); - let wb = - Array2::from_shape_vec(wb0.raw_dim(), vec).unwrap(); - workbuffer_edges.1.east = Some(wb); + let mut vec = boundary.into_raw_vec(); + vec.resize(wb.0.len(), 0.0); + let boundary = + Array2::from_shape_vec(wb.0.raw_dim(), vec).unwrap(); + wb_push.push(boundary); } _ => unreachable!(), } - euler::SAT_east(sbp.deref(), k.east_mut(), y, metrics, wb0.view()); + euler::SAT_east(sbp.deref(), k.east_mut(), y, metrics, wb.0.view()); }; - if let Some(boundary) = boundaries.west { + if let Some(mut boundary) = boundaries.west { boundaries_remaining.west = false; - let wb0 = workbuffer_edges.0.west_mut(); - let wb1 = workbuffer_edges.1.west.insert(boundary); + let wb = workbuffer_edges.west_mut(); + let wb_push = workbuffer_free.west_mut(); match boundary_conditions.west() { DistributedBoundaryConditions::Channel => { - std::mem::swap(wb0, wb1); + std::mem::swap(&mut wb.0, &mut boundary); + wb_push.push(boundary); } DistributedBoundaryConditions::Interpolate(int_op) => { - let is_fine2coarse = wb1.shape()[1] > wb0.shape()[2]; - for (to, from) in wb0.outer_iter_mut().zip(wb1.outer_iter()) + let is_fine2coarse = boundary.shape()[1] > wb.0.shape()[2]; + for (to, from) in + wb.0.outer_iter_mut().zip(boundary.outer_iter()) { if is_fine2coarse { int_op.fine2coarse(from, to); @@ -1148,16 +1145,15 @@ impl DistributedSystemPart { } } // Reshape edge buffer to correct size - let wb = workbuffer_edges.1.west.take().unwrap(); - let mut vec = wb.into_raw_vec(); - vec.resize(wb0.len(), 0.0); - let wb = - Array2::from_shape_vec(wb0.raw_dim(), vec).unwrap(); - workbuffer_edges.1.west = Some(wb); + let mut vec = boundary.into_raw_vec(); + vec.resize(wb.0.len(), 0.0); + let boundary = + Array2::from_shape_vec(wb.0.raw_dim(), vec).unwrap(); + wb_push.push(boundary); } _ => unreachable!(), } - euler::SAT_west(sbp.deref(), k.west_mut(), y, metrics, wb0.view()); + euler::SAT_west(sbp.deref(), k.west_mut(), y, metrics, wb.0.view()); }; }); } From 0ec3e16566da914cf185bba1648f656463d8eba4 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Thu, 30 Sep 2021 05:15:31 +0000 Subject: [PATCH 25/25] Align portable-simd with master --- sbp/Cargo.toml | 2 +- sbp/src/operators/algos.rs | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/sbp/Cargo.toml b/sbp/Cargo.toml index e642625..aac3e27 100644 --- a/sbp/Cargo.toml +++ b/sbp/Cargo.toml @@ -12,7 +12,7 @@ serde = { version = "1.0.115", optional = true, default-features = false, featur num-traits = "0.2.14" float = { path = "../utils/float" } constmatrix = { path = "../utils/constmatrix" } -core_simd = { git = "https://github.com/rust-lang/stdsimd" } +core_simd = { git = "https://github.com/rust-lang/portable-simd" } [features] # Use f32 as precision, default is f64 diff --git a/sbp/src/operators/algos.rs b/sbp/src/operators/algos.rs index 74c351f..cb5dee3 100644 --- a/sbp/src/operators/algos.rs +++ b/sbp/src/operators/algos.rs @@ -393,7 +393,6 @@ pub(crate) fn diff_op_2d_sliceable_y_simd