SummationByParts/multigrid/src/main.rs

use argh::FromArgs;
use rayon::prelude::*;

use euler::eval::Evaluator;
use sbp::operators::SbpOperator2d;
use sbp::*;

mod file;
mod parsing;
use file::*;
mod eval;

struct System {
    fnow: Vec<euler::Field>,
    fnext: Vec<euler::Field>,
    wb: Vec<euler::WorkBuffers>,
    k: [Vec<euler::Diff>; 4],
    grids: Vec<grid::Grid>,
    metrics: Vec<grid::Metrics>,
    bt: Vec<euler::BoundaryCharacteristics>,
    eb: Vec<euler::BoundaryStorage>,
    time: Float,
    operators: Vec<Box<dyn SbpOperator2d>>,
}

use std::sync::atomic::{AtomicBool, Ordering};
pub(crate) static MULTITHREAD: AtomicBool = AtomicBool::new(false);

impl integrate::Integrable for System {
    type State = Vec<euler::Field>;
    type Diff = Vec<euler::Diff>;
    fn assign(s: &mut Self::State, o: &Self::State) {
        if MULTITHREAD.load(Ordering::Acquire) {
            s.par_iter_mut()
                .zip(o.par_iter())
                .for_each(|(s, o)| euler::Field::assign(s, o))
        } else {
            s.iter_mut()
                .zip(o.iter())
                .for_each(|(s, o)| euler::Field::assign(s, o))
        }
    }
    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<grid::Grid>,
        bt: Vec<euler::BoundaryCharacteristics>,
        operators: Vec<Box<dyn SbpOperator2d>>,
    ) -> Self {
        let fnow = grids
            .iter()
            .map(|g| euler::Field::new(g.ny(), g.nx()))
            .collect::<Vec<_>>();
        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::<Vec<_>>();
        let k = [k.clone(), k.clone(), k.clone(), k];
        let metrics = grids
            .iter()
            .zip(&operators)
            .map(|(g, op)| g.metrics(&**op).unwrap())
            .collect::<Vec<_>>();

        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<euler::Diff>, prev: &Vec<euler::Field>, 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::<integrate::Rk4, System, _>(
            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
    }
}

#[derive(Debug, FromArgs)]
/// Options for configuring and running the solver
struct CliOptions {
    #[argh(positional)]
    json: std::path::PathBuf,
    /// number of simultaneous threads
    #[argh(option, short = 'j')]
    jobs: Option<usize>,
    /// name of output file
    #[argh(
        option,
        short = 'o',
        default = "std::path::PathBuf::from(\"output.hdf\")"
    )]
    output: std::path::PathBuf,
    /// number of outputs to save
    #[argh(option, short = 'n')]
    number_of_outputs: Option<u64>,
    /// print the time to complete, taken in the compute loop
    #[argh(switch)]
    timings: bool,
    /// print error at the end of the run
    #[argh(switch)]
    error: bool,
    /// disable the progressbar
    #[argh(switch)]
    no_progressbar: bool,
    /// output information regarding time elapsed and error
    /// in json format
    #[argh(switch)]
    output_json: bool,
}

#[derive(Default, serde::Serialize)]
struct OutputInformation {
    filename: std::path::PathBuf,
    #[serde(skip_serializing_if = "Option::is_none")]
    time_elapsed: Option<std::time::Duration>,
    #[serde(skip_serializing_if = "Option::is_none")]
    error: Option<Float>,
}

fn main() {
    let opt: CliOptions = argh::from_env();
    let filecontents = std::fs::read_to_string(&opt.json).unwrap();

    let config: parsing::Configuration = match json5::from_str(&filecontents) {
        Ok(config) => config,
        Err(e) => {
            eprintln!("Configuration could not be read: {}", e);
            if let json5::Error::Message {
                location: Some(location),
                ..
            } = e
            {
                eprintln!("\t{:?}", location);
            }
            return;
        }
    };

    let parsing::RuntimeConfiguration {
        names,
        grids,
        grid_connections,
        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 dt = sys.max_dt();

    let ntime = (integration_time / dt).round() as u64;

    {
        let nthreads = opt.jobs.unwrap_or(1);
        if nthreads > 1 {
            MULTITHREAD.store(true, Ordering::Release);
            rayon::ThreadPoolBuilder::new()
                .num_threads(nthreads)
                .build_global()
                .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);
    }
    progressbar.finish_and_clear();

    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 {
        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);
    }

    if opt.output_json {
        println!("{}", json5::to_string(&outinfo).unwrap());
    } else {
        if let Some(duration) = outinfo.time_elapsed {
            println!("Time elapsed: {} seconds", duration.as_secs_f64());
        }
        if let Some(error) = outinfo.error {
            println!("Total error: {:e}", error);
        }
    }
}

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})"),
        )
    }
}