adaptive integration option

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<SBP: SbpOperator2d>(universe: &mut System<SBP>, n: usize) {
     for _ in 0..n {
@@ -14,6 +15,19 @@ fn advance_system_upwind<UO: UpwindOperator2d>(universe: &mut System<UO>, n: usi
     }
 }
 
+fn advance_embedded<UO: UpwindOperator2d>(universe: &mut System<UO>, 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();
 }
 

sbp/src/euler.rs

@@ -99,6 +99,13 @@ impl<SBP: SbpOperator2d> System<SBP> {
     pub fn y(&self) -> ArrayView2<Float> {
         self.grid.0.y.view()
     }
+
+    pub fn nx(&self) -> usize {
+        self.grid.0.nx()
+    }
+    pub fn ny(&self) -> usize {
+        self.grid.0.ny()
+    }
 }
 
 impl<UO: UpwindOperator2d> System<UO> {
@@ -127,6 +134,43 @@ impl<UO: UpwindOperator2d> System<UO> {
         );
         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::<integrate::BogackiShampine, _, _>(
+                &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)]

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<BTableau: ButcherTableau, F, RHS>(
     mut rhs: RHS,
@@ -140,6 +189,30 @@ pub fn integrate<BTableau: ButcherTableau, F, RHS>(
     }
 }
 
+#[allow(clippy::too_many_arguments)]
+pub fn integrate_embedded_rk<BTableau: EmbeddedButcherTableau, F, RHS>(
+    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<ArrayView3<'r, Float>>,
+    for<'r> &'r mut F: std::convert::Into<ArrayViewMut3<'r, Float>>,
+    RHS: FnMut(&mut F, &F, Float),
+{
+    integrate::<BTableau, F, RHS>(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<BTableau: ButcherTableau, F, RHS>(