2017-04-04

OpenFOAM-2.4.0をmacのgccでmakeする

OpenFOAM

Motivation

Macはclangをつかってコンパイルをしている。どうやらコンパイラによって微妙に挙動が違うらしくMacでgccが使いたいと先生に言われて色々やった。

OpenFOAM-2.4.0向け

概要

流れ的には

gccを入れる
openmpiを入れる
flexを入れる
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に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

Motivation
概要
gccを入れる
openmpiを入れる
flexを入れる
OpenFOAMを入れる

2017-03-26

decomposeParでmanual切りをする

OpenFOAM

注意: バージョンによってファイル形式が異なる可能性があります.(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

ソースコードの解説

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 プログラミングメモ

2017-02-27

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

2017-02-27

pvpythonでアニメーション作成した

OpenFOAM

追記

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

可視化対象

今回は適当にキャビティ流れを対象とします．

tutorials/incompressible/icoFoam/cavity

実行

blockMesh
icoFoam

解説

このコードは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]を押すって意味です意味不明だとお思いますが，これをすると後が楽です．

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)

カメラの座標の設定

カメラの座標はparaviewの3Dの右にあるカメラボタンを押せば，現在の情報が表示されます． xmlで吐くこともできるので，それを読み取るとスムーズです f:id:kurenaif:20171019182544p:plain

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 moduleｄで外部コマンドである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

2017-02-26

gdbがちょっと便利になるcgdb

OpenFOAMをデバッグしようとしてClionやらEclipseやら試してたんだけどどうもうまく動かない

できればIDEでやりたいけどgdbを便利にしたものは無いものかと調べた結果がこれ

インストール方法

ubuntuならapt-getで入った． macならbrewで入るらしい

使い方

適当にicoFoamで使ってみます tutorialから適当なサンプルを持ってきて…

cp -r ~/OpenFOAM/OpenFOAM-2.4.0/tutorials/incompressible/icoFoam/cavity .

cgdbをかまして実行します

cgdb icoFoam

するとこんな感じ

f:id:kurenaif:20170226184436p:plain

上にソースコード，下にgdbが表示されています．いいですね

