Swift-Hohenberg equation

README.md

@@ -25,6 +25,8 @@ Contents
 - PDE
   - Kuramoto-Sivashinsky equation
     - [notebook](KSE.ipynb)
+  - Swift-Hohenberg equation
+    - [notebook](Swift-Hohenberg.ipynb)
 
 Lyapunov analysis
 -----------------
@@ -43,6 +45,9 @@ Gallery
 ### Kuramoto-Sivashinsky equation
 ![KSE](kse.png)
 
+### Swift-Hohenberg equation
+![SWE](swe.png)
+
 License
 -------
 MIT-License, see [LICENSE](LICENSE) file.

examples/swe.rs

@@ -0,0 +1,40 @@
+extern crate eom;
+extern crate ndarray;
+extern crate ndarray_linalg;
+
+use eom::traits::*;
+use eom::*;
+use ndarray::*;
+use ndarray_linalg::*;
+
+fn main() {
+    let n = 128;
+    let l = 100.0;
+    let dt = 1e-3;
+    let interval = 100;
+    let step = 1000;
+
+    let eom = pde::SWE::new(n, l, 1.0, 6.0);
+    let mut pair = pde::Pair::new(n);
+    let n_coef = eom.model_size();
+    let teo = semi_implicit::DiagRK4::new(eom, dt);
+    let mut teo = adaptor::nstep(teo, interval);
+
+    let x: Array1<c64> = c64::new(0.1, 0.0) * random(n_coef);
+    let ts = adaptor::time_series(x, &mut teo);
+
+    // output CSV header
+    print!("time");
+    for i in 0..n {
+        print!(",x{}", i);
+    }
+    println!("");
+
+    for (t, v) in ts.take(step).enumerate() {
+        let time = dt * t as f64;
+        print!("{:e},", time);
+        let u = pair.to_r(v.as_slice().unwrap());
+        let nums: Vec<_> = u.iter().map(|x| format!("{:e}", x)).collect();
+        println!("{}", nums.join(","));
+    }
+}

src/pde/kse.rs

@@ -57,15 +57,15 @@ impl SemiImplicit for KSE {
     where
         S: DataMut<Elem = Self::Scalar>,
     {
-        for i in 0..self.nf {
-            self.u.c[i] = uf[i];
-            self.ux.c[i] = self.k[i] * uf[i];
-        }
+        azip!(mut u(self.u.coeff_view_mut()), mut ux(self.ux.coeff_view_mut()), k(&self.k), uf(&*uf) in {
+            *u = uf;
+            *ux = k * uf;
+        });
         self.u.c2r();
         self.ux.c2r();
-        for (u, ux) in self.u.r.iter_mut().zip(self.ux.r.iter()) {
-            *u = -*u * *ux;
-        }
+        azip!(mut u(self.u.real_view_mut()), ux(self.ux.real_view()) in {
+            *u = -*u * ux;
+        });
         self.u.r2c();
         uf.as_slice_mut().unwrap().copy_from_slice(&self.u.c);
         uf

src/pde/mod.rs

@@ -1,11 +1,15 @@
 //! Example nonlinear PDEs
 
 pub mod kse;
+pub mod swe;
+
 pub use self::kse::KSE;
+pub use self::swe::SWE;
 
 use fftw::array::*;
 use fftw::plan::*;
 use fftw::types::*;
+use ndarray::*;
 
 /// Pair of one-dimensional Real/Complex aligned arrays
 pub struct Pair {
@@ -36,7 +40,6 @@ impl Pair {
     pub fn c2r(&mut self) {
         self.c2r.c2r(&mut self.c, &mut self.r).unwrap();
     }
-
     pub fn to_r<'a, 'b>(&'a mut self, c: &'b [c64]) -> &'a [f64] {
         self.c.copy_from_slice(c);
         self.c2r();
@@ -48,6 +51,22 @@ impl Pair {
         self.r2c();
         &self.c
     }
+
+    pub fn real_view(&self) -> ArrayView1<f64> {
+        ArrayView::from_shape(self.r.len(), &self.r).unwrap()
+    }
+
+    pub fn coeff_view(&self) -> ArrayView1<c64> {
+        ArrayView::from_shape(self.c.len(), &self.c).unwrap()
+    }
+
+    pub fn real_view_mut(&mut self) -> ArrayViewMut1<f64> {
+        ArrayViewMut::from_shape(self.r.len(), &mut self.r).unwrap()
+    }
+
+    pub fn coeff_view_mut(&mut self) -> ArrayViewMut1<c64> {
+        ArrayViewMut::from_shape(self.c.len(), &mut self.c).unwrap()
+    }
 }
 
 impl Clone for Pair {

src/pde/swe.rs

@@ -0,0 +1,80 @@
+use fftw::types::c64;
+use ndarray::*;
+use std::f64::consts::PI;
+
+use super::Pair;
+use crate::traits::*;
+
+/// One-dimensional Swift-Hohenberg equation with spectral-method
+#[derive(Clone)]
+pub struct SWE {
+    n: usize,
+    nf: usize,
+
+    /// System size (length of entire periodic domain)
+    length: f64,
+    /// Parameter for linear stablity
+    r: f64,
+    /// Length scale of instablity
+    qc: f64,
+
+    k: Array1<c64>,
+    u: Pair,
+}
+
+impl ModelSpec for SWE {
+    type Scalar = c64;
+    type Dim = Ix1;
+    fn model_size(&self) -> usize {
+        self.nf
+    }
+}
+
+impl SWE {
+    pub fn new(n: usize, length: f64, r: f64, lc: f64) -> Self {
+        let nf = n / 2 + 1;
+        let k0 = 2.0 * PI / length;
+        let qc = 2.0 * PI / lc;
+        SWE {
+            n,
+            nf,
+            length,
+            r,
+            qc,
+            k: Array::from_iter((0..nf).map(|i| c64::new(0.0, k0 * i as f64))),
+            u: Pair::new(n),
+        }
+    }
+}
+
+impl SemiImplicit for SWE {
+    fn nlin<'a, S>(
+        &mut self,
+        uf: &'a mut ArrayBase<S, Self::Dim>,
+    ) -> &'a mut ArrayBase<S, Self::Dim>
+    where
+        S: DataMut<Elem = Self::Scalar>,
+    {
+        let c2 = 1.64;
+        let c3 = 1.0;
+
+        azip!(mut u(self.u.coeff_view_mut()), uf(&*uf) in {
+            *u = uf;
+        });
+        self.u.c2r();
+        azip!(mut u(self.u.real_view_mut()) in {
+            *u = c2 * *u * *u - c3 * *u * *u * *u;
+        });
+        self.u.r2c();
+        uf.as_slice_mut().unwrap().copy_from_slice(&self.u.c);
+        uf
+    }
+
+    fn diag(&self) -> Array1<c64> {
+        let r = c64::new(self.r, 0.0);
+        let qc2 = c64::new(self.qc.powi(2), 0.0);
+        let k2 = &self.k * &self.k;
+        let d = qc2 + &k2;
+        r - &d * &d
+    }
+}