いろんな言語で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
decomposeParでmanual切りをする
注意: バージョンによってファイル形式が異なる可能性があります.(OpenFOAM-2.4.0)
Motivation
需要は少ないですが,manualでdecomposeParをする必要がありました. その手法について資料がないので解説します
ケースファイル
tutorials/incompressible/icoFoam/cavity
メッシュ作成
blockMesh
decomposePar(自動領域分割)
tutorials/incompressible/pisoFoam/les/motorBike/motorBike/system
を少しいじったものです
以下を使用します
system
以下に入れておくと良いです.
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object decomposeParDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // numberOfSubdomains 4; // numberOfSubdomains 3; method simple; simpleCoeffs { n ( 2 2 1 ); delta 0.001; } manualCoeffs { dataFile "cellDecomposition"; } // ************************************************************************* //
領域分割を実行します cellDist
を指定すると後でどういうふうに領域分割したのかを見れます.
decomposePar -cellDist
OpenFOAMで領域分割結果を確認します.(cellDistにチェックを入れて見ると良いです)
使用するソースコード
ここに書いてるソースコードをwmakeしてcavityのディレクトリで実行したら動きます.
https://github.com/kurenaif/manualDecomposeTutorial
これを実行すると,constant/cellDecomposition
が生まれます.
その中にどのセルがどのprocessorに該当するのかが一セルずつ記述されています.
また,manualでdecomposeをするためには,system/decomopsePar
を弄る必要があります.
具体的には,
- 今回は3分割のため,numberOfSubdomainsを3に
- methodをmanual
に変更します.変更したら以下のようになります
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 2.4.0 | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object decomposeParDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // numberOfSubdomains 3; method manual; simpleCoeffs { n ( 2 2 1 ); delta 0.001; } manualCoeffs { dataFile "cellDecomposition"; } // ************************************************************************* //
ソースコードの解説
argsの準備
コマンドライン引数をOpenFOAMのものに変換する.
Foam::argList args(argc,argv); if(not args.checkRootCase()) Foam::FatalError.exit();
runTimeの準備
meshを用意するために必要. argsはここで必要
Foam::Time runTime(Foam::Time::controlDictName, args);
meshの準備
Foam::fvMesh mesh ( Foam::IOobject ( Foam::fvMesh::defaultRegion, runTime.timeName(), runTime, Foam::IOobject::MUST_READ ) );
各セルがどのprocessorに該当するのかを記述するための変数の用意
OpenFOAMのListはconstructorで渡してやるとそのサイズを用意してくれる.
mesh.cells().size()
はメッシュのセルのサイズを返してくれるので,この容量があればそれぞれのセルがどのprocessorに該当するのかを収容できる.
labelList procIds(mesh.cells().size());
各セルをループで回しそれぞれのセルにprocessorIDを振り分ける
回すだけ 斜めと円形に切ってます
forAll(mesh.cells(), cid){ const vector &C = mesh.C()[cid]; if(C[0]*C[0]+C[1]*C[1] < 0.01){ procIds[cid] = 1; } else{ procIds[cid] = 0; } if(C[0] + C[1] < 0.1){ ++procIds[cid]; } }
出力
IOobjectを使います. ちょっと詳細はよくわかりませんがprocIdsを渡してやるとどうやら出力してくれるらしいです.
labelIOList cellDecomposition ( IOobject ( "cellDecomposition", mesh.facesInstance(), mesh, IOobject::NO_READ, IOobject::NO_WRITE, false ), procIds ); cellDecomposition.write();
References
decomopseParのソースコード
OpenFOAMで計算結果以外をコピーするスクリプト
constantとおかsystemとかを残してコピーします 数字のファイルの判別方法が雑なのでミスってるかも…
一応これによって生じた不利益等は責任は負いません.
#!/bin/bash if [ $# -ne 2 ]; then echo "usage: $0 source dist [first step]" exit 1 fi mkdir $2 find $1 -maxdepth 1 | grep -vx $1 | grep -v [0-9]$ | grep -v log | grep -v postProcessing | xargs cp -r -t $2 if [ $# -ne 3 ]; then cp -r $1/0 $2/ else cp -r $1/$3 $2/ fi
pvpythonでアニメーション作成した
追記
2017-03-26 追記 --- pvbatchの並列実行
可視化対象
今回は適当にキャビティ流れを対象とします.
tutorials/incompressible/icoFoam/cavity
実行
blockMesh icoFoam
可視化をするソースコード
とりあえず実行
OpenFOAMの可視化はpvpython
を使って行います.
以下のファイルを適当に保存して,
pvbatch filename.py
とかで実行してみましょう するとaniフォルダが生成されて,連番のpngが生成されているはずです. それを適当に結合すると下の実行結果みたいな感じになります(writeInterbalをいじっているのでちょっと細かく出てます.)
from paraview.simple import * cavity_OpenFOAM = PV4FoamReader( FileName='./cavity.OpenFOAM' ) cavity_OpenFOAM.UiRefresh = 0 cavity_OpenFOAM.VolumeFields = ['p', 'U'] cavity_OpenFOAM.MeshParts = ['internalMesh'] RenderView1 = GetRenderView() RenderView1.CenterAxesVisibility = 0 RenderView1.OrientationAxesVisibility = 0 RenderView1.Background = [1.0, 1.0, 1.0] DataRepresentation6 = Show() DataRepresentation6.EdgeColor = [0.0, 0.0, 0.5000076295109483] DataRepresentation6.SelectionPointFieldDataArrayName = 'p' DataRepresentation6.SelectionCellFieldDataArrayName = 'p' DataRepresentation6.ScalarOpacityUnitDistance = 0.01924175606617764 DataRepresentation6.ExtractedBlockIndex = 2 DataRepresentation6.ScaleFactor = 0.010000000149011612 RenderView1.CameraClippingRange = [0.2611983736150351, 0.2905205542589584] DataRepresentation6.ScalarOpacityFunction = [] DataRepresentation6.ColorArrayName = ('CELL_DATA', 'U') DataRepresentation6.ColorAttributeType = 'CELL_DATA' a3_U_PVLookupTable = GetLookupTableForArray( "U", 3, RGBPoints=[0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0] ) at = AnnotateTime() annotateShow = Show(at) annotateShow.Color = [0.0,0.0,0.0] view = GetActiveView() size0 = view.ViewSize view.ViewSize = [1024, 768] l = servermanager.rendering.ScalarBarWidgetRepresentation() l.LookupTable = a3_U_PVLookupTable l.Title = 'U' l.Enabled = 0 l.TitleFontSize = 12 l.LabelFontSize = 10 l.LabelColor = [0.0,0.0,0.0] l.TitleColor = [0.0,0.0,0.0] l.AutomaticLabelFormat = 0 l.LabelFormat = '%1.1f' view.Representations.append(l) AnimationScene1 = GetAnimationScene() AnimationScene1.EndTime = 0.5 AnimationScene1.PlayMode = 'Snap To TimeSteps' import commands import os if not os.path.exists("./ani"): os.mkdir("./ani") times = map(float, commands.getoutput("foamListTimes").split("\n")) for time in times: AnimationScene1.AnimationTime = time Render() WriteImage('./ani/U'+str('%.04f' % time)+'.png')
output
解説
このコードは1から書いたものではありません. Paraviewで生成されたものをどうにかするとこんな感じになります.
Paraviewでおおまかなソースコードの作成
- paraFoamを起動する
- cavity.OpenFOAMをDeleteする
- [MenuBar]->[Tools]->[Start Trace]を起動する
- cavity.OpenFOAMを再び読み込む(MenuBar->File->Recent Filesから読み込むと早いです)
- Uの可視化を行う.
- [MenuBar]->[Tools]->[Stop Trace]を押す
今こんな感じだと思います. 2番目はPipeline Borwserの[☓ Delete]を押すって意味です意味不明だとお思いますが,これをすると後が楽です.
で,同時にScript Editorからなんか出てくると思います. そのスクリプトを回すと今回起きた内容を再現できます.
簡易的なソースコードの解説といろいろ付け足し
これをベースにいろいろ付け加えていくと,上記のソースコードのようになるのですが,一部わかっておいたほうがいいものがあるので少し解説します
caseの読み込み
cavity_OpenFOAM = PV4FoamReader( FileName='./cavity.OpenFOAM' )
軸の非表示,背景色変更
RenderView1 = GetRenderView() RenderView1.CenterAxesVisibility = 0 RenderView1.OrientationAxesVisibility = 0 RenderView1.Background = [1.0, 1.0, 1.0]
AnnotateTime追加
at = AnnotateTime() annotateShow = Show(at) annotateShow.Color = [0.0,0.0,0.0]
出力画像サイズの変更
view = GetActiveView() size0 = view.ViewSize view.ViewSize = [1024, 768]
色の設定,range設定
Blue -> Redの例 [min, r, g,b, max, r, g, b]の順
a3_U_PVLookupTable = GetLookupTableForArray( "U", 3, RGBPoints=[0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0] )
ScalarBar(右についてるカラーバー)の表示,設定
LookupTableはGetLookupTableForArrayで受け取ったものを表示する
l = servermanager.rendering.ScalarBarWidgetRepresentation() l.LookupTable = a3_U_PVLookupTable l.Title = 'U' l.Enabled = 0 l.TitleFontSize = 12 l.LabelFontSize = 10 l.LabelColor = [0.0,0.0,0.0] l.TitleColor = [0.0,0.0,0.0] l.AutomaticLabelFormat = 0 l.LabelFormat = '%1.1f' view.Representations.append(l)
カメラの座標の設定
カメラの座標はparaviewの3Dの右にあるカメラボタンを押せば,現在の情報が表示されます. xmlで吐くこともできるので,それを読み取るとスムーズです
view = GetActiveView() view.CameraViewUp = [0, 1, 0] view.CameraFocalPoint = [0.25, 0.0199999995529652, -0.00499999988824129] view.CameraViewAngle = 45 view.CameraPosition = [0.25,0.0199999995529652 , 0.224204409914699]
時間ステップの変更,再描画
commands moduledで外部コマンドであるfoamListTimesを呼び,描く時間ステップを取得
AnimationTimeを設定し,Render()
WriteImage()
を呼ぶことで,描く時間を取得し,animationを出力することができます.
最終的にこの連番pngをimagemagickなりffmpegなりで動画化すれば,動くようになります.
AnimationScene1 = GetAnimationScene() AnimationScene1.EndTime = 0.5 AnimationScene1.PlayMode = 'Snap To TimeSteps' import commands import os if not os.path.exists("./ani"): os.mkdir("./ani") times = map(float, commands.getoutput("foamListTimes").split("\n")) for time in times: AnimationScene1.AnimationTime = time Render() WriteImage('./ani/U'+str('%.04f' % time)+'.png')
References
[Paraview] Slices generation with pvpython
PENGUINITIS - ParaView Python Script によるポスト処理
ParaView’s Python documentation! — ParaView/Python 5.4.1-766-g40689f7 documentation
(2017-03-26追記) pvbatchの並列実行
mpirun -np n pvbatch filename.py
で並列実行できるらしいです
pvbatchで描画せずに画像だけ表示
pvbatch --use-offscreen-rendering filename.py
http://www.opencae.or.jp/wp-content/uploads/2015/06/course20111201handout.pdf