読者です 読者をやめる 読者になる 読者になる

くれなゐの雑記

身の回りの人や自分が困っていたことを記事にしています

いろんな言語で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

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?

環境

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

www.rustup.rs

qiita.com

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

qiita.com

dependenciesの設定

行列を扱いたいけどFortranでもないしrustにはそんなライブラリは標準搭載されていない. rustで行列が扱えるライブラリを探したら以下の3つが見つかった

一番使いやすそうで,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

qiita.com

連立方程式を解いてみる

行列の表示

以下のソースコード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()をすると良い.

詳しくはここを参照する:

qiita.com

これの結果があってるかわからないので元の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をつかってコンパイルをしている。 どうやらコンパイラによって微妙に挙動が違うらしくMacgccが使いたいと先生に言われて色々やった。

OpenFOAM-2.4.0向け

概要

流れ的には

  1. gccを入れる
  2. openmpiを入れる
  3. flexを入れる
  4. OpenFOAMを入れる

といった順番になる。

一つづつやっていこう

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にgccln -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にチェックを入れて見ると良いです) f:id:kurenaif:20170326220734j:plain

使用するソースコード

ここに書いてるソースコードを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";
}


// ************************************************************************* //

f:id:kurenaif:20170326220641p:plain

ソースコードの解説

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のソースコード

PENGUINITIS - 領域分割

PENGUINITIS - OpenFOAM プログラミングメモ

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

f:id:kurenaif:20170227012459g:plain


解説

このコードは1から書いたものではありません. Paraviewで生成されたものをどうにかするとこんな感じになります.

Paraviewでおおまかなソースコードの作成

  1. paraFoamを起動する
  2. cavity.OpenFOAMをDeleteする
  3. [MenuBar]->[Tools]->[Start Trace]を起動する
  4. cavity.OpenFOAMを再び読み込む(MenuBar->File->Recent Filesから読み込むと早いです)
  5. Uの可視化を行う.
  6. [MenuBar]->[Tools]->[Stop Trace]を押す

今こんな感じだと思います. 2番目はPipeline Borwserの[ Delete]を押すって意味です意味不明だとお思いますが,これをすると後が楽です.

f:id:kurenaif:20170227170820p:plain

で,同時に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)
時間ステップの変更,再描画

commands moduledで外部コマンドであるfoamListTimesを呼び,描く時間ステップを取得 AnimationTimeを設定し,Render() WriteImage()を呼ぶことで,描く時間を取得し,animationを出力することができます. 最終的にこの連番pngimagemagickなり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.2.0-RC4-23-g404b27b documentation

(2017-03-26追記) pvbatchの並列実行

mpirun -np n pvbatch filename.py

で並列実行できるらしいです

http://www.opencae.or.jp/wp-content/uploads/2015/06/course20111201handout.pdf