PHITS Multi-Purpose Particle and Heavy Ion Transport code System PHITS 講習会 基礎実習(III): 計算条件の設定 2015年9月改訂 title 1 はじめに • PHITS計算は,[Parameters]セクションで定義され る様々なパラメータ(物理モデルの選択など)により コントロールされています。 • 全てのパラメータには初期設定値があり,基本的に は変更する必要はありません。 • ただし,最適な計算結果を得るためには,いくつか のパラメータを変更する必要があります。 本実習では,パラメータの設定方法について学習します 2 本実習の目標 PHITSの便利な機能を活用し,統計数や核反応モデル を変化させ,最適な計算結果が得られるようになる 宿題体系における初期設定での 陽子(上)・中性子(下)フルエンス ヒストリー数,切断エネルギー,核データを適切 に設定した陽子(上)・中性子(下)フルエンス 3 実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 まとめ 4 計算モードの選択 通常輸送計算 反応なしモード 体系確認モード 5 計算モードの選択 体系確認モード まずは、基礎実習(II)で 行ったジオメトリを確認 するタリーを用いて、サ ンプルインプットの体系 を確認してみましょう。 6 体系の確認([t-3dshow]) lec03.inp file = lec03.inp [Title] ・・・・・・ [Parameters] icntl = 11 [t-3dshow]を実行 maxcas = 100 maxbch = 10 file(6) = phits.out [T-3Dshow] output = 3 material = -1 6 x0 = 0 y0 = 0 z0 = 0 e-the = 70 $ eye e-phi = 20 e-dst = 80 l-the = 20 $ light l-phi = 0 l-dst = 100 w-wdt = 50 $ window w-hgt = 50 w-dst = 25 heaven = z line = 1 shadow = 2 1層5cm幅の玉ねぎ構造 file = 3dshow.out title = Check onion structure using [T-3dshow] tally epsout = 1 ・・・・・・ 3dshow.eps 7 体系の確認(gshowオプション) lec03.inp file = lec03.inp [Title] ・・・・・・ [Parameters] icntl = 11 8 修正して実行 maxcas = 100 maxbch = 10 file(6) = phits.out [T-Track] mesh = xyz 各領域のセル番号と満たさ x-type = 2 れている物質名をチェック! nx = 100 xmin = -50. xmax = 50. y-type = 1 ny = 1 -5.0 5.0 z-type = 2 nz = 100 zmin = -50. zmax = 50. ・・・・・・ axis = xz file = track_xz.out title = Check source direction using [T-track] tally gshow = 3 epsout = 1 track_xz.eps 8 線源の確認([t-track]) lec03.inp file = lec03.inp [Title] ・・・・・・ [Parameters] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out 修正して実行 [T-Track] mesh = xyz x-type = 2 nx = 100 xmin = -50. xmax = 50. y-type = 1 ny = 1 -5.0 5.0 z-type = 2 nz = 100 zmin = -50. zmax = 50. e-type = 1 ne = 1 0. 200. unit = 1 axis = xz file = track_xz.out title = Check source direction using [T-track] tally epsout = 1 9 線源の飛跡確認 track-xz.eps 全ての物質がvoidとなるため,まっすぐに飛んでいく (線源の発生位置・方向が確認できる) 10 実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 まとめ 11 外部ファイルの挿入 lec03.inp file = lec03.inp [Title] ・・・・・・ [Parameters] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out set: c1[10] [Source] ・・・・・・ • onion.inp [Material] ・・・・・・ [Mat Name Color] ・・・・・・ [Surface] 10 so 500. 11 so 5. 12 so 10. 13 so 15. 14 so 20. 15 so 25. [Cell] 100 -1 10 infl: {onion.inp}[1-33] 101 1 -19.32 -11 102 2 -1. 11 -12 103 3 -8.93 12 -13 インプットファイルとは別のファイルを取り込む場合は”infl:”コマンドを使う 104 4 -1. 13 -14 {}でファイル名,[n1-n2]で取り込む行の範囲を指定する 105 5 -0.9 14 -15 106 6 -1.20e-3 15 -10 インプットファイルの一行目に“file = インプットファイル名”を書く(重要) • • “off”したセクションの下にある場合、無視されるので注意 12 ユーザー定義定数と数式の利用 lec03.inp file = lec03.inp [Title] ・・・・・・ [Parameters] ・・・・・・ set: c1[10] [Source] s-type = 1 proj = proton dir = 1.0 r0 = c1 z0 = -30. z1 = -30. e0 = 150 totfact = pi*c1**2 ユーザー定義定数の利用: 同じ値を使う場合は,予め定義しておくと便利! PHITSでは定数を定義して、インプットファイル内でその定数を参照 することが可能。 “set: ci[x]” •名前はc1からc99まで使用可能 •書いた箇所より下の部分で有効となる。また、何度も同じ名前を 定義できる •“off”で飛ばしたセクションの下にある場合、無視されるので注意 数式の利用: PHITS入力ファイルでは,FORTRAN形式の数式が 利用可能 • 円周率(pi)は定義しなくても使える • べき乗は「**」なので注意 • expやcosなどの関数も使える 13 規格化について lec03.inp file = lec03.inp [Title] ・・・・・・ [Parameters] ・・・・・・ set: c1[10] [Source] s-type = 1 proj = proton dir = 1.0 r0 = c1 z0 = -30. z1 = -30. e0 = 150 totfact = pi*c1**2 • PHITSのタリー結果は,発生する線源数*(試行 回数)当たりで出力される (*厳密には,線源粒子のウェイト当たり) • 測定値などと比較するには,Bqや/cm2/sで表され る線源の発生頻度で規格化する必要がある • タリー結果を定数倍(totfactで定義)することによ り,直接, /cm2/sやmGy/hなどの単位で出力でき るようになる • 今回の場合,規格化しないと線源のフルエンスは 1/πr2 (/cm2/source)となるため,タリー結果にπr2を 乗じることでフルエンスが1(/cm2)となるように規 格化している s-type=1でz0=z1とすると、半径 r0 r0の円の面が線源領域となる 14 実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 まとめ 15 体積、面積計算の原理 モンテカルロ積分: 乱数を用いて定積分の近似解を求める 方法。積分範囲の内か外かを判定することにより、積分値を 確率的に評価する。計算精度は試行回数に依存する。 体積計算にこの考え方を適用すると… 線の長さの合計 [cm] = 線の数(試行回数) [source] × 平均track length [cm/source] 求める体積 [cm3] = 線の長さの合計 ÷ 線の面密度 = 平均track length ÷ 平均flux = [t-track]の結果 × 線源面積 試行回数を増やすこ とで精度が向上する 線の面密度(間隔) [1/cm2] = 線の数(試行回数) [source] × 平均flux [1/cm2/source] あらかじめ 予測可能 (今回はπr2 ) 16 課題1: 体系の体積を計算してみよう lec03.inp [Parameters] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out set: c1[30] [Source] s-type = 1 proj = proton dir = 1.0 r0 = c1 z0 = -30. z1 = -30. e0 = 150 totfact = pi*c1**2 [T-Track] mesh = reg reg = 101 102 103 104 105 ・・・・・・ ・・・・・・ file = volume.out 修正して実行 体系の体積や面積を計算する場合: • icntl=5, • r0を体系が全部入る大きさにする 17 体積計算の結果 モンテカルロ積分の要領で体系の体積を求めることが可能 volume.out ・・・・ #num 1 2 3 4 5 ・・ reg 101 102 103 104 105 volume 1.0000E+00 1.0000E+00 1.0000E+00 1.0000E+00 1.0000E+00 flux 3.6914E+02 3.5499E+03 1.0295E+04 1.9822E+04 3.1834E+04 r.err 0.2362 0.0965 0.0577 0.0380 0.0247 各層の体積は? •4π(5)3/3=524 cm3 •4π(10)3/3 - 524=3665 cm3 •4π(15)3/3 - 4189=9948 cm3 •4π(20)3/3 - 14137=19373 cm3 •4π(25)3/3 - 33510=31940 cm3 内側の結果ほど一致していない → 統計量が十分ではない → 統計量を増やすには? 18 実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 まとめ 19 試行回数(ヒストリ数)の設定 モンテカルロ計算の統計精度は、シミュレーションの試行回 数に依存しています。したがって、信頼できる計算結果を得る ためには、十分な数の試行回数を設定する必要があります。 maxcas (D=10) 1バッチのヒストリー数 maxbch (D=10) バッチ数 rseed (D=0.0) rseed < 0 rseed = 0 rseed > 0 irskip (D=0) 乱数のコントロール irskip > 0 irskip回ヒストリーをスキップする irskip < 0 Irskip回乱数をスキップする maxcas ×maxbch =全試行回数 • 初期設定では,常に 初期乱数オプション 同じ乱数を使う 計算開始時間から初期乱数を決定 • 毎回違う結果を得た デフォルトの初期乱数 い場合は,rseed < 0 Rseedの値を初期乱数とする とする 20 ヒストリーとバッチの関係 Q. ヒストリー数とは? [source]セクションで指定した線源を発生させる回数→ぶどうの粒の数 Q. バッチ数とは? ある一定のヒストリー数(maxcas)を束(バッチ)にして,その束を実行する回数(maxbch) → maxcas: 1房当たりのぶどうの粒の数 maxbch:ぶどうの房の数 Q. 各バッチの終了時には何をするの? タリーの結果を取りまとめて平均値や統計誤差を導出し,途中経過を出力する (itall=1) →収穫したぶどうの房の味見をする Q. なぜ束(バッチ)に分けるの? 一気に全てのヒストリー数を実行すると,何か間違っていたときに再計算が大変! →味見もせずに全てのぶどうを収穫すると,収穫時期を間違えたときに大損害! 各ヒストリー終了時に統計処理をすると,計算時間が掛かりすぎてしまう →1粒1粒,味見をしながら収穫すると,収穫に時間が掛かりすぎてしまう 統計処理の時間が気にならない程度に,maxcasとmaxbchを調整する 例)1バッチ当たりの計算時間を2~3分にする 21 試行回数と統計誤差の関係 試行回数を増やすと統計誤差が小さくなる lec03.inp [Parameters] icntl = 5 maxcas = 100 1000 10000 maxbch = 10 file(6) = phits.out set: c1[30] [Source] s-type = 1 proj = proton dir = 1.0 r0 = c1 z0 = -30. z1 = -30. e0 = 150 totfact = pi*c1**2 相対統計誤差 volume.out ・・・・ #num 1 2 3 4 5 ・・ reg 101 102 103 104 105 volume 1.0000E+00 1.0000E+00 1.0000E+00 1.0000E+00 1.0000E+00 flux 3.6914E+02 4.8308E+02 5.2339E+02 3.5499E+03 3.6685E+03 3.6356E+03 1.0295E+04 1.0019E+04 9.9282E+03 1.9822E+04 1.9150E+04 1.9319E+04 3.1834E+04 3.1346E+04 3.1999E+04 r.err 0.0638 0.2362 0.0199 0.0965 0.0299 0.0095 0.0577 0.0184 0.0058 0.0380 0.0123 0.0039 0.0247 0.0080 0.0025 各層の体積は? •4π(5)3/3=524 cm3 •4π(10)3/3 - 524=3665 cm3 •4π(15)3/3 - 4189=9948 cm3 •4π(20)3/3 - 14137=19373 cm3 •4π(25)3/3 - 33510=31940 cm3 22 再開始計算 lec03.inp 統計が足りない過去のタリー結果を読み込んで計算の再開が可能 file = lec03.inp [Title] ・・・・・・ [Parameters] icntl = 5 maxcas = 10000 maxbch = 10 file(6) = phits.out istdev = -1 再開始計算が実行された場合、コンソール画面に メッセージが表示されます。 set: c1[30] [Source] s-type = 1 proj = proton dir = 1.0 r0 = c1 z0 = -30. z1 = -30. e0 = 150 totfact = pi*c1**2 計算が終了した各バッチの情報も出力されます。 23 バッチ計算の活用 PHITSでは、タリーの結果を各バッチが終了する度に出力さ せることができ、それを見ながらモンテカルロ計算を中断する ことができます。 itall (D=0) タリーの出力をバッチ毎に行うオプション = 0 出力なし = 1 同じファイルに出力 batch.now 1 <--- 1:continue, 0:stop 計算の途中で1を0に書き換え保存すると その次のバッチで計算をやめる。 ------------------------------------------------------------------------------bat[ 560] ncas = 560. : rijk = 151264979546685. low neutron = 0. : ncall/s = 4.000000000E+00 cpu time = 0.288 s. date = 2012-05-02 time = 15h 08m 25 24 lec03.inp 課題2: batch.nowを使って 計算を中断させてみましょう file = lec03.inp [Title] ・・・・・・ [Parameters] icntl = 5 itall = 1 maxcas = 10000 maxbch = 100 file(6) = phits.out $ istdev = -1 • maxbchを増やして長時間の計算を実行 • itallを使って途中結果を出力させる (結果が刻々と変化するのを確認する) • batch.nowを開いて中断 track_xz.eps batch.now 01 <--- 1:continue, 0:stop ------------------------------------------------------------------------------- 途中のバッチで終わっているか phits.outで確認してみましょう。 25 実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 まとめ 26 切断エネルギーとは? • 注目していない、或いは重要でないと思われる粒子の輸送や核 反応の発生をストップし、数値計算の量を減らすことが可能です。 • 初期設定では輸送しない粒子(低エネルギー中性子,光子,電子 など)を輸送するように変更できます。 粒子 ID 切断エネルギー 陽子 emin(1) 1 MeV 中性子 emin(2) 1 MeV π,μ emin(3~8) 1 MeV 電子・陽電子 emin(12,13) 109 MeV 光子 emin(14) 109 MeV 原子核 emin(15~19) 1 MeV/u 27 通常の輸送計算を実行 lec03.inp [Parameters] icntl = 5 0 itall = 1 maxcas = 10000 1000 maxbch = 100 10 file(6) = phits.out $ istdev = -1 set: c1[30] [Source] s-type = 1 proj = proton ・・・・・・ ・・・・・・ ・・・・・・ e0 = 150 輸送計算を実行する 修正して実行 [ T - C r o s s ] off mesh = reg reg = 5 r-in r-out area 102 101 1.0 103 102 1.0 104 103 1.0 105 104 1.0 106 105 1.0 e-type = 2 ne = 100 emin = 0. emax = 200. axis = eng unit = 1 output = flux file = cross.out epsout = 1 玉ねぎ構造の各層を通過する粒子 をタリーする。 28 切断エネルギーの設定 lec03.inp 50MeV以下の中性子を輸送させない。 [Parameters] icntl = 0 itall = 1 maxcas = 1000 maxbch = 10 emin(2) = 50 file(6) = phits.out $ istdev = -1 修正して実行 set: c1[30] cross.eps 50MeV未満の中性子はカット → 輸送されず各層を通過していない 29 実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 まとめ 30 断面積データライブラリとは? • 光子・電子・陽電子並びに低エネルギー中性子(20MeV以下)の 断面積は,ターゲットやエネルギーに複雑に依存するため,画一 的なモデルでは表現できない → 断面積データライブラリが必要 • 断面積ライブラリを使った計算は,計算時間が掛かる場合がある ため,デフォルトでは,それらの粒子の切断エネルギーは高く設 定されている 正しい輸送計算結果を得るため には,断面積データライブラリを 使うように[parameters]セクショ ンでeminなどを設定する必要が ある 113Cdの中性子反応断面積(JENDL-4.0より) 31 データライブラリの設定 1. データライブラリファイルの確認 c:/phits/XS/neu 中性子用核データライブラリ 2. アドレス(xsdir)ファイルの確認 c:/phits/data/xsdir.jnd アドレスファイル 3. アドレスファイルの1行目の確認 datapath=c:/phits/XS データライブラリを含むフォルダ名 4. [Parameters]セクションにおいて“emin(i)” “dmax(i)” ,“file(7)”パラメータを設定 以上はWindowsの場合です。Macの場合は、/Users/ユーザ名/phits/***の ようになりますので、***の部分をそれぞれの内容に置き換えてください。 32 課題3: 中性子用核データライブラリの利用 lec03.inp 修正して実行 [Parameters] icntl = 0 itall = 1 maxcas = 1000 中性子のデータを使う maxbch = 10 (emin < dmaxを確認) emin(2)= 1.0e-10 dmax(2) = 20 file(6) = phits.out file(7) = c:/phits/data/xsdir.jnd $ istdev = -1 アドレスファイル名 [T-Track] mesh = reg reg = (101 102 103 104 105) e-type = 2 ne = 100 emin = 0 emax = 20. axis = eng unit = 1 part = proton neutron photon alpha file = track_eng.out epsout = 1 • 各領域の合計値を出力したい場 合は,合計したい領域を()で囲む • 中性子の輸送を10-10(MeV)=0.1 (meV)まで計算する • track_eng.eps から、20MeV以下の中性子が輸送されていることを確認する 33 課題3: 中性子用核データライブラリの利用 track_eng.eps phits.out --------------------------------------------------------------CPU time and number of event called in PHITS --------------------------------------------------------------・・・・・・ ・・・・・・ ・・・・・・ dklos = 0. hydro = 255. n-data = 46556. h-data = 0. p-data = 0. e-data = 0. p-egs5 = 0. e-egs5 = 0. ・・・・・・ ・・・・・・ track_eng.epsで、中性子のスペクトルが見えた! phits.outで、核データを使ったことが確認できる 課題4: 光子・電子用データライブラリの利用 • 電子・陽電子・光子の切断エネルギーemin(12~14) を設定する →計算時間が長くなるので,ここでは1MeVに設定 • 電子・陽電子・光子のライブラリ上限エネルギーdmax(12~14) を設定する → 全て1GeVに設定 PHITSを実行 • track_eng.epsで、光子スペクトルを確認 • phits.outにあるe-data,p-dataを確認 演習 [t-track](track_xz.out)のpartを変えて電子・光 子・陽電子の飛跡を確認してみよう 注)中性子の設定emin(2)&dmax(2)は,そのまま残しておいて下さい 35 課題5: EGS5の利用 lec03.inp phits.out [Parameters] icntl = 0 itall = 1 maxcas = 1000 maxbch = 10 emin(2)= 1.0e-10 dmax(2) = 20 emin(12) = 1.0 emin(13) = 1.0 emin(14) = 1.0 dmax(12) = 1000.0 dmax(13) = 1000.0 EGS5用ライブラリ dmax(14) = 1000.0 格納フォルダ名 file(6) = phits.out file(7) = c:/phits/data/xsdir.jnd file(20) = c:/phits/XS/egs negs = 1 $ istdev = -1 光子・電子の輸送に --------------------------------------------------------------CPU time and number of event called in PHITS --------------------------------------------------------------・・・・・・ ・・・・・・ ・・・・・・ dklos = 112. hydro = 252. n-data = 43719. h-data = 0. p-data = 0. e-data = 0. p-egs5 = 1382. e-egs5 = 29835. ・・・・・・ ・・・・・・ EGS5を使う (eminからdmaxまで) 36 実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 まとめ 37 γ脱励起オプション igamma: 核反応で生成される放射性核種からのγ線放出 を考慮するためのオプション。デフォルトはγ線放出を考慮 しないので注意が必要 igamma (D=0) =0 =1 =2 =3 残留核のγ崩壊オプション γ崩壊を考慮しない γ崩壊を考慮する γ崩壊をEBITEMモデルを用いて考慮する γ崩壊とアイソマー生成をEBITEMモデルを 用いて考慮する Version 2.64以降の奨励値は igamma = 2 38 荷電粒子のビームライン設計オプション nspred & nedisp: 荷電粒子の物質中での角度・エネルギー 分散を考慮するためのオプション。加速器ビームのビーム ライン設計の際,不可欠となる nspred (D=0) =0 =1 =2 = 10 荷電粒子のクーロン散乱(角度分散)オプション クーロン散乱を考慮しない クーロン散乱をNMTCモデルを用いて考慮する クーロン散乱をLynchの式を用いて考慮する クーロン散乱をATIMAの式を用いて考慮する nedisp (D=0) =0 =1 = 10 荷電粒子のエネルギー分散オプション エネルギー分散を考慮しない エネルギー分散をLandau Vavilovの式を用いて考慮する エネルギー分散をATIMAの式を用いて考慮する ビームライン設計時の奨励値は nspred = 2 & nedisp = 1 39 核反応モデルの変更 PHITSには、Bertini, JAM, JQMD, INCL, INC-ELFという 核反応モデルがあり、状況に応じてユーザーが使い分け ることができます。 inclg (D=1) ejamnu (D=20.) JAMへの切り替えエネルギー(MeV) ejampi (D=20.) パイオン入射反応のJAMへの切り替えエネル ギー (MeV) eqmdnu (D=20.) JQMDへの切り替えエネルギー (MeV) eqmdmin (D=10.) JQMD適用の下限エネルギー (MeV/u) ejamqmd (D=3500.) 原子核反応のJQMDからJAMQMDへの切り替 えエネルギー (MeV/u) incelf (D=0) INC-ELFを使用する場合の切り替えオプション dmax(i) (D=emin(i)) i-th粒子のライブラリー利用の上限エネルギー INCLを使用する場合の切り替えオプション 40 核反応モデルの切替エネルギー (1MeV) emin(i) Nucleon (=emin) (3.0GeV) dmax(i) 核データ einclmax INCL (inclg=1) (1MeV) (3.0GeV) emin(i) Pion einclmax INCL (inclg=1) ejamqmd eqmdmin Kaon, Hyperon JAM (3.5GeV/u) (10MeV/u) Nucleus (d, t, 3He, α) JAM JQMD INCL (inclg=1) JAMQMD JAM 41 イベントジェネレーターモード 核データ+統計崩壊モデルの利用により、残留核からの 放出粒子まで考慮し、かつエネルギーと運動量の保存則 を満たした核反応イベントを生成するモード。 イベントジェネレータモードを使った方がよい場合 • 低エネルギー中性子が放出する陽子やα粒子スペクトルが見たい • 残留核の情報(反跳エネルギーなど)が知りたい • イベント毎の情報が知りたい(検出器応答関数の計算など) イベントジェネレータモードを使わない方がよい場合 • 中性子と光子の情報だけ知りたい(遮へい計算など) • 中性子の透過率を計算したい • 熱中性子の正確な挙動が見たい 42 イベントジェネレーターモードの設定 • 使用する方法は、核データを使用する設定を施した上で、 [parameters]セクションにおいて“e-mode=2”とするだけです (指定しない場合は“e-mode=0”となっています) • この場合“igamma”パラメータは自動的に2になります e-mode (D=0) =0 =1 =2 Event generator modeオプション 通常のモード Event generator mode version 1 Event generator mode version 2 43 課題6: イベントジェネレーターモードの利用 lec03.inp [Parameters] icntl = 0 itall = 1 maxcas = 1000 maxbch = 10 emin(2)= 1.0e-10 dmax(2) = 20 emin(12) = 1.0 emin(13) = 1.0 emin(14) = 1.0 dmax(12) = 1000.00000 dmax(13) = 1000.00000 dmax(14) = 1000.00000 file(6) = phits.out … 修正して実行 [Source] s-type = 1 proj = neutron dir = 1.0 r0 = c1 z0 = -30. z1 = -30. e0 = 20.0 totfact = pi*c1**2 e-mode = 0 • まずはイベントジェネレーターモードを使わずに計算してみる。(デフォルト設定) • 中性子の影響を見るため,線源は20MeV中性子とする • 球全体の粒子フルエンスを出力して,陽子・α粒子が検出されたか確認する 44 課題6: イベントジェネレーターモードの利用 lec03.inp 修正して実行 [Parameters] icntl = 0 itall = 1 maxcas = 1000 maxbch = 10 emin(2)= 1.0e-10 dmax(2) = 20 emin(12) = 1.0 emin(13) = 1.0 emin(14) = 1.0 dmax(12) = 1000.00000 dmax(13) = 1000.00000 dmax(14) = 1000.00000 file(6) = phits.out file(7) = c:/phits/data/xsdir.jnd ... e-mode = 2 track_eng.eps 陽子やα粒子のスペクトルが見えた! 45 入出力ファイルの設定 サマリーの出力ファイル名。指定しない時は、 標準出力 File(6) (D=phits.out) File(7) (D=c:/phits/data/xsdir.jnd) File(15) (D=dumpall.dat) File(18) (D=voxel.bin) Ivoxel=1,2を指定した際のバイナリデータ用 ファイル名 File(20) (D=c:/phits/XS/egs/) EGS用断面積データを格納したフォルダ名 (negs=1の際に必要) File(21) 断面積のアドレスファイル名 Dumpall=1を指定した時の出力ファイル名 (D=c:/phits/dchain-sp/data/) DCHAIN-SP用データを格納したフォルダ名 46 実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 まとめ 47 まとめ • PHITSでは、[Parameters]セクションの設定を変更することで、 様々な計算モードや物理モデルの選択を行うことができる。 • 各計算モード(icntl)を使い分けることにより、ジオメトリの確 認,ソースのチェック、粒子の輸送計算を行う。 • 統計精度と試行回数(maxcas, maxbch)は密接に関連して おり、適切な数を設定する必要がある。 • 低エネルギーの中性子や光子,電子,陽電子等は,切断エ ネルギー(emin)や核データの上限エネルギー(dmax)を変 更することで精度の高い計算結果を得ることができる。その 際,file(7)でアドレスファイルを指定する必要がある • 状況に応じてイベントジェネレーターモード(e-mode)など物理 モデルの設定を変更する必要がある。 最適な設定が分からない場合は,奨励設定ファイル(recommendation)をご参照ください 48 宿題 • 宿題ファイルで,熱中性子まで考慮し, 20MeV以下の中性子に対して核データラ イブラリを使う • イベントジェネレータモードを使ってみる • 線量-深さ分布の相対標準偏差2%以内を 目指す(maxcas, istdev, batch.nowな どを活用する) 49 宿題(解答例) 陽子(上)・中性子(下)フルエンス 円柱中心(上)と端側(下)における 付与エネルギーの深さ分布 50
© Copyright 2025 ExpyDoc