fortran中的射击方法（中子星振荡） [英] Shooting method in fortran (neutron star oscillation)

问题描述

我一直在用fortran 90写剧本，用射击方法解决中子星的径向振荡问题。但由于不明原因，我的程序从未解决。如果没有拍摄方法组件，程序会顺利运行，因为它成功构建了恒星。但是一旦射门进入，一切都会消失。

  PROGRAM ROSCILLATION2 
 USE eos_parameters 
 IMPLICIT NONE 
 INTEGER :: i，j，k ，l 
 INTEGER，PARAMETER :: N_ode = 5 
 REAL，DIMENSION（N_OD）:: y 
 REAL（8）:: rho0_cgs，rho0，P0，r0，phi0，pi 
 REAL（8）:: r，rend，mass，P，phi，delta，xi，eta 
 REAL（8）:: step，omega，omegastep，tiny，rho_print，Radius，B，a2，s0 ，lamda，E0，E 
 EXTERNAL :: fcn 
 
 !!!!用户输入
 
 rho0_cgs = 2.D + 15！以cgs为单位的中央密度
 step = 1.D-4！步长dr 
 omegastep = 1.D-2！步长d（Ω）
 tiny = 1.D-8！小数字P（R）/ P（0）来定义星表面
 !!!!!!!!! 
 open（unit = 15，file =data.dat，status =new）
 pi = ACOS（-1.D0）
 
 a2 =（（（ （1.6022D-13）** 4）*（6.674D-11）*（（2.997D8）**  -  7）*（（1.0546D-34）**  -  3）*（1.D6））** （0.5D0））* a2_MeV！转换为编码单位（km ^ -1）
 
 B =（（1.6022D-13）** 4）*（6.674D-11）*（（2.997 D8）**  -  7）*（（1.0546D-34）**  -  3）*（1.D6）* B_MeV！转换为代码单位（km ^ -2）
 
 s0 =（ 1.D0 / 3.D0） - （1 /（6 * pi ** 2））* a2 *（（1 /（16 * pi ** 2）* a2 ** 2 +（pi **  -  2）* a1 *（rho0-B））**  -  0.5）！在r = 0时声速的平方
 
 lamda = -0.5D0 * log（1-2 * y（1）/ r ）
 
 E0 =（r0 **  -  2）* s0 * exp（lamda + 3 * phi0）
 
 
 
 rho0 = rho0_cgs * 6.67 D-18 / 9.D0！将rho0转换为代码单位（km ^ -2）
 
！计算中心压力P0 
 
 P0 =（1.D0 / 3.D0）* rho0  - （4.D0 / 3.D0）* B  - （1.D0 /（a4 *（12.D0 ）*（pi ** 2）））* a2 ** 2  - & 
&（a2 /（（3.D0）* a4））*（（（1.D0 /（16.D0 * pi ** 4））* a2 ** 2 +（1.D0 /（ pi ** 2））* a4 *（rho0-B））** 0.5D0）
 
 
 !!度量函数phi 
 phi0 = 0.1D0的初始值！任意（需要稍后调整）
 r0 = 1.D-30！整合起点
 
！设置初始条件
 !!!!!!!!!!!!!!!!! 
 !!开始整合循环
 !!!!!!!!!!!!!!!!! （3）= 1 /（3 * E0）（2）= P0 $ b $由（3）= phi0 $ b $（1）= 0.D0 $ b $ 
y（5）= 1 
 omega = 2 * pi * 1000 /（2.997D5）！1kHz以代码单位
 DO l = 1,100 
 omega = omega + omegastep！拍摄方法部分
 DO i = 1，1000000000 
 
 rend = r0 + REAL（i）* step 
 call oderk（r，rend，y，N_ode，fcn） （4）
 $ br = rend 
 mass = y（1）
 P = y（2）
 phi = y（3）
 xi = y 
 eta = y（5）
 
 IF（P  WRITE（*，*）Central density（10 ^ 14 cgs）= rho0_cgs / 1.D14 
 WRITE（*，*）质量（太阳质量）=，质量/1.477D0 
 WRITE（*，*）半径（km）=，r 
 WRITE（*，*）紧凑性M / R，质量/ r 
 WRITE（15，*）（Ω* 2.997D5 /（2 * pi）），y（5）
 
 GOTO 21 
 ENDIF 
 
 ENDDO 
 ENDDO 
 21继续
 
结束程序roscillation2 
 
！ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
 $ b $ SUBROUTINE fcn（r，y，yprime）
 USE eos_parameters 
 IMPLICIT NONE 
 REAL（8），DIMENSION（5）:: y，yprime 
 REAL（8）:: r，m，P，phi，rho，pi，B，a2，xi，eta，W，Q，E，s，lamda，omega 
 INTEGER :: j 
 $（b）= ACOS（-1.D0）
 a2 =（（（（1.6022D-13）** 4）*（6.674D-11）*（（2.997D8）**  -  7） *（（1.0546D-34）**  -  3）*（1.D6））**（0.5D0））* a2_MeV！转换为代码单位（km ^ -1）
 B =（（1.6022D -13）** 4）*（6.674D-11）*（（2.997D8）** -7）*（（1.0546D-34）** -3）*（1.D6）* B_MeV！转换为编码单位（km ^ -2）
m = y（1）
 P = y（2）
 phi = y（3）
 xi = y（4）
 eta = y（5）
 
 rho = 3.D0 * P + 4.D0 * B +（（3.D0）/（4.D0 * a4 *（pi ** 2））） * A2 ** 2 +（A2 / A4）*&安培; （（3.D0 /（pi ** 2））* a4 * b2 +（（（（D.D0）*（pi ** 4）））* a2 ** 2 + （P + B）））** 0.5D0）
 
s =（1.D0 / 3.D0） - （1 /（6 * pi ** 2））* a2 *（（1 / （16 * pi ** 2）* a2 ** 2 +（pi **  -  2）* a4 *（rho-B））**  -  0.5）！声速的平方
 
 W =（r **  -  2）*（rho + P）* exp（3 * lamda + phi）
 
 E =（r **  -  2）* s * exp（lamda + 3 * phi ）$（$）
 Q =（r **  -  2）* exp（lamda + 3 * phi）*（rho + P）*（（yprime（3）** 2）+ 4 *（r * * -1）* yprime（3）-8 * pi * P * exp（2 * lamda））
 
 yprime（1）= 4.D0 * pi * rho * r ** 2 
 $ b $ yprime（2）=  - （rho + P）*（m + 4.D0 * pi * P * r ** 3）/（r *（r-2.D0 * m））
 
 yprime（3）=（m + 4.D0 * pi * P * r ** 3）/（r *（r-2.D0 * m））
 
 yprime （4）= y（5）/（3 * E）
 
 yprime（5）=  - （W * omega ** 2 + Q）* y（4）
 
 END SUBROUTINE fcn 
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!! 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
！ 
！ Runge-Kutta方法（来自Numerical Recipes）
！ 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
 
子程序oderk（ri，re，y，n，derivs）
 INTEGER，PARAMETER :: NMAX = 16 
 REAL（8）:: ri，re，step 
 REAL（8），DIMENSION（NMAX）:: y，dydx，yout 
 EXTERNAL :: derivs，rk4 
 
调用derivs（ri，y，dydx）
 step = re-ri 
 CALL rk4（y，dydx，n，ri，step，yout，derivs）
 do i = 1，n 
y（i）= yout（i）
 enddo 
返回
结束子程序oderk 
 
 SUBROUTINE RK4（Y，DYDX，N，X，H，YOUT，DERIVS）
 INTEGER，PARAMETER :: NMAX = 16 
 REAL（8）:: H，HH，XH，X，H6 
 REAL（8），DIMENSION（N）:: Y，DYDX，YOUT 
 REAL（8），DIMENSION （NMAX）:: YT，DYT，DYM 
 EXTERNAL :: derivs 
 
 
 HH = H * 0.5D0 
 H6 = H / 6D0 
 XH = X + HH 
 
 DO I = 1，N 
 YT（I）= Y（I）+ HH * DYDX（I）
 ENDDO 
 
 CALL DERIVS（XH，YT，DYT）
 
 DO I = 1，N 
 YT（I）= Y（I）+ HH * DYT（I）
 ENDDO 
 
 CALL DERIVS（XH，YT，DYM）
 
 DO I = 1，N 
 YT（I）= Y（I）+ H * DYM（I） 
 DYM（I）= DYT（I）+ DYM（I）
 ENDDO 
 
 CALL DERIVS（X + H，YT，DYT）
 
请问我= 1，N 
 YOUT（I）= Y（I）+ H6 *（DYDX（I）+ DYT（I）+ 2 * DYM（I））
 ENDDO 
 
 
 END SUBROUTINE RK4 
  

