chiba-aplab / cansplus

3D MHD simulation code based on HLLD approximate Riemann solver and MP5 interpolation
http://www.astro.phys.s.chiba-u.ac.jp/cans/doc/index.html
14 stars 1 forks source link

phiの時間更新 #2

Open ymatumot opened 5 years ago

ymatumot commented 5 years ago

(reported by asahina)

integrate__TVDRK3内のphiの時間更新に関して、現在は phi(i,j,k) = k1phi1(i,j,k)+k2(+phi(i,j,k) & +dtodx(fphix(i-1,j,k)-fphix(i,j,k)) & +dtody(fphiy(i,j-1,k)-fphiy(i,j,k)) & +dtodz(fphiz(i,j,k-1)-fphiz(i,j,k)) & )exp(-dtch2/cp2) となっていますが、exp(..)は全体にかかって、 phi(i,j,k) = (k1phi1(i,j,k)+k2(+phi(i,j,k) & +dtodx(fphix(i-1,j,k)-fphix(i,j,k)) & +dtody(fphiy(i,j-1,k)-fphiy(i,j,k)) & +dtodz(fphiz(i,j,k-1)-fphiz(i,j,k)) & ))exp(-dtch2/cp2) ではないでしょうか?

phi1はfluxがなければdt後、減衰してphi1exp(-dtch^2/cp^2)となると思うのですが、 実際に計算してみたら、そうはなっていなかったので疑問に思いました。

以下、fluxがゼロの時の時間発展の計算です。 phi(n=1) = 0 + phi1exp phi(n=2) = 3/4phi1 + 1/4phi(n=1)exp = 3/4phi1 + 1/4phi1exp^2 phi(n=3) = 1/3phi1 + 2/3phi(n=2)exp = 1/3phi1 + 1/2phi1exp + 1/6phi1*exp^3

全体にexpをかけた場合、 phi(n=1) = 0 + phi1exp phi(n=2) = [3/4phi1 + 1/4phi(n=1)]exp = 3/4phi1exp + 1/4phi1exp^2 phi(n=3) = [1/3phi1 + 2/3phi(n=2)]exp = [1/3phi1 + 1/2phi1exp + 1/6phi1exp^2]*exp

計算してみるとどちらもphi1*expになっていませんでした。