DEV Community: Toshiki Teramura

stencil; abstract stencil calculation

Toshiki Teramura — Sat, 17 Mar 2018 09:57:12 +0000

I am developing a library for stencil calculation in Rust.

https://github.com/termoshtt/stencil
Abstract stencil calculation

stencil v0.2 is published on crates.io, and you can use it through cargo:

[dependencies]
stencil = "0.2"

The full document is available on docs.rs.

StencilArray trait

This library abstract a stencil calculation through StencilArray trait:

/// Array with stencil calculations
pub trait StencilArray<St>: NdArray
where
    St: Stencil<Elem = Self::Elem, Dim = Self::Dim>,
{
    /// Execute a stencil calculation
    fn stencil_map<Output, Func>(&self, out: &mut Output, Func)
    where
        Output: NdArray<Dim = Self::Dim>,
        Func: Fn(St) -> Output::Elem;
}

stencil_map function is the core of this crate. Basic example is here:

let mut a = torus::Torus::<f64, Ix1>::zeros(n);
let mut b = a.clone();
a.coordinate_fill(|x| x.sin());
let dx = a.dx();
a.stencil_map(&mut b, |n: N1D1<f64>| (n.r - n.l) / (2.0 * dx));

This crate defines Torus<A, D>, which represents N-dimensional torus. a denotes one-dimensional torus, we can regard it as a region with periodic boundary condition. The second argument of stencil_map, a lambda function |n: N1D1<f64>| (n.r - n.l) / (2.0 * dx) denotes a central difference. N1D1 means 1-neighbor on 1-dimensional space. So, this example describes a central difference of sin(x) function, i.e. cos(x) is set b.

N1D1<A> is a small struct:

/// one-neighbor, one-dimensional stencil
#[derive(Clone, Copy)]
pub struct N1D1<A: LinalgScalar> {
    /// left
    pub l: A,
    /// right
    pub r: A,
    /// center
    pub c: A,
}

and implements the trait Stencil. It has values around "what we focused on"; when we update b[i] using a, n has a[i] as n.c and its neighbors a[i-1] as n.l and a[i+1] as n.r. This crate abstracts the stencil calculation in this way.
The implementation of StencilArray for the Torus<A, D> describes how to get N1D1 from the array a. Currently, there are three Stencils: N1D1<A>, N2D1<A>, and N1D2<A>, and two StencilArrays: Torus<A, D> and Line<A, P: Padding, E: Edge>. The Line has a padding at the each edge of line, and it can be regarded as a fixed boundary condition.

Manifold trait

stencil crate introduces one more useful trait, Manifold:

/// Uniformly coordinated array
pub trait Manifold: NdArray {
    /// Type of coordinate
    type Coordinate;

    /// Increment of coordinate
    fn dx(&self) -> Self::Coordinate;

    /// Fill manifold by a function
    fn coordinate_fill<F>(&mut self, F)
    where
        F: Fn(Self::Coordinate) -> Self::Elem;

    /// Map values on manifold using a function
    fn coordinate_map<F>(&mut self, F)
    where
        F: Fn(Self::Coordinate, Self::Elem) -> Self::Elem;
}

This trait gives array a coordinate. This is useful to set value of array by mathematical function as you see above example

a.coordinate_fill(|x| x.sin());

x is set by the implementation of Manifold for each objects.

Example1: Diffusion equation

Let us solve the diffusion equation for both periodic and fixed boundary conditions:

fn periodic(r: Sender<Vec<f64>>) {
    let dt = 1e-3;
    let mut t = torus::Torus::<f64, Ix1>::zeros(32);
    t.coordinate_fill(|x| x.sin());
    let dx = t.dx();
    let mut s = t.clone();
    for _step in 0..3000 {
        t.stencil_map(&mut s, |n: N1D1<f64>| {
            let d2 = (n.l + n.r - 2.0 * n.c) / (dx.powi(2));
            n.c + dt * d2
        });
        r.send(s.as_view().to_vec()).unwrap();
        swap(&mut t, &mut s);
    }
}

