7.16. [ T-WWG ] セクション

[t-wwg] は、[weight window] セクションのパラメータを出力するタリーです。 想定している3次元体系に対して、weight window 機能が効果的に動作するパラメータを自動で得ることができる Weight Window Generator (WWG) です。 基本的には、まず [t-wwg] を設定したインプットで仮計算を行い、次に、得られた [weight window] のパラメータを加えたインプットで本計算を行う、という流れで本タリーを利用します。 mesh=reg の場合は [volume] セクション等により各セルの体積を与えておく必要があります。 また、再開始計算機能や sumtally 機能を利用したい場合は、wwg 以外の axis を定義した上で、[parameters] セクションで ireschk=1 とする必要があります。 ただし、粒子誘導機能を利用して chwei(i) パラメータを指定した場合は、再開始計算機能や sumtally 機能は使えません。

出力する [weight window] の下限値 \(W_{klm,{\rm low}}\) は、仮計算で得られた各メッシュのフルエンス \(\Sigma_j\Phi_{jklm}\) に基づいて次式で決定します。ここで、\(j\) はヒストリー、\(k\) は粒子、\(l\) はエネルギーメッシュ、\(m\) は幾何学メッシュを表します。

(7.16.1)\[W_{klm,{\rm low}} = C_{kl} \left[ \left( \Sigma_j \Phi_{jklm} \right)_{k,l,x{\rm \%ile}} + \Sigma_j \Phi_{jklm} \right].\]

ここで \(\left( \Sigma_j \Phi_{jklm} \right)_{k,l,x{\rm \%ile}}\) は、粒子 \(k\)、エネルギーメッシュ \(l\) に対するフルエンスの下限 \(x\) パーセンタイル値を表し、\(x\) の値は pedestal パラメータより指定可能です。初期値は 10 パーセンタイルです。 Version 3.34 までは自動的に最小値を選ぶようにしていましたが、\(\Sigma_j\Phi_{jklm}\) が 0 の場合が多いと粒子の無限増殖が頻発してしまうため、ある程度の大きさの値を設定するように変更しました。 式 (7.16.1)\(C_{kl}\) は規格化定数で、normww パラメータに応じて次式で決定します。

(7.16.2)\[\begin{split}C_{kl} = \begin{cases} \frac{W_{\rm max}}{ \left( \Sigma_j \Phi_{jklm} \right)_{m,{\rm max}} } & {\rm for \, normww=0} \\ \frac{W_{\rm max}}{ \left( \Sigma_j \Phi_{jklm} \right)_{lm,{\rm max}} } & {\rm for \, normww=1} \\ \frac{W_{\rm max}}{ \left( \Sigma_j \Phi_{jklm} \right)_{klm,{\rm max}} } & {\rm for \, normww=2} \\ \frac{W_{\rm max}}{ \left( \Sigma_j \Phi_{jklm} \right)_{m,{\rm max}} } \frac{ \left( \Sigma_{jm} \Phi_{jklm} \right) } {\left( \Sigma_{jm} \Phi_{jklm} \right)_{l,{\rm min}}} & {\rm for \, normww=-1} \\ \frac{W_{\rm max}}{ \left( \Sigma_j \Phi_{jklm} \right)_{m,{\rm max}} } \frac{ \left( \Sigma_{jm} \Phi_{jklm} \right) } {\left( \Sigma_{jm} \Phi_{jklm} \right)_{kl,{\rm min}}} & {\rm for \, normww=-2} \end{cases}\end{split}\]

ここで添字の \(m,max, lm,max, klm,max\) は、それぞれの添字に対する最大値、 \(l,min, lm,min\) はそれぞれの添字に対する最小値を示します。 例えば、 \(\left( \Sigma_j \Phi_{jklm} \right)_{lm,{\rm max}}\) は全てのエネルギー群l及び空間メッシュmの中で最大となる \(\left( \Sigma_j \Phi_{jklm} \right)\) という意味で、粒子kごとに異なる値となります。 normww が 0, 1, 2 の場合は、それぞれ m、lm、klm のメッシュ点において最大となるフルエンスに基づいて規格化が行われます。 一方、 normww が -1, -2 の場合は、m に関する和、すなわち空間全体にわたる総フルエンス量が最小となるエネルギー群(-1 の場合)あるいは粒子・エネルギー群(-2 の場合)を基準とした補正項が加えられます。 これにより、エネルギー群ごとの総フルエンス量の偏り(-1 の場合)や、粒子・エネルギー群ごとの総フルエンス量の偏り(-2 の場合)が補正され、空間的にモンテカルロ粒子の分布がより均衡するようなウェイトウインドウ値が得られます。 \(W_{\rm max}\)wwmax パラメータで指定する \(W_{klm,{\rm low}}\) の最大値で、デフォルトは 0.99 です。 \(W_{\rm max}\) のデフォルト値が 1 でない理由は、桁落ちの関係でウェイト 1 の線源が発生した際にいきなり分割されてしまうのを防ぐためです。