上に移動するのはescかctrl+[ 下に移動するのはi でできます vimlikeなキーバインドになってます

あとはソースコードをhjklで移動できたりspaceでbreakpointを設置できたりできます．キーバインドに関してはこちらが詳しいです

miettal.hatenablog.com

今日から楽しくLet’s gdb

Eclipseで使えるようにならないかなぁ…

2017-02-24

epslatexなgnuplotで保存するたびにpdfを更新してそれを見ながら編集する(gnuplot, omake)

デモ

ちょっとタイムラグと画質がきになりますがこんな感じです
今回はtexのフォントを使いたかったのでplatexを１回挟んでます

makeplt_tex demo

.plt(gnuplotのファイル) -> .eps+.tex -> .pdfファイル

以下のスクリプトを使います makeplt_texと名付けています
このコマンドはあくまでもpdfファイルを生成するものであって自動更新は別です
filename.pdfが生成されます
個人的趣味でA4サイズで出力して，全体のサイズを確認しながら作業してます．

 makeplt_tex filename.plt

#!/bin/sh
set -e
PLTNAME=$(basename $1)
DIR=$(dirname $1)
gnuplot<<EOF
cd "$DIR"
set terminal epslatex color 
set output "${PLTNAME%plt}tex"
load "$PLTNAME"
EOF

OUTDIR=source_${PLTNAME%.plt}

mkdir -p $OUTDIR

cat << EOS > $OUTDIR/main.tex
\documentclass[a4j,10pt]{jsarticle}
 
\usepackage[dvipdfmx]{graphicx}
 
\begin{document}
\begin{figure}[htb]
	\centering
	\input{${PLTNAME%plt}tex}
	\caption{caption}
	\label{fig:${PLTNAME%.plt}}
\end{figure}
 
\end{document}
EOS

cat << EOS > $OUTDIR/Makefile
main.pdf:
	platex main.tex
	dvipdfmx main.dvi
EOS

mv ${PLTNAME%plt}eps ${PLTNAME%plt}tex $OUTDIR
cd $OUTDIR

make

mv main.pdf ../${PLTNAME%plt}pdf

.pdfファイルの自動更新

ファイル監視ツールみたいなの作ったのですがいいやつがありました
omake(1. ガイド — OMakeマニュアル日本語訳)を使います
まず，omakeをapt-getなりmakeするなりしてinstallしたら該当ディレクトリに移動して，

omake --install

でOMakefileを生成します．
OMakefileを以下のように編集します

filename.pdf:
    makeplt_tex filename.plt

.DEFAULT: filename.pdf

あとは

omake -P

とするとファイルの監視が起動して，.pltファイルを編集するたびに.pdfが更新されるので，okularなりSumatraPDFなりskimなりで開いて更新していく様子を観察すればOKです

pltを各際の注意

動画で一度しくじってますが，バックスラッシュは二回必要です

2017-01-29

OpenFOAMでgdbを使ってデバッグっぽいことをしてみる

OpenFOAM

OpenFOAMでgdbを使いたい

Intro

Motivation

先日 OpenCAE勉強会でgdbが便利だという話だという話をしました。詳細について教えてほしいということなのでブログの記事を書くことにしました。

この記事では,

概要
デバッグオプションをいれたOpenFOAMのmake
実際にエラーの原因を探してみる
gdbOF(WIP まだ描いてないよ)

の4つを紹介していきます結構長いですが, 冷静に見ると内容薄いのでささっと読めると思います(出力結果おおいし)

Target

OpenFOAMのcavity流れのtutorialが回せるくらいの人
gdbが何かわからない人
OpenFOAMでgdbを使えることがしらなかった人

環境

ubutnu16.04, OpenFOAM-2.4.0

どんなことができるの

OpenFOAMで適当なものを実行してみる. デバッグをしたいので, なにか計算が落ちるケースがほしい. 今回は OpenFOAM/OpenFOAM-2.4.0/tutorials/incompressible/icoFoam/cavity のcavity流れの流速を $U=100m/s$ とかにして計算を飛ばす

通常のエラーメッセージ

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::symGaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ??:?
#4  Foam::symGaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
#5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#6  ? at ??:?
#7  ? at ??:?
#8  ? at ??:?
#9  ? at ??:?
#10  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11  ? at ??:?

デバッグオプション

#0  Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-2.4.0/src/OSspecific/POSIX/printStack.C:219
#1  Foam::sigFpe::sigHandler(int) at ~/OpenFOAM/OpenFOAM-2.4.0/src/OSspecific/POSIX/signals/sigFpe.C:108
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::symGaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ~/OpenFOAM/OpenFOAM-2.4.0/src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:167
#4  Foam::symGaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ~/OpenFOAM/OpenFOAM-2.4.0/src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:237
#5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ~/OpenFOAM/OpenFOAM-2.4.0/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C:169良いかと思います
#6  ? at ~/OpenFOAM/OpenFOAM-2.4.0/src/finiteVolume/lnInclude/fvMatrixSolve.C:193 (discriminator 9)
#7  ? at ~/OpenFOAM/OpenFOAM-2.4.0/src/finiteVolume/lnInclude/fvMatrixSolve.C:82
#8  ? at ~/OpenFOAM/OpenFOAM-2.4.0/src/finiteVolume/lnInclude/fvMatrixSolve.C:331
#9  ? at ~/OpenFOAM/OpenFOAM-2.4.0/src/finiteVolume/lnInclude/fvMatrix.C:1390
#10  ? at ~/OpenFOAM/OpenFOAM-2.4.0/applications/solvers/incompressible/icoFoam/icoFoam.C:63 (discriminator 7)
#11  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12  ? at ??:?

パット見で両者に違いがわかります.

通常のエラーメッセージ: どこでエラーが起きてるのかよくわからない. ギリギリメソッドはわかるけどどのファイル名かはわからない.
デバッグオプション込: どこのファイルで呼ばれたのかわかる. 通常のエラーメッセージでは#6くらいから完全に?になってるが, デバッグオプションでは, icoFoam.Cまで遡ることができる.

しかし, デバッグオプション込では実行速度が遅くなる欠点があります. 体感実行時間2倍は硬いです. なので基本的にはデバッグオプション無しで実行し, エラーが起きた時にデバッグオプションをつけてコンパイル, 実行が良いかと思います.

また, gdbを使えば

プログラムを途中で止めてその状態で変数等確認
1行ずつ実行
backtrace(今いる関数はどこから呼ばれたのかなどを表示)

などができます.

使ってみる

Debugオプションを入れた状態でのmake

OpenFoam-x.x.x/etc/bashrcの

#- Optimised, debug, profiling:
#    WM_COMPILE_OPTION = Opt | Debug | Prof
export WM_COMPILE_OPTION=Opt

を(Debugという文字はここにしかでてこないので, Debugという文字列で検索したらここがすぐ出てくると思います.)

#- Optimised, debug, profiling:
#    WM_COMPILE_OPTION = Opt | Debug | Prof
export WM_COMPILE_OPTION=Debug

に変更します.(Opt->Debug) この状態で, bashrcを読みなおします bashrcの読みなおしは,

source etc/bashrc

などで行います. そのあと

cd ~/OpenFOAM/OpenFOAM-x.x.x
./Allwmake

などで再コンパイルするとDebug用実行ファイルが生成されます. 例えば, OpenFOAM-2.4.0/src/finiteVolume/Make などのファイルを確認すると linux64GccDPDebug みたいなDebugっぽい感じのディレクトリが生成されます. OptとDebugでフォルダがわかれているので, bashrcのオプションをoptに戻すとまた全部コンパイルせずに済みますね.

その後適当に実行時エラーを吐くようになコードを書いたり, ケースを作ったりしたら先ほど述べたようなエラーメッセージが出ます.

DebugとOptを切り替える

bashrcをDebugにしたり, Optにしたり, 変更があった場合はmakeをすればOKです.

エラー落ちさせてみる

コードをいじってもいいですが, めんどくさいのでエラー落ちするようなケースを作ります.

cp ~/OpenFOAM/OpenFOAM-2.4.0/tutorials/incompressible/icoFoam/cavity ~/Desktop # 一応避難
cd ~/Desktop/cavity
流速を100m/sなどにして計算を発散させる

以下のようなメッセージが出る

#0  Foam::error::printStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-2.4.0/src/OSspecific/POSIX/printStack.C:219
#1  Foam::sigFpe::sigHandler(int) at ~/OpenFOAM/OpenFOAM-2.4.0/src/OSspecific/POSIX/signals/sigFpe.C:108
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::symGaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ~/OpenFOAM/OpenFOAM-2.4.0/src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:167
#4  Foam::symGaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ~/OpenFOAM/OpenFOAM-2.4.0/src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:237
#5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ~/OpenFOAM/OpenFOAM-2.4.0/src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.C:169良いかと思います
#6  ? at ~/OpenFOAM/OpenFOAM-2.4.0/src/finiteVolume/lnInclude/fvMatrixSolve.C:193 (discriminator 9)
#7  ? at ~/OpenFOAM/OpenFOAM-2.4.0/src/finiteVolume/lnInclude/fvMatrixSolve.C:82
#8  ? at ~/OpenFOAM/OpenFOAM-2.4.0/src/finiteVolume/lnInclude/fvMatrixSolve.C:331
#9  ? at ~/OpenFOAM/OpenFOAM-2.4.0/src/finiteVolume/lnInclude/fvMatrix.C:1390
#10  ? at ~/OpenFOAM/OpenFOAM-2.4.0/applications/solvers/incompressible/icoFoam/icoFoam.C:63 (discriminator 7)
#11  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12  ? at ??:?

sigFpeより上はエラーをお知らせするだけなので, #3でエラーおちしてそうです

gdbを使ってみる

macだとlldbになります. 私はlldbは使ったことがないのでよくわかりませんが, 多分同様にやればOKです.

gdbのコマンド一覧

説明がめんどくさいのでこちらの記事に任せます. gdb の使い方・デバッグ方法まとめ

エラー落ちの原因を探る

当然流速が早すぎるからですソースコードではどこで落ちているのでしょうか

エラーメッセージを見てみると, ~/OpenFOAM/OpenFOAM-2.4.0/src/OpenFOAM/matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.Cの167行目 Foam::symGaussSeidelSmoother::smooth(...) メソッド内で死んでいるっぽいです. ソースコードはこれみたいですね github.com

167行目は

psii /= diagPtr[celli];

とかかいてるのでゼロ除算とかで死んでるのかもしれないですね見てみましょう

gdbを実際に使う

通常は, icoFoamと打つところですが, 今回は少しコマンドを増やして

gdb icoFoam

と実行します.

(gdb)

みたいなのが出るので, runと打ち込んでやりましょう

(gdb) run

出力結果

Starting program: /home/kurenaif/OpenFOAM/OpenFOAM-2.4.0/platforms/linux64GccDPDebug/bin/icoFoam
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
(gdb) /*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.4.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 2.4.0-bf7c62d394a7
Exec   : /home/kurenaif/OpenFOAM/OpenFOAM-2.4.0/platforms/linux64GccDPDebug/bin/icoFoam
Date   : Jan 29 2017
Time   : 01:55:16
Host   : "kurenaif-server"
PID    : 22737
Case   : /home/kurenaif/Desktop/cavity
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading transportProperties

Reading field p

Reading field U

Reading/calculating face flux field phi


Starting time loop

Time = 0.005

Courant Number mean: 0 max: 0
smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 8.90511e-06, No Iterations 19
smoothSolver:  Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG:  Solving for p, Initial residual = 1, Final residual = 7.55423e-07, No Iterations 35
time step continuity errors : sum local = 5.03808e-07, global = -4.99749e-18, cumulative = -4.99749e-18
DICPCG:  Solving for p, Initial residual = 0.523588, Final residual = 9.72371e-07, No Iterations 34
time step continuity errors : sum local = 1.07766e-06, global = -2.2097e-17, cumulative = -2.70945e-17
ExecutionTime = 0.02 s  ClockTime = 0 s

Time = 0.01

Courant Number mean: 9.76803 max: 58.5723

Program received signal SIGFPE, Arithmetic exception.
0x00007ffff59e5d14 in Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:167
167             psii /= diagPtr[celli];

167行目でArithmetic exceptionですねゼロ除算だとかオーバーフローだとかででるエラーメッセージです. 例外をはいただけでまだプログラムは終了してませんこのまま 変数を覗いてみましょう. pコマンドで変数の中身を見ることができます.

(gdb) p psii

出力結果

$1 = 1.7457973770058812e+305

あ、やっぱり発散してますね除算する値はどうでしょうか

(gdb) p diagPtr[celli]
value has been optimized out

困りました。見れません “このエラーは通常コンパイラの最適化によって変数が無いよ” みたいなエラーです元をたどりましょう 85行目にこんなのがありました

 register const scalar* const __restrict__ diagPtr = matrix_.diag().begin();

どうやらdiagPtrはmatrix_.diag().begin()と同じ意味みたいです. 同じ意味なのでコンパイラが省略しちゃったんですねではdiagPtrを見る代わりに, matrix_.diag().begin()で確認してみましょう

(gdb) p matrix_.diag().begin()[celli]
$2 = 0.00055000001416059699

無事見れましたねどうやら0除算ではなかったみたいです. ただの値の発散ですね

次に値が発散していく様子を観察してみましょう いちどgdbを閉じて (Ctrl+d) 再度起動します

gdb icoFoam

とりあえず問題の symGaussSeidelSmoother.Cの167行目で止めてみましょう途中のやつはyでOKです

(gdb) b symGaussSeidelSmoother.C:167
No source file named symGaussSeidelSmoother.C.
Make breakpoint pending on future shared library load? (y or [n]) y
Breakpoint 1 (symGaussSeidelSmoother.C:167) pending.

これでrunをしてみます

Breakpoint 1, Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:167
167             psii /= diagPtr[celli];
(gdb)

止まりました. sやnを押すと1行次に進みたいと思います.

関数があったときに, 関数の中身に飛んでいくほうが s
関数があってもそれを1行とみなして中身をスルーするのがnです。

gdbは, 前のコマンドと同じ入力をするときに, 何も入力せずにEnterを押すと繰り返しができます.

nを押してみて, Enterを連打してみましょう

167              psii /= diagPtr[celli];
(gdb) n
170             for (register label facei=fStart; facei<fEnd; facei++)
(gdb)
172                 bPrimePtr[uPtr[facei]] -= lowerPtr[facei]*psii;
(gdb)
170             for (register label facei=fStart; facei<fEnd; facei++)
(gdb)
172                 bPrimePtr[uPtr[facei]] -= lowerPtr[facei]*psii;
(gdb)
170             for (register label facei=fStart; facei<fEnd; facei++)
(gdb)
175             psiPtr[celli] = psii;
(gdb)
151         for (register label celli=0; celli<nCells; celli++)
(gdb)
154             fStart = fEnd;
(gdb)
155             fEnd = ownStartPtr[celli + 1];

forループが観測できますね psii の発散する様子が見たいので, psiiをdisplayしてみましょう

display: 値を監視して変更等があったときに, その都度printしてくれる

(gdb) display psii
1: psii = 0

はい監視対象に入りました nコマンドではforに捕まった時に面倒です. breakpoint 1まで飛ばしてしまってよいでしょうそういうときはcを使います

(gdb) c
Continuing.

Breakpoint 1, Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:167
167             psii /= diagPtr[celli];
1: psii = 0

はい breakpointまで飛びました.　ついでに監視対象のpsiiも見れますね続けて見ていきましょう Enter連打してcを何度も回してみます

Breakpoint 1, Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:167
167             psii /= diagPtr[celli];
1: psii = 0
(gdb)
Continuing.

Breakpoint 1, Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:167
167             psii /= diagPtr[celli];
1: psii = 0
(gdb)
Continuing.

Breakpoint 1, Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:167
167             psii /= diagPtr[celli];
1: psii = 0
(gdb)
Continuing.

Breakpoint 1, Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:167
167             psii /= diagPtr[celli];
1: psii = 0
(gdb)
Continuing.

Breakpoint 1, Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:167
167             psii /= diagPtr[celli];
1: psii = 0

うーんなかなか発散してくれませんねこういう時はウォッチポイントを使います. psiiが1以上になった時に止まるようにしてみましょう breakpointで止まるようになるのが邪魔なので消します.

(gdb) delete
Delete all breakpoints? (y or n) y
(gdb) watch psii >= 1
Watchpoint 2: psii >= 1

消えました. ついでにwatchpointも増えました continueしてみましょう

(gdb) c
Continuing.

Watchpoint 2: psii >= 1

Old value = false
New value = true
Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:170
170             for (register label facei=fStart; facei<fEnd; facei++)
1: psii = 30.769230769230489

おっ psiiが30になりましたこれは発散しそうですでは, 次はpsiiが変化した時に止まるようにしてみましょう

(gdb) watch psii
Watchpoint 3: psii
(gdb) c
Continuing.

Watchpoint 3: psii

Old value = 30.769230769230489
New value = 0.023076923076922929
Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:161
161             for (register label facei=fStart; facei<fEnd; facei++)
1: psii = 0.023076923076922929
(gdb)
Continuing.

Watchpoint 3: psii

Old value = 0.023076923076922929
New value = 41.958041958041726
Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:170
170             for (register label facei=fStart; facei<fEnd; facei++)
1: psii = 41.958041958041726
(gdb)
Continuing.

### Watchpoint 3: psii

Old value = 41.958041958041726
New value = 0.024195804195803933
Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:161
161             for (register label facei=fStart; facei<fEnd; facei++)
1: psii = 0.024195804195803933
(gdb)
Continuing.

Watchpoint 3: psii

Old value = 0.024195804195803933
New value = 43.9923712650982
Foam::symGaussSeidelSmoother::smooth (fieldName_=..., psi=..., matrix_=..., source=..., interfaceBouCoeffs_=..., interfaces_=..., cmpt=0 '\000', nSweeps=1) at matrices/lduMatrix/smoothers/symGaussSeidel/symGaussSeidelSmoother.C:170
170             for (register label facei=fStart; facei<fEnd; facei++)
1: psii = 43.9923712650982

徐々に発散してますね今回はこれくらいにしておきましょう

backtrace書き忘れましたこんな感じで関数の呼び出し元を見てくれます

(gdb) bt
#0  0x00007ffff5ed496d in Foam::symGaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) ()
   from /home/kurenaif/OpenFOAM/OpenFOAM-2.4.0/platforms/linux64GccDPOpt/lib/libOpenFOAM.so
#1  0x00007ffff5ed4e88 in Foam::symGaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const ()
   from /home/kurenaif/OpenFOAM/OpenFOAM-2.4.0/platforms/linux64GccDPOpt/lib/libOpenFOAM.so
#2  0x00007ffff5ecc07a in Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const ()
   from /home/kurenaif/OpenFOAM/OpenFOAM-2.4.0/platforms/linux64GccDPOpt/lib/libOpenFOAM.so
#3  0x0000000000437d9f in Foam::fvMatrix<Foam::Vector<double> >::solveSegregated(Foam::dictionary const&) ()
#4  0x0000000000457bc9 in Foam::fvMatrix<Foam::Vector<double> >::solve(Foam::dictionary const&) ()
#5  0x0000000000457e24 in Foam::fvMatrix<Foam::Vector<double> >::solve() ()
#6  0x0000000000419e59 in main ()

まとめ

今回はエラーを起こしてそのエラーの発生箇所の特定, 及び原因の一部を見てみました. お役に立てれば幸いです

References

gdb の使い方・デバッグ方法まとめ

(pp.44~)

OpenFOAM －空間の離散化と係数行列の取り扱い（Spatial Discretization and Coefficient Matrix）－ from Fumiya Nozaki

Motivation

概要

gccを入れる

openmpiを入れる

flexを入れる

OpenFOAMを入れる

Motivation

ケースファイル

メッシュ作成

decomposePar(自動領域分割)

使用するソースコード

ソースコードの解説

argsの準備

runTimeの準備

meshの準備

各セルがどのprocessorに該当するのかを記述するための変数の用意

各セルをループで回しそれぞれのセルにprocessorIDを振り分ける

出力

References

追記

可視化対象

実行

可視化をするソースコード

とりあえず実行

output

解説

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

簡易的なソースコードの解説といろいろ付け足し

caseの読み込み

軸の非表示，背景色変更

AnnotateTime追加

出力画像サイズの変更

色の設定，range設定

ScalarBar(右についてるカラーバー)の表示，設定

カメラの座標の設定

時間ステップの変更，再描画

References

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

pvbatchで描画せずに画像だけ表示

インストール方法

使い方

デモ

.plt(gnuplotのファイル) -> .eps+.tex -> .pdfファイル

.pdfファイルの自動更新

pltを各際の注意

OpenFOAMでgdbを使いたい

Intro

Motivation

Target

環境

どんなことができるの

通常のエラーメッセージ

デバッグオプション

使ってみる

Debugオプションを入れた状態でのmake

DebugとOptを切り替える

エラー落ちさせてみる

gdbを使ってみる

gdbのコマンド一覧

エラー落ちの原因を探る

gdbを実際に使う

まとめ

References