任何答复都会很棒我真的很郁闷调试。

解决方案

你的程序因为这一行而炸毁：

< pre $ yprime（5）= - （W * omega ** 2 + Q）* y（4）


子程序 fcn 中的

。在这个子程序中， omega 完全独立于主程序中声明的那个。如果你的编译器足够好（或者被告知）初始化变量，那么这个未初始化并用于表达式，该表达式将包含随机值或零。

如果您希望主程序中的变量 omega 与您在<$ c $中使用的变量相同c> fcn ，那么你需要以某种方式将该变量传递给 fcn 。由于你设计这个程序的方式，传递它将需要修改你的所有程序来传递 omega ，以便它可以被提供给你所有的 DERIVS （这是您与 fcn 关联的虚拟参数）。

另一种方法是将 omega 放入模块中，使用该模块，您需要访问 omega ，例如在 eos_parameters 中声明它，而不是在 fcn 的作用域单元和你的主程序中声明它。

I have been writing a script in fortran 90 for solving the radial oscillation problem of a neutron star with the use of shooting method. But for unknown reason, my program never works out. Without the shooting method component, the program runs smoothly as it successfully constructed the star. But once the shooting comes in, everything dies.

   PROGRAM ROSCILLATION2
USE eos_parameters
IMPLICIT NONE
INTEGER ::i, j, k, l
INTEGER, PARAMETER :: N_ode = 5
REAL, DIMENSION(N_ode) :: y
REAL(8) :: rho0_cgs, rho0, P0, r0, phi0, pi
REAL(8) :: r, rend, mass, P, phi, delta, xi, eta
REAL(8) :: step, omega, omegastep, tiny, rho_print, Radius, B, a2, s0, lamda, E0, E
EXTERNAL :: fcn

!!!! User input 

rho0_cgs = 2.D+15           !central density in cgs unit
step = 1.D-4                ! step size dr
omegastep = 1.D-2           ! step size d(omega)
tiny = 1.D-8                ! small number P(R)/P(0) to define star surface
!!!!!!!!!
open(unit=15, file="data.dat", status="new")
pi = ACOS(-1.D0)