通常の設定では、指定した領域内のモンテカルロ粒子ができる限り均等になるような [weight window] を生成します。 一方、遮へい体後方における特定領域の線量を短時間で精度良く評価したい場合などは、ヒストリーカウンターと chwei(i) パラメータで指定する誘導強度 \(B_{\rm max}\) を組み合わせることにより、その領域にモンテカルロ粒子が届きやすくなるようバイアスを掛けた粒子誘導型の [weight window] を生成することができます。 各領域のバイアスファクター \(B_{klm}\) は、その領域のモンテカルロ粒子フルエンス \(\varphi_{jklm}\) より次式で決定します。

(7.16.3)\[B_{klm} = \frac{\Sigma_j H_{ij} \varphi_{jklm}} {\left( \Sigma_j H_{ij} \varphi_{jklm} \right)_{k,l,{\rm max}} } \left[ B_{\rm max} -1 \right] +1.\]

ここで \(H_{ij}\) は、\(j\) 番目のヒストリーにおけるカウンター番号 \(i\) のヒストリーカウンター値を表します。 \(B_{klm}\) は、\(\Sigma_j H_{ij} \varphi_{jklm}\) が最大となる領域で \(B_{\rm max}=\textbf{chwei(i)}\) となり最も大きくなります。 本機能を用いた場合、式 (7.16.1) 及び (7.16.2)\(\Sigma_j \Phi_{jklm}\)\(\Sigma_j \Phi_{jklm}/B_{klm}\) に置換され、\(B_{klm}\) が大きい領域ほど \(W_{klm,{\rm low}}\) が小さくなります。 詳細は、文献 [1] や講習会資料 [2] を参照してください。

[t-wwg] の入力書式は、以下の通りです。

表 7.16.1 mesh

value

explanation

reg, xyz, tet

メッシュ型。メッシュ型サブセクションが必要です。 reg, xyztet のみ指定可能です。

表 7.16.2 part

value

explanation

all (省略時), 粒子名

タリーする粒子。

表 7.16.3 material

value

explanation

(省略可)

タリーする物質を指定します。

all

all が default で、未定義と同じです。

物質数

物質数を指定し、次の行で物質番号を定義します。 この値を負値にすると、指定した物質をタリー対象から除外します。

(次の行)

2 5 8 のように物質番号を指定します。 1行に書ききれない場合は、継続行の記号( \ )を使わずに複数行にわたって記述できます。

表 7.16.4 e-type

value

explanation

1, 2, 3, 4, 5

エネルギーメッシュ。 メッシュサブセクションの定義方法は 6.6.1 章 を参照。

表 7.16.5 t-type

value

explanation

1, 2, 3, 4, 5 (省略可)

時間メッシュ。 メッシュサブセクションの定義方法は 6.6.1 章 を参照。

表 7.16.6 unit

value

explanation

1 (省略可)

1: [1/cm^2/source]。1 のみ指定可。

表 7.16.7 axis

value

explanation

eng, reg, tet, xy, yz, xz, t, wwg

出力データの x 軸。 wwg 以外は [t-track] と同じです。 再開始計算が必要な場合は、2つ以上の axis を定義し、最初の1つを wwg 以外で指定します。

このタリーでは、meshregxyztet のみ指定可能で、r-z は指定できません。これは、[weight window] セクションのパラメータが r-z の形式に対応していないためです。 mesh=xyz を指定した場合は、x, y, z メッシュの面と cell の面が一致しないよう注意してください。特に、パラメータを使ってジオメトリと xyz メッシュ範囲を連動させている場合は注意が必要です。 なお、Version 3.34 より、[weight window] 及び [t-wwg]xyz メッシュ最大値・最小値が整数で定義された場合は、その整数から最小値は -1*deltxyz、最大値は +2*deltxyz ずらすようにしています。Q2.9: Lost particleが頻発します。どう対処したらよいでしょうか? も参照してください。 mesh=tet は連続四面体形状を含む体系でしか使用できません。連続四面体形状を構成する各四面体要素の値を出力することができます。 出力の axis には、基本的に [weight window] セクションのパラメータを出力する wwg を指定してください。 この他に engregtetxyyzxzt も指定できますが、これらは [t-track] で指定した場合と同じ結果しか得られず、[weight window] に関係する値は得られません。 また、xyyzxz を選択した場合、誤差ファイルが出力されないため再開始計算を行うことができません。

