くれなゐの雑記

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

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)
カメラの座標の設定

カメラの座標は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 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.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

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が表示されています. いいですね

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

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

miettal.hatenablog.com

今日から楽しくLet’s gdb

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

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

デモ

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



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を各際の注意

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

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

OpenFOAMでgdbを使いたい

Intro

Motivation

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

この記事では,

  1. 概要
  2. デバッグオプションをいれたOpenFOAMのmake
  3. 実際にエラーの原因を探してみる
  4. 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)

止まりました. snを押すと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 の発散する様子が見たいので, psiidisplayしてみましょう

  • 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~)

www.slideshare.net

文書作成ソフト使い方まとめ

Kobe University Advent Calendar 2016 - Adventar の 16日目を担当させていただきました。

さて, そろそろ卒論シーズンですねもうかきはじめている人もいるかと思いますが、皆様何を使っていますか?
文書作成ソフトですとWordが簡単で使いやすいみたいなのをよく聞きますが、本当に使いやすいですかね…? 結構上級者向けソフトだと思ってます。
TeXを使ったことがないみたいな人も知り合いにいましたので, Word, TeXとついでにmarkdownで使ってるTipsをもりもり描いていこうかなと思います。
エディタとかは色々使ってみて使いやすいのを選ぶといいと思います。いまだにしっくり来るのがないのでオススメを教えてください

TeXの使い方

概要

TeXコンパイルをしていい感じの見た目に自動的に揃えてくれます。個人的にはWordを使う前に、TeXを触ってからWordにいくといいと思ってます。
あとTeXLive使ってる人はアップデートしておきましょう。(外部コマンドの実行 - TeX Wiki)

エディタについて

TeXは色々なエディタで使用できます. 最初はVimを使っていたのですが、なんかプラグインがバグりまくってまともに動作しないのでIDE的なやつを入れることにしました。結果、以下の3つのエディタに行き着いた.
IDe的なやつを使うと表の挿入が簡単にできたりするので好き

本当はそれぞれのエディタを説明したいところだけど, 説明することがおおいのでざっくり特徴を抑えておく。

  • Sublimetext: デフォルトのスニペットが優秀, 普段良く使ってる人も多いので初期の学習コストも少ない人も多いハズ
  • TeXstudio: WindowsではしんどいTeXとエディタの連携がこれを使うと簡単にできる。 内蔵のビューワーと補完が優秀 Vimキーバインドが使えなくてこれは使ってなかったけどそれさえ使えればって感じ
  • kile: 一応クロスプラットフォームでも動くと書いているけどLinux以外ではいい感じではない. 最初Tab幅とか補完とかの設定が少し面倒 使い勝手はTeXstudioと同じくらい。 最初は設定を上からみてポチポチボタンを押していこう あとひとつの数式だけコンパイルしてプレビューみたいなのもできる 便利

バージョン管理

文章を推敲していった時、他人に何処を変更したか、あるいは変更してもらったか等を見る、また複数人で一つの文章を書く と言ったときにgitなどを使ってバージョン管理が出来ると非常に便利。
適当に章でbranchを切って複数人で作業をしたりできる. Githubのプライベートリポジトリとかに上げておけばどこでも編集できるし、先生に変更を通知できたりして便利。

一つの文で改行

特に習ったわけではないが、こうしている。
TeXでは、1回の改行は特に意味はない。
1文ずつは短いほうが良いので, こうすることで長い文を簡単に見つけ出すことができる.

大きな文章(卒論等)を書く場合、それぞれ章を分けて別々にコンパイルできるようにする

なっがい文章を描く時、いちいちコンパイルが長くてしんどくなってくる。章ごとにコンパイルして確認するとコンパイル時間が短くなって便利.
これを見たらできる
qiita.com

Wordの使い方

概要

WYSIWYG(What You Get Is What You See)なエディタ(?)で非常に直感的なUIをしていると思います。
それ故に機能をほぼほぼ使わずにある程度形にできてしまい、Wordを使った気分になった人も多い印象です。(酷いケースだとスペース(' ')で無理やり中央寄せをしていたケースを見た)
TeXと比較しながらよく聞かれる機能を紹介する。

見出し機能(Section相当)

[ctrl+alt+1,2,3]でsection, subsection, subsubsection相当の挙動ができる.
このように追加した見出しは, wordのナビゲーションウィンドウから一覧を表示することができて、後の検索にも便利。
もちろん, 前に章番号を置いたり、それぞれフォントを変えることもできる。
積極的に使っていこう
またこの機能を使うと自動目次挿入機能でしっかり目次を作成してくれる。