a2 =((((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6))**(0.5D0))*a2_MeV !convert to code unit (km^-1)

B = ((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6)*B_MeV !convert to code unit (km^-2)

s0 = (1.D0/3.D0) - (1/(6*pi**2))*a2*((1/(16*pi**2)*a2**2 + (pi**-2)*a4*(rho0 - B))**-0.5) !square of the spped of sound at r=0

lamda = -0.5D0*log(1-2*y(1)/r)

E0 = (r0**-2)*s0*exp(lamda + 3*phi0)



rho0 = rho0_cgs*6.67D-18 / 9.D0  !convert rho0 to code unit (km^-2)

!! Calculate central pressure  P0

P0 = (1.D0/3.D0)*rho0 - (4.D0/3.D0)*B - (1.D0/(a4*(12.D0)*(pi**2)))*a2**2 - &
&(a2/((3.D0)*a4))*(((1.D0/(16.D0*pi**4))*a2**2+(1.D0/(pi**2))*a4*(rho0-B))**0.5D0)


!! initial value for metric function phi 
phi0 = 0.1D0     ! arbitrary (needed to be adjusted later)
r0 = 1.D-30      ! integration starting point

!! Set initial conditions
!!!!!!!!!!!!!!!!!
!!Start integration loop
!!!!!!!!!!!!!!!!!
r = r0
y(1) = 0.D0
y(2) = P0
y(3) = phi0
y(4) = 1/(3*E0)
y(5) = 1
omega = 2*pi*1000/(2.997D5) !omega of 1kHz in code unit 
DO l = 1, 1000
    omega = omega + omegastep !shooting method part
DO i = 1, 1000000000

   rend = r0 + REAL(i)*step
   call oderk(r,rend,y,N_ode,fcn)

   r = rend
   mass = y(1)
   P = y(2)
   phi = y(3)
   xi = y(4)
   eta = y(5)

   IF (P < tiny*P0) THEN
      WRITE(*,*) "Central density (10^14 cgs) = ",  rho0_cgs/1.D14
      WRITE(*,*) " Mass (solar mass) = ", mass/1.477D0
      WRITE(*,*) " Radius (km) = ", r
      WRITE(*,*) " Compactness M/R ", mass/r
      WRITE(15,*) (omega*2.997D5/(2*pi)), y(5)

      GOTO 21
   ENDIF

ENDDO
ENDDO
21 CONTINUE

END PROGRAM roscillation2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE fcn(r,y,yprime)
USE eos_parameters
IMPLICIT NONE 
REAL(8), DIMENSION(5) :: y, yprime
REAL(8) :: r, m, P, phi, rho, pi, B, a2, xi, eta, W, Q, E, s, lamda, omega
INTEGER :: j 

pi = ACOS(-1.D0)
a2 =((((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6))**(0.5D0))*a2_MeV !convert to code unit (km^-1)
B = ((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6)*B_MeV !convert to code unit (km^-2)
m = y(1)
P = y(2)
phi = y(3)
xi = y(4)
eta = y(5)

rho = 3.D0*P + 4.D0*B +((3.D0)/(4.D0*a4*(pi**2)))*a2**2+(a2/a4)*&
&(((9.D0/((16.D0)*(pi**4)))*a2**2+((3.D0/(pi**2))*a4*(P+B)))**0.5D0)

s = (1.D0/3.D0) - (1/(6*pi**2))*a2*((1/(16*pi**2)*a2**2 + (pi**-2)*a4*(rho - B))**-0.5) !square of speed of sound

W = (r**-2)*(rho + P)*exp(3*lamda + phi)

E = (r**-2)*s*exp(lamda + 3*phi)

Q = (r**-2)*exp(lamda + 3*phi)*(rho + P)*((yprime(3)**2) + 4*(r**-1)*yprime(3)- 8*pi*P*exp(2*lamda))

yprime(1) = 4.D0*pi*rho*r**2

yprime(2) = - (rho + P)*(m + 4.D0*pi*P*r**3)/(r*(r-2.D0*m))

yprime(3) = (m + 4.D0*pi*P*r**3)/(r*(r-2.D0*m)) 

yprime(4) = y(5)/(3*E)

yprime(5) = -(W*omega**2 + Q)*y(4)

END SUBROUTINE fcn

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!! Runge-Kutta method (from Numerical Recipes)
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine oderk(ri,re,y,n,derivs) 
INTEGER, PARAMETER :: NMAX=16
REAL(8) :: ri, re, step
REAL(8), DIMENSION(NMAX) :: y, dydx, yout
EXTERNAL :: derivs,rk4 

call derivs(ri,y,dydx) 
step=re-ri 
CALL rk4(y,dydx,n,ri,step,yout,derivs) 
do i=1,n 
   y(i)=yout(i) 
enddo
return 
end subroutine oderk

SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,DERIVS) 
INTEGER, PARAMETER :: NMAX=16 
REAL(8) :: H,HH,XH,X,H6 
REAL(8), DIMENSION(N) :: Y, DYDX, YOUT 
REAL(8), DIMENSION(NMAX) :: YT, DYT, DYM 
EXTERNAL :: derivs


HH=H*0.5D0 
H6=H/6D0 
XH=X+HH

DO I=1,N 
   YT(I)=Y(I)+HH*DYDX(I) 
ENDDO

CALL DERIVS(XH,YT,DYT) 

DO I=1,N 
   YT(I)=Y(I)+HH*DYT(I) 
ENDDO

CALL DERIVS(XH,YT,DYM) 

DO I=1,N 
   YT(I)=Y(I)+H*DYM(I) 
   DYM(I)=DYT(I)+DYM(I) 
ENDDO

CALL DERIVS(X+H,YT,DYT) 

DO I=1,N 
   YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2*DYM(I)) 
ENDDO


END SUBROUTINE RK4

Any reply would be great i am just really depressed for the long debugging.

解决方案

Your program is blowing up because of this line:

yprime(5) = -(W*omega**2 + Q)*y(4)

in subroutine fcn. In this subroutine, omega is completely independent of the one declared in your main program. This one is uninitialized and used in an expression, which will either contain random values or zero, if your compiler is nice enough (or told) to initialize variables.

If you want the variable omega from your main program to be the same variable you use in fcn then you need to pass that variable to fcn somehow. Due to the way you've architected this program, passing it would require modifying all of your procedures to pass omega so that it can be provided to all of your calls to DERIVS (which is the dummy argument you are associating with fcn).

An alternative would be to put omega into a module and use that module where you need access to omega, e.g. declare it in eos_parameters instead of declaring it in the scoping units of fcn and your main program.

这篇关于fortran中的射击方法（中子星振荡）的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！