SAD部屋 - Ticket-45 Diff

  • Added parts are displayed like this.
  • Deleted parts are displayed like this.

! PhotonListにNaNが入る

:Priority:Normal
:Reporter:amorita
:Status:New
:Assigned to:?
:Version:5322
:Milestone:?
:Created:2019-02-28

!! Description

'''PHTONS'''・'''RAD'''フラグ付きのTrackParticleが生成するPhotonListに'''NaN'''が含まれるケースが存在するとの報告が菊池氏からよせられた

!! Attachement files
{{attach_map('Ticket-45')}}
{{attach_form('','')}}

!! Changelog
!!!amorita (2019-02-28 (木) 14:24:19)
* NaNが入る位置は、NPARAやSeedRandomに依存していることから、PRNGの内部状態に依存した現象である
* 適切に乱数種を調整すると、NPARA=1の状態で、PhotonListの最初の要素にNaNが入る状態を選べる
* src/tbrad.fの調査から当該事例では
** tphotonstore関数の引数 dpgx, dpgy, dpgzが NaNで汚染されている
** tphotonconv関数の変数 dpzが NaNで汚染されている
*** dpz = dp * sqrt( 1 - dpx^2 - dpy^2 )
*** 補足事例 dp = 1.4493605746231978E-021, dpx = -8.8016648617202851, dpy = -2.1949146672237063
**** dpzがNaNになる直接の原因は、sqrt(1 - dpx^2 - dpy^2)の domain errorである
**** dpx, dpyは、tphotonconv関数の引数
!!!amorita (2019-02-28 (木) 14:34:02)
* tphotonconv関数に渡す dpx = pxi - thx1, dpy = pyi - thy1
** 当該事例では、thx1, thy1のスケールが他のサンプルから大きく外れる(O(1))
** thx1, thy1のsource termのうちO(1)異常値を示しているのは、thu, thvである
** thu, thvは、SYNRADCLからの戻り値
!!!amorita (2019-02-28 (木) 14:44:09)
* 当該事例では、SYNRADCL関数の変数'''TH'''が O(1)の異常値になっている
** thu = TH * COS(PHI), thv = TH * SIN(PHI)
** PHIは、[0, 2π)の乱数(境界は多少怪しいが、SIN, COSの中なので結論は同じ)
!!!amorita (2019-02-28 (木) 14:52:42)
* 当該事例で変数'''TH'''は、下記のフローに依存する
P1 = tran()
V=2*SQRT(2D0)*SIN(ASIN(5*P1/(4*SQRT(2D0)))/3)
X=1.5D0*UPSILON*U*V**3
EG=X/(1+X)*E0
PHI=2*PI*tran()
TH=SQRT(MAX(0D0,1-V**2))/(V*GAM)
* '''TH'''は、'''V'''→0の極限で発散する
* '''P1'''→0の極限で、'''V'''→0
!!!amorita (2019-02-28 (木) 14:53:29)
* 当該事例では、'''P1'''に'''9.6571166068315506e-6'''を拾っている
!!!amorita (2019-02-28 (木) 15:02:35)
以上から、光子放出モデル(物理モデル)の不具合と思われる
当該モデルの原典・来歴の調査が必要
!!!amorita (2019-03-04 (月) 17:59:52)
k64で再現しない原因は、'''tphotonconv'''関数内の式'''dpz = dp * sqrt( 1 - dpx^2 - dpy^2 )'''中の'''sqrt'''が、実質的に'''sqrt(max(1e-20, 1 - dpx^2 - dpy^2 ))'''に置き換わっており、dpzが'''NaN'''に汚染されないため

k64コードの問題点は二つ
* '''dpx^2 + dpy^2 > 1''' を実質的に '''dpx^2 + dpy^2 == 1'''と解釈する物理的な妥当性
** 当該事象では、演算誤差を越えて '''dpx^2 + dpy^2 >> 1'''なので、物理モデルの意味論が問われる
* 下限値 '''1e-20'''(dpz >= 1e-10 * dp)の妥当性


{{its_edit_ticket_form}}