MMS: 2次元定常熱伝導
考え方:答えを先に決める
概要
自作の熱伝導ソルバーが正しいか確かめたいのですが、複雑な形状や境界条件では解析解がありません。どうすればいいですか?
答えを自分で作ってしまえばいい。それが製造解の方法(MMS)だ。たとえば T = cos(πx/3)·sin(πy) のような、なめらかで好きな関数を『答え』と決める。それを熱伝導の方程式に代入すると、その答えになるために必要な熱源 q と境界条件が逆算できる。その熱源と境界条件をソルバーに与えて解かせ、出てきた温度を決めておいた答えと比べる。格子を細かくしたとき、誤差が理論どおりの速さで減っていけば、コードは正しく作られていると言える。
方程式と製造解
熱伝導率 $k = 1$ の単位正方形。$x = 0$、$y = 0$、$y = 1$ は製造解の値を温度として与え、$x = 1$ には製造解から求めた熱流(ノイマン条件)を与える。$a$ は製造解の形を変えるパラメータ。
収束次数
$e_h$ は格子幅 $h$ のときの誤差。中心差分の5点差分法は2次精度なので、$p$ は2に近づくはず。境界条件や熱源の実装に誤りがあると、$p$ が理論値より小さくなる。
境界条件の2つの書き方
左は境界の外に仮想の点を置く2次精度の書き方、右は片側の差分による1次精度の書き方。右を『バグ』として入れて、MMSで見つけられるかを試す。
『答えから問題を作る』発想
製造解の考え方は、1980年代から数値計算の研究者のあいだで使われていたが、コードの検証の方法として体系的にまとめたのはRoache(ローチ)らで、1998年の本や2000年代初めの論文で広まった。『解析解がないから検証できない』という壁を、逆向きの発想で取り払った点が画期的だった。今では多くの商用ソフトや研究用のコードが、MMSによる検証の結果を公開しており、ASME V&V 10や20でもコード検証の中心的な方法として位置づけられている。
計算例
例1:正しい実装の収束
| 分割数 | 最大誤差(a = 2) | 次数 | 最大誤差(a = 3) | 次数 |
|---|---|---|---|---|
| 8 | 5.93×10⁻³ | — | 8.27×10⁻³ | — |
| 16 | 1.49×10⁻³ | 1.99 | 2.07×10⁻³ | 2.00 |
| 32 | 3.74×10⁻⁴ | 2.00 | 5.16×10⁻⁴ | 2.00 |
| 64 | 9.34×10⁻⁵ | 2.00 | 1.29×10⁻⁴ | 2.00 |
どちらの製造解でも、格子を半分にするたびに誤差は1/4になり、理論どおりの2次収束が確認できる。
例2:境界条件のバグ
| 分割数 | a = 2(境界で曲率0) | 次数 | a = 3(境界で曲率あり) | 次数 |
|---|---|---|---|---|
| 8 | 5.97×10⁻³ | — | 5.77×10⁻³ | — |
| 16 | 1.50×10⁻³ | 2.00 | 4.13×10⁻³ | 0.48 |
| 32 | 3.74×10⁻⁴ | 2.00 | 2.39×10⁻³ | 0.79 |
| 64 | 9.35×10⁻⁵ | 2.00 | 1.28×10⁻³ | 0.90 |
a = 2だと、バグがあるのに2次で収束して見えるんですね!
そこがMMSの落とし穴だ。片側差分の誤差は、境界での温度の2階微分(曲率)に比例する。cos(πx/2)は x = 1 でちょうど曲率が0になるから、バグの誤差も0になって隠れてしまう。a = 3にすると境界で曲率が残るので、次数が0.48→0.79→0.90と1に近づき、1次精度のバグがはっきり見える。製造解は、すべての項が0でなく効くように、しかも境界で特別な値にならないように選ぶ。三角関数に指数関数や多項式を混ぜた、少し『いびつな』関数を使うのが定石なんだ。
MMSの進め方
手順
- なめらかで、すべての項(拡散、移流、非線形項、境界の微分)が0にならない製造解を選ぶ。
- 記号計算ソフトで熱源と境界条件を導く。手計算の誤りを避ける。
- 3〜4段階の格子で解き、誤差と次数を求める。
- 次数が理論値に近づくかを確かめる。合わなければ、項や境界条件を1つずつ切り分けて原因を探す。
「誤差は小さいのに次数が合わない」
MMSで誤差が0.1%以下と十分小さいので『合格』にしたが、あとで調べると次数は1だった——という例がある。誤差の大きさは格子の細かさでいくらでも小さくできるので、それだけでは正しさの証拠にならない。1次精度のバグがあっても、格子を細かくすれば誤差は小さくなる。ただし理論より遅く。実際の問題では、そうした格子の細かさを確保できないことが多く、バグが結果に効いてくる。MMSで見るべきは誤差の大きさではなく、格子を細かくしたときの誤差の減り方、つまり次数なんだ。
よくある間違い
間違いと対策
| 間違い | 起きること | 対策 |
|---|---|---|
| 単純すぎる製造解 | バグが隠れる | すべての項が効く関数 |
| 誤差の大きさで合否判定 | 低次のバグを見逃す | 次数で判定 |
| 熱源を手で導く | 導出の誤りをバグと誤認 | 記号計算で導く |
| 粗い格子だけで次数を見る | 漸近域に入っていない | 3〜4段階で傾向を見る |
| 境界条件を検証しない | 境界のバグが残る | 境界の種類ごとに試す |
MMSや収束次数についてもっと知りたいです。
全体の考え方はMMSの概要、次数の評価は収束次数の評価、熱源の導き方は製造解の熱源項、コード検証全体の位置づけはコード検証、実問題での数値誤差の見積もりはRoacheのGCIの記事にある。リチャードソン外挿ツールで次数の計算、2次元熱伝導で2次元の熱伝導も試せるよ。
関連トピック
なった
詳しく
報告