OpenFOAMのcase fileをgitで管理する
OpenFOAMのcase fileをgitで管理する
お久しぶりです。くれなゐです
gitをソースコードの管理以外で使ったことがない方っていませんか?
gitは非常に便利なツールで、様々な履歴を残すことができます。
今日はOpenFOAMのケースファイルを例に、ソースコードだけじゃないんだぞと言うところをお見せしたいと思います。
このブログの対象
git
の使い方をまだ良くわかっていない人- ソースコードの管理以外に
git
の使用価値を見いだせない人 - OpenFOAMの計算条件の履歴を残したい人
バージョンを管理する対象
今回も、いつものcavity
を使用しています。
OpenFOAMのバージョンはちょっと古いですが2.4.0を使用しています。
gitのバージョン管理をする本質ではバージョンはあまり関係ないと考えています(blockMesh等のパスは少し変わりますが…)
OpenFOAM/OpenFOAM-2.4.0/tutorials/incompressible/icoFoam/cavity
あと忘れないうちに
blockMesh
しておきましょう
まずは init, first commit
init
gitでバージョン管理していくために、まずはバージョン管理の環境を用意してあげる必要があります。 それをするコマンドが
git init
です。
こんなメッセージが出てればOKです
Initialized empty Git repository in /home/kurenaif/Desktop/cavity/.git/
add
とりあえず今の状態をcommit
しておきましょう
commitとは簡単に言うとセーブポイントを設置するような感じです。こうしておくことあとでこの状態に復元できますね。
commitをする前に、どのファイルをコミットするか選択しなければなりません。それがadd
になります。すべてのファイルをadd
する場合は、
git add --all
などとしましょう。
status
add
をしたら、 statusを確認しましょう
git status
add
されたファイルが確認できます。
On branch master Initial commit Changes to be committed: (use "git rm --cached <file>..." to unstage) new file: 0/U new file: 0/p new file: constant/polyMesh/blockMeshDict new file: constant/polyMesh/boundary new file: constant/polyMesh/faces new file: constant/polyMesh/neighbour new file: constant/polyMesh/owner new file: constant/polyMesh/points new file: constant/transportProperties new file: system/controlDict new file: system/fvSchemes new file: system/fvSoluti
commit
commitメッセージは必ず付けなければなりません。 commitメッセージは、その変更に対する理由付けだとかそういうものです。あとでこのコミットメッセージを参照して、必要なときに戻ってきたりするわけです。 たとえば、今回の場合、「計算が回ってくれた」「良い計算結果が出た」、「前つまっていたTimeが突破できた」 などなどといったケースで使うことができます。
最初のcommitはfirst commit
とすることが通例です。(これも流派があるので、一概にfirst commitではないです。)
git commit -m "first commit"
計算結果をバージョン履歴に残さないようにする
これで最初の準備が整いました。それではicoFoamを実行しましょう。
icoFoam
blockMesh忘れてなければ実行できると思います。忘れてたらblockMeshしてからicoFoamしましょう。
計算結果は出ましたが、これは私はバージョン管理対象に含めるべきではないと考えています。(というのは、ファイルが膨大になることがおおいため、.gitディレクトリが大変になることが予想される)
では、計算結果がgit add --all
の対象に含まれないようにしてあげましょう
そのために、.gitignore
を設定します。
現在、
git status
を見てみると、
On branch master Untracked files: (use "git add <file>..." to include in what will be committed) 0.1/ 0.2/ 0.3/ 0.4/ 0.5/ nothing added to commit but untracked files present (use "git add" to track)
このようになっています。このメッセージは、add
されてないよーという趣旨ですので、これからの目標はstatus
をした時にこいつらを消してやることにあります。
消すためには、以下のように.gitignoreを書きます。 .gitignoreファイルがない場合は、ファイルを作ってディレクトリに配置しましょう。
[0-9]*/ !0
そしてふたたび
git status
をすると、
On branch master Untracked files: (use "git add <file>..." to include in what will be committed) .gitignore nothing added to commit but untracked files present (use "git add" to track)
というメッセージが出てきます。 .gitignoreがaddされてないよーというメッセージですね。 .gitignoreはバージョン管理したほうがいいので、追加しましょう
git add .gitignore # .gitignoreを選択 git commit -m "[add] git ignore" # コミット
ではここでgit logを見てみましょう
git log
commit 6a3b160619f88ec23816b707e32e61b2ade96848 Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:08:56 2017 +0900 [add] git ignore commit 30351f767bce8d5665d0a1dc95f3d1a35fa1d927 Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:00:40 2017 +0900 first commit
2つのコミットが確認できますね
では、計算条件を変えてみましょう。
変更を加えてそれをcommitする
流速を100m/sにしてみましょうかね
0/Uの該当する部分を100m/sに変更してみましょう。
そして、
git diff
コマンドを入力します。 すると前回の計算と今のどこが違うのかわかるようになります。
diff --git a/0/U b/0/U index a12d708..5a2a2a0 100644 --- a/0/U +++ b/0/U @@ -23,7 +23,7 @@ boundaryField movingWall { type fixedValue; - value uniform (1 0 0); + value uniform (100 0 0); } fixedWalls
解説不要ですね。では
icoFoam
計算は止まります。流速が早すぎるからですね では、その旨をコミットしてみましょう
git add --all git commit -m "流速を100m/sに変更 計算が止まりました。" git log
commit b699f012c914037fb499dad001737234cbe2478f Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:18:28 2017 +0900 流速を100m/sに変更 計算が止まりました。 commit 6a3b160619f88ec23816b707e32e61b2ade96848 Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:08:56 2017 +0900 [add] git ignore commit 30351f767bce8d5665d0a1dc95f3d1a35fa1d927 Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:00:40 2017 +0900 first commit
追加されましたね。
計算が失敗したので、計算が動いていた頃に戻したいです。そういう場合はgit reset
を使用します。
"[add] git ignore"のcommit id が "commit b699f012c914037fb499dad001737234cbe2478f" となっていますので、最初数文字をとって
git reset --hard 6a3b160
としてやります。このidは人によって違うと思うので、自分でlogをして確認してください。 0/Uのディレクトリを見ていると、確かに流速が1m/sに戻ってますね。
ではもう一度icoFoamを実行してみましょう
icoFoam
計算が回りましたね。
ケースファイルをコピーして編集してみる
こういう事例、ありませんか?
ある計算Aの計算結果がイマイチでした。 原因の候補はいくつか思い当たるのですが、どれが正しいのかはわかりません。
よくある話だと思います。 この際によくやる手続きとして、
- 計算結果を"含めずに"コピーを行う
- 編集を行う。
- 良かった結果を今後採用する
があると思います。
では順番にやっていきましょう。
計算結果を含めずにコピーする
計算結果をgitignoreしてますので、実は
cd .. # gitで管理してるディレクトリの一個上へ git clone cavity cavity-A git clone cavity cavity-B
こんな感じで複製できます。計算結果は入っていないはずです。
編集を行う。
それでは、cavity-A
を編集しましょう。
cavity-Aでは、出力タイムステップを短くしてみましょう.
このように、cloneしてファイルを編集する場合は、branchを切ります。利点は後でわかります。
バージョンが枝分かれするイメージですね
git checkout -b fine-delta
こうするとfine-delta
というブランチを作ってそのブランチに移動することができます。
master
に戻りたいときは
git checkout master
となります。移動したらまたfine-deltaに戻っておきましょう
git checkout fine-delta
では、controlDictを以下のように編集して
diff --git a/system/controlDict b/system/controlDict index 084dc33..b8f6631 100644 --- a/system/controlDict +++ b/system/controlDict @@ -25,7 +25,7 @@ stopAt endTime; endTime 0.5; -deltaT 0.0005; +deltaT 0.005; writeControl timeStep;
icoFoam # ねんのため確認 git add --all git commit -m "出力タイムステップを短くしました"
良かった結果を今後採用する
よいな!と思ったら元のディレクトリに変更に還元しましょう
git push origin fine-delta
これはfine-deltaブランチを元のディレクトリにアップロードするという意味です。
元のディレクトリはmasterブランチ(デフォルト)で作成してますので、新しくブランチが生成されます。(即座にmasterに反映されるわけではありません。)
cavity-Bでも変更する
cavity-Bでも同湯に変更を加えます。 今回はmesh幅を変えましょう
diff --git a/constant/polyMesh/blockMeshDict b/constant/polyMesh/blockMeshDict index aac305a..5fdcfd2 100644 --- a/constant/polyMesh/blockMeshDict +++ b/constant/polyMesh/blockMeshDict @@ -30,7 +30,7 @@ vertices blocks ( - hex (0 1 2 3 4 5 6 7) (20 20 1) simpleGrading (1 1 1) + hex (0 1 2 3 4 5 6 7) (40 40 1) simpleGrading (1 1 1) ); edges
倍細かくなりましたね
git checkout -b fine-mesh blockMesh icoFoam git add --all git commit -m "メッシュを細かくしました" git push origin fine-mesh
元のディレクトリで変更を反映する
cd .. cd cavity # もとのディレクトリ git br
とすると、
fine-delta fine-mesh * master
となっていることがわかります。 時間間隔を短くしたやつとメッシュを細かくしたやつがいますね。 ここから
git checkout fine-mesh
とすると、メッシュが細かい計算のケースファイルに移行しますし、
git checkout fine-delta
とるすと出力の時間間隔が短い計算のケースファイルに移行することができます。 もちろん
git checkout master
とすると、大本の状態にもどることができます
さて、それでは変更をmasterに反映します。反映するためには、merge
を使います
git checkout master # masterに移動 git merge fine-delta
これでmasterに時間間隔が短くなる変更が適応されました。
git log
commit 39611f0253452a1165b26b62fa358435e41cac2c Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:31:28 2017 +0900 出力タイムステップを短くしました commit 6a3b160619f88ec23816b707e32e61b2ade96848 Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:08:56 2017 +0900 [add] git ignore commit 30351f767bce8d5665d0a1dc95f3d1a35fa1d927 Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:00:40 2017 +0900 first commit
logを見てもいい感じですね。 実際に実行してみると
icoFoam
確かに時間間隔が短くなってます。
mergeができたものは削除してOKです
git br -d fine-delta
それでは、メッシュの方も反映しちゃいましょう
git merge fine-mesh
先ほどとは違い、なんかテキストエディタが出てきますが今回は気にせず閉じちゃってOKです
これは"mergeが走る"という現象ですが、説明するのはそろそろ手がつかれたので割愛します。 変更が衝突したりした時に、いろいろしなければならくなるやつです。(例えばcavity-Aとcavity-Bで同じUを編集した場合、どちらを採用するのか等)
git log
をみてみると
commit 24e3b70e11e407086533ed4ae9d42b096d71e6e3 Merge: 39611f0 35c342c Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:49:32 2017 +0900 Merge branch 'fine-mesh' commit 35c342c111ecd706f10b0bf6c0fdbc8a28800578 Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:43:17 2017 +0900 メッシュを細かくしました commit 39611f0253452a1165b26b62fa358435e41cac2c Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:31:28 2017 +0900 出力タイムステップを短くしました commit 6a3b160619f88ec23816b707e32e61b2ade96848 Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:08:56 2017 +0900 [add] git ignore commit 30351f767bce8d5665d0a1dc95f3d1a35fa1d927 Author: kurenaif <antyobido@gmail.com> Date: Thu Nov 2 14:00:40 2017 +0900 first commit
のように、なんか変なのが混じってますがメッシュを細かくした内容が反映されています。
icoFoam paraFoam
とすると、確かに時間ステップも短いし、メッシュも細かいケースができました。
おわりに
研究する上で、研究ノートやログを取る。という試行錯誤を残すことは非常に重要なことです。 今回のように、gitを使えば、直接ファイルを見ながらバージョン管理できますし、コメントアウトしながら計算条件をバックアップする必要はなくなります。
commitごとに計算結果を残しておきたい場合は, commit idに対応したディレクトリを用意して、gitignoreすると良いでしょう。 そのcommit idに対応したcommitに移動した時にそれを戻してくればいいです。
このようなgitの使い方はまだまだ初歩の初歩で、git-flowやgitHubを有効にあつかったGitHub-Flow. 編集の衝突などちょくちょく問題が起きるのでそれの対応が必要になります。
もしこの記事がgitの使用価値について気づける機会になれば幸いです。
SnackDown Elimination 2017 PLUSMUL
問題要約(相方担当)
- N個の数列が与えられる
- それぞれの数の間に"*" か "+"を入れた式を作る(通りの式ができる)
- 個の式の総和 MOD を出力せよ.
概要
なんとなくDPで解けそうだったので適当にやったら遷移できた
実験
を,数列を左から個読み込んだ時の式の総和とする.
とする
1つの要素の数列が与えられ,その間に"*"か"+"を挿入すると,
またこの総和をとする.
2つの要素の数列が与えられ,その間に"*"か"+"を挿入すると,
またこの総和をとする.
3つの要素の数列が与えられ,その間に"*"か"+"を挿入すると,
またこの総和をとする.
0->1の遷移に注目すると
1->2の遷移に注目すると
2->3の遷移に注目すると
との部分とそれ以外の部分()に分けることができた.
それ以外の部分もDPで求めることが出来る 具体的な証明や導出は省略するが,今回の場合は
と新たに増えた要素を使って簡単に表記できる.に係数2が付いているが,これは注目している数列の個数をとすると,となる.
大体法則性がわかったので一般化する.N個の数列をではなく, と表記する.
に関するDP遷移は以下のようになる.
では
に関するDP遷移は以下のようになる.
では
ここで,にかかってしまうので,計算しながらメモすることによってにする.
総計算オーダーはとなり,これで間に合う.
ソースコード
実は実装ミスって上のやつと添字がちょっと違う
void solve() { int N; cin >> N; vector<int> v; REP(i, N) v.push_back(in()); vector<int> twoPow(N); vector<int> m(N); twoPow[0] = 1; FOR(i, 1, N) twoPow[i] = twoPow[i - 1] * 2 % MOD; vector<int> dp(N), dpsum(N); dp[0] = v[0]; dpsum[0] = v[0]; m[0] = v[0]; FOR(i,0,N-1){ m[i + 1] = ((m[i] * v[i + 1])%MOD + (twoPow[i] * v[i + 1])%MOD)%MOD; dp[i + 1] = (dpsum[i] + m[i+1])%MOD; dpsum[i + 1] = (dpsum[i] + dp[i + 1])%MOD; } cout << dp[N - 1] << endl; }
いろんな言語でSI単位系の計算をやってみよう(Cpp, python, rust, Julia)(WIP)
C++
C++なのでコンパイル時になんやかんやできるやつを探してたらBoostであった
Boost.Units を使えば良い.(http://www.boost.org/doc/libs/1_64_0/doc/html/boost_units.html)
sample
使用例として,簡単な問題を解いてみる.
h=634m地点から初速v0=2.4m/sでm=5kgのボールを落とす. ボールにかかる力は重力のみとすると,何秒で地面に到達するか.
#include <iostream> #include <boost/units/systems/si.hpp> #include <boost/units/io.hpp> #include <boost/units/pow.hpp> #include <boost/units/cmath.hpp> #include <boost/mpl/arithmetic.hpp> using namespace boost::units; // create new unit!! typedef boost::mpl::times<length_dimension,mass_dimension>::type LM_type; typedef unit<LM_type, si::system> kurenaif; int main() { const quantity<si::length> height = 634.0 * (si::meter); //x0 = 634m const quantity<si::mass> m = 5 * (si::kilogramme); const quantity<si::velocity> v0 = 2.4 * (si::meter / si::second); // const quantity<si::force> F = m*v0; compile error!!! because of invalid unit type const quantity<si::acceleration> g = 9.8 * (si::meter/pow<2>(si::second)); //create new unit!!! const quantity<kurenaif> test = height*m; const quantity<si::force> F = m*g; const quantity<si::time> t = -(v0 - sqrt(v0*v0+2.0*g*height))/g; std::cout << "height(x0) = " << height << std::endl; std::cout << "weight = " << m << std::endl; std::cout << "init. velocity = " << v0 << std::endl; std::cout << "g: " << g << std::endl; std::cout << "kurenaif:" << test << std::endl; std::cout << "force:" << F << std::endl; std::cout << "time:" << t << std::endl; return 0; } /* output height(x0) = 634 m weight = 5 kg init. velocity = 2.4 m s^-1 g: 9.8 m s^-2 kurenaif:3170 m kg force:49 m kg s^-2 time:11.1326 s */
memo
- 次元間違えはコンパイル時エラー
- auto型使えます
Juliaで競技プログラミング(もうしない, ABC061)
Motivation
研究でちょっと触っているのでABCで使ってみた
解説は色々みたらわかるので割愛
A問題
http://abc061.contest.atcoder.jp/tasks/abc061_a
A,B,C = parse.(split(readline(STDIN))) if A <= C && C <= B println("Yes") else println("No") end
http://abc061.contest.atcoder.jp/submissions/1283186
B問題
http://abc061.contest.atcoder.jp/tasks/abc061_b
N,M = parse.(split(readline(STDIN))) # ブロードキャスト 参考: http://bicycle1885.hatenablog.com/entry/2016/12/13/205646 count = zeros(Int64, N) for i in 1:M a,b = parse.(split(readline(STDIN))) count[a] += 1 count[b] += 1 end for value in count println("$value") end
http://abc061.contest.atcoder.jp/submissions/1284568
C問題 (TLE)
mapのかわりにdictを使ってみたがTLEしてしまった.
http://abc061.contest.atcoder.jp/tasks/abc061_c
N,M = split(readline(STDIN)) N = parse(N) M = parse(M) dict = Dict{Int64,Int64}() for i in 1:N a,b = split(readline(STDIN)) a = parse(a) b = parse(b) dict[a] = get(dict, a, 0) dict[a] += b end for key in sort(collect(keys(dict))) M -= dict[key] if M <= 0 println("$key") break end end
http://abc061.contest.atcoder.jp/submissions/1285715
所感
A問題すらすごい時間がかかって謎だった C問題は解法はO(NlogN)だしWAだとしてもTLEはしないはずなんだけど… 原因が分かるまで競技プログラミングでjuliaは封印します.
rustをインストールから行列で連立方程式を解くところまで
Motivation
なぜrust?
- 最近C++がつらい
- rustが流行っているらしい
- https://imoz.jp/note/rust-functions.html
環境
- ubuntu 16.04 64bit
rustのインストール
このコマンドで簡単にインストールする
curl https://sh.rustup.rs -sSf | sh
以下のようなメッセージができる. 初心者なのでdefaultのstableをinsatllする.
default host triple: x86_64-unknown-linux-gnu default toolchain: stable modify PATH variable: yes 1) Proceed with installation (default) 2) Customize installation 3) Cancel installation
インストールが終わると以下のメッセージが出てくるので
To configure your current shell run source $HOME/.cargo/env
.bashrc
かなんかに
source $HOME/.cargo/env
を追加する
rustc --version
を入力して,それらしいものがでてきたらとりあえずOK.
References
crateの作成
Rust
では create
(木箱) 単位でプログラムを開発するらしい
今回は簡単な行列計算をやってみるので,mat_calc
という名前でcrateを作る
cargo new mat_calc --bin
なんかできた
output: Created binary (application) `mat_calc` project
ディレクトリを覗いてみると
$ ls mat_calc . .. Cargo.toml .git .gitignore src
どうやら.git
も生成されているらしいぱっと見た感じ.gitignore
もいい感じに定義されてた
実行テスト
src/main.rs
にはすでにHello Worldを出力するプログラムが入っているのでとりあえず実行してみる
とりあえずディレクトリに移動して
$ cd mat_calc $ cargo run Compiling mat_calc v0.1.0 (file:///home/kurenaif/Desktop/mat_calc) Finished dev [unoptimized + debuginfo] target(s) in 0.53 secs Running `target/debug/mat_calc` Hello, world!
無事出力に成功した
References
dependenciesの設定
行列を扱いたいけどFortranでもないしrustにはそんなライブラリは標準搭載されていない. rustで行列が扱えるライブラリを探したら以下の3つが見つかった
- https://github.com/sebcrozet/nalgebra
- https://github.com/termoshtt/ndarray-linalg
- https://github.com/masonium/linxal
一番使いやすそうで,commitが盛んなnalgebra
を使うことにした.
速度面はよくわからなかった.
ここで, Cargo.toml
を編集し,nalgebraを使用できるようにする(pip install
に近いイメージ)
[dependencies] nalgebra = "0.11.0"
versionは2017-05-05の時点での最新バージョンを使用した.
githubなりなんなりをみてうまいこと指定する.
githubのURLでの指定も可能だけど https://crates.io
に登録されているのでこんな感じでURLなしでもcrateを入れることが出来る.
References
連立方程式を解いてみる
行列の表示
以下のソースコードをsrc/main.rs
に実装し,実行してみる
今後の流体計算を考慮して普通のMatrix
よりも実行時に長さを決められるDMatrixのほうが都合がよさそうなのでDMatrix
を使用する.
使い方その他もろもろはソースコードみて察していただきたい.
個々で気をつけなければならないのが行列の数字はすべて実数で表記しないとちゃんとprintlnしてくれない.
つまり 1,2,3… のような整数を書く場合でも, 1., 2., 3. と書いたほうが後々楽になる.
rustは整数と小数,また同じ整数どうしや小数どうしてもbit幅が違えばコンパイルエラーを出すほど厳密.
extern crate nalgebra as na; fn main(){ let mat = na::DMatrix::from_iterator(3,3,[ 1.,2.,3., 4.,5.,6., 7.,8.,9. ].iter().cloned()); println!("{}",mat); let b = na::DVector::from_iterator(3,[ 1.,2.,3. ].iter().cloned()); println!("{}",b); }
$ cargo run ┌ ┐ │ 1.000 4.000 7.000 │ │ 2.000 5.000 8.000 │ │ 3.000 6.000 9.000 │ └ ┘ ┌ ┐ │ 1.000 │ │ 2.000 │ │ 3.000 │ └ ┘
- めっちゃいい感じの出力をしてくれる.(DMatrixじゃなくて普通のMatrixはこうはいかない)
- どうやらソースコードに書いたものに比べて転置したものが出てくるらしい 縦ベクトルが3つ横に並んでるイメージ.
逆行列
方程式を解くためには逆行列が必要不可欠 さっき用意したmat
は実は逆行列が存在しない.
面白そうなのでさっきのmat
の逆行列を計算してみる.
nalgebraでは,.try_inverse()
を用いて逆行列を計算する.
無理だったら無理なりの挙動をしてくれる.
try_inverseを出力するプログラム.
println!("{}")
と println!("{:?}")
では呼びされるメソッドが微妙に異なっており,後者の方がちょっと詳しいメッセージが出てくることが多い(デバッグ用出力?).
今回,前者を使うとエラーを吐くので後者を使っている.
extern crate nalgebra as na; fn main(){ let mat = na::DMatrix::from_iterator(3,3,[ 1.,2.,3., 4.,5.,6., 7.,8.,9. ].iter().cloned()); println!("{:?}",mat.try_inverse()); }
$ cargo run None
Noneが出力された.
リファレンスに書いているが,try_inverse()
をしてなかった時はNone
が返ってくるらしい
では逆行列が存在する時はどんな出力をするのだろうか?
今回はtry_inverseの結果を"{}"
と"{:?}"
の二つを出力している.
最初のいもすさんの記事でも書いてある通り,rustでは基本所有権がmoveしちゃうので最初のprintln!
でmat
の所有権が切れる.
そこで,clone()
することで所有権を継続させている.
extern crate nalgebra as na; fn main(){ let mat = na::DMatrix::from_iterator(3,3,[ 1.,3.,5., 2.,3.,8., 2.,1.,5. ].iter().cloned()); println!("{:?}",mat.clone().try_inverse()); println!("{}",mat.try_inverse().unwrap()); }
Some(Matrix { data: MatrixVec { data: [1.4, -2, 1.8, 1.2, -1, 0.4, -0.8, 1, -0.6], nrows: Dynamic { value: 3 }, ncols: Dynamic { value: 3 } }, _phantoms: PhantomData }) ┌ ┐ │ 1.400 1.200 -0.800 │ │ -2.000 -1.000 1.000 │ │ 1.800 0.400 -0.600 │ └ ┘
一つ目の出力がtry_inverse()
の中身を詳細に出力したものである.
Some()
がついていて,実はこのままでは{}
の出力ができない.
ここでSome()
は値がNullになるかもしれない的なやつらしい.
Some()
を取るためには,unwrap()をすると良い.
詳しくはここを参照する:
これの結果があってるかわからないので元のmatにかけて単位行列が出ること確認する
rustでは演算子のオーバーロードができる.
このライブラリでは行列の掛け算もちゃんと定義されているので,以下のソースコードみたいな感じで表記できる.
extern crate nalgebra as na; fn main(){ let mat = na::DMatrix::from_iterator(3,3,[ 1.,3.,5., 2.,3.,8., 2.,1.,5. ].iter().cloned()); println!("mat\n{}",mat.clone()); println!("inv(mat)\n{}",mat.clone().try_inverse().unwrap()); println!("mat*inv(mat)\n{}",mat.clone().try_inverse().unwrap()*mat); }
$ cargo run mat ┌ ┐ │ 1.000 2.000 2.000 │ │ 3.000 3.000 1.000 │ │ 5.000 8.000 5.000 │ └ ┘ inv(mat) ┌ ┐ │ 1.400 1.200 -0.800 │ │ -2.000 -1.000 1.000 │ │ 1.800 0.400 -0.600 │ └ ┘ mat*inv(mat) ┌ ┐ │ 1.000 -0.000 0.000 │ │ 0.000 1.000 0.000 │ │ 0.000 0.000 1.000 │ └ ┘
ちゃんと逆行列が計算されているっぽい
連立方程式を解いてみる
いよいよ本題
とりあつかう問題はこれの例1: 連立方程式の解き方
連立方程式を $$ x = A^{-1} b$$ の形にしてinvして解くだけ
extern crate nalgebra as na; fn main(){ let a = na::DMatrix::from_iterator(3,3,[ 1.,1.,0., 2.,-1.,1., 3.,2.,1. ].iter().cloned()); println!("A\n{}",a); let b = na::DVector::from_iterator(3,[ 0.,3.,-1. ].iter().cloned()); println!("B\n{}",b); let invA = a.try_inverse().unwrap(); println!("invA\n{}",invA); println!("x\n{}", invA*b) }
A ┌ ┐ │ 1.000 2.000 3.000 │ │ 1.000 -1.000 2.000 │ │ 0.000 1.000 1.000 │ └ ┘ B ┌ ┐ │ 0.000 │ │ 3.000 │ │ -1.000 │ └ ┘ invA ┌ ┐ │ 1.500 -0.500 -3.500 │ │ 0.500 -0.500 -0.500 │ │ -0.500 0.500 1.500 │ └ ┘ x ┌ ┐ │ 2.000 │ │ -1.000 │ │ 0.000 │ └ ┘
ページの答えと一致してるしちゃんと解けてるっぽい.
おまけ
今回色々やったリポジトリをここにおいておきます github.com
次回は簡単な熱拡散問題とか解いてみるかも
OpenFOAM-2.4.0をmacのgccでmakeする
Motivation
Macはclangをつかってコンパイルをしている。 どうやらコンパイラによって微妙に挙動が違うらしくMacでgccが使いたいと先生に言われて色々やった。
OpenFOAM-2.4.0向け
概要
流れ的には
といった順番になる。
一つづつやっていこう
gccを入れる
自分の環境と合わせたかったので最新のものではなくgcc-5.4.0をぶっこんだ
ソースコードからbuildする
まずはソースコードを入手し、解凍
$ wget http://ftp.tsukuba.wide.ad.jp/software/gcc/releases/gcc-5.4.0/gcc-5.4.0.tar.gz $ tar zxvf gcc-5.4.0.tar.gz
解凍したディレクトリに移動して依存ファイル取得
$ cd gcc-5.4.0 $ contrib/download_prerequisites
build
ディレクトリを作成し、configure, build
$ mkdir build $ cd build $ ../configure --enable-languages=c,c++,fortran --prefix=/usr/local/gcc-5.4.0 --disable-bootstrap --disableultilib
$ make -j4 $ make install
/usr/local/binにgccをln -s
$ ln -s /usr/local/gcc-5.4.0/bin/gcc /usr/local/bin/gcc $ ln -s /usr/local/gcc-5.4.0/bin/g++ /usr/local/bin/g++
以下をbashrcに追記
export PATH=/usr/local/bin:$PATH
確認
$ gcc --version gcc (GCC) 5.4.0 Copyright (C) 2015 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE
参考: http://qiita.com/liveralmask/items/6ed4a98ebb3bf6b7f707
openmpiを入れる
新しすぎるopenmpiはなんかエラー吐いてめんどくさいのでちょっと古いの使う
$ wget https://www.open-mpi.org/software/ompi/v1.10/downloads/openmpi-1.10.2.tar.gz $ tar zxvf openmpi-1.10.2.tar.gz
参考: https://www.open-mpi.org/software/ompi/v1.10/
buildディレクトリ作ってconfigureしてmake
$ mkdir build $ ../configure --prefix=/usr/local/openmpi $ make $ make install
以下を.bashrcに追加
export MANPATH=/usr/local/openmpi/share/man:$MANPATH export LD_LIBRARY_PATH=/usr/local/openmpi/lib:$LD_LIBRARY_PATH export PATH=/usr/local/openmpi/bin:$PATH
チェック
$ mpicc --version
参考: http://qiita.com/himo/items/c30d83d0f7642fb3af57
flexを入れる
これはmake – installでもbrewでもなんか無理だった macportを使う
sudo port install flex
OpenFOAMを入れる
なんかいろいろする
$ hdiutil create -size 8.3g -type SPARSEBUNDLE -fs HFSX -volname OpenFOAM -fsargs -s OpenFOAM.sparsebundle $ mkdir -p OpenFOAM $ hdiutil attach -mountpoint $HOME/OpenFOAM OpenFOAM.sparsebundle $ cd OpenFOAM
OpenFOAM2.4.0をダウンロード、解凍する
$ wget http://downloads.sourceforge.net/foam/OpenFOAM-2.4.0.tgz?use_mirror=mesh $ tar zxvf OpenFOAM-2.4.0.tar.gz
参考: https://openfoam.org/download/2-4-0-source/ パッチを当てる
$ cd OpenFOAM-2.4.0 これをダウンロード -> https://sourceforge.net/p/openfoam-extend/svn/HEAD/tree/trunk/Breeder_2.4/distroPatches/MacPatch/OpenFOAM-2.4.x-Mac.patch $ patch -p1<OpenFOAM-2.4.x_Mac.patch
make
./Allwmake
エラーが出るのでsrc/OpenFOAM/lnIncludeにFlexLexer.hを無理やりリンクさせて通す
$ ln -s /opt/local/include/FlexLexer.h $FOAM_SRC/OpenFOAM/lnInclude
もう一度make
./Allwmake
チェック
which icoFoam