Scheme::Core

examples/l63_accuracy.rs

@@ -11,11 +11,11 @@ use ndarray_linalg::*;
 use eom::*;
 use eom::traits::*;
 
-fn check_accuracy<A, D, F, Sc>(teo: Sc, init: Array<A, D>, fname: &str)
+fn check_accuracy<A, D, Sc>(teo: Sc, init: Array<A, D>, fname: &str)
 where
     A: Scalar<Real = f64>,
     D: Dimension,
-    Sc: Scheme<F, Scalar = A, Dim = D, Time = f64>,
+    Sc: Scheme<Scalar = A, Dim = D, Time = f64>,
 {
     let acc = adaptor::accuracy(teo, init, 0.01, 1000, 12);
     let mut f = File::create(fname).unwrap();

src/adaptor.rs

@@ -5,7 +5,7 @@ use ndarray_linalg::*;
 use super::traits::*;
 
 /// Test time accuracy of equation of motion
-pub fn accuracy<A, D, F, Sc>(
+pub fn accuracy<A, D, Sc>(
     mut teo: Sc,
     init: Array<A, D>,
     dt_base: A::Real,
@@ -15,7 +15,7 @@ pub fn accuracy<A, D, F, Sc>(
 where
     A: Scalar,
     D: Dimension,
-    Sc: Scheme<F, Scalar = A, Dim = D, Time = A::Real>,
+    Sc: Scheme<Scalar = A, Dim = D, Time = A::Real>,
 {
     let data: Vec<_> = (0..num_scale)
         .map(|n| {
@@ -38,15 +38,15 @@ where
 /// Iterate equation of motion by `step` at once
 pub fn iterate<S, TEO>(
     teo: &mut TEO,
-    x0: ArrayBase<S, TEO::Dim>,
+    mut x0: ArrayBase<S, TEO::Dim>,
     step: usize,
 ) -> ArrayBase<S, TEO::Dim>
 where
     S: DataMut<Elem = TEO::Scalar> + DataClone,
     TEO: TimeEvolution,
 {
-    let mut ts = time_series(x0, teo);
-    ts.nth(step - 1).unwrap()
+    teo.iterate_n(&mut x0, step);
+    x0
 }
 
 /// An iterator generated by [time_series] for the time-series of the given EoM.

src/explicit.rs

@@ -23,11 +23,18 @@ impl<A: Scalar, F: Explicit<Scalar = A>> TimeStep for Euler<F> {
     }
 }
 
-impl<F: Explicit> Scheme<F> for Euler<F> {
+impl<F: Explicit> Scheme for Euler<F> {
+    type Core = F;
     fn new(f: F, dt: Self::Time) -> Self {
         let x = Array::zeros(f.model_size());
         Self { f, dt, x }
     }
+    fn core(&self) -> &Self::Core {
+        &self.f
+    }
+    fn core_mut(&mut self) -> &mut Self::Core {
+        &mut self.f
+    }
 }
 
 impl<F: Explicit> ModelSpec for Euler<F> {
@@ -72,12 +79,19 @@ impl<A: Scalar, F: Explicit<Scalar = A>> TimeStep for Heun<F> {
     }
 }
 
-impl<F: Explicit> Scheme<F> for Heun<F> {
+impl<F: Explicit> Scheme for Heun<F> {
+    type Core = F;
     fn new(f: F, dt: Self::Time) -> Self {
         let x = Array::zeros(f.model_size());
         let k1 = Array::zeros(f.model_size());
         Self { f, dt, x, k1 }
     }
+    fn core(&self) -> &Self::Core {
+        &self.f
+    }
+    fn core_mut(&mut self) -> &mut Self::Core {
+        &mut self.f
+    }
 }
 
 impl<F: Explicit> ModelSpec for Heun<F> {
@@ -135,7 +149,8 @@ impl<A: Scalar, F: Explicit<Scalar = A>> TimeStep for RK4<F> {
     }
 }
 
-impl<F: Explicit> Scheme<F> for RK4<F> {
+impl<F: Explicit> Scheme for RK4<F> {
+    type Core = F;
     fn new(f: F, dt: Self::Time) -> Self {
         let x = Array::zeros(f.model_size());
         let k1 = Array::zeros(f.model_size());
@@ -150,6 +165,12 @@ impl<F: Explicit> Scheme<F> for RK4<F> {
             k3,
         }
     }
+    fn core(&self) -> &Self::Core {
+        &self.f
+    }
+    fn core_mut(&mut self) -> &mut Self::Core {
+        &mut self.f
+    }
 }
 
 impl<F: Explicit> ModelSpec for RK4<F> {

src/semi_implicit.rs

@@ -52,8 +52,8 @@ impl<F: SemiImplicit> TimeEvolution for Diagonal<F> {
     }
 }
 
-impl<F: SemiImplicit> Scheme<F> for Diagonal<F> {
-    fn new(f: F, dt: Self::Time) -> Self {
+impl<F: SemiImplicit> Diagonal<F> {
+    fn new(f: F, dt: <Self as TimeStep>::Time) -> Self {
         let diag = f.diag();
         let mut exp_diag = diag.to_owned();
         for v in exp_diag.iter_mut() {
@@ -79,7 +79,8 @@ pub struct DiagRK4<F: SemiImplicit> {
     k3: Array<F::Scalar, F::Dim>,
 }
 
-impl<F: SemiImplicit> Scheme<F> for DiagRK4<F> {
+impl<F: SemiImplicit> Scheme for DiagRK4<F> {
+    type Core = F;
     fn new(nlin: F, dt: Self::Time) -> Self {
         let lin = Diagonal::new(nlin.clone(), dt / into_scalar(2.0));
         let x = Array::zeros(lin.model_size());
@@ -98,6 +99,12 @@ impl<F: SemiImplicit> Scheme<F> for DiagRK4<F> {
             k3,
         }
     }
+    fn core(&self) -> &Self::Core {
+        &self.nlin
+    }
+    fn core_mut(&mut self) -> &mut Self::Core {
+        &mut self.nlin
+    }
 }
 
 impl<F: SemiImplicit> TimeStep for DiagRK4<F> {

src/traits.rs

@@ -43,9 +43,27 @@ pub trait TimeEvolution: ModelSpec + TimeStep {
     ) -> &'a mut ArrayBase<S, Self::Dim>
     where
         S: DataMut<Elem = Self::Scalar>;
+
+    /// calculate next step
+    fn iterate_n<'a, S>(
+        &mut self,
+        a: &'a mut ArrayBase<S, Self::Dim>,
+        n: usize,
+    ) -> &'a mut ArrayBase<S, Self::Dim>
+    where
+        S: DataMut<Elem = Self::Scalar>,
+    {
+        for _ in 0..n {
+            self.iterate(a);
+        }
+        a
+    }
 }
 
 /// Time evolution schemes
-pub trait Scheme<F>: TimeEvolution {
-    fn new(f: F, dt: Self::Time) -> Self;
+pub trait Scheme: TimeEvolution {
+    type Core: ModelSpec<Scalar = Self::Scalar, Dim = Self::Dim>;
+    fn new(f: Self::Core, dt: Self::Time) -> Self;
+    fn core(&self) -> &Self::Core;
+    fn core_mut(&mut self) -> &mut Self::Core;
 }