炉物理の研究

第 54 号 (2002 年 9 月)

高速・高精度モンテカルロ計算の実現をめざして―

MVP の開発－
日本原子力研究所
森貴正

mori@mike.tokai.jaeri.go.jp
1.

開発の経緯

私とモンテカルロ法との係わりは修士課程の学生時代まで遡る。当時、私は、断面積データの積分的評価を目的として、京大炉の電子線ライナックを用いて原子炉材料中の keV領域の中性子スペクトル測定を行っているグループに属していた。この実験は、直径 60cmあまりの体系の中央で光中性子を発生させ、体系中の設けたリエントラントホールと呼ぶ孔の底から中性子を取り出し、time-of-flight 法によって体系中の角度中性子束スペクトルを測定するもので、その特長として 1 次元球体系モデルで解析できることを謳っていた。
リエントラントホールは 1 次元球体系モデルでは考慮できないので、その影響を評価するというのが私の課題であった。そこでモンテカルロコードを作成し、解析を試みたのであるが、当時の計算機は今のパソコンと比較しても能力が劣る上に、私のモンテカルロ法に対する知識が乏しかったこともあり、十分な評価を行うことができず、結局実験的にその影響が小さいことを確認するのみとなってしまった。そのため、モンテカルロ法に対する最初の印象は決して良いものではなかった。
ところが、原研に入所した私に与えられたテーマは「核融合炉物理解析法の研究」、「FNSを用いた日米協力核融合炉ブランケット工学実験の解析」であり、複雑な実験体系を取り扱うために再びモンテカルロ法と付き合うこととなった。ここでは、日本側の主力解析コードとして多群二重微分型断面積とモンテカルロコード MORSE-DD を開発し、それなりの貢献ができたと思っている。しかし、その実験解析では 1 ケース 5 時間以上の計算時間を要したことから、モンテカルロ法を気楽に用いるにはその高速化が不可欠と感じた。また、米国側が協力実験の後期から使用し始めた MCNP コードの結果や原研が ANL より入手した VIM コードの使用経験などから、我が国独自の連続エネルギー法に基づく高速・高精度モンテカルロコードの戦略的重要性を認識し、1980 年代末頃からその開発に着手した。開発スタッフは、中川正幸さん（現原子力発電㈱）、佐々木誠さん（日本総研㈱、当時原研外来研究員）と私の 3 人であった。
1980 年代後半当時は、FACOM－VPP シリーズなどのベクトルスーパーコンピュータが出現し、原研計算センターを中心に、既存の原子力コードのベクトル化が盛んに進められていた。モンテカルロコードも例外ではなく、KENO、MORSE などのベクトル化が試みられていたが、その高速化は 2 倍程度であり十分ではなかった。一方、ミシガン大学のF. Brown などからは、全く新たな計算アルゴリズムを採用することによってモンテカルロ法でも十分高い高速化が可能であることが報告されていた。こうした状況から、我々3 人

39

炉物理の研究第 54 号 (2002 年 9 月)

は、ベクトルスーパーコンピュータの有効利用を想定した連続エネルギー法に基づく汎用モンテカルロコードを全く新たに開発するとの方針を固めた。
2.

ベクトル化モンテカルロコード MVP の開発

ベクトル化モンテカルロ計算が、従来の個々の中性子について生成から消滅するまでを順次追跡する history-base のアルゴリズムのかわりに、多数の中性子を同時に追跡する
event-base のアルゴリズムによって可能となることは広く知られていた。しかし、モンテカルロ計算全体をどのような計算タスクに分割しベクトル計算を行うか、中性子の属性（位置、飛行方向、エネルギー等）、幾何形状データ、核データなどの多量のデータをどのように取り扱うかには任意性があり、効率的な方法の開発が必要であった。また、これらのデータを用いて行うベクトル計算に適した各種サンプリング法が不可欠であった。
我々は、ベクトル計算の特性把握からスタートし、種々の試行錯誤の後に、広範囲の問題に対して平均 15 倍程度の高速化を達成できる方法として、「スタック駆動ゾーン選択方式」と呼ぶアルゴリズムに到達した。
「ゾーン選択方式」とは、飛行解析において、この処理を待っている全中性子を対象とはしないで、待っている中性子の最も多い空間領域（ゾーン）を選んで、そのゾーンに存在する中性子に対してのみ処理を行う方法である。直感
的には、全ての中性子を対象とした方が同時に処理する粒子数が多く（ベクトル長が長く）なり効率的と思われるが、この方法では、必要となる前処理やそれに伴うメモリーアクセスの競合のために、どうしても高速計算が実現できなかった。我々のアルゴリズムの詳細はこれまでに学会誌欧文誌などに発表した論文を参照されたい。ここでは、ベクトル計算のために考案した２つの方法を紹介する。
(イ) MVP における衝突核種の決定法

