Rust-ndarray自主練⑦ 湧出し渦輪供養

  • 供養① 十字衝突
    4方向から中央にめがけて十字衝突

    中央でぶつかって

    横に逃げる
  • 供養② 19方向衝突
    19方向から中央にめがけて衝突
    中央でぶつかって
    逃げられずにグダグダ
  • 供養③ 湧出し直角クロス
    端が相手の円中心の位置になるようにし湧出し口を直角にクロス 変形して ちえの輪みたいになった
  • 供養④ 渦輪リング
    ちいさい渦輪を連ねてリングにする 渦輪リングをクロス グダグダ
  • 供養⑤ 渦輪横並び
    横に並べて なんやかんやとなって くちびるみたいになって 鎖みたいになった
  • 供養⑥ 渦輪壁衝突
    渦輪で壁を作って そこに渦輪を向かわせる ぶつかって グニャグニャとなって 猫ちゃんに見えなくもない

    出来上がり

    コード

    手入れ中なのでまた今度

Rust-ndarray自主練⑥ 渦輪

まだやっています何周回遅れのRust-ndarray自主練。
前回やったrayon化はわたくしの低級マシンでは何の効果もないことがその後の計測で判明したため一旦なかったこととし、記述の縮約化を目指すこととします。
配られたカードでやるしかないのでございます。 前回は速度場u,v,wを3つの変数に割当てていましたが、ひとつのarrayにまとめた方が記述がすくなくすみます。 前回なぜそうしなかったかというと、2次元と4次元の行列積(テンソル積)をいい感じに書く方法を知らなかったからです。

ndarray-einsum(アインシュタイン和)

現代物理学の御輿に縛り付けられた悲劇のチート師匠アインシュタイン様が考案なされたチート技の1つ、アインシュタイン縮約をndarrayで使えるクレートです。

Rust-ndarrayでeinsumを扱うクレートはいくつかあるのですが、numpyのeinsumと同様の記法で縮約指定ができる当クレートを使うこととしました。

これを使うと例えば3x19のarrayと19x75x25x25とのを"dq,qijk->dijk"と指定するだけで積とq軸の縮約をして3x75x25x25のarrayにしてくれる超便利なクレートです。 Result対応と後続の静的化がうざいですが暴力的にunwrapで押さえ込みます。

今回の題材:渦輪

2種類作ります。
前回に引き続きLBMでBGKでSRTで周期境界です。

①湧出し型

円平面から流体がある量、時間湧き出すと周辺の流体が引き摺られて渦輪ができるパターン。
でんじろう先生がやってるアレですね。

②あたかも最初からそこにあった型

ラム・オーセン分布と言う数式(ビオサバールを元にしたような)を使って渦輪の速度場を作り、あたかも最初から空間内にいたかのように初期分布に染み込ませるパターン。

画像

  • ①-1:湧出し渦輪1つ
    右側の枠は密度変化の表示。0.95〜0.995くらいの範囲をフィルタして表示。

  • ①-2:渦輪正面衝突
    相向って湧出し
    ぶつかって
    大きくなった

  • ①-3:渦輪上下ずらし正面衝突
    上下ずらして衝突
    斜めになって
    六角形になった。土星の極みたい。

  • ①-4:斜め衝突
    片方45度傾けて衝突
    中間の角度になった

  • ①-5:三連星
    同方向、同軸に同間隔で3つ湧出し
    最初の2つがさきに合体
    ドムの顔みたいな形になった
    Xが現れた

  • ②-1:ラム・オーセン渦輪1つ
    スマートで大人しい印象

  • ②-2:ラム・オーセン渦輪正面衝突
    自走しないので近くに置く
    衝突 というより融合

出来上がり

コード(GitHub)

ToroidalVortex

|--simulation
|  |--sdl2-draw-collision_lamb_oseen
|  | ※②あたかも最初からそこにあった型
|  |--sdl2-draw-collision_vortex_wheel_an
|  | ※①湧出し型
|  |--sdl2-draw-kalman_vortex
|  | ※カルマン渦(無影響確認)
|--vortex_ring
| ※格子ボルツマン法(アインシュタイン和バージョン)

