From 5003c0d55c1e8f854aa23fc3d6f3a0103970d83f Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Sun, 19 Apr 2020 20:45:06 +0200 Subject: [PATCH] adaptive integration option --- sbp/benches/euler.rs | 30 ++++++++++++++++++ sbp/src/euler.rs | 44 ++++++++++++++++++++++++++ sbp/src/integrate.rs | 73 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 147 insertions(+) diff --git a/sbp/benches/euler.rs b/sbp/benches/euler.rs index d912573..69611c0 100644 --- a/sbp/benches/euler.rs +++ b/sbp/benches/euler.rs @@ -1,6 +1,7 @@ use criterion::{black_box, criterion_group, criterion_main, Criterion}; use sbp::euler::System; use sbp::operators::{SbpOperator2d, Upwind4, UpwindOperator2d, SBP4}; +use sbp::Float; fn advance_system(universe: &mut System, n: usize) { for _ in 0..n { @@ -14,6 +15,19 @@ fn advance_system_upwind(universe: &mut System, n: usi } } +fn advance_embedded(universe: &mut System, embedded: bool) { + let dt = 0.2 / std::cmp::max(universe.nx(), universe.ny()) as Float; + let t = 1.0; + if embedded { + let mut dt = dt; + universe.advance_adaptive(t, &mut dt, 1e-2); + } else { + for _ in 0..(t / dt).round() as isize { + universe.advance_upwind(dt); + } + } +} + fn performance_benchmark(c: &mut Criterion) { let mut group = c.benchmark_group("EulerSystem"); group.sample_size(25); @@ -48,7 +62,23 @@ fn performance_benchmark(c: &mut Criterion) { advance_system(&mut universe, black_box(20)) }) }); + group.finish(); + let mut group = c.benchmark_group("adaptive integration"); + group.sample_size(10); + let mut universe = System::new(x.into_owned(), y.into_owned(), Upwind4); + group.bench_function("static dt", |b| { + b.iter(|| { + universe.init_with_vortex(0.0, 0.0); + advance_embedded(&mut universe, false); + }) + }); + group.bench_function("adaptive dt", |b| { + b.iter(|| { + universe.init_with_vortex(0.0, 0.0); + advance_embedded(&mut universe, true); + }) + }); group.finish(); } diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs index ac55c35..57c3a23 100644 --- a/sbp/src/euler.rs +++ b/sbp/src/euler.rs @@ -99,6 +99,13 @@ impl System { pub fn y(&self) -> ArrayView2 { self.grid.0.y.view() } + + pub fn nx(&self) -> usize { + self.grid.0.nx() + } + pub fn ny(&self) -> usize { + self.grid.0.ny() + } } impl System { @@ -127,6 +134,43 @@ impl System { ); std::mem::swap(&mut self.sys.0, &mut self.sys.1); } + pub fn advance_adaptive(&mut self, dt: Float, guess_dt: &mut Float, maxerr: Float) { + let bc = BoundaryCharacteristics { + north: BoundaryCharacteristic::This, + south: BoundaryCharacteristic::This, + east: BoundaryCharacteristic::This, + west: BoundaryCharacteristic::This, + }; + let op = &self.op; + let grid = &self.grid; + let wb = &mut self.wb.0; + let mut rhs_upwind = |k: &mut Field, y: &Field, _time: Float| { + let (grid, metrics) = grid; + let boundaries = boundary_extractor(y, grid, &bc); + RHS_upwind(op, k, y, metrics, &boundaries, wb) + }; + let mut time = 0.0; + let mut sys2 = self.sys.0.clone(); + while time < dt { + integrate::integrate_embedded_rk::( + &mut rhs_upwind, + &self.sys.0, + &mut self.sys.1, + &mut sys2, + &mut time, + *guess_dt, + &mut self.k, + ); + let err = self.sys.0.h2_err(&sys2, &self.op); + if err < maxerr { + time += *guess_dt; + std::mem::swap(&mut self.sys.0, &mut self.sys.1); + *guess_dt *= 1.05; + } else { + *guess_dt *= 0.8; + } + } + } } #[derive(Clone, Debug)] diff --git a/sbp/src/integrate.rs b/sbp/src/integrate.rs index 9886e5e..2256340 100644 --- a/sbp/src/integrate.rs +++ b/sbp/src/integrate.rs @@ -8,6 +8,10 @@ pub trait ButcherTableau { const C: &'static [Float]; } +pub trait EmbeddedButcherTableau: ButcherTableau { + const BSTAR: &'static [Float]; +} + pub struct Rk4; impl ButcherTableau for Rk4 { const A: &'static [&'static [Float]] = &[&[0.5], &[0.0, 0.5], &[0.0, 0.0, 1.0]]; @@ -87,6 +91,51 @@ impl ButcherTableau for Rk6 { ]; } +pub struct Fehlberg; + +impl ButcherTableau for Fehlberg { + #[rustfmt::skip] + const A: &'static [&'static [Float]] = &[ + &[1.0 / 4.0], + &[3.0 / 32.0, 9.0 / 32.0], + &[1932.0 / 2197.0, -7200.0 / 2197.0, 7296.0 / 2197.0], + &[439.0 / 216.0, -8.0, 3680.0 / 513.0, -845.0 / 4104.0], + &[-8.0 / 27.0, 2.0, -3544.0 / 2565.0, 1859.0 / 4104.0, -11.0 / 40.0], + ]; + #[rustfmt::skip] + const B: &'static [Float] = &[ + 16.0 / 135.0, 0.0, 6656.0 / 12825.0, 28561.0 / 56430.0, -9.0 / 50.0, 2.0 / 55.0, + ]; + const C: &'static [Float] = &[0.0, 1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0]; +} + +impl EmbeddedButcherTableau for Fehlberg { + const BSTAR: &'static [Float] = &[ + 25.0 / 216.0, + 0.0, + 1408.0 / 2565.0, + 2197.0 / 4104.0, + -1.0 / 5.0, + 0.0, + ]; +} + +pub struct BogackiShampine; + +impl ButcherTableau for BogackiShampine { + const A: &'static [&'static [Float]] = &[ + &[1.0 / 2.0], + &[0.0, 3.0 / 4.0], + &[2.0 / 9.0, 1.0 / 3.0, 4.0 / 9.0], + ]; + const B: &'static [Float] = &[2.0 / 9.0, 1.0 / 3.0, 4.0 / 9.0, 0.0]; + const C: &'static [Float] = &[0.0, 1.0 / 2.0, 3.0 / 4.0, 1.0]; +} + +impl EmbeddedButcherTableau for BogackiShampine { + const BSTAR: &'static [Float] = &[7.0 / 24.0, 1.0 / 4.0, 1.0 / 3.0, 1.0 / 8.0]; +} + #[allow(clippy::too_many_arguments)] pub fn integrate( mut rhs: RHS, @@ -140,6 +189,30 @@ pub fn integrate( } } +#[allow(clippy::too_many_arguments)] +pub fn integrate_embedded_rk( + rhs: RHS, + prev: &F, + fut: &mut F, + fut2: &mut F, + time: &mut Float, + dt: Float, + k: &mut [F], +) where + for<'r> &'r F: std::convert::Into>, + for<'r> &'r mut F: std::convert::Into>, + RHS: FnMut(&mut F, &F, Float), +{ + integrate::(rhs, prev, fut, time, dt, k); + fut2.into().assign(&prev.into()); + for (&b, k) in BTableau::BSTAR.iter().zip(k.iter()) { + if b == 0.0 { + continue; + } + fut2.into().scaled_add(b * dt, &k.into()); + } +} + #[cfg(feature = "rayon")] #[allow(clippy::too_many_arguments)] pub fn integrate_multigrid(