k 番目の核種と衝突する確率は Pk = ρ k σ tk / Σ t であり、通常のモンテカルロコードでは、k に対する積算確率を求め、それを用いて衝突する核種を決定している。この方法は、ベクトル計算で行うことも可能であるが、多数の中性子を取り扱うために粒子毎の積算確率が同時に必要となり、膨大なメモリーが必要となる。一方、飛行解析に必要となる巨視的断面積の計算では、その過程で、対象とする全ての中性子に対して、kΣ tk = ∑ ρ j σ t jj =1が計算されている。また、確率 Pｋは、次式のように書き直すことができる。

Pk =ρ kσ tkΣt=ρ kσ tk Σtk Σtk +1×L×k +2Σtk Σtk +1 Σt

ΣtK −1ΣtK

この式は、巨視的断面積の計算が k 番目の核種まで進んだ時点で、確率 ρ σ t / Σ t で k 番kkk目の核種を衝突核種の候補とし（既に候補があればそれに置き換える）、全核種 K までチェックが終了した時点で候補とされている核種が衝突核種として選ばれることを意味する。

40

第 54 号 (2002 年 9 月)

炉物理の研究この方法だと、同時に必要となるのは、Σ t と ρ σ t のみであり、必要となるメモリーは大kkk幅に低減される。さらに、この方法はベクトル計算効率も高いことが分かった。
（ロ）六方格子体系内の格子セル位置の決定法
幾何形状表現法の工夫のjm3

例として、六方格子体系内の格子セル位置の決定法を32紹介する。VIM コードなど

1

では、セル位置を右図に示した 2 つのインデックス iと j で識別しており、六角柱の各側面を区別していない。そのために、(i, j)の決

m24m10l =m=3i21定には複雑な IF 文による0012345= k判定を必要とし、ベクトル計算には適さなかった。そこで新たに、それぞれの面に垂直な 3 つの方向に、k、l、m の3 つのインデックスを導入した。これによって、六角柱内部の三角柱、すなわち六角柱の各側面を区別することとなり、IF 文による判定が不要となり、単純でベクトル計算に適するものとなった。
その他にも、核データ、幾何形状の表現法などに種々の工夫をした。中には、高速計算を重視して採用した等確率ビンによる非弾性散乱中性子のエネルギー分布表現など、後になって計算精度の観点から表現法を改めたものもある。また、断面積ライブラリー作成の簡便化のために、任意温度で計算できる機能を追加した。さらに、多くのユーザーの使用計算機環境がワークステーションやパソコンに移っていることを考慮し、ベクトル計算のために必要であった多大なメモリー使用量の低減や、ワークステーション・クラスターなどを利用した並列計算機能などを追加し、現在に至っている。
3.

おわりに原研では、現在も MVP コードの機能拡張を続けている。このようにして開発してきたMVP コードが、炉物理解析分野のかなり広い範囲で使用されていることには感慨深いものがある。しかし、モンテカルロ法による固有値（keff）計算には、核分裂源の収束加速法、その収束判定法、微小反応度評価などのための摂動計算法など多くの課題が残されている。また、国内でも、MCNP などの外国産コードを使用する人も多い。若い研究者の方がモンテカルロ法に興味を持って、手法の高度化とともに我が国独自のコードである MVPの発展に寄与して下さることを期待する。
MVP コードの開発に当っては、学生時代からご指導頂いた中川正幸さんと MORSE-DDのころから一緒にコード開発を進めた佐々木誠さんに負うところが大きい。ここに感謝いたします。
41