バグの原因は・・・思わぬところに。

　プログラミングの世界において避けて通れないのがバグ。そう、あのお目々がクリクリッとして、小さなブルドックを髣髴とさせる可愛らしいワンチャン・・・って、それはパグ犬だよーーー！！・・・というギャグは置いといて、バグと言うのは検出するのも難しく、原因を特定するのも非常に難しく、そしていざ判明してみると、たった１文字のスペルミスだったり、ほんの些細な事が諸悪の根源だった・・・という事もザラだったりします。

　前回のブログでお話したように、円筒座標系の流体解析にチャレンジしているのですが、参考にできるＨＰも少なく、検索エンジンで引っかかる研究論文を見ても、基礎方程式がテンソル表示で記述されているだけで、こちらが本当に知りたい事（数値解法で用いる途中の方程式など）については日本語でちょろちょろっと書いてある程度で・・・。本当にこの●●●●方程式の形であってるんだろうか・・・と不安になって３～４回導出を試みて同じ結果になる事を確認して(←これはもちろん手計算！)、コーディングに入ったわけです。

　しかし・・・案の定と言うべきかウンともスンとも言わず・・・撃沈。確かに、はじめの内は式の形がおかしなところもありました。座標変換に伴って余分な項が新たに発生しますが、その取り扱いが抜けていた所がありました。でも、この類のミスは理論とコードを付き合わせれば早期に解決できます。これで、コーディングされた方程式と手計算の理論は完全に一致しました。

　しかし・・・何度やってもウンともスンとも言わず・・・。今回テストケースに選んだのはサーマルグラディエントを与えた円筒容器内の自然対流なので、要は温度の差が重力(浮力)に反映されていないのは一目瞭然でした。でも・・・この項の記述には間違いはありませんでした。まさか「●●●●方程式の形」が、間違っている！？・・・

　疑うべき所はこの方程式のルーチンに絞られました。この方程式に組み込まれる全てのパラメータの値を逐次表示させるように修正して再び走らせました。

　なんじゃこりゃーー！！。画面には目を疑う光景が・・・本来物性値であるはずプラントル数Ｐｒが開始直後は当初から設定されている値を示していましたが、突然ゼロに切り替わりました。・・・そんなカバな！いやバカな！・・・プラントル数はずっと一定のはず（おいおい、これって世に言う怪奇現象か）。冗談はさておき、温度の差が重力(浮力)に反映されるためにはプラントル数は当初の設定値のままである必要があります。それがゼロになってしまったので、温度差が重力(浮力)に反映される事は無く、流れはいつまでたっても駆動される事はない・・・従って、ウンともスンとも動かない。ようやく繋がった。

　このプラントル数が勝手に書き換わる、という事はどこかで上書きする記述があるという事ですね。プラントル数の名前(ＰＲ)をプログラムの一番上からガンガン検索していくと・・・ありました！・・・何と何と半径方向の圧力傾度の変数名をＰＲにしてました！。初期の段階では半径方向の圧力傾度は特に与えていない(初期化された状態)ので、確かにゼロになります。

　結論としては同じ変数名をプラントル数と圧力傾度の２つで使っていたのです。これを修正したら・・・アレヨアレヨと動いちゃったではありませんか！？しかも、計算結果も間違ってはなさそう。気になる●●●●方程式の形も問題なかったようです。

　と言う訳で、変数名には気をつけましょう。