Transcript 中性子深層透過計算

PHITS
Multi-Purpose Particle and Heavy Ion Transport code System
PHITS 応用実習②：
計算効率を上げるためのvariance
reduction（粒子のウェイト）の利用
2015年08月改訂
title
1
実習内容
１．はじめに
２．中性子深層透過計算
[importance]
Geometry splitting and Russian Roulette (セルインポータンス法）
[weight window]
Weight windows (ウェイト・ウインドウ法)
３．薄膜からの生成粒子計算
[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.50GHｚ 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 タリー
ウェイト： モンテカルロ計算における粒子の重要度、通常の計算は１*。
ウェイトを操作して、統計上発生しにくいイベントを人為的に発生させる。
ただし，ウェイトを操作した場合，ヒストリー毎の頻度分布（[t-deposit], output
= depositなど）の計算ができなくなるので注意が必要。
*ただし、核データを使用した中性子輸送計算では、ウェイトが変化する
5
実習内容
１．はじめに
２．中性子深層透過計算
[importance]
Geometry splitting and Russian Roulette (セルインポータンス法）
[weight window]
Weight windows (ウェイト・ウインドウ法)
３．薄膜からの生成粒子計算
[forced collision]
強制衝突法
6
セルインポータンス法
•領域(セル)にインポータンスIを設定
•中性子がセル1とセル2を横切る時、I2/I1を設定
I2/I1 > 1: Particle Split （粒子の分割）
例 I2/I1 が整数の場合
I2/I1 が非整数の場合
I2/I1 = 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%の確率で３つの粒子に分割
25%の確率で２つの粒子に分割
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
標準では、インポータンスは全て１
imp-trajectory.eps
8
インポータンスを大きくしたときの中性子挙動を確認
imp-hist.inp
I10=1
[importance]
part = neutron
reg imp
10 1
1 2
2 4
I1=2
ここで反応が起きた
I2=4
この時点で既に２つ
に分かれている
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
セル間のインポータンスの比は最大２～３程度が良い
“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
悪い例
粒子
セル間のインポータンスの比は最大２～３程度が良い
参考文献
“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
実習内容
１．はじめに
２．中性子深層透過計算
[importance]
Geometry splitting and Russian Roulette (セルインポータンス法）
[weight window]
Weight windows (ウエイト・ウインドウ法)
３．薄膜からの生成粒子計算
[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の１つ目の[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
ウェイトウィンドウを下げたときの中性子挙動を確認
領域２のウェイトウインドウの下限値を
下げて実行してみよう
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
領域２で粒子ウェイトが領域ウェイトの上限値を上回るので２つに分裂した
22
ウェイトウィンドウを極端に下げたときの中性子挙動を確認
領域２のウェイトウインドウの下限値を
極端に下げて実行してみよう
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粒子まで）
セル間のウェイトウィンドウの比も最大２～３程度が良い
23
エネルギー群毎にウェイトウィンドウを設定
weight-hist.inp
１つ目の[weight window]を無効にして、
２つ目の[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ページ前と同じ）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
実習内容
１．はじめに
２．中性子深層透過計算
[importance]
Geometry splitting and Russian Roulette (セルインポータンス法）
[weight window]
Weight windows (ウェイト・ウインドウ法)
３．薄膜からの生成粒子計算
[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 = -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法が有効
隣り合うセル間のインポータンスやウェイトウィンドウ値の比は、大き
すぎないようにすること。 ファクター２～３以下が望ましい。
ウェイトウィンドウを活用することで、さらに効率的な中性子輸送の計
算が可能。
薄膜における生成粒子情報を計算するには、forced collisionsが有効
35