表 7.16.8 file

value

explanation

file name

出力ファイル名を定義します。 これは axis の各設定に対して必要です。

表 7.16.9 resfile

value

explanation

(省略可, D=file

restart 計算で過去の tally を読み込むファイル名を定義します。 axis を複数定義していても、指定する resfile は1つだけです。

表 7.16.10 factor

value

explanation

(省略可, D=1.0

規格化定数。

表 7.16.11 title

value

explanation

(省略可)

タイトル。

表 7.16.12 angel

value

explanation

(省略可)

ANGEL パラメータ。

表 7.16.13 rshow

value

explanation

0 (省略時), 1, 2, 3

mesh=reg, axis=xy,yz,xz の時、領域境界(1)、物質番号(2)、領域番号(3)を表示します。 この下に xyz のメッシュ型サブセクションが必要です。

表 7.16.14 resol

value

explanation

1 (default)

gshow または rshow オプションの時、領域境界を求める分解能を各辺 resol 倍します。

表 7.16.15 width

value

explanation

0.5 (default)

線の太さを定義します。

表 7.16.16 volume

value

explanation

(省略可)

reg メッシュの時に、各領域の体積を定義します。 この行の下に volume 定義文が必要です。 省略した場合、[volume] セクションで定義された体積の値が使用されます。 また、reg= において ( ) 等により特殊な領域指定を行った場合は、volume を省略すると、 インプットエコーに PHITS 内部で与えられた領域番号が出力されます。

reg vol

volume 定義文。書式は [volume] を参照。

表 7.16.17 epsout

value

explanation

0 (default), 1, 2

epsout=1 では、結果を eps ファイルとして出力します。 ファイル名は拡張子を .eps に置き換えたものになります。 epsout=2 では、2次元型の axis=xy, yz, xz, rz を除き、eps ファイルにエラーバーも表示します。

表 7.16.18 ph5out

value

explanation

0 (default), 1

1でタリー出力結果をPHITS用HDF5ファイルフォーマットph5でタリーを出力する。ファイル名は出力ファイルの拡張子をph5に変えたファイル名。

表 7.16.19 trcl

value

explanation

(省略可)

r-z または xyz メッシュに対する座標変換番号、またはその定義。

表 7.16.20 gslat

value

explanation

2 (default)

gshow または rshow を指定した際の lattice および連続四面体の境界線に関するオプション。

0

描画しない。

1

描画する。 ただし、ボクセルファントムの線は明瞭に描画されないことがあります。

2

同じセル内の線は描画しない。

3

同じ物質内の線は描画しない。 この場合、同じ物質で満たされた隣接セルの境界は、lattice や連続四面体でなくても描画されません。 また、 gshow>=3 の場合でもセル番号や LAT 番号は表示されなくなります。

表 7.16.21 stdcut

value

explanation

(省略可, D=-1

STD cut off のしきい値。

stdcut を指定すると、PHITS は STD、すなわち標準偏差の値に応じて自動的 に計算を停止します。 この機能は stdcut が正で、 [parameters] セクションにおいて itall=0,1 を指定した場合に利用できます。 1バッチの最後に、タリー結果の STD の相対値がすべて 0 より大きく stdcut より小さい場合、計算は停止します。 2つ以上のタリーセクションで stdcut を設定した場合は、それらすべてが 条件を満たした時にこの機能が動作します。

表 7.16.22 pedestal

value

explanation

(省略可、D=0.1)

フラックスが 0 の領域における weight window の下限値を決定するパーセンタイル値 (式 (7.16.1)\(x\))÷100。

表 7.16.23 normww

value

explanation

(省略可、D=-2)

粒子、エネルギー群間でモンテカルロ粒子数の調整の有無を決めるオプション。 0: 調整を行わない。 +/-1: 各エネルギー群のモンテカルロ粒子数が等しくなるように調整する。 +/-2: 各粒子・エネルギー群のモンテカルロ粒子数が等しくなるように調整する。

normww を正値で指定した場合は、全体的にウェイトウィンドウ値が低くなるため、主に粒子分割により各粒子のウェイトをコントロールします。 一方、負値で指定した場合は、全体的にウェイトウィンドウ値が高くなるため、ロシアンルーレットが頻繁に行われるようになります。 幅広いエネルギー分布を持つ線源の場合など、ヒストリ毎の分散が大きい計算条件に対しては、1 ヒストリあたりの計算時間が短くなる負値の方が単位計算時間当たりの統計が良くなります。したがって、normww のデフォルト値は -2 としています。 ただし、normww=-2 かつ各粒子・各エネルギー群の中で 1 つでも粒子フルエンスが極端に低い群がある場合、それ以外の群に属するほとんどの粒子がロシアンルーレットによりカットされてしまいますので、粒子やエネルギー群の設定には十分注意してください。 part に関しては、all もしくは neutron photon を設定することを推奨します。エネルギー群に関しては、最大でも 2 から 3 群程度とし、タリーへの寄与の大きいエネルギー帯の群幅を少し狭く設定してください。 どのような設定が最適か不明な場合は、part=allne=1 を推奨します。その場合、計算結果は normww には依存しません。

表 7.16.24 wwmax

value

explanation

(省略可、D=0.99)

規格化定数。 \(W_{klm,{\rm low}}\) の最大値に相当します。( (7.16.2) 参照)

表 7.16.25 updateww

value

explanation

1 (default)

計算中に weight window の下限値を自動更新するかを指定します。 インプットに [weight window] セクションが既に存在する場合は、自動的に 0 に設定されます。

0

[t-wwg] は weight window ファイルの出力のみを行い、 計算中の weight window は変化しません。

1

各バッチの集計結果から weight window を自動的に更新します(単一計算での weight window 最適化)。

表 7.16.26 wwgcut

value

explanation

(省略可、D=0.1

weight window の収束による計算停止のしきい値。他のタリーにおけるstdcutと似たような役割を果たします。 updateww=1 (各バッチ間の終了時にWW値を更新する)の場合、計算条件が途中で変更されるため 原理的にその統計誤差を出すことができません。その代わり、 バッチ間でのWW値の変化率が wwgband 以上の条件割合(change fraction)を出力します。 その値が2バッチ連続で wwgcut より小さくなった場合、PHITSは計算結果が収束したと判定し計算を終了します。 2 つ以上の [t-wwg] セクションで wwgcut を設定した場合は、それらすべてが 条件を満たした時にこの機能が動作します。

表 7.16.27 wwgband

value

explanation

(省略可、D=0.1

wwgcut による収束判定で、バッチ間で更新前後の WW 下限値の比がこの値以上となった条件を「変化した条件」として数えます。

表 7.16.28 ctmin(i)

value

explanation

(省略可, D=-9999

i 番目の counter の最小値。

表 7.16.29 ctmax(i)

value

explanation

(省略可, D=9999

i 番目の counter の最大値。

表 7.16.30 chmin(i)

value

explanation

(省略可, D=-9999

i 番目の history-counter の最小値。 このパラメータは batch variance mode の istdev=1 では指定できません。

表 7.16.31 chmax(i)

value

explanation

(省略可, D=9999

i 番目の history-counter の最大値。 このパラメータは batch variance mode の istdev=1 では指定できません。

表 7.16.32 chwei(i)

value

explanation

(省略可、D=0)

粒子誘導型 [t-wwg] を利用する際の誘導強度(式 (7.16.3)\(B_{\rm max}\) )。 i はカウンター番号です。 このパラメータは、カウンター番号に関わらず 1 つの [t-wwg] に対して 1 つしか 定義することはできません。 このパラメータを定義した場合、再開始計算や sumtally による結果の合算ができなくなります。 奨励される値は、おおよそ 20~200 で、値が大きいほどより強く粒子誘導するようになります。

表 7.16.33 chplane

value

explanation

(省略可、D=all

(7.16.3) の最大値を決定する際に最大値を検索する面。 デフォルト値は all で、全ての領域から最大値を検索します。 xy, yz, xz を指定した場合は、それらの面毎に最大値を決定して規格化します。 このパラメータは、mesh=xyz かつ chwei が定義された場合のみ有効です。 なお、トレースバック機能(tracemem > 0)を使えば、より効果的に粒子誘導が可能であり、 現在は利用が奨励されていません。

表 7.16.34 tracemem

value

explanation

(省略可、D=0

粒子誘導機能の精度を向上させるトレースバック機能で、1 スレッドあたりに確保する トレースバック用メモリの要素数(real(8) = 8 バイト単位)。 0 の場合はトレースバック機能を使用しません。 ただし、chwei が定義されており、かつ tracemem を定義した [t-wwg] が他に 存在しない場合は、自動的に 1000000 が設定されます。

粒子誘導機能( chwei )を用いた weight window 生成では、[counter] で対象とした領域に到達した粒子の履歴に重みを付けて weight window を生成します。しかし、version 3.36までの手法では、1ヒストリー内に対象領域に粒子が1つでも到達すれば、そのヒストリーで発生した全ての粒子が通った飛跡上にも粒子を誘導してしまうため、2次粒子の数が多い場合に効果的に機能しない問題がありました。そこで、version 3.37から、全ての粒子飛跡を一時メモリ(バッファ)に記録してからヒストリー終了後に再解析(トレースバック)することにより、対象領域に到達した粒子の親粒子の飛跡上に限定して粒子誘導できるようにしました。

tracemem (> 0)は、以下をすべて満たす場合にのみ使用できます。満たさない場合はエラーで停止します。

  • メッシュタイプが xyz または tet であること(reg は不可)。

  • chwei と併用すること(トレースバックは chwei モードを前提とします)。

  • tracemem > 0 を指定できる [t-wwg] タリーは、1 つのインプット内で 1 つのみです。

トレースバック用バッファは、tracemem で指定したサイズ × 並列数ぶん確保されます(OpenMP の各スレッドが独立したバッファを持つため)。したがって並列数が大きいほど必要メモリが増えます。 記録中にバッファが不足した場合は、適宜メモリを補充(再確保して拡張)しながら記録を継続します。 ただし、メモリ不足により補充(再確保)に失敗した場合は、警告を出力したうえで、そのヒストリにおける以降のトレースバック記録を中止します(中止したヒストリの飛跡は weight window 生成に反映されません)。 tracemem を十分大きく設定しておくと、補充の発生やメモリ不足による記録中止を避けられます。

トレースバックで使用したメモリの情報は、計算終了時に標準出力およびコンソール出力ファイルへ、次の形式で出力されます。

tracemem summary for [t-wwg]: max used =       <値>, allocated =       <値>
         reallocations =   <値>, allocation failures =   <値>
  • max used:いずれかのスレッドで使用されたバッファ要素数の最大値(ピーク使用量)。

  • allocated:最終的に確保されているバッファサイズ(1 スレッドあたりの要素数。補充により拡張された場合はその値)。

  • reallocations:バッファを補充(再確保で拡張)した回数。

  • allocation failures:補充(再確保)に失敗した回数。

max usedallocated**(初期の **tracemem)を大きく下回っていれば tracemem を小さくでき、reallocationsallocation failures が出ている場合は tracemem を大きくすることが推奨されます。

表 7.16.35 elowthre

value

explanation

(省略可、D=0

低エネルギー・アンバイアス法のしきい値エネルギー [MeV]。 このエネルギー以下の weight window 下限値を定数倍することにより、人為的に 低エネルギー側のモンテカルロ粒子密度を下げて計算時間を短縮することができます。 0 の場合は使いません。 また、この機能を使う場合は ne=1 としてください。

表 7.16.36 elowbias

value

explanation

(省略可、D=5.0

低エネルギー・アンバイアス法を利用する場合の低エネルギー群 weight window 下限値の倍率。

wwg の出力例を以下に示します。

リスト 7.16.1 wwg の出力例
[ Weight Window ]
     mesh =  reg

set:  c101[0.0]  c102[c101+2.13496E-08]  c103[1.25478E-04]
set:  c104[0.0]  c105[c104+9.29569E-09]  c106[2.08987E-04]

    part = neutron
     eng =  2
          1.00000E-03   1.00000E+03

    reg          ww1                    ww2
    1   (8.94092E-05+c101)/c103  (2.08987E-04+c104)/c106
    2   (1.25478E-04+c101)/c103  (1.26817E-04+c104)/c106
    3   (8.53835E-05+c101)/c103  (5.78131E-05+c104)/c106

ここでは、中性子に対して、2 群のエネルギービンの weight window のパラメータが定義されています。 パラメータの定義には c** で表される定数が使われています。 この場合、c101-c103 が第 1 エネルギー群、c104-c106 が第 2 エネルギー群に対するパラメータです。 これらのパラメータの中でユーザー自身が変更する可能性があるのは、各群に対して左端で定義される初期設定値が 0 のパラメータ c101c104 などです。 これらのパラメータは、各群の weight window にある一定値を加える際に利用します。