計算効率を上げるためのvariance reduction（粒子のウェイト）の利用

Presentation on theme: "計算効率を上げるためのvariance reduction（粒子のウェイト）の利用"— Presentation transcript:

1 計算効率を上げるためのvariance reduction（粒子のウェイト）の利用
PHITS Multi-Purpose Particle and Heavy Ion Transport code System 計算効率を上げるためのvariance reduction（粒子のウェイト）の利用 PHITS講習会　入門実習 2016年5月改訂 title

2 実習内容 １．はじめに ２．中性子深層透過計算 [importance]
　Geometry splitting and Russian Roulette (セルインポータンス法） [weight window] 　Weight windows (ウェイト・ウインドウ法) ３．薄膜からの生成粒子計算 [forced collisions] 　強制衝突法

3 中性子深層透過計算 厚い遮蔽体の中性子輸送を計算し、深さ毎の実効線量率分布を計算 まずはimp.inpを走らせて，計算時間を計ってみよう
[ P a r a m e t e r s ] icntl = 0 maxcas = maxbch = 1 ... [ T - T r a c k ] (3つ目) mesh = xyz … z-txt = (uSv/h)/(source/sec) file = imp-dose-xz.out part = neutron epsout = 1 multiplier = all emax = mat mset1 all ( ) ヒストリ数は1,500 計算時間はphits.outの最後の方（total CPU time）に出力される 半径50cm、深さ180cmのコンクリート円柱 空気 コンクリート 14MeV 中性子 実効線量は，[t-track]で計算した中性子フルエンスに内蔵線量換算係数を乗じて導出 詳細はマニュアル（4.25[multiplier]＆6.1[t-track]）もしくはlecture\exercise1を参照

4 中性子深層透過計算 シングルコア（1.50GHｚ AMD)の計算 効率良く計算したい。 分散低減法 ウェイトの概念を利用
ヒストリー数 1,500 imp-dose-xz.eps total cpu time = sec imp-dose-xz.eps ヒストリー数 270,000 total cpu time = sec 効率良く計算したい。 分散低減法　ウェイトの概念を利用

5 モンテカルロ計算におけるWeightの概念
例：track length タリー ウェイト：　モンテカルロ計算における粒子の重要度、通常の計算は１*。 ウェイトを操作して、統計上発生しにくいイベントを人為的に発生させる。 ただし，ウェイトを操作した場合，ヒストリー毎の頻度分布（[t-deposit], output = depositなど）の計算ができなくなるので注意が必要。 *ただし、核データを使用した中性子輸送計算では、ウェイトが変化する

6 実習内容 １．はじめに ２．中性子深層透過計算 [importance]
　Geometry splitting and Russian Roulette (セルインポータンス法） [weight window] 　Weight windows (ウェイト・ウインドウ法) ３．薄膜からの生成粒子計算 [forced collision] 　強制衝突法

7 中性子がセル1とセル2を横切る時、I2/I1を設定
セルインポータンス法 領域(セル)にインポータンスIを設定 中性子がセル1とセル2を横切る時、I2/I1を設定 I2/I1 > 1: Particle Split （粒子の分割） 例 I2/I1 が整数の場合 　I2/I1 = 3: ３つの粒子にsplit(分割) I2/I1 が非整数の場合 I1=1 I2=3 I2/I1 = 2.75: 75%の確率で３つの粒子に分割 25%の確率で２つの粒子に分割 W=1/3 W=1/3 W=1/3 W=1 I1=1 I2=1/3 W=3 I2/I1 < 1: Russian Roulette （消滅か生存) 1/3の確率 W=1 消滅 2/3の確率

8 1ヒストリーのときの中性子飛跡を確認 imp-hist.inpを走らせて，中性子の飛跡を確認してみよう 空気 コンクリート （半径20cm，
5MeV 中性子 I10=1.0 I1=1.0 I2=1.0 標準では、インポータンスは全て１ imp-trajectory.eps

9 インポータンスを大きくしたときの中性子挙動を確認
[importance] part = neutron reg imp 10 1 1 2 2 4 imp-hist.inp ここで反応が起きた I10=1 I1=2 I2=4 この時点で既に２つに分かれている imp-trajectory.eps

10 インポータンスを小さくしたときの中性子挙動を確認
[importance] part = neutron reg imp 10 1 1 1/2 2 1/4 imp-hist.inp 1/2の確率で生存 1/2の確率で消滅 I10=1 I1=1/2 I2=1/4 imp-trajectory.eps

