› 製造解法 (MMS) コード検証シミュレーター 戻る
検証・妥当性確認 (V&V)

製造解法 (MMS) コード検証シミュレーター

CFD/FEM コードの離散化が正しく実装されているかを、「製造解法 (Method of Manufactured Solutions)」で検証するツールです。任意の解析解を「製造」してソース項を逆算し、4段階の格子で L2 誤差と観測収束次数を測ります。観測次数が設計次数と一致すれば、コード検証 PASS です。

パラメータ設定
領域長さ L
1次元計算領域 x ∈ [0, L] の長さ
粗グリッド分割数 N
最も粗い格子のセル数。以下、4段階に細分化
リファインメント比 r
隣接する格子間で h を 1/r に縮める比率
期待収束次数 p_design
スキームの理論次数（中央差分なら 2）
拡散係数 D
−D·u'' = f(x) の物性係数
製造解 u_mms(x)
本ツールでは u_mms(x) = sin(πx/L) を固定使用。f(x) = D·(π/L)²·sin(πx/L) を逆算してコードに与えます。
計算結果
—
粗格子誤差 E₁
—
細格子誤差 E₄
—
観測収束次数 p_obs（平均）
—
期待次数 p_design
—
MMS 検証結果
—
誤差比 E₁/E₄
—
格子細分化と数値解 vs 製造解

4段階の格子（粗→細）で u_mms(x) = sin(πx/L) を解いた数値解を重ね描き。誤差バーが格子細分化とともに縮んでいく様子を可視化します。

収束プロット — log E vs log h
格子間の観測収束次数 p_obs
理論・主要公式

$$f_{mms}(x) = L[u_{mms}(x)],\quad E_k = \|u_h^{(k)} - u_{mms}\|_2,\quad p_{obs} = \frac{\log(E_1/E_2)}{\log r}$$

L は支配方程式の差分演算子、u_mms は任意の滑らかな関数、p_obs ≈ p_design なら検証成功。

$$-D\,u''(x) = f(x),\qquad u_{mms}(x) = \sin\!\left(\tfrac{\pi x}{L}\right),\qquad f(x) = D\!\left(\tfrac{\pi}{L}\right)^{2}\sin\!\left(\tfrac{\pi x}{L}\right)$$

本ツールが解く 1次元定常拡散方程式と、製造解からソース項を逆算した f(x)。境界条件は u(0)=u(L)=0。

$$E_k \;\propto\; h_k^{\,p_{design}},\qquad \frac{E_k}{E_{k+1}} = r^{\,p_{obs}}$$

2次中央差分なら E ∝ h²、r=2 で隣接格子の誤差比が 2² = 4 になるのが理論値。本ツールはこの比から p_obs を逆算します。

製造解法 (Method of Manufactured Solutions, MMS) によるコード検証

🙋
「製造解法」って、なんだか不思議な名前ですね。普通は方程式を解いて解を「求める」のに、なんで「製造する」んですか？
🎓
いい着眼点だね。考え方が逆なんだ。普通の検証は「ベンチマーク問題の解析解と比較する」やり方だけど、複雑な方程式だと解析解が手に入る問題が限られる。MMS はその逆で、まず好きな滑らかな関数 u_mms(x) を勝手に決めてしまう。例えば u_mms = sin(πx/L) とかね。それを支配方程式 -D·u''(x) = f(x) に代入すると、左辺は計算できるから「この u を成り立たせるためにはソース項 f がいくらでなきゃいけないか」が逆算できる。今回なら f(x) = D·(π/L)²·sin(πx/L) になる。
🙋
なるほど、最初に答えを決めて、そこから問題を作るんですね。それでコードのバグが分かるんですか？
🎓
そう、その f をテストコードに入力して数値解 u_h を計算させる。本来なら u_h は u_mms と一致するはずだ。実際には離散化の打ち切り誤差で少しズレる。そのズレ E = ‖u_h − u_mms‖₂ を、格子をどんどん細かくしながら測っていく。2次中央差分なら理論上 E ∝ h² で減るから、h を半分にすれば誤差は 1/4 になるはず。本ツールの右上にある E₁/E₄ が 64（=2⁶）になっているのが、それだよ（h を3回半分にしたから 2² × 2² × 2² = 64）。
🙋
じゃあ「観測収束次数」が設計次数と一致したら、コードが正しいって言えるんですね？
🎓
かなり強いお墨付きになる。観測次数 p_obs = log(E_k/E_{k+1})/log(r) が p_design とピタリ一致するのは「離散化が想定通り」「境界条件も正しく組まれている」「ソース項の符号も合っている」が全部揃わないと起きないからね。実務で多いハマりどころは、コードを書き換えて「動くけど精度が落ちた」ケース。テストが既存ベンチマークだけだと「許容範囲内」で見逃すんだけど、MMS だと p_obs が 2 から 1 にハッキリ落ちるから、回帰バグを定量的に検出できる。Roy・Oberkampf や ASME V&V 10/40 が MMS を gold standard として推奨しているのは、この感度の高さが理由なんだ。
🙋
逆に p_obs が設計次数からズレたときは、何を疑えばいいですか？
🎓
パラメータの「期待収束次数 p_design」を 1 に変えてみて、収束プロットの傾きを観察してごらん。観測次数が想定より低くなる典型パターンは4つ。まず境界条件の実装ミスで、境界1点だけ1次精度になると全体が1次に引きずられる。次にソース項の符号ミス——f の符号を反転すると誤差が格子に依存せず一定値に張り付く。三つ目が反復ソルバの収束不足——線形方程式の残差が打ち切り誤差を上回ると p_obs が 0 に近くなる。最後がasymptotic 域に入っていない——粗格子側で誤差が頭打ちになる場合は、N をもっと大きく取り直す必要があるんだ。