図番号, 表番号の相互参照(\label, \ref相当)

図番号を描くときに 図xxをyyyして… のxxの部分は手打ちしていないだろうか? 前共同で描いていた人がこれを手打ちで描いていて反映が非常に苦だった。
相互参照という機能がある。図を追加し、その図を右クリックして[図表番号の挿入]を選択し、図を挿入。
[alt+n r f]と入力すると、相互参照のメニューが出てきて適当にガチャガチャやると挿入ができる。
これをすることで、参照先の図番号が変わったりしても反映してくれる。
しかし全自動で反映してくれるわけではなく、全選択→右クリック→フィールドの更新をする必要がある。
ちなみに参照先がなければエラー! 参照元が見つかりません。と出る。

参考文献リスト(bibtex相当)

[alt s m] を入力すると資料文献の管理ダイアログが出てくる。これを使用すると、参考文献のリストを一括管理することができ、 [alt s c] で登録したリストを参照して\cite相当の挙動をすることができる。
スタイルも変更可能で、必要があれば作ることも出来る。
最後にリストをまとめてReferencesとしてのっけたいときも文献目録機能で乗っけることが出来るので使っていこう。

図を並べたい

Wordで図を横に並べる時は、表機能を使うと便利。
表を並べて、その中でうまいこと調整してあげよう。

数式番号を挿入したい

この記事を参考にしよう
toomva.blog60.fc2.com

数式のフォントを変えたい

note.chiebukuro.yahoo.co.jp

バージョン管理

変更履歴機能を使おう 印刷したときも赤字で変更箇所を明示してくれて非常に便利 推敲していく時に使っていこう
support.office.com

wordの数式で等号揃えをしたい

kenkitagawa.cocolog-nifty.com

この文字を行末に持ってきてほしくない

禁則文字機能を使いましょう

、。→, . にしたい

オートコレクト機能を使えばよい。
michisugara.jp

markdown

概要

ささっと描くときとgithubのreadmeを書くときに使う

エディタ

  • Sublimetext, Atom, Marp(プレゼン用)

.pdfにしたい

こちらの記事を参考にする
kurenaif.hatenablog.com
kurenaif.hatenablog.com

.htmlにしたい

sublimetextもAtomもhtmlにする機能があるのでそれを使えば良い。

あとがき

全体的にババババッってまとめる感じにしたかった
それぞれ気になる機能があればおそらく自分でググればすぐ出てくるので調べていただければと思います。

ubuntuにsourcecodeProを入れる

タイトル通り こちらの記事をコピペするだけで行けました
ubuntu 16.04

askubuntu.com

