#![feature(str_strip)] use either::*; use structopt::StructOpt; use sbp::operators::{SbpOperator2d, UpwindOperator2d}; use sbp::*; mod parsing; use parsing::{json_to_grids, json_to_vortex}; mod file; use file::*; pub(crate) type DiffOp = Either, Box>; struct System { fnow: Vec, fnext: Vec, wb: Vec, k: [Vec; 4], grids: Vec, metrics: Vec, bt: Vec, eb: Vec, time: Float, operators: Vec, } 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 = [fnow.clone(), fnow.clone(), fnow.clone(), fnow.clone()]; let metrics = grids .iter() .zip(&operators) .map(|(g, op)| { let sbpop: &dyn SbpOperator2d = op.as_ref().either(|op| &**op, |uo| uo.as_sbp()); g.metrics(sbpop).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, pool: &rayon::ThreadPool) { 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 [euler::Field], prev: &[euler::Field], time: Float| { let prev_all = &prev; pool.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); match op.as_ref() { Left(sbp) => { euler::RHS_trad(&**sbp, fut, prev, metrics, &bc, &mut wb.0); } Right(uo) => { euler::RHS_upwind(&**uo, fut, prev, metrics, &bc, &mut wb.0); } } }) } }); }; let mut k = self .k .iter_mut() .map(|k| k.as_mut_slice()) .collect::>(); sbp::integrate::integrate_multigrid::( rhs, &self.fnow, &mut self.fnext, &mut self.time, dt, &mut k, pool, ); std::mem::swap(&mut self.fnow, &mut self.fnext); } } #[derive(Debug, StructOpt)] struct Options { json: std::path::PathBuf, /// Disable the progressbar #[structopt(long)] no_progressbar: bool, /// Number of simultaneous threads #[structopt(short, long)] jobs: Option>, /// Name of output file #[structopt(default_value = "output.hdf", long, short)] output: std::path::PathBuf, /// Number of outputs to save #[structopt(long, short)] number_of_outputs: Option, /// Print the time to complete, taken in the compute loop #[structopt(long)] timings: bool, /// Print error at the end of the run #[structopt(long)] error: bool, } fn main() { let opt = Options::from_args(); let filecontents = std::fs::read_to_string(&opt.json).unwrap(); let json = json::parse(&filecontents).unwrap(); let vortexparams = json_to_vortex(json["vortex"].clone()); let (names, grids, bt, operators) = json_to_grids(json["grids"].clone(), vortexparams); let integration_time: Float = json["integration_time"].as_number().unwrap().into(); let mut sys = System::new(grids, bt, operators); sys.vortex(0.0, vortexparams); let max_n = { let max_nx = sys.grids.iter().map(|g| g.nx()).max().unwrap(); let max_ny = sys.grids.iter().map(|g| g.ny()).max().unwrap(); std::cmp::max(max_nx, max_ny) }; let dt = 0.2 / (max_n as Float); let ntime = (integration_time / dt).round() as u64; let pool = { let builder = rayon::ThreadPoolBuilder::new(); if let Some(j) = opt.jobs { if let Some(j) = j { builder.num_threads(j) } else { builder } } else { builder.num_threads(1) } .build() .unwrap() }; 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 output = File::create(&opt.output, sys.grids.as_slice(), names).unwrap(); let mut output = OutputThread::new(output); let progressbar = progressbar(opt.no_progressbar, ntime); let timer = if opt.timings { Some(std::time::Instant::now()) } else { None }; for itime in 0..ntime { if should_output(itime) { output.add_timestep(itime, &sys.fnow); } progressbar.inc(1); sys.advance(dt, &pool); } progressbar.finish_and_clear(); if let Some(timer) = timer { let duration = timer.elapsed(); println!("Time elapsed: {} seconds", duration.as_secs_f64()); } 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(); fvort.vortex(grid.x(), grid.y(), time, vortexparams); let sbpop: &dyn SbpOperator2d = op.as_ref().either(|op| &**op, |uo| uo.as_sbp()); e += fmod.h2_err(&fvort, sbpop); } println!("Total error: {:e}", e); } } 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})"), ) } }