adapt multigrid to dynamic dispatch

This commit is contained in:
Magnus Ulimoen 2020-04-15 19:49:59 +02:00
parent d4831b4336
commit a37d9f4bce
2 changed files with 33 additions and 80 deletions

View File

@ -13,3 +13,4 @@ indicatif = "0.14.0"
structopt = "0.3.13" structopt = "0.3.13"
ndarray = "0.13.0" ndarray = "0.13.0"
json = "0.12.4" json = "0.12.4"
either = "1.5.3"

View File

@ -1,4 +1,6 @@
#![feature(str_strip)] #![feature(str_strip)]
use either::*;
use sbp::operators::{SbpOperator2d, UpwindOperator2d};
use sbp::utils::json_to_grids; use sbp::utils::json_to_grids;
use sbp::*; use sbp::*;
use structopt::StructOpt; use structopt::StructOpt;
@ -12,31 +14,20 @@ struct System {
wb: Vec<euler::WorkBuffers>, wb: Vec<euler::WorkBuffers>,
k: [Vec<euler::Field>; 4], k: [Vec<euler::Field>; 4],
grids: Vec<grid::Grid>, grids: Vec<grid::Grid>,
metrics: Vec<Metrics>, metrics: Vec<grid::Metrics>,
bt: Vec<euler::BoundaryCharacteristics>, bt: Vec<euler::BoundaryCharacteristics>,
eb: Vec<euler::BoundaryStorage>, eb: Vec<euler::BoundaryStorage>,
time: Float, time: Float,
operators: Vec<Either<Box<dyn SbpOperator2d>, Box<dyn UpwindOperator2d>>>,
interpolation_operators: Vec<euler::InterpolationOperators>, interpolation_operators: Vec<euler::InterpolationOperators>,
} }
enum Metrics {
Upwind4(grid::Metrics<operators::Upwind4, operators::Upwind4>),
Upwind9(grid::Metrics<operators::Upwind9, operators::Upwind9>),
Upwind4h2(grid::Metrics<operators::Upwind4h2, operators::Upwind4h2>),
Trad4(grid::Metrics<operators::SBP4, operators::SBP4>),
Trad8(grid::Metrics<operators::SBP8, operators::SBP8>),
Upwind4Upwind4h2(grid::Metrics<operators::Upwind4, operators::Upwind4h2>),
Upwind4h2Upwind4(grid::Metrics<operators::Upwind4h2, operators::Upwind4>),
}
impl System { impl System {
fn new( fn new(
grids: Vec<grid::Grid>, grids: Vec<grid::Grid>,
bt: Vec<euler::BoundaryCharacteristics>, bt: Vec<euler::BoundaryCharacteristics>,
interpolation_operators: Vec<euler::InterpolationOperators>, interpolation_operators: Vec<euler::InterpolationOperators>,
operatorx: &str, operators: Vec<Either<Box<dyn SbpOperator2d>, Box<dyn UpwindOperator2d>>>,
operatory: &str,
) -> Self { ) -> Self {
let fnow = grids let fnow = grids
.iter() .iter()
@ -50,34 +41,10 @@ impl System {
let k = [fnow.clone(), fnow.clone(), fnow.clone(), fnow.clone()]; let k = [fnow.clone(), fnow.clone(), fnow.clone(), fnow.clone()];
let metrics = grids let metrics = grids
.iter() .iter()
.map(|g| match (operatorx, operatory) { .zip(&operators)
("upwind4", "upwind4") => Metrics::Upwind4( .map(|(g, op)| {
g.metrics::<operators::Upwind4, operators::Upwind4>() let sbpop: &dyn SbpOperator2d = op.as_ref().either(|op| &**op, |uo| uo.as_sbp());
.unwrap(), g.metrics(sbpop).unwrap()
),
("upwind9", "upwind9") => Metrics::Upwind9(
g.metrics::<operators::Upwind9, operators::Upwind9>()
.unwrap(),
),
("upwind4h2", "upwind4h2") => Metrics::Upwind4h2(
g.metrics::<operators::Upwind4h2, operators::Upwind4h2>()
.unwrap(),
),
("trad4", "trad4") => {
Metrics::Trad4(g.metrics::<operators::SBP4, operators::SBP4>().unwrap())
}
("trad8", "trad8") => {
Metrics::Trad8(g.metrics::<operators::SBP8, operators::SBP8>().unwrap())
}
("upwind4", "upwind4h2") => Metrics::Upwind4Upwind4h2(
g.metrics::<operators::Upwind4, operators::Upwind4h2>()
.unwrap(),
),
("upwind4h2", "upwind4") => Metrics::Upwind4h2Upwind4(
g.metrics::<operators::Upwind4h2, operators::Upwind4>()
.unwrap(),
),
(opx, opy) => panic!("operator combination {}x{} not known", opx, opy),
}) })
.collect::<Vec<_>>(); .collect::<Vec<_>>();
@ -98,6 +65,7 @@ impl System {
eb, eb,
time: 0.0, time: 0.0,
interpolation_operators, interpolation_operators,
operators,
} }
} }
@ -114,6 +82,7 @@ impl System {
let wb = &mut self.wb; let wb = &mut self.wb;
let mut eb = &mut self.eb; let mut eb = &mut self.eb;
let intops = &self.interpolation_operators; let intops = &self.interpolation_operators;
let operators = &self.operators;
let rhs = move |fut: &mut [euler::Field], let rhs = move |fut: &mut [euler::Field],
prev: &[euler::Field], prev: &[euler::Field],
@ -122,36 +91,22 @@ impl System {
_mt: &mut ()| { _mt: &mut ()| {
let bc = euler::extract_boundaries(prev, &bt, &mut eb, &grids, time, Some(intops)); let bc = euler::extract_boundaries(prev, &bt, &mut eb, &grids, time, Some(intops));
pool.scope(|s| { pool.scope(|s| {
for ((((fut, prev), bc), wb), metrics) in fut for (((((fut, prev), bc), wb), metrics), op) in fut
.iter_mut() .iter_mut()
.zip(prev.iter()) .zip(prev.iter())
.zip(bc) .zip(bc)
.zip(wb.iter_mut()) .zip(wb.iter_mut())
.zip(metrics.iter()) .zip(metrics.iter())
.zip(operators.iter())
{ {
s.spawn(move |_| match metrics { s.spawn(move |_| match op.as_ref() {
Metrics::Upwind4(metrics) => { Left(sbp) => {
euler::RHS_upwind(fut, prev, metrics, &bc, &mut wb.0) euler::RHS_trad(&**sbp, fut, prev, metrics, &bc, &mut wb.0);
} }
Metrics::Upwind9(metrics) => { Right(uo) => {
euler::RHS_upwind(fut, prev, metrics, &bc, &mut wb.0) euler::RHS_upwind(&**uo, fut, prev, metrics, &bc, &mut wb.0);
} }
Metrics::Upwind4h2(metrics) => { })
euler::RHS_upwind(fut, prev, metrics, &bc, &mut wb.0)
}
Metrics::Trad4(metrics) => {
euler::RHS_trad(fut, prev, metrics, &bc, &mut wb.0)
}
Metrics::Trad8(metrics) => {
euler::RHS_trad(fut, prev, metrics, &bc, &mut wb.0)
}
Metrics::Upwind4Upwind4h2(metrics) => {
euler::RHS_trad(fut, prev, metrics, &bc, &mut wb.0)
}
Metrics::Upwind4h2Upwind4(metrics) => {
euler::RHS_trad(fut, prev, metrics, &bc, &mut wb.0)
}
});
} }
}); });
}; };
@ -201,7 +156,6 @@ struct Options {
} }
fn main() { fn main() {
type SBP = operators::Upwind4;
let opt = Options::from_args(); let opt = Options::from_args();
let filecontents = std::fs::read_to_string(&opt.json).unwrap(); let filecontents = std::fs::read_to_string(&opt.json).unwrap();
@ -249,23 +203,20 @@ fn main() {
west: Some(Box::new(operators::Interpolation4)), west: Some(Box::new(operators::Interpolation4)),
}) })
.collect::<Vec<_>>(); .collect::<Vec<_>>();
let grids = jgrids.into_iter().map(|egrid| egrid.grid).collect();
let grids = jgrids
.into_iter()
.map(|egrid| egrid.grid)
.collect::<Vec<grid::Grid>>();
let integration_time: Float = json["integration_time"].as_number().unwrap().into(); let integration_time: Float = json["integration_time"].as_number().unwrap().into();
let (operatorx, operatory) = { let operators = grids
if json["operator"].is_object() { .iter()
( .map(|_| Right(Box::new(operators::Upwind4) as Box<dyn UpwindOperator2d>))
json["operator"]["x"].as_str().unwrap(), .collect::<Vec<_>>();
json["operator"]["y"].as_str().unwrap(),
)
} else {
let op = json["operator"].as_str().unwrap_or("upwind4");
(op, op)
}
};
let mut sys = System::new(grids, bt, interpolation_operators, operatorx, operatory); let mut sys = System::new(grids, bt, interpolation_operators, operators);
sys.vortex(0.0, vortexparams); sys.vortex(0.0, vortexparams);
let max_n = { let max_n = {
@ -331,10 +282,11 @@ fn main() {
if opt.error { if opt.error {
let time = ntime as Float * dt; let time = ntime as Float * dt;
let mut e = 0.0; let mut e = 0.0;
for (fmod, grid) in sys.fnow.iter().zip(&sys.grids) { for ((fmod, grid), op) in sys.fnow.iter().zip(&sys.grids).zip(&sys.operators) {
let mut fvort = fmod.clone(); let mut fvort = fmod.clone();
fvort.vortex(grid.x(), grid.y(), time, vortexparams); fvort.vortex(grid.x(), grid.y(), time, vortexparams);
e += fmod.h2_err::<SBP>(&fvort); 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); println!("Total error: {:e}", e);
} }