#!/bin/bash
mkdir /tmp/adodefont
cd /tmp/adodefont
wget https://github.com/adobe-fonts/source-code-pro/archive/2.010R-ro/1.030R-it.zip
unzip 1.030R-it.zip
mkdir -p ~/.fonts
cp source-code-pro-2.010R-ro-1.030R-it/OTF/*.otf ~/.fonts/
fc-cache -f -v

CODE FESTIVAL 2016 qual A. D - マス目と整数 / Grid and Integers

解説に乗ってるけどちょっとわかりにくかったので自分の言葉に置き換えて整理します。 あとイラレの練習(
http://code-festival-2016-quala.contest.atcoder.jp/data/other/code-festival-2016-quala/editorial.pdf
)

要約

行列A = [a_{i,j}], (i \in [1,R], j \in [1,C]) が与えられる.
N個の要素には, すでに値が記入されている.
すべての値は, 以下の条件を満たさなければならない.

条件1. すべての範囲内のi,jに対して,  a_{i,j}+a_{i+1,j+1} = a_{i+1,j} + a_{i,j+1}
条件2. a_{i,j} >= 0

問題の変形

この問題では, 変形が重要です.

条件1の変形

条件1を変形すると,
 a_{i,j}-a_{i,j+1} = a_{i+1,j} - a_{i+1,j+1}
となり, 左右同士の差は上下で等しいことがわかります。図の矢印方向の差が等しいということですね
f:id:kurenaif:20160925221650j:plain
上下で等しい 上下で等しい… と下へ下へつなげていくと, 結局同じ列どうしの比較だと, 行にかかわらず差が等しいことがわかります。
f:id:kurenaif:20160925222232j:plain

さすがに, 列が違うと差は変わってきます.(同じだと思ってた時期が私にもありました…)
ここで, それぞれの行の一番左の値(a_{1,1} ,\cdots,a_{R,1}) を 簡単のためにx_1, \cdots, x_R とおき,
a_{i,0}-a_{i,1}, a_{i,1}-a_{i,2}, \cdots, a_{i,C-1}-a_{i,C}y_1, \cdots, y_C とおきます.
f:id:kurenaif:20160925223126j:plain

すると, a_{i,j} = x_i + y_1 + y_2 + \cdots + y_C と表記できることがわかります.

yy_j = y_1 + y_2 + \cdots + y_j = \sum_{i=1}^{j} y_i と置き換えると,

a_{i,j} = x_i + y_j と簡単にすることができます。(画像つくり忘れたんですけど, 矢印が左端からそのマスまでたどり着いてるイメージですね)

詳しい証明は解説参照で, ざっくりいうと, 条件1が成立するためには,
すべてのa_{i,j} = x_i + y_j を満たすような, (x_1,\cdots,x_R), (y_1,\cdots,y_C) が存在しなければならない.
ということになります.

条件2の変形

条件2は, 条件1と組み合わせて,
すべての i,jで,
a_{i,j} = x_i + y_j >= 0
が成立する必要があります. 全通り試すともちろん時間が足りないので, 工夫が必要です.
実は, x_i, y_j のとり方はuniqueではありません. 例えば, x_iを+1したら, y_jを-1したら良いです.
一番小さいところからスタートして(x), 一番小さい値を足す(負の数)と, a_{i,j} の最小値が求まり, これが0以上になればOKっぽそうです

実装方針

a_{i,j} = x_i + y_j を求める.

グラフを使って求めます. x_iy_iを頂点と捉えます.
その間を入力r,c,aを使って
x_ry_c をつなぐ重さaの辺を貼ります.
この時辺と頂点の関係を, x_r + y_c = a と定義します.

二変数なので, なにか一つ値を決めないといけません. 一つ値を適当に0(なんでもいいです)と決めると, あとは随時決まっていくのがわかると思います.
(例えば, x_0 + y_1 = 3, x_0+y_2 = 4 という風に辺が張られていたら, x_0=0と定めると, y_1, y_2は自ずと決まってきます)

こんな感じでDFSで求めていきます.

条件2

DFSと一緒に, xの最小値とyの最小値を更新しながら探索していきます.
xの最小値+yの最小値 < 0 ならば, "No" となります.

実装する上での注意点

xの最小値とyの最小値の初期値に注意しなければなりません
処理をするときに, x_i + y_j のように描いていると, オーバーフローするのでこのような実装をするならば, LLONG_MAX/2-1程度を上限にすると良いでしょう.

SourceCode

const ll INF = LLONG_MAX / 2 - 1;


struct Edge {
	ll t;
	ll c;
	Edge(ll to, ll cost):t(to),c(cost){}
	Edge(){}
};

vector<ll> weigh;
ll R, C;

bool DFS(ll node, vector<vector<Edge> >& edge, ll& xmin, ll& ymin) {
	if (node < R) xmin = min(xmin, weigh[node]);
	else ymin = min(ymin, weigh[node]);
	for (auto &to : edge[node]) {
		if (weigh[to.t] == INF) {
			weigh[to.t] = to.c - weigh[node];
			if (not DFS(to.t, edge, xmin, ymin)) return false;
		}
		else {
			if (weigh[to.t] + weigh[node] != to.c) return false;
		}
	}
	return true;
}

void solve() {
	cin >> R >> C;
	ll N = in();
	weigh.resize(R + C, INF);
	vector<vector<Edge> > edge(R+C);
	REP(i, N) {
		ll r, c, a; cin >> r >> c >> a;
		--r; --c;
		edge[r].emplace_back(R+c,a);
		edge[R+c].emplace_back(r, a);
	}

	bool res = true;
	ll xmin, ymin;
	REP(i, R + C) if(weigh[i] == INF and edge[i].size() > 0) {
		xmin = ymin = numeric_limits<int>::max();
		weigh[i] = 0;
		if (not DFS(i,edge,xmin,ymin) or xmin+ymin < 0) {
			cout << "No" << endl;
			return;
		}
	}
	cout << "Yes" << endl;