内部依存クレート

emb_arrowgraph
emb_bargraph
emb_linetrim
emb_shapegraph
linspacef32
wio_sbcamera_ortho
wio_sbeye

RushJobRust(14) ndarray自主練⑤格子ボルツマン法 レーヨン化

毎度おなじみ周回遅れのrust-ndarray自主練の5回目
前回改悪した格子ボルツマン法をrayonの力で格上げします。
前回はこちら
ボルツマン方程式の理解に至らないまま(というか本の入手が遅れてるので何も進んでいないまま)プログラムの改竄は続いてゆきます。

ndarrayに絡んだrayonの使用方法は以下のお方の記事を参考にしました。

rayonはイテレータの処理を並列化してくれるという楽ちんなクレートです。
いにしえのクレートなので皆さん百も承知でしょうね。

ndarrayとの組み合わせとなる今回のプログラムでは、以下2つのやり方で切り抜けました。
①方向軸の単位で配列を並列化して更新したい場合は軸イテレータミュータブルにする。軸イテレータはインデックスアクセス可能。

    // 流入 ----------------------------------
    // 両端で右向きの安定した流れを強制する
    f.axis_iter_mut(Axis(0))
    .into_par_iter()
    .enumerate()
    .filter(|(i, _)| ci[(*i, 2)] != 0)
    .for_each(|(i, mut fi)|
      fi.slice_mut(s![.., .., 0]).fill(
        w[i] * 
        (1.+3.*v0*c[(i,2)]+4.5*v0p2-1.5*v0p2)
      )
    );

②reduceはstdのとは違い、単位元を用意する。
単位元は整数なら0、配列ならzeros。

    // 速度 ---------------------------------
    // vx(x,t)=1/p(Σᵢ cᵢxfᵢ(u,t))
    let vx =
      f.axis_iter(Axis(0))
      .into_par_iter()
      .enumerate()
      .filter(|(i, _)| ci[(*i, 2)] != 0)
      .map(|(i, fi)| &fi * c[(i, 2)])
      .reduce(
        || Array3::<T>::zeros((nz,ny,nx)),
        |sum, fc| &sum + &fc
      )
      / &p
    ;

もうパッと見ただけでは何をやろうとしていたのか解りづらくなってしまうところが難点ですかね。

バウンスバックに関しては並列化の目処が立っていません。

で、結果どうだったかというと 残念ながらわたくしの非力なマシンではちょっと速くなった程度です。
16年前のiMac 2009late core2duoでcpu使用率140%(撮影ではvncかましているのでもっと落ちる)ではそんなものかもしれません。
せめて4コアあれば全然違うんだろうなあと夢見ております。
別にこのマシンをこだわって使っているわけてはなく。完全に買い替え時を逸してしまっていたわけでございます。

出来上がり

コード(GitHub)

CFDPython

|--CFDPython
|  |--cavityflow2d
|  |--cavityflow3d
|  |--simulation
|  |  |--.DS_Store
|  |  |--draw-cavityflow2d
|  |  |--draw-cavityflow3d
|  |  |--graph_supply_quiver2d
|  |  |--graph_supply_quiver3d
|  |  |--graph_supply_vortex_lbm
|  |  |--graph_supply_vortex_lbm3d
|  |  |--graph_supply_vortex_lbm_arrow
|  |  |--sdl2-draw-cavityflow2d
|  |  |--sdl2-draw-cavityflow3d
|  |  |--sdl2-draw-vortex_lbm
|  |  |--sdl2-draw-vortex_lbm2dmix
|  |  |--sdl2-draw-vortex_lbm2dmix_par
|  |  | ※SDL2での表示(2Dパラレル)
|  |  |--sdl2-draw-vortex_lbm2dmix_sgl
|  |  | ※SDL2での表示(2Dシングル)
|  |  |--sdl2-draw-vortex_lbm3d
|  |  |--sdl2-draw-vortex_lbm3d_par
|  |  | ※SDL2での表示(3Dパラレル)
|  |  |--sdl2-draw-vortex_lbm3d_sgl
|  |  | ※SDL2での表示(3Dシングル)
|  |  |--sdl2-draw-vortex_lbm_arrow
|  |--vortex_lbm
|  |--vortex_lbm2dmix
|  |--vortex_lbm2dmix_par
|  | ※格子ボルツマン法 2Dパラレル
|  |--vortex_lbm2dmix_sgl
|  | ※格子ボルツマン法 2Dシングル
|  |--vortex_lbm3d
|  |--vortex_lbm3d_par
|  | ※格子ボルツマン法 3Dパラレル
|  |--vortex_lbm3d_sgl
|  | ※格子ボルツマン法 3Dシングル
|--ndarray_aulait

