中性子深層透過計算 - PHITS Homepage

PHITS
Multi-Purpose Particle and Heavy Ion Transport code System
PHITS 応用実習②:
計算効率を上げるためのvariance
reduction(粒子のウェイト)の利用
2015年08月改訂
title
1
実習内容
1.はじめに
2.中性子深層透過計算
[importance]
Geometry splitting and Russian Roulette (セルインポータンス法)
[weight window]
Weight windows (ウェイト・ウインドウ法)
3.薄膜からの生成粒子計算
[forced collisions]
強制衝突法
2
中性子深層透過計算
厚い遮蔽体の中性子輸送を計算し、深さ毎の実効線量率分布を計算
まずはimp.inpを走らせて,計算時間を計ってみよう imp.inp
 ヒストリ数は1,500
 計算時間はphits.outの最後の方
(total CPU time)に出力される
半径50cm、深さ180cmのコンクリート円柱
空気
14MeV
中性子
コンクリート
[Parameters]
icntl = 0
maxcas = 1500
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
part = neutron
emax = 1000.0
mat
mset1
all ( 1.0 -102 )
実効線量は,[t-track]で計算した中性子フルエンスに内蔵線量換算係数を乗じて導出3
詳細はマニュアル(4.25[multiplier]&6.1[t-track])もしくはlecture\exercise1を参照
中性子深層透過計算
シングルコア(1.50GHz AMD)の計算
imp-dose-xz.eps
ヒストリー数 1,500
total cpu time = 19.53 sec
imp-dose-xz.eps
効率良く計算したい。
ヒストリー数 270,000
total cpu time = 660.05 sec
分散低減法 ウェイトの概念を利用
4
モンテカルロ計算におけるWeightの概念
例:track length タリー
ウェイト: モンテカルロ計算における粒子の重要度、通常の計算は1*。
ウェイトを操作して、統計上発生しにくいイベントを人為的に発生させる。
ただし,ウェイトを操作した場合,ヒストリー毎の頻度分布([t-deposit], output
= depositなど)の計算ができなくなるので注意が必要。
*ただし、核データを使用した中性子輸送計算では、ウェイトが変化する
5
実習内容
1.はじめに
2.中性子深層透過計算
[importance]
Geometry splitting and Russian Roulette (セルインポータンス法)
[weight window]
Weight windows (ウェイト・ウインドウ法)
3.薄膜からの生成粒子計算
[forced collision]
強制衝突法
6
セルインポータンス法
•領域(セル)にインポータンスIを設定
•中性子がセル1とセル2を横切る時、I2/I1を設定
I2/I1 > 1: Particle Split (粒子の分割)
例 I2/I1 が整数の場合
I2/I1 が非整数の場合
I2/I1 = 3: 3つの粒子にsplit(分割)
I1=1
W=1
I2/I1 = 2.75:
I2=3
W=1/3
W=1/3
W=1/3
I2/I1 < 1: Russian Roulette
(消滅か生存)
75%の確率で3つの粒子に分割
25%の確率で2つの粒子に分割
I1=1
I2=1/3
W=3
1/3の確率
W=1
消滅 2/3の確率
7
1ヒストリーのときの中性子飛跡を確認
imp-hist.inpを走らせて,中性子の飛跡を確認してみよう
空気
5MeV I =1.0
中性子 10
コンクリート
(半径20cm,
深さ8cm×2)
I1=1.0
I2=1.0
標準では、インポータンスは全て1
imp-trajectory.eps
8
インポータンスを大きくしたときの中性子挙動を確認
imp-hist.inp
I10=1
[importance]
part = neutron
reg imp
10 1
1 2
2 4
I1=2
ここで反応が起きた
I2=4
この時点で既に2つ
に分かれている
imp-trajectory.eps
9
インポータンスを小さくしたときの中性子挙動を確認
imp-hist.inp
I10=1
[importance]
part = neutron
reg imp
10 1
1 1/2
2 1/4
I1=1/2
I2=1/4
1/2の確率で生存
1/2の確率で消滅
imp-trajectory.eps
10
インポータンスを極端に大きくしたときの中性子挙動を確認
imp-hist.inp
I10=1
[importance]
part = neutron
reg imp
10 1
1 2
2 100
I1=2
粒子が分割しすぎて計算時間が無駄
になってしまう
(統計はあまり良くならない)
I2=100
imp-trajectory.eps
セル間のインポータンスの比は最大2~3程度が良い
“A Sample Problem for Variance Reduction in MCNP” LA-10363-MS DE86 004380
11
インポータンスを使用した中性子深層透過計算例
imp.inpの [importance]セクションを有効にして、
実行してみよう
半径50cm、深さ180cmのコンクリート円柱
セルの厚さは15cm
半径1cmの14MeV中性子をz=0からz軸に沿って入射
14MeV
中性子
Ii+1/Ii = 2.5 と設定
[importance]
part = neutron
reg imp
1 2.5**0
2 2.5**1
3 2.5**2
4 2.5**3
5 2.5**4
6 2.5**5
7 2.5**6
8 2.5**7
9 2.5**8
10 2.5**9
11 2.5**10
12 2.5**11
12
インポータンスを使用した中性子深層透過計算例
1.50GHz, single
[importance] 未使用
ヒストリー数 1,500
total cpu time = 19.53 sec
1
Dose
imp-dose-xz.eps
[importance] 使用
total cpu time = 50.49 sec
Relative
error
imp-dose-xz_err.eps
相対誤差の確認重要(赤1.0、黄 約0.1、緑 約0.01 )
13
インポータンスを使用した中性子深層透過計算例
元データ:imp-dose-reg.out
元データ:imp-energy-reg.out
領域(15cm厚)毎の線量変化
・各領域のスペクトルが良い一致
セル3の中性子エネルギースペクトル
妥当なインポータンスの設定
14
インポータンス利用の注意点
セル間の大きなインポータンスの比は良くない。多くの粒子分割を一度に引き起こす。
良い例
粒子
1
2
4
8
16
32
1
1
8
8
8
32
悪い例
粒子
セル間のインポータンスの比は最大2~3程度が良い
参考文献
“A Sample Problem for Variance Reduction in MCNP” LA-10363-MS DE86 004380 15
インポータンスの誤った使用例
大きなインポータンスのギャップを作って実行してみよう。
imp-dose-xz.eps
imp-dose-xz_err.eps
悪い例:
Cell 5と6の間に急激
なImportanceギャップ
[[importance]
part = neutron
reg imp
1 2.5**0
2 2.5**0
3 2.5**0
4 2.5**0
5 2.5**0
6 2.5**5
7 2.5**6
8 2.5**7
9 2.5**8
10 2.5**9
11 2.5**10
12 2.5**11
先ほどの結果と比較して相対誤差が非常に大きい
16
実習内容
1.はじめに
2.中性子深層透過計算
[importance]
Geometry splitting and Russian Roulette (セルインポータンス法)
[weight window]
Weight windows (ウエイト・ウインドウ法)
3.薄膜からの生成粒子計算
[forced collisions]
強制衝突法
17
セルインポータンス法とウェイトウィンドウ法の違い
•
•
セルインポータンス法は,領域ごとにウェイトが取りえる値を設定
ウェイトウィンドウ法は,領域かつエネルギー群ごとに粒子が取
りえるウェイトの幅(ウィンドウ)を設定
重要となるエネルギー領域(例えば高速中性子など)に,
より焦点を絞ったシミュレーションが可能
セルインポータンス法
W2=1/3
領域2
粒子数 1 → I2/I1=3
W=1
WU=2.5
WU=1.25
W1=1
ウェイト
I2=3
W1=1
領域1
WU=2.5
ウェイト
ウェイト
I1=1
ウェイトウィンドウ法
WL=0.5
WL=0.25
領域1
領域2
エネルギー群1 エネルギー群2
粒子数 1 → 1
WL=0.5
WU=0.83
W2=0.5
WL=0.17
粒子数 1→W1/W2=2
18
ウェイトウィンドウ法に関連するパラメータ
 wupn: ウェイトウィンドウの上限値の下限値の比(D=5.0)
 mxspln: 1回の粒子分割で分割する粒子の最大数(D=5.0)
 mwhere: ウェイトウィンドウが考慮されるタイミング(D=0)
-1:核反応時,0:両方,1:領域横断時
詳細はマニュアル「4.2.4 時間カット,ウェイトカット,ウェイトウィンドウ」を参照
19
1ヒストリーのときの中性子挙動を確認
weight-hist.inpを実行してみよう
5MeV
中性子
空気
水
(半径10cm,
深さ5cm×2)
weight-trajectory.eps
20
ウェイトウィンドウを設定した場合の中性子挙動を確認
weight-hist.inpの1つ目の[weight window]
を有効にして実行してみよう
 ウ ェ イ ト の 下 限 値 WL を [weight
window]セクションで設定する
 上限値は,自動的にその5倍に設定
される(wupnパラメータ)
weight-hist.inp
[weight window]
part = neutron
reg ww1
1 0.25
2 0.25
領域毎の指定のみ
(エネルギー毎に分割しない)
WU=5 x WL =1.25
ウェイト
W1=1
WL=0.25 W =0.25
L
領域1
領域2
weight-trajectory.eps
中性子のウェイトがウェイトの幅の中にあるので、粒子分割や粒子消滅は起きない
21
ウェイトウィンドウを下げたときの中性子挙動を確認
領域2のウェイトウインドウの下限値を
下げて実行してみよう
weight-hist.inp
[weight window]
part = neutron
reg ww1
1 0.25
2 0.125
WU=1.25
ウェイト
W1=1
WU=0.625
WL=0.25
領域1
WL=0.125
領域2
粒子分割
weight-trajectory.eps
領域2で粒子ウェイトが領域ウェイトの上限値を上回るので2つに分裂した
22
ウェイトウィンドウを極端に下げたときの中性子挙動を確認
領域2のウェイトウインドウの下限値を
極端に下げて実行してみよう
weight-hist.inp
[weight window]
part = neutron
reg ww1
1 0.25
2 0.005
WU=1.25
ウェイト
W1=1
WL=0.25
WU=0.025
WL=0.005
領域1
領域2
weight-trajectory.eps
40個になるまで粒子分割(ただし1度に分割するのは最大5粒子まで)
セル間のウェイトウィンドウの比も最大2~3程度が良い
23
エネルギー群毎にウェイトウィンドウを設定
weight-hist.inp
1つ目の[weight window]を無効にして、
2つ目の[weight window]を有効にしよう
 エネルギー毎に異なるウェイトウィンドウをかける
 ただし,今回は,別々に設定するが値は同じ
WU=1.25
WU=1.25
W1=1
W1=1
WL=0.25
WL=0.125
領域1
領域2
粒子分割
E < 0.01 MeV
WU=0.625
ウェイト
WU=0.625
ウェイト
[weight window]
part = neutron 各エネルギー群
の最大値
eng = 2
0.01
1.0e5
reg ww1
ww2
1
0.25
0.25
2
0.125 0.125
WL=0.25
WL=0.125
領域1
領域2
粒子分割
0.01 MeV < E < 105MeV
weight-trajectory.eps
領域2で粒子ウェイトが領域ウェイトの上限値を上回るので2つに分裂した(2ページ前と同じ)24
低エネルギー中性子のウェイトウィンドウを上げる
weight-hist.inp
0.01MeV以下の中性子のウェイトウインドウ
の下限値のみ10倍に上げてみよう
WU=12.5
WU=6.65
WL=2.5
W1=1
WU=1.25
WL=1.33
W1=1
ウェイト
ウェイト
WU=0.67
領域1
領域2
粒子消滅
E < 0.01 MeV
[weight window]
part = neutron
eng = 2
0.01
1.0e5
reg ww1
ww2
1
0.25*10 0.25
2
0.133*10 0.133
低エネルギー
中性子が消滅
WL=0.25
WL=0.133
領域1
領域2
粒子分割
0.01 MeV < E < 105MeV
weight-trajectory.eps
 低エネルギー中性子の発生を抑え,高エネルギー中性子を中心に解析する
 浅い場所での統計が下がり,より深い場所での統計が溜まりやすくなる 25
ウェイトウィンドウを使用した中性子深層透過計算例
weight.inp
weight.inpを走らせて,計算時間を確認してみよう
[weight window]
set: c10[1.0]
part = neutron
eng = 2
1.0e-3
reg ww1
1 c10*2.5**(-1)
2 c10*2.5**(-2)
3 c10*2.5**(-3)
・・・
体系,ヒストリー数などはimp.inpと同じ
半径50cm、深さ1.8mの
コンクリート円柱
1e5
ww2
2.5**(-1)
2.5**(-2)
2.5**(-3)
14MeV
中性子
セルの厚さは15cm
total cpu time = 25.89 sec
 領域毎にウェイトウィンドウを1/2.5に下げて粒子分割する
 全エネルギーに同じウェイトウィンドウをかける(c10=1.0のため)
weight-dose-xz.eps
weight-dose-xz_err.eps
26
ウェイトウィンドウを使用した中性子深層透過計算例
weight.inp
weight.inpの低エネルギー中性子のウェイトウィンドウを
10倍にして実行し,計算時間を計ってみよう
[weight window]
set: c10[10.0]
part = neutron
eng = 2
1.0e-3
reg ww1
1 c10*2.5**(-1)
2 c10*2.5**(-2)
3 c10*2.5**(-3)
・・・
低エネルギー中性子は,生成された瞬間にウェイトウィン
ドウ幅から外れるので,ロシアンルーレットにかけられ,ほ
とんどが殺されてしまうため,計算時間が短縮される
1e5
ww2
2.5**(-1)
2.5**(-2)
2.5**(-3)
weight-dose-xz.eps
ただし,低エネルギー中性子は飛程が短いので,低エネ
ルギー中性子をカットしても深部での統計誤差はそれほ
ど変わらない
total cpu time = 25.89 -> 14.53 sec
weight-dose-xz_err.eps
27
計算時間の比較
普通の計算
total cpu time = 19.53 sec
[importance] 使用
total cpu time = 50.49 sec
[weight window] 領域毎に設定
total cpu time = 25.89 sec
[weight window] 領域・エネルギー毎に設定
total cpu time = 14.53 sec
28
コンクリート中の線量率減衰
90 cm
180 cm
通常計算の深い場所を除いて良く一致
エネルギー群で異なるウェイトウィンドウを設定すれば,
より効率的な粒子輸送の計算が可能となる
29
実習内容
1.はじめに
2.中性子深層透過計算
[importance]
Geometry splitting and Russian Roulette (セルインポータンス法)
[weight window]
Weight windows (ウェイト・ウインドウ法)
3.薄膜からの生成粒子計算
[forced collisions]
強制衝突法
30
薄膜を用いた断面積測定実験を模擬しよう
• [t-product]: ターゲット内で生成するα粒子とSiの生成エネルギー分布 → 真の断面積
• [t-cross]:検出器に入射したα粒子とSiのエネルギー分布 → 測定上の断面積
強制衝突無し:殆どの入射粒子は衝突しない
force.inp
force.inpを実行してみよう
maxcas =
maxbch =
product.eps
5000
5
半径0.5cm, 0.01mm厚さSi
cross.eps
100MeV陽子
真空
検出器
31
強制衝突法
衝突の確率の低いターゲット(薄膜)から生成する粒子を解析するときに利用
二つに分割
非衝突粒子 l:セルを横切る長さ
入射粒子 Wi
衝突粒子
非衝突のウェイト:
Wi×exp(-σl)
衝突のウェイト:
Wi×{1-exp(-σl)}
強制衝突領域
σ: 巨視的断面積
衝突位置は断面積に従って確率的に決定
強制衝突法はヒストリー分散を減らすが、ヒストリー当たりの計算時間が増える。
セルでは一つ以上の衝突が起こる。
fcl:強制衝突コントロールパラメータ (通常 fcl = -1 を使用)
• fcl = -1:強制衝突による生成粒子は通常の衝突をさせる(ウェイトカットオフは考慮しない)
• fcl = 1:強制衝突による生成粒子は,ウェイトカットオフになるまで強制衝突させる
32
強制衝突させてみよう
ターゲット内での粒子生成時点の情報をタリー
force.inp
 強制衝突の設定を行おう
[Forced Collisions] off
part = proton
reg fcl
1 -1.0
track-xz.eps
α
product.eps
offを削除して実行
強制衝突による生成粒子のウェイトが、ウェイトカットオフwc2を下回る。
⇒ Russian Rouletteが働くため、粒子輸送されにくい。
cross.eps
33
強制衝突で発生させた2次粒子を輸送するには?
force.inp
 強制衝突で生成する粒子を輸送させるため、
ウェイトカットオフwc2を極力低く設定してみよう。
[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
34
多くの生成Siがターゲット内で止まるため,この実験からは断面積が正しく測定できない
まとめ
深層透過計算の場合は、importance法やweight window法が有効
隣り合うセル間のインポータンスやウェイトウィンドウ値の比は、大き
すぎないようにすること。 ファクター2~3以下が望ましい。
ウェイトウィンドウを活用することで、さらに効率的な中性子輸送の計
算が可能。
薄膜における生成粒子情報を計算するには、forced collisionsが有効
35