fn fixed(r: Sender<Vec<f64>>) {
    let dt = 1e-3;
    let mut t = Line::<f64, P1, Closed>::new(64, 0.0, 1.0, 2.5 * PI);
    t.coordinate_fill(|x| x.sin());
    let dx = t.dx();
    let mut s = t.clone();
    for _step in 0..10000 {
        t.stencil_map(&mut s, |n: N1D1<f64>| {
            let d2 = (n.l + n.r - 2.0 * n.c) / (dx.powi(2));
            n.c + dt * d2
        });
        r.send(s.as_view().to_vec()).unwrap();
        swap(&mut t, &mut s);
    }
}

This example uses asink for storing data. Please see the last post.

Example2: KdV equation

At last, we solve KdV equation using the ZK-scheme:

fn kdv(r: Sender<Vec<f64>>) {
    let dt = 1e-5;
    let mut t = torus::Torus::<f64, Ix1>::zeros(128);
    t.coordinate_fill(|x| x.sin());
    let dx = t.dx();
    let mut b = t.clone();
    let mut n = t.clone();
    for step in 0..1_000_000 {
        t.stencil_map(&mut n, |n: N2D1<f64>| {
            // ZK scheme
            let uux = (n.l + n.c + n.r) * (n.r - n.l) / (3.0 * dx);
            let u3 = (n.rr - 2.0 * n.r + 2.0 * n.l - n.ll) / dx.powi(3);
            uux + 0.01 * u3
        });
        azip!(mut n, b in { *n = b - dt * *n });
        if step % 1000 == 0 {
            r.send(n.as_view().to_vec()).unwrap();
        }
        swap(&mut t, &mut b);
        swap(&mut t, &mut n);
    }
}

Leap-frog algorithm can also be written easily.

Please see the examples/ directory for full examples.
I'd like to get your feedback. Thanks.

asink; async sink for time-series data

Toshiki Teramura — Thu, 15 Mar 2018 18:14:56 +0000

I'd like to share my small library for storing time-series data yielded from numerical simulations.

https://github.com/termoshtt/asink

asink is published to crates.io:

[dependencies]
asink = "0.1.1"

and full document is placed on docs.rs

Sink and mpsc::channel

This library defines Sink, like most logging library, which stores data into the storage.

let sink = msgpack::MsgpackSink::from_str("test.msg");
let (s, th) = sink.run();
// execute simulation...
th.join().unwrap();

s generated by sink.run() is a Sender<T> belongs to std::sync::mpsc. As its name, "Multi-Producer, Single-Consumer", mpsc::channel create a pair of copyable Sender<T> and uncopyable Receiver<T>. Sink starts a thread for writing received data from the main thread through this channel onto the storage.
The second value th is a thread handler to join this thread. When all copied Sender are droped, the channel is closed, and a finalize the Sink thread.

In this case, the sink is a msgpack file. The data received from main thread is encoded into msgpack and written down to this file. Thanks to the serde library, we can serialize any structure implementing Serialize into JSON/msgpack/BSON. Or, you can automatically implement Serialize trait for your struct using serde-derive

extern crate serde;
#[macro_use]
extern crate serde_derive;

#[derive(Serialize)]
struct Doc {
    time: f64,
    data: Vec<f64>,
}

Working Example

extern crate asink;
extern crate serde;
#[macro_use]
extern crate serde_derive;

use asink::*;
use std::sync::mpsc::Sender;

/// This will be serialized into msgpack
#[derive(Serialize)]
struct Doc {
    id: usize,
    data: Vec<f64>,
}

fn experiment(s: Sender<Doc>) {
    for i in 0..5 {
        let doc = Doc {
            id: i,
            data: vec![i as f64],
        };
        s.send(doc).unwrap(); // Send data to sink
    }
}

fn main() {
    let sink = msgpack::MsgpackSink::from_str("test.msg");
    let (s, th) = sink.run(); // Sink start thread to write recieved data into msgpack
    experiment(s);
    th.join().unwrap();
}

We can switch sink by replacing one line:

let sink = json::JsonSink::from_str("test.json"); // output to JSON file
let sink = mongo::MongoSink:::local("mydb", "collection1"); // output to MongoDB on localhost

This library will be useful both storing large binary data (e.g. velocity/megnetic field) and logging datas (e.g. a residual of iterative method like BiCGStab).

I'd like to get your feedback. Thanks.