# Blog: Meeting Update June 10th - Erini: LE_module_update.f90

File LE_module_update.f90, 4.8 KB (added by , 11 years ago) |
---|

Line | |
---|---|

1 | Module LE_module |

2 | |

3 | USE GlobalDeclarations |

4 | USE PhysicalDeclarations |

5 | USE Profiles |

6 | Implicit None |

7 | |

8 | Private |

9 | |

10 | PUBLIC ConstructPolyTropeProfile |

11 | |

12 | CONTAINS |

13 | |

14 | SUBROUTINE ConstructPolyTropeProfile(Profile, Mass_Profile,P_profile,rho_profile,R_outer, poly_index, rho_cent, M_star, xi_last,a_scale,p_cent) |

15 | TYPE(ProfileDef), POINTER :: Profile |

16 | INTEGER :: nPoints |

17 | REAL :: R_outer, poly_index, rho_cent, xi_last |

18 | REAL :: h, a_scale, K, p_cent, dxi, xi, theta, dtheta, M_star |

19 | PARAMETER(pi=3.14159265359) |

20 | |

21 | a_scale=(-M_star/(4*pi*rho_cent*(xi_last**2)*dtheta_last))**(1/3) |

22 | R_outer=a_scale*xi_last |

23 | nPoints=max(R_outer/levels(maxLevel)%dx, 1000) |

24 | h=R_outer/nPoints !Computational step size |

25 | |

26 | IF (poly_index >= 5) THEN |

27 | IF (xi_last <= 0) THEN |

28 | IF (MPI_ID == 0) write(*,*) 'Error - must specify dimensionless truncation cutoff xi_last if the polytropic index > 5 in ConstructPolyTropeProfile' |

29 | STOP |

30 | END IF |

31 | ELSE |

32 | !Initilize solution at some small xi |

33 | xi=1e-4 |

34 | theta=1-(xi**2.0)/6.0+(poly_index/120.0)*(xi**4.0)-(8.0*poly_index**2.0-5.0*poly_index)/(15120.0)*(xi**5.0) |

35 | dtheta=-xi/3.0+(poly_index/30.0)*(xi**3.0)-(8.0*poly_index**2.0-5.0*poly_index)/(2520.0)*(xi**5.0) |

36 | |

37 | !Initial step size |

38 | dxi=1e-2 |

39 | |

40 | DO |

41 | theta_last=theta |

42 | xi_last=xi |

43 | dtheta_last=dtheta |

44 | CALL LE_Advance(xi, theta, dtheta, dxi) |

45 | IF (xi > 10d0) dxi=.01*theta/dtheta !adjust step size once were beyond xi=10 in case poly_index is just under 5... |

46 | IF (theta < 0d0) EXIT |

47 | END DO |

48 | dtheta = dtheta - (xi-xi_last)*(dtheta-dtheta_last)/dxi |

49 | xi = xi + (theta/(theta_last-theta))*dxi |

50 | theta=0d0 |

51 | xi_last=xi |

52 | theta_last=theta |

53 | dtheta_last=dtheta |

54 | END IF |

55 | |

56 | !Construct Profile Object |

57 | CALL CreateProfile(Profile, nPoints, (/Mass_Field, Press_Field/), RADIAL) |

58 | |

59 | |

60 | !Calculate various other scales |

61 | a_scale=(-M_star/(4*pi*rho_cent*(xi_last**2)*dtheta_last))**(1/3) |

62 | p_cent=4*pi*ScaleGrav*(a_scale**2)*(rho_cent**2)*(theta_last**(poly_index+1d0))/(poly_index+1) |

63 | |

64 | !Start solution at small positive xi |

65 | xi=1e-4 |

66 | theta=1-(xi**2.0)/6.0+(poly_index/120.0)*(xi**4.0)-(8.0*poly_index**2.0-5.0*poly_index)/(15120.0)*(xi**5.0) |