11 インポータンスを極端に大きくしたときの中性子挙動を確認
[importance] part = neutron reg imp 10 1 1 2 粒子が分割しすぎて計算時間が無駄になってしまう （統計はあまり良くならない） imp-hist.inp I10=1 I1=2 I2=100 imp-trajectory.eps セル間のインポータンスの比は最大２～３程度が良い “A Sample Problem for Variance Reduction in MCNP” LA MS DE

12 インポータンスを使用した中性子深層透過計算例
Ii+1/Ii　= 2.5 と設定 imp.inpの [importance]セクションを有効にして、実行してみよう [importance] part = neutron reg imp **0 **1 **2 **3 **4 **5 **6 **7 **8 **9 **10 **11 半径50cm、深さ180cmのコンクリート円柱 セルの厚さは15cm 半径1cmの14MeV中性子をz=0からz軸に沿って入射 14MeV 中性子

13 インポータンスを使用した中性子深層透過計算例
1.50GHz, single ヒストリー数 1,500 [importance] 未使用 total cpu time = sec 1 Relative error Dose imp-dose-xz.eps imp-dose-xz_err.eps [importance] 使用 total cpu time = sec 1 相対誤差の確認重要（赤1.0、黄 約0.1、緑 約0.01 ）

14 インポータンスを使用した中性子深層透過計算例
元データ：imp-dose-reg.out 領域（15cm厚）毎の線量変化 14

15 インポータンス利用の注意点 セル間の大きなインポータンスの比は良くない。多くの粒子分割を一度に引き起こす。 良い例 粒子 1 2 4 8 16 32 悪い例 粒子 1 1 8 8 8 32 セル間のインポータンスの比は最大２～３程度が良い 参考文献 “A Sample Problem for Variance Reduction in MCNP” LA MS DE