よくある質問

通常の検証は「ベンチマーク問題の解析解と比較する」方式で、扱える問題が限られます。MMS は逆で、まず任意の滑らかな解析解 u_mms(x,t) を自由に「製造」し、それを支配方程式に代入してソース項 f を逆算します。そのソース項をコードに入力すれば、どんな複雑な方程式・境界条件でも厳密解付きのテストケースが作れます。これにより理論収束次数（設計次数 p）と観測次数 p_obs を比較でき、コード実装の正しさを定量的に証明できる、V&V の gold standard です。
格子 k の代表サイズ h_k で離散化したときの L2 誤差を E_k = sqrt(Σ(u_h − u_mms)² · h) で測ります。隣接する2格子について、リファインメント比 r で h を縮めたときの誤差比から p_obs = log(E_k / E_{k+1}) / log(r) を計算します。本ツールでは粗・中・細・最細の4段階で h を半分（r=2）に縮めて、3組の p_obs を出し、その平均を表示します。p_obs が設計次数 p（例: 2次差分なら 2.0）と±0.15 以内で一致すれば、コードに離散化次数を落とす実装バグはないと判断できます。
代表的な原因は4つです。(1) 境界条件の実装ミス — 境界での次数低下が全領域の誤差を支配し、p_obs が 1 次に落ちる。(2) ソース項の符号ミス — 多くは f の正負を間違えており、誤差が格子に依らず一定になる。(3) 反復解法の収束不足 — 残差が打ち切り誤差より大きいと、p_obs が見かけ上 0 に近くなる。(4) 格子がまだ asymptotic 域に入っていない — 粗格子側で誤差が頭打ちになる。本ツールで p_obs が設計次数と大きく外れた場合、まず境界条件と f の符号を再確認してください。
いいえ、非線形・非定常・多次元・連成問題のすべてで使えるのが MMS の最大の強みです。Navier-Stokes、移流拡散、構造非線形、流体構造連成、輻射伝熱 — いずれも u_mms を選び、その演算子 L[u_mms] をシンボリック計算（SymPy 等）で展開すれば f が得られます。製造解は物理的に意味がある必要はなく、コードの全項を活性化する滑らかな関数（多項式・三角関数・指数関数の組み合わせ）を選ぶのがコツです。本ツールは線形1次元拡散を題材に MMS の基本ロジックを可視化していますが、同じ手順がそのまま大規模 CFD コードの検証にも適用できます。

実世界での応用

商用 CFD・FEM ソルバの検証：ANSYS Fluent、OpenFOAM、CFD++、Code_Saturne、ABAQUS、LS-DYNA など主要ソルバの開発元は、リリース前に必ず MMS でコード検証を行っています。新しい乱流モデル、輸送方程式、境界条件タイプを追加するたびに、対応する u_mms と f を生成して観測次数をチェックする回帰テストが自動化されています。NASA の CFL3D・FUN3D、Sandia 国立研の SIERRA など、安全クリティカル分野で使われるコードは MMS テストスイートが品質保証の根幹です。