内部依存クレート

emb_arrowgraph
emb_bargraph
emb_linetrim
emb_shapegraph
linspacef32
wio_sbcamera_ortho
wio_sbeye

RushJobRust(14) ndarray自主練④格子ボルツマン法 改悪

周回遅れのrust-ndarray自主練4回目
前回丸パクりした格子ボルツマン法のPGを改悪します。
前回はこちら ボルツマン方程式の理解に至らないわたくしは前回丸パクりしたものを再構成してちょっとだけでも理解に近づこうとしますが、実行速度が落ちるので「改悪」となってしまいます。

D2Q9の速度ベクトル、重み、逆方向インデックスを以下のように定義します。

// D2Q9 モデルの速度ベクトル 
fn d2q9() -> Array2<T> {
  d2q9index().mapv(|q| q as T)
}
// D2Q9 インデックス
fn d2q9index() -> Array2<I> {
  arr2(&[
  //  y ,  x
    [ 0 ,  0 ], // 0
    [ 0 ,  1 ], // E
    [ 1 ,  0 ], // N
    [ 0 , -1 ], // W
    [-1 ,  0 ], // S
    [ 1 ,  1 ], // NE
    [ 1 , -1 ], // NW
    [-1 , -1 ], // SW
    [-1 ,  1 ], // SE
  ])
}
// 重み w[i]
fn w() -> Array1<T> {
  Array1::<T>::from_iter(
    [4./9. as T].into_iter()
    .chain(
      [1./9.].into_iter().cycle().take(4)
    ).chain(
      [1./36.].into_iter().cycle().take(4)
    )
  )
}
// 逆方向インデックス opp[i]
fn opp() -> [U; 9] {
  [0, 3,  4,  1,  2,  7,  8,  5,  6]
}

粒子の取りうる速度の方向が9つあり、それぞれy,x方向への速度を持ちます。

平均分布関数を求めるための重みを9つの方向に以下のように割当てます。
・方向0に4/9
・上下左右方向の4つに1/9
・斜め方向の4つに1/36

逆方向インデックスにはバウンスバック時に180度跳ね返す先のインデックスを設定します。
E(1)の障壁に入ってきたものをW(3)へコピーしてE-W=0とするためのものですmaybe(知らんけど)。

まずは初期状態を作ります。