16 先ほどの結果と比較して相対誤差が非常に大きい
インポータンスの誤った使用例 悪い例： Cell 5と6の間に急激なImportanceギャップ 大きなインポータンスのギャップを作って実行してみよう。 imp-dose-xz.eps [[importance] part = neutron reg imp **0 **0 **0 **0 **0 **5 **6 **7 **8 **9 **10 **11 imp-dose-xz_err.eps 先ほどの結果と比較して相対誤差が非常に大きい

17 実習内容 １．はじめに ２．中性子深層透過計算 [importance]
　Geometry splitting and Russian Roulette (セルインポータンス法） [weight window] 　Weight windows (ウエイト・ウインドウ法) ３．薄膜からの生成粒子計算 [forced collisions] 　強制衝突法

18 セルインポータンス法とウェイトウィンドウ法の違い
セルインポータンス法は，領域ごとにウェイトが取りえる値を設定 ウェイトウィンドウ法は，領域かつエネルギー群ごとに粒子が取りえるウェイトの幅（ウィンドウ）を設定 重要となるエネルギー領域（例えば高速中性子など）に， より焦点を絞ったシミュレーションが可能 セルインポータンス法 ウェイトウィンドウ法 WU=2.5 WU=2.5 I1=1 I2=3 WU=1.25 W1=1 W=1 W1=1 ウェイト ウェイト ウェイト WU=0.83 WL=0.5 WL=0.5 W2=1/3 WL=0.25 W2=0.5 領域1 エネルギー群1 領域2 エネルギー群2 領域1 領域2 WL=0.17 粒子数 1 　→ I2/I1=3 粒子数 1 → 1 粒子数 1→ W1/W2=2

19 ウェイトウィンドウ法に関連するパラメータ
wupn: ウェイトウィンドウの上限値の下限値の比（D=5.0） mxspln: 1回の粒子分割で分割する粒子の最大数（D=5.0） mwhere: ウェイトウィンドウが考慮されるタイミング（D=0） 　　　　　　　　-1：核反応時，0:両方，1:領域横断時 詳細はマニュアル「4.2.4　時間カット，ウェイトカット，ウェイトウィンドウ」を参照

20 1ヒストリーのときの中性子挙動を確認 weight-hist.inpを実行してみよう 水 空気 （半径10cm， 深さ5cm×2） 5MeV
weight-trajectory.eps

21 ウェイトウィンドウを設定した場合の中性子挙動を確認
weight-hist.inpの１つ目の[weight window]を有効にして実行してみよう weight-hist.inp [weight window] part = neutron reg ww1 ウェイトの下限値WLを[weight window]セクションで設定する 上限値は，自動的にその5倍に設定される（wupnパラメータ） 領域毎の指定のみ （エネルギー毎に分割しない） WU=5 x WL =1.25 W1=1 ウェイト WL=0.25 WL=0.25 領域1 領域2 weight-trajectory.eps 中性子のウェイトがウェイトの幅の中にあるので、粒子分割や粒子消滅は起きない

22 ウェイトウィンドウを下げたときの中性子挙動を確認
weight-hist.inp [weight window] part = neutron reg ww1 領域２のウェイトウインドウの下限値を下げて実行してみよう WU=1.25 W1=1 WU=0.625 ウェイト WL=0.25 WL=0.125 領域1 領域2 粒子分割 weight-trajectory.eps 領域２で粒子ウェイトが領域ウェイトの上限値を上回るので２つに分裂した

23 ウェイトウィンドウを極端に下げたときの中性子挙動を確認
weight-hist.inp [weight window] part = neutron reg ww1 領域２のウェイトウインドウの下限値を極端に下げて実行してみよう WU=1.25 W1=1 ウェイト WL=0.25 WU=0.025 WL=0.005 領域1 領域2 weight-trajectory.eps 40個になるまで粒子分割（ただし1度に分割するのは最大5粒子まで） セル間のウェイトウィンドウの比も最大２～３程度が良い

24 エネルギー群毎にウェイトウィンドウを設定
weight-hist.inp [weight window] part = neutron eng = 2 e5 reg ww ww2 １つ目の[weight window]を無効にして、 ２つ目の[weight window]を有効にしよう 各エネルギー群の最大値 エネルギー毎に異なるウェイトウィンドウをかける ただし，今回は，別々に設定するが値は同じ WU=1.25 WU=1.25 W1=1 W1=1 WU=0.625 WU=0.625 ウェイト ウェイト WL=0.25 WL=0.25 WL=0.125 WL=0.125 領域1 領域2 領域1 領域2 粒子分割 粒子分割 E < 0.01 MeV 0.01 MeV < E < 105MeV weight-trajectory.eps 領域２で粒子ウェイトが領域ウェイトの上限値を上回るので２つに分裂した（2ページ前と同じ）

25 低エネルギー中性子のウェイトウィンドウを上げる
weight-hist.inp [weight window] part = neutron eng = 2 e5 reg ww ww2 * * 0.01MeV以下の中性子のウェイトウインドウの下限値のみ10倍に上げてみよう WU=12.5 WU=6.25 WU=1.25 WL=2.5 低エネルギー中性子が消滅 WL=1.25 W1=1 W1=1 WU=0.625 ウェイト ウェイト WL=0.25 WL=0.125 領域1 領域2 領域1 領域2 粒子消滅 粒子分割 E < 0.01 MeV 0.01 MeV < E < 105MeV weight-trajectory.eps 低エネルギー中性子の発生を抑え，高エネルギー中性子を中心に解析する 浅い場所での統計が下がり，より深い場所での統計が溜まりやすくなる

26 ウェイトウィンドウを使用した中性子深層透過計算例
weight.inpを走らせて，計算時間を確認してみよう weight.inp 体系，ヒストリー数などはimp.inpと同じ [weight window] set: c10[1.0] part = neutron eng = 2 1.0e e5 reg ww ww2 1 c10*2.5**(-1) 2.5**(-1) 2 c10*2.5**(-2) 2.5**(-2) 3 c10*2.5**(-3) 2.5**(-3) ・・・ 半径50cm、深さ1.8mの コンクリート円柱 14MeV 中性子 セルの厚さは15cm total cpu time = sec 領域毎にウェイトウィンドウを1/2.5に下げて粒子分割する 全エネルギーに同じウェイトウィンドウをかける（c10=1.0のため） weight-dose-xz.eps weight-dose-xz_err.eps

27 ウェイトウィンドウを使用した中性子深層透過計算例
weight.inpの低エネルギー中性子のウェイトウィンドウを10倍にして実行し，計算時間を計ってみよう weight.inp [weight window] set: c10[10.0] part = neutron eng = 2 1.0e e5 reg ww ww2 1 c10*2.5**(-1) 2.5**(-1) 2 c10*2.5**(-2) 2.5**(-2) 3 c10*2.5**(-3) 2.5**(-3) ・・・ 低エネルギー中性子は，生成された瞬間にウェイトウィンドウ幅から外れるので，ロシアンルーレットにかけられ，ほとんどが殺されてしまうため，計算時間が短縮される ただし，低エネルギー中性子は飛程が短いので，低エネルギー中性子をカットしても深部での統計誤差はそれほど変わらない total cpu time = > 14.53 sec weight-dose-xz.eps weight-dose-xz_err.eps

28 計算時間の比較 total cpu time = 19.53 sec total cpu time = 50.49 sec
普通の計算 total cpu time = sec [importance] 使用 total cpu time = sec [weight window] 領域毎に設定 total cpu time = sec [weight window] 領域・エネルギー毎に設定 total cpu time = sec

29 コンクリート中の線量率減衰 通常計算の深い場所を除いて良く一致
90 cm 180 cm 通常計算の深い場所を除いて良く一致 エネルギー群で異なるウェイトウィンドウを設定すれば，より効率的な粒子輸送の計算が可能となる

30 実習内容 １．はじめに ２．中性子深層透過計算 [importance]
　Geometry splitting and Russian Roulette (セルインポータンス法） [weight window] 　Weight windows (ウェイト・ウインドウ法) ３．薄膜からの生成粒子計算 [forced collisions] 　強制衝突法

31 強制衝突とは？ 放射線を衝突の確率の低いターゲット（薄膜）に強制的に衝突させる。 応用例： 薄膜を用いた断面積測定実験を模擬するときに利用
10μm厚さSi 10μm厚さSi 100MeV陽子 強制衝突なし 陽子は薄膜をそのまま通過。 強制衝突あり 陽子は薄膜と必ず衝突。 応用例：　薄膜を用いた断面積測定実験を模擬するときに利用

32 薄膜を用いた断面積測定実験を模擬しよう force.inpを実行してみよう
[t-product]: ターゲット内で生成するα粒子とSiの生成エネルギー分布　→　真の断面積 [t-cross]：検出器に入射したα粒子とSiのエネルギー分布　→　測定上の断面積 force.inp 強制衝突無し：殆どの入射粒子は衝突しない force.inpを実行してみよう product.eps maxcas = maxbch = 半径0.5cm, 10μm厚さSi cross.eps 100MeV陽子 真空 検出器

33 強制衝突法 薄膜から生成する粒子を解析するときに利用 二つに分割 非衝突粒子 非衝突のウェイト：Wi×exp(-σl) l:セルを横切る長さ
強制衝突領域 σ: 巨視的断面積 衝突位置は断面積に従って確率的に決定 強制衝突法はヒストリー分散を減らすが、ヒストリー当たりの計算時間が増える。 セルでは一つ以上の衝突が起こる。 ｆｃｌ：強制衝突コントロールパラメータ　 （通常　fcl = -1 を使用） fcl = -1：強制衝突による生成粒子は通常の衝突をさせる（ウェイトカットオフは考慮しない） fcl = 1：強制衝突による生成粒子は，ウェイトカットオフになるまで強制衝突させる

34 強制衝突させてみよう 強制衝突の設定を行おう force.inp ターゲット内での粒子生成時点の情報をタリー product.eps
[Forced Collisions] off part = proton reg fcl offを削除して実行 強制衝突による生成粒子のウェイトが、ウェイトカットオフwc2を下回る。 ⇒　Russian Rouletteが働くため、粒子輸送されにくい。 track-xz.eps cross.eps α

35 強制衝突で発生させた2次粒子を輸送するには？
強制衝突で生成する粒子を輸送させるため、 ウェイトカットオフwc2を極力低く設定してみよう。 force.inp [parameters] ………. wc2(1) = 1.0E-12 # weight cutoff of proton wc2(2) = 1.0E-12 # weight cutoff of neutron wc2(14) = 1.0E-12 # weight cutoff of photon wc2(15) = 1.0E-12 # weight cutoff of deuteron wc2(16) = 1.0E-12 # weight cutoff of triton wc2(17) = 1.0E-12 # weight cutoff of 3He wc2(18) = 1.0E-12 # weight cutoff of Alpha wc2(19) = 1.0E-12 # weight cutoff of Nucleus product.eps cross.eps track-xz.eps α Si 多くの生成Siがターゲット内で止まるため，この実験からは断面積が正しく測定できない

36 まとめ 深層透過計算の場合は、importance法やweight window法が有効
隣り合うセル間のインポータンスやウェイトウィンドウ値の比は、大きすぎないようにすること。　ファクター２～３以下が望ましい。 ウェイトウィンドウを活用することで、さらに効率的な中性子輸送の計算が可能。 薄膜における生成粒子情報を計算するには、forced collisionsが有効