研究コードの査読・公開審査：近年は学術誌（JCP, CMAME, IJNMF など）で新しい数値手法を提案する論文に、MMS による収束次数の証明を求めるケースが増えています。「2次精度のスキームを提案した」と書くなら、観測次数が 2.00±0.05 程度に揃っているプロットを添えるのが事実上の要件です。査読者にとっても、解析解との単純比較より MMS のほうがバグ検出能力が高く、信頼性の高い検証手段とされます。

原子力・航空宇宙の安全評価：原子炉熱水力（CFD/CFD-NS）、ロケット燃焼室、再突入機の空力加熱など、解析結果が安全規制に直結する分野では、コード検証 (Code Verification) → 計算検証 (Solution Verification) → 妥当性確認 (Validation) の三段階が ASME V&V 10/20/40 や ANS-2.27 で標準化されています。コード検証フェーズの中核技術が MMS で、製造解で観測次数 PASS が取れていない実装は、その後の Validation 結果も信頼できないと見なされます。

社内開発スクリプトの品質ゲート：大企業のエンジニアリング部門では、社内 Python/MATLAB で書いた解析スクリプトを実プロジェクトに使う前の品質ゲートとして MMS を採用する例が増えています。少数の格子で観測次数を測るだけなので計算コストが小さく、CI/CD パイプラインに組み込みやすいのが理由です。本ツールのインタラクションは、その「最小限の MMS テストランナー」のミニチュアです。

よくある誤解と注意点

まず最大の誤解が、「MMS が PASS したらコードは正しい＝結果も信用できる」と思い込むこと。MMS が証明するのはあくまで「コード検証 (Code Verification)」——プログラムが意図した数学モデルを正しく解いている、という一点だけです。実際の物理現象との一致（Validation）は別途、実験データとの比較が必要です。例えば乱流の RANS 解析で MMS が PASS でも、モデル自体が実流れと合っていなければ無意味です。MMS は「方程式の解法が正しい」ことを保証するが、「方程式が現実に正しい」かは保証しないという区別を、常に意識してください。

次に、製造解を「物理的に妥当」にしようとして失敗すること。u_mms は計算しやすく滑らかであれば何でもよく、物理的に意味がある必要はまったくありません。むしろ単純な多項式や sin だけだと、コードの一部の項（移流項、非線形項、境界条件の特殊な分岐など）が活性化されず、バグを見逃す危険があります。実務では「方程式の全項にゼロでない寄与が出る」よう、多項式・三角関数・指数関数を組み合わせた人工的な u_mms を意図的に選びます。本ツールでは sin(πx/L) を採用しており、境界条件 u(0)=u(L)=0 を自動的に満たす最もシンプルな例です。

最後に、格子が asymptotic 域に入っていないのに観測次数を出してしまう失敗。p_obs が意味を持つのは、誤差が「主項の打ち切り誤差」に支配されている領域だけです。格子が粗すぎると高次項や non-asymptotic 効果が混じり、観測次数が設計次数からズレます。判別法は単純で、4段階の格子の p_obs がほぼ一定であれば asymptotic、ばらついていれば非 asymptotic です。本ツールで N（粗格子分割数）を 4 まで下げると、p_obs にわずかなばらつきが出ることがあります。実コードでも、粗側の格子だけで判断せず、必ず最低3比較（4格子）以上で p_obs の収束性も確認してください。

使い方ガイド

  1. 拡散方程式のパラメータ（熱伝導率λ、領域長さL）と解析解の多項式次数（1次～5次）を入力します
  2. 格子分割数を4段階設定し、各段階でMMSソースターム q(x) を製造して数値計算を実行します
  3. 粗格子（N₁）から細格子（N₄）までのL2ノルム誤差 E₁, E₂, E₃, E₄ を取得します
  4. 隣接格子間の誤差比から観測収束次数 p_obs = log₂(Eᵢ/Eᵢ₊₁) を3組計算し、平均値と設計次数を比較判定します

具体的な計算例

アルミニウム導熱棒（λ=237 W/m·K、L=0.5m）で2次多項式解析解 u(x)=x²+2x を製造します。格子N₁=8→16→32→64で数値計算時、誤差は E₁=3.2×10⁻³、E₂=7.8×10⁻⁴、E₃=1.9×10⁻⁴、E₄=4.8×10⁻⁵となり、p_obs=(3.04+3.03+3.02)÷3≈3.03で設計次数p_design=3に一致。観測次数が設計次数の±0.1内なら検証成功と判定します。

実務での注意点