67 | dtheta=-xi/3.0+(poly_index/30.0)*(xi**3.0)-(8.0*poly_index**2.0-5.0*poly_index)/(2520.0)*(xi**5.0) |

68 | |

69 | !Integrate to first profile zone |

70 | dxi=.5*(h/a_scale)-xi |

71 | nsteps=ceiling(dxi/.01) |

72 | ddxi=dxi/nsteps |

73 | DO j=1,nsteps |

74 | CALL LE_Advance(xi, theta, dtheta, ddxi) |

75 | END DO |

76 | |

77 | |

78 | !Calculate step sizes for rest of profile |

79 | dxi=h/a_scale |

80 | nsteps=ceiling(h/a_scale/dxi) |

81 | ddxi=dxi/nsteps |

82 | |

83 | data(:,1)=h*(/(i,i=1,npoints)/)-.5*h |

84 | data(:,2:3:4)=0 |

85 | |

86 | !Initialize profile |

87 | DO i=1, size(data, 1) |

88 | ! Convert from dimenionless to computational units |

89 | data(i,2)=rho_cent*theta**poly_index |

90 | data(i,3)=4*pi*ScaleGrav*(a_scale**2)*(rho_cent**2)*(theta**(poly_index+1d0))/(poly_index+1) |

91 | data(i,4)=-4*pi*(a_scale**3)*(xi**2)*rho_cent*dtheta/(Msun) |

92 | DO j=1,nsteps |

93 | !backup values for theta, xi, dtheta in case we overshoot |

94 | theta_last=theta |

95 | xi_last=xi |

96 | dtheta_last=dtheta |

97 | |

98 | !Advance to next position |

99 | CALL LE_Advance(xi, theta, dtheta, ddxi) |

100 | IF (theta < 0d0) EXIT |

101 | END DO |

102 | IF (theta < 0d0) EXIT |

103 | END DO |

104 | |

105 | ! |

106 | dtheta = dtheta - (xi-xi_last)*(dtheta-dtheta_last)/ddxi |

107 | xi = xi + (theta/(theta_last-theta))*ddxi |

108 | theta=0d0 |

109 | |

110 | M_star=4d0*pi*rho_cent*(R_star**3)*(-1d0/xi_last)*(dtheta) |

111 | |

112 | IF (MPI_ID == 0) THEN |

113 | print*,'final value for xi',xi |

114 | print*,'final value for theta',theta |

115 | print*,'final value for dtheta/dphi',dtheta |

116 | |

117 | print*, 'first root of phi at xi=',xi |

118 | print*, 'at this point dtheta/dt = ',dtheta |

119 | END IF |

120 | END SUBROUTINE ConstructPolyTropeProfile |

121 | |

122 | SUBROUTINE LE_Advance(xi,dtheta,poly_index,theta,h) |

123 | Real:: xi, dtheta, poly_index, theta, h |

124 | REAL :: k1(2), k2(2),k3(2),k4(2),g(2) |

125 | |

126 | g=(/theta, dtheta/) |

127 | k1=h*f(poly_index,xi,g) |

128 | k2=h*f(poly_index,xi+h/2.0,g+k1/2d0) |

129 | k3=h*f(poly_index,xi+h/2.0,g+k2/2d0) |

130 | k4=h*f(poly_index,xi+h,g+k3) |

131 | |

132 | g=g+(k1+2.0*k2+2.0*k3+k4)/6.0 |

133 | theta=g(1) |

134 | dtheta=g(2) |

135 | xi=xi+h |

136 | end SUBROUTINE LE_Advance |

137 | |

138 | |

139 | !function f calculate the second order dtheta |

140 | |

141 | function f(poly_index,xi,g) |

142 | implicit none |

143 | real poly_index,xi,g(2), f(2) |

144 | f(1)=g(2) |

145 | f(2)=-g(1)**poly_index-2.0*g(2)/xi |

146 | end function f |