type T = f32;
type U = usize;
type I = isize;
// 格子ボルツマン法
pub fn vortex_lbm2dmix(
  ny       : U, // 格子行
  nx       : U, // 格子列
  viscosity: T, // 動粘性係数
  v0       : T, // 初速と流入速度 
                // 障壁
  barrier  : impl Iterator<Item=(U, U)>
             + Clone
) -> impl Iterator<Item=Array2<T>> {
  // 緩和パラメータ 1/τ
  // vi=1/3(τ-1/2)→τ=3vi+1/2
  // Re=代表流れ*代表長さ/vi
  let omega = 1. / (3.*viscosity + 0.5);
  // D2Q9 モデルの速度ベクトル cᵢ
  let (c, ci) = (d2q9(), d2q9index());
  // 重み、逆方向インデックス
  let (w, opp) = (w(), opp());
  // v0の二乗
  let v0p2 = v0*v0;
  // すべての配列を右肩上がりの流れに初期化
  // 9方向の粒子密度
  let mut f = Array3::<T>::zeros((9,ny,nx));
 {let arr = Array2::<T>::ones((ny,nx));
  (0..9)
  .for_each(|i| {
    let cv = c[(i, 1)] * v0;
    f.slice_mut(s![i, .., ..]).assign(&(
      w[i] * 
      (&arr + 3.0*cv + 4.5*cv*cv - 1.5*v0p2)
    ))
  });
 }// arrをdrop

配列のメモリレイアウトの選択は結構大事で体感でも違いがわかります。
今回はインデックスが1番大きいxを最後の軸に割当てます。

また、タイプエイリアスは本来こういう使い方でないのは重々承知しておりますゆえお気遣いなくお願い申し上げます。

バウンスバックは単純なものを使うので外部から障壁をイテレータで受け取りメインループ前に対象インデックスと反対方向インデックスの対を作っておきます。

  // バウンスバック境界条件 -----------------
  let bounceback = 
    Array1::<((U, U, U), (U, U, U))>
          ::from_iter(
      iproduct!(
        opp.iter().enumerate().skip(1),
        barrier,
      )
      .map(|((i, oi), (y, x))| (
        (*oi,
         (y as I + ci[(*oi, 0)]) as U,
         (x as I + ci[(*oi, 1)]) as U,
        ), 
        (i, y, x)
      ))
      .filter(|((_,oy,ox), (_,y,x))|
        *y  < ny && *x  < nx &&
        *oy < ny && *ox < nx
      )
    ) 
  ;

並進…
この意味が理解しかねています。
説明によると「格子点上の粒子(流体密度分布関数)が、それぞれの速度方向に従って隣接する格子点へ移動するプロセス」とのことですが、何で?ってなります。
これはボルツマン方程式の理解に至っていないゆえでしょう…

  // 無限イテレータ ------------------------
  (0..).scan(0 , move |_, _| {
    // stream ------------------------------
    // すべてのパーティクルを進行方向に沿って
    // 一歩ずつ移動させる(pbc):
    // fᵢ(u+cᵢΔt,t+Δt)-f*ᵢ(u,t)
    iproduct!(1..9, 0..2)
    .filter(|(i, j)| ci[(*i, *j)] != 0)
    .for_each(|(i, j)| {
      rolla2v(
        f.slice(s![i, .., ..]),
        ci[(i, j)],
        j,
      )
      .assign_to(f.slice_mut(s![i, .., ..]));
    });

方向軸(1..9)と次元軸(0..2)をそれぞれ網羅し並進する方向へロールします。前回作成したオレオレroll関数が活きてきます。
速度が0のものはあらかじめフィルタではじいておきます。

並進後、バウンスバックして密度を求めます。

    // バウンスバック -----------------------
    bounceback
    .iter()
    .for_each(|(oidx, idx)| f[*oidx]=f[*idx]); 
    // 密度 ---------------------------------
    // p(u,t) = Σᵢ fᵢ(u,t)
    let p = f.sum_axis(Axis(0));

次に速度を求めます。

    // 速度 ---------------------------------
    // vx(u,t)=1/p[(f₁+f₅+f₈)-(f₃+f₆+f₇)]
    let vx = 
      (0..9)
      .filter(|i| ci[(*i, 1)] != 0)
      .map(|i| 
        &f.slice(s![i, .., ..]) * c[(i, 1)]
      )
      .reduce(|sum, fc| &sum + &fc).unwrap()
      / &p
    ;
    // vy(u,t)=1/p[(f₂+f₅+f₆)-(f₄+f₇+f₈)]
    let vy =
      (0..9)
      .filter(|i| ci[(*i, 0)] != 0)
      .map(|i| 
        &f.slice(s![i, .., ..]) * c[(i, 0)]
      )
      .reduce(|sum, fc| &sum + &fc).unwrap()
      / &p
    ;

てな感じであとは衝突、流入、渦度の出力と続きます。

こうして見てみると、方向軸の演算は独立しているので並列化がしやすいという評判もうなずけます。

あとは形だけならD2Q9をD3Q19などに拡張するのはたやすく行えます。

D3Q19の速度ベクトル、重み、逆方向インデックスは以下になります。

// D3Q19 モデルの速度ベクトル 
fn d3q19() -> Array2<T> {
  d3q19index().mapv(|q| q as T)
}
// D3Q19 インデックス
fn d3q19index() -> Array2<I> {
  arr2(&[
  //  z,   y ,  x 
    [ 0 ,  0 ,  0 ], // 0
    [ 0 ,  0 ,  1 ], // 1  W
    [ 0 ,  0 , -1 ], // 2  E
    [ 1 ,  0 ,  0 ], // 3  B
    [-1 ,  0 ,  0 ], // 4  F
    [ 0 ,  1 ,  0 ], // 5  N
    [ 0 , -1 ,  0 ], // 6  S
    [-1 ,  0 , -1 ], // 7  EB
    [-1 ,  0 ,  1 ], // 8  WB
    [ 1 ,  0 , -1 ], // 9  EF
    [ 1 ,  0 ,  1 ], // 10 WF
    [ 0 ,  1 ,  1 ], // 11 NE
    [ 0 ,  1 , -1 ], // 12 NW
    [ 0 , -1 ,  1 ], // 13 SE
    [ 0 , -1 , -1 ], // 14 SW
    [-1 , -1 ,  0 ], // 15 SB
    [ 1 ,  1 ,  0 ], // 16 NF
    [ 1 , -1 ,  0 ], // 17 SF
    [-1 ,  1 ,  0 ], // 18 NB
  ])
}
// 重み w[i]
fn w() -> Array1<T> {
  Array1::<T>::from_iter(
    [1./3. as T].into_iter()
    .chain(
      [1./18.].into_iter().cycle().take(6)
    ).chain(
      [1./36.].into_iter().cycle().take(12)
    )
  )
}
// 逆方向インデックス opp[i]
fn opp() -> [U; 19] {
  [0, 2,  1,  4,  3,  6,  5,  8,  7, 10, 9,
     12, 11, 14, 13, 16, 15, 18, 17]
}

プログラム的には4次元にするのも難しくなさそうですが今度はお手本となるモデルが見つかりません。

出来上がり

コード(GitHub)

CFDPractise

|--CFDPython
|  |--cavityflow2d
|  |--cavityflow3d
|  |--simulation
|  |  |--draw-cavityflow2d
|  |  |--draw-cavityflow3d
|  |  |--graph_supply_quiver2d
|  |  |--graph_supply_quiver3d
|  |  |--graph_supply_vortex_lbm
|  |  |--graph_supply_vortex_lbm3d
|  |  | ※グラフ供給
|  |  |--graph_supply_vortex_lbm_arrow
|  |  |--sdl2-draw-cavityflow2d
|  |  |--sdl2-draw-cavityflow3d
|  |  |--sdl2-draw-vortex_lbm
|  |  |--sdl2-draw-vortex_lbm2dmix
|  |  | ※SDL2での表示(改悪2D)
|  |  |--sdl2-draw-vortex_lbm3d
|  |  | ※SDL2での表示(改悪3D)
|  |  |--sdl2-draw-vortex_lbm_arrow
|  |--vortex_lbm
|  |--vortex_lbm2dmix
|  | ※格子ボルツマン法 改悪2D
|  |--vortex_lbm3d
|  | ※格子ボルツマン法 改悪3D
|--ndarray_aulait

内部依存クレート

emb_arrowgraph
emb_bargraph
emb_linetrim
emb_shapegraph
linspacef32
wio_sbcamera_ortho
wio_sbeye

RushJobRust(14) ndarray自主練③格子ボルツマン法

周回遅れのndarray自主練3回目
今回のお題は格子ボルツマン法

以下のサイトで説明されているそのまたソース元のnumpyで作られたものを、今回は特に手を入れずほぼパクりでrust-ndarrayで書く練習となります。 格子ボルツマン法(LBM)は、流体の連続体方程式を直接解くのではなく、離散化されたボルツマン方程式を格子上で解く数値解析手法とのことです。

ボルツマン方程式とは「気体など多数の粒子(分子)の集団が、時間とともにどのように振る舞い、変化していくかを記述する方程式」とのことです。

ボルツマン方程式を理解するにはエントロピーまで理解しなきゃいけないと思ってはいるのですが、現時点ではそこまで理解に至っておりません。

理解が至らないならせめて少しでもボルツマンを感じようと思い講義録や言い伝えを読んだりしたのですが、とてもこの人がみずから命を絶つような人とは思えなく「これは◯られたな...」と思わず呟いてしまったものでした。

プログラムの詳しい説明は上記のサイトで行われていますのでここで付け足すことはありません。というか出来ません。

で、実行して表示されたこれは何を見ているのか。
というと渦度というものを見ているそうです。

渦度...
渦そのものではなく回転の強さを見ているとのことです。

今回のオレオレ:rollのオレオレ関数作成

おそらくみんなが隠し持っているであろうnumpyで言うところのroll関数。 みなさん隠蔽が完璧で、いろいろ探したけれども見つかりません。
なので作らざるを得ない状況でございます。
ついでにmaskcopyの関数も作成した次第です。

SDL2での実行

今回はメモリ使用量がでかいのでwio-terminalでの実行は断念です。
それでもメモリさあればembedで動く...(はず)

出来上がり

コード(GitHub)

CFDPractise

|--CFDPython
|  |--cavityflow2d
|  |--cavityflow3d
|  |--simulation
|  |  |--draw-cavityflow2d
|  |  |--draw-cavityflow3d
|  |  |--graph_supply_quiver2d
|  |  |--graph_supply_quiver3d
|  |  |--graph_supply_vortex_lbm
|  |  | ※グラフ供給
|  |  |--graph_supply_vortex_lbm_arrow
|  |  | ※グラフ供給(ベクトル表示)
|  |  |--sdl2-draw-cavityflow2d
|  |  |--sdl2-draw-cavityflow3d
|  |  |--sdl2-draw-vortex_lbm
|  |  | ※SDL2での表示
|  |  |--sdl2-draw-vortex_lbm_arrow
|  |  | ※SDL2での表示(ベクトル表示)
|  |--vortex_lbm
|  | ※格子ボルツマン法
|--ndarray_aulait
| ※オレオレ関数

内部依存クレート

emb_arrowgraph
emb_bargraph
emb_linetrim
emb_shapegraph
linspacef32
wio_sbcamera_ortho
wio_sbeye

RushJobRust(14) ndarray自主練②3次元キャビティフロー

周回遅れのrust-ndarray自主練2回目になります。
今回は前回の2次元キャビティフローを3次元化します。
前回はこちら 毎度そうなのですが先達のnumpyで作られたものをパクって練習の題材とさせていただいています。

今回のお題

上記本の中の 7.1:3次元立方体キャビティ内の流れ

※ちなみに本ブログはアフィなどは一切行っておりませんので心置きなくリンクを踏んでいただいて構わないと思われるのですが、かと言ってその先でのトラブルには一切責任を取りませんことを謹んでお詫び申し上げます。
※なおこのブログ内に出てくる広告は、はてなブログのショバ代でして私の手元には何も入らないことを表明いたしたく存じます。

あと、おそらくは上記本の内容に沿った以下のサイト
流れ解析の練習-1 キャビティ流れ 3D – おっさんの備考録

の2つにインスパイアしてオマージュとなる次第でございます。

今回の自己肯定

参考元のPGは配列をインデックスアクセスする方式となっているので逆手をとってスライスで攻めてみたところ冗長でダサいコードが出来上がりました。
これも次のステップのためであると強引にポジティブになっていますが、本来の「自己肯定」とは異なり、無理にポジティブな自分になろうとすることは、真の自己肯定から遠ざかる可能性があります。
真の自己肯定には、自己肯定感が低い自分を受け入れることや、「肯定しなければいけない」と思わないことが大切です。

wio-terminalでの実行

毎度おなじみembedded_allocを使用し、11x11x11のグリッドを使用する見積もりで80KBヒープに割当てます。
関数化やスコープを付けるなりして適宜メモリdropするようにすれば60KBまで減らせます。

また、グリッド数の多い場合をembedded_graphics_simulatorを使用してPC上のSDL2で表示します。
これはつまりメモリさえあればwio-terminal上でも動くということを意味します。

出来上がり

コード(GitHub)

CFDPython

|--cavityflow2d
|--cavityflow3d
| ※3次元キャビティフロー
|--simulation
|  |--draw-cavityflow2d
|  |--draw-cavityflow3d
|  | ※wio-termimalでの表示
|  |--graph_supply_quiver2d
|  |--graph_supply_quiver3d
|  | ※グラフ供給
|  |--sdl2-draw-cavityflow2d
|  |--sdl2-draw-cavityflow3d
|  | ※SDL2での表示

内部依存クレート

emb_arrowgraph
emb_bargraph
emb_linetrim
linspacef32
wio_sbcamera_ortho
wio_sbeye

RushJobRust(14) ndarray自主練①2次元キャビティフロー

唐突ではありますが遅ればせながらndarrayを始めたいと思い立ち、自主練を開始させていただきたく存じます。
情報源は以下のお方の記事を辿って行けば大体のものはカバー出来ると思います。
nalgebraを始めた時もこのお方の記事を起点としていました。

自主練の目標

さて、rust-ndarrayのいい所は何といっても人様がnumpyで書いたものをサクッとパクれるところにあると思います。
とはいえ、不足している機能があったりndarray自身の書き方のクセがあったりとか、やり方が公式のndarray_for_numpy_usersに載っていなかったりとかあります。
そういうポイントを一つずつシノギながら高効率でnumpyをパクれることを目標として自主練を行っていこうと思った次第でございます。

ndarrayの実行速度

結論はnumpyに「劣る」と思っていいような結果がいろいろなサイトやredditなどで語られています。
やる意味あるのかと思われるでしょうがやりたくなったものは致し方ない状況とご理解いただければ幸いです。

今回のお題

最初は以下のお方が解説してくださっている2Dキャビティフローを題材としたく勝手にリンクさせていただいています。 ナビエストークス方程式を丁寧に解説されている素晴らしいサイトです。
しかしこれがゴールではなくスタート地点に立つようなhello world的な扱いであることが流体解析の恐ろしいところであります。

今回のシノギ脳死で配列演算に&付け

スライスを定義するs![]マクロの公式解説ページに以下の使用例があります。

fn laplacian(v: &ArrayView2<f32>) -> Array2<f32> {
    -4. * &v.slice(s![1..-1, 1..-1])
    + v.slice(s![ ..-2, 1..-1])
    + v.slice(s![1..-1,  ..-2])
    + v.slice(s![1..-1, 2..  ])
    + v.slice(s![2..  , 1..-1])
}

そうそうこういうのがやりたいんですよ。

なんですけど、ぱっと見違和感を感じます。
なんで最初のスライスだけ&がついてるのか?

これはスライスの戻り値であるArrayView が演算のためのトレイトを &ArrayView に実装しているからとのことで起点となるスライスには&を付けないといけないようです。
しかし、実際流体解析では上記の様なシンプルな形ではなくカッコでくるまれまくった演算になり、そのひとつずつどれに&を付けるか考えないといけないようです。

理解が間違ってるかもしれませんが、それでもそんなことに頭を使いたくない...

ということで、スライスに限らず配列演算には脳死で全て&を付けまして、不要な&はClippyが教えてくれるので後で消していくという方法を使用し効率化をはかります。
(ちなみに以前にも書きましたがわたくしスマホテキストエディタで作ってるものですから上記のようなやり方をしています)

wio-terminalでの実行

ndarrayはno_std対応ですがヒープアロケータ必須になります。
毎度おなじみembedded_allocを使用し、今回は70KBヒープに割当てます。

2次元キャビティフロー

出来上がり

コード(GitHub)

CFDPractise

|--CFDPython
|  |--cavityflow2d
|  |  ※2Dキャビティフロー
|  |--simulation
|  |  |--draw-cavityflow2d
|  |  |  ※wio-terminalでの実行
|  |  |--graph_supply_quiver2d
|  |  |  ※グラフ供給
|  |  |--sdl2-draw-cavityflow2d
|  |  |  ※sdl2での表示

内部依存クレート

emb_arrowgraph
※矢印グラフ
emb_bargraph
emb_linetrim
linspacef32