Changes between Version 1 and Version 2 of LE_Module


Ignore:
Timestamp:
06/29/13 18:35:59 (12 years ago)
Author:
Erini Lambrides
Comment:

—

Legend:

Unmodified
Added
Removed
Modified

  • LE_Module

    v1 v2  
    88== Requirements ==
    99|| poly_index ||  This is the index,n, of the polytrope equation of state  ||
    10 || rho_cent ||  The central denisty of the clump object  ||
    11 || xi_last ||  The root of the solve lame-emden equation - this represents the polytrope radius  ||
    12 || R_outer ||  The maximum domain size, needed to determine computational step size ||
     10|| rho_cent ||  The central density of the clump object (cgs) ||
     11|| M_star  ||  The mass of the object (cgs) ||
    1312
    1413If the polytrope index n that is given is greater than or equal to 5, we will get a solution with infinite radius - if the user explicitly provides a xi_last then the profile can be constructed, otherwise it is stopped.
    … …  
    1615== Output ==
    1716
    18 || K || K is the free parameter from the polytropic equation of state ||
     17|| a_scale || scaling factor from dimensionless to physical values ||
    1918|| p_cent || Central pressure of clump ||
    20 || M_star || Mass of clump ||
    21 || data(*,1:3) || profiles - data(*,2): density profile, data(*,3): pressure profile ||
     19|| t_cent  || Central temp ||
     20|| R_outer || Radius of object ||
     21|| data(*,1:3) || profiles - data(*,2): density profile, data(*,3): pressure profile, data(*,4): temp profile ||
    2222
    2323== Code ==
    24 
     24{{{
    2525Module LE_module
    26 
    2726   USE GlobalDeclarations
    28    USE PhysicalDeclarations
     27   USE PhysicsDeclarations
    2928   USE Profiles
    3029   Implicit None
    … …  
    3635CONTAINS
    3736
    38    SUBROUTINE ConstructPolyTropeProfile(Profile, R_outer, poly_index, rho_cent, M_star, xi_last)
     37   SUBROUTINE ConstructPolyTropeProfile(Profile,R_outer, poly_index, rho_cent, M_star, xi_last,a_scale,p_cent,temp_cent,rho_outer)
    3938      TYPE(ProfileDef), POINTER :: Profile
    40       INTEGER :: nPoints
    41       REAL :: R_outer, poly_index, rho_cent, xi_last
    42       REAL :: h, a_scale, K, p_cent, dxi, xi, theta, dtheta, M_star
     39      INTEGER(KIND=8) :: nPoints
     40      REAL(KIND=QPREC) :: nPointsf
     41      REAL(KIND=QPREC) :: h, poly_index,press_profile,rho_profile
     42      REAL(KIND=QPREC) :: M_star,rho_cent  ! INPUT PARAMETERS
     43      REAL(KIND=QPREC) ::  xi, xi_last, dxi, dtheta, dtheta_last, theta, theta_last,ddxi,nsteps !calculated dimensionless variables
     44      REAL(KIND=QPREC) :: a_scale, R_outer, p_cent,temp_cent !calculated scaling values and other useful relations
     45      REAL(KIND=qPREC) :: theta_target, rho_outer
     46      INTEGER :: i,j
     47      !Profile data stored in Profile object
    4348     
    44       nPoints=max(R_outer/levels(maxLevel)%dx, 1000)
    45       h=R_outer/nPoints !Computational step size
    46 
    4749      IF (poly_index >= 5) THEN
     50         IF (MPI_ID == 0) THEN
     51            write(*,*) 'not yet fully implemented for poly_index >= 5'
     52            STOP
     53         END IF
    4854         IF (xi_last <= 0) THEN
    4955            IF (MPI_ID == 0) write(*,*) 'Error - must specify dimensionless truncation cutoff xi_last if the polytropic index > 5 in ConstructPolyTropeProfile'
    … …  
    5157         END IF
    5258      ELSE
     59         
     60         theta_target=(rho_outer/rho_cent)**(1d0/poly_index)
     61
    5362         !Initilize solution at some small xi
    54          xi=1e-4
    55          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)
    56          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)
     63         xi=1d-4
     64         theta=1d0-(xi**2d0)/6d0+(poly_index/120d0)*(xi**4d0)-(8d0*poly_index**2d0-5d0*poly_index)/(15120d0)*(xi**5d0)
     65         dtheta=-xi/3d0+(poly_index/30d0)*(xi**3d0)-(8d0*poly_index**2d0-5d0*poly_index)/(2520d0)*(xi**5d0)
    5766
    5867         !Initial step size
    59          dxi=1e-2
    60 
     68         dxi=1d-2
    6169         DO
    62             theta_last=theta
    63             xi_last=xi
    64             dtheta_last=dtheta
    65             CALL LE_Advance(xi, theta, dtheta, dxi)
     70            CALL LE_Advance(xi,dtheta,poly_index,theta,dxi,1,0d0)
    6671            IF (xi > 10d0) dxi=.01*theta/dtheta !adjust step size once were beyond xi=10 in case poly_index is just under 5...
    67             IF (theta < 0d0) EXIT
     72            IF (theta <= 0d0) EXIT
     73!            write(*,*) xi, theta
    6874         END DO
    69          dtheta = dtheta - (xi-xi_last)*(dtheta-dtheta_last)/dxi
    70          xi = xi + (theta/(theta_last-theta))*dxi
    71          theta=0d0
    7275         xi_last=xi
    7376         theta_last=theta
    74       END IF
     77         dtheta_last=dtheta
     78         write(*,*) "xi_last :" ,xi_last   
     79         write(*,*) "dtheta_last: ", dtheta_last
     80       END IF
    7581
    76       !Construct Profile Object
    77       CALL CreateProfile(Profile, nPoints, (/Mass_Field, Press_Field/), RADIAL)
    78      
     82       !Calculate length scale
     83       a_scale=(-M_star/(4d0*Pi*rho_cent*(xi_last**2)*dtheta_last))**(1d0/3d0)
     84       write(*,*) "a_scale", a_scale, M_star, rho_cent, xi_last, dtheta_last
    7985
    80       !Calculate various other scales
    81       a_scale=xi_last/R_outer
    82       K=4d0*pi*ScaleGrav/((poly_index+1)*a_scale**2)*rho_cent**((poly_index-1)/poly_index)
    83       p_cent=K*rho_cent**(1d0+1d0/poly_index)
     86       !Calculate polytrope physical radius
     87       R_outer=a_scale*xi_last
    8488
    85       !Start solution at small positive xi
    86       xi=1e-4
    87       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)
     89       !Determine number of points based on resolution with a minimum of 1000
     90       nPoints=max(ceiling(R_outer/(levels(maxLevel)%dx*lScale)),999)+1
     91
     92       !Physical step size between profile points
     93       h=R_outer/real(nPoints-1)
     94       write(*,*) "step size h: ",h
     95
     96      !Construct Profile Object and initialize default values for no polytrope
     97      CALL CreateProfile(Profile, nPoints, (/Mass_Field, P_Field, Temp_Field/), RADIAL)         
     98      Profile%data(:,1)=h*(/(i-1,i=1,nPoints)/)
     99      Profile%data(:,2:4)=0d0
     100
     101      !Calculate central pressure
     102      p_cent=4d0*Pi*G*(a_scale**2)*(rho_cent**2)/(poly_index+1d0)
     103      write(*,*) "p_cent: ", p_cent
     104
     105      !Initialize center cell
     106      profile%data(1,1:3)=(/0d0, rho_cent, p_cent/)
     107     
     108
     109      !Start solution at small positive xi - but don't go past 1st point
     110      xi=min(h, 1d-4)
     111      theta=1d0-(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)
    88112      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)
    89      
    90       !Integrate to first profile zone
    91       dxi=.5*(h/a_scale)-xi
    92       nsteps=ceiling(dxi/.01)
    93       ddxi=dxi/nsteps
    94       DO j=1,nsteps
    95          CALL LE_Advance(xi, theta, dtheta, ddxi)
    96       END DO
     113      dxi=h/a_scale-xi
    97114
     115      !Advance to second profile point
     116      write(*,*), "dxi:  ",dxi
     117      CALL LE_Advance(xi,dtheta,poly_index,theta,dxi,100,theta_target)
     118      write(*,*),"xi...WORK",xi 
     119      profile%data(2,1:3)=(/xi*a_scale,rho_cent*theta**poly_index, p_cent*theta**(poly_index+1d0)/)
    98120     
    99121      !Calculate step sizes for rest of profile
    100122      dxi=h/a_scale
    101       nsteps=ceiling(h/a_scale/dxi)
    102       ddxi=dxi/nsteps
    103 
    104       data(:,1)=h*(/(i,i=1,npoints)/)-.5*h
    105       data(:,2:3)=0
     123      write(*,*) "dxi for rest of profile: ", dxi     
    106124
    107125      !Initialize profile
    108       DO i=1, size(data, 1)
    109          ! Convert from dimenionless to computational units
    110          data(i,2)=rho_cent*theta**poly_index
    111          data(i,3)=k*data(i,2)**(1d0+1d0/poly_index)
    112 
    113          DO j=1,nsteps
    114             !backup values for theta, xi, dtheta in case we overshoot
    115             theta_last=theta
    116             xi_last=xi
    117             dtheta_last=dtheta
    118            
    119             !Advance to next position
    120             CALL LE_Advance(xi, theta, dtheta, ddxi)
    121             IF (theta < 0d0) EXIT
    122          END DO
    123          IF (theta < 0d0) EXIT
     126      DO i=2, size(Profile%data, 1)
     127         CALL LE_ADvance(xi,dtheta,poly_index,theta,dxi,100,theta_target)
     128         Profile%data(i,1:3)=(/xi*a_scale,rho_cent*theta**poly_index,p_cent*theta**(poly_index+1d0)/)
     129!         Print *, "Density: ", Profile%data(i,2)
     130!         IF (theta <= theta_target) EXIT
    124131      END DO
    125 
    126       !
    127       dtheta = dtheta - (xi-xi_last)*(dtheta-dtheta_last)/ddxi
    128       xi = xi + (theta/(theta_last-theta))*ddxi
    129       theta=0d0
    130 
    131       M_star=4d0*pi*rho_cent*(R_star**3)*(-1d0/xi_last)*(dtheta)
    132 
     132      theta=theta_target
     133      Print *, xi
     134      profile%data(:,4)=profile%data(:,3)/profile%data(:,2) * Xmu*Hmass/Boltzmann
     135     
    133136      IF (MPI_ID == 0) THEN
    134137         print*,'final value for xi',xi
    135138         print*,'final value for theta',theta
    136139         print*,'final value for dtheta/dphi',dtheta
    137          
    138          print*, 'first root of phi at xi=',xi
    139          print*, 'at this point dtheta/dt = ',dtheta
    140140      END IF
     141
     142
    141143   END SUBROUTINE ConstructPolyTropeProfile
    142144
    143    SUBROUTINE LE_Advance(xi,dtheta,poly_index,theta,h)
    144       Real:: xi, dtheta, poly_index, theta, h
    145       REAL :: k1(2), k2(2),k3(2),k4(2),g(2)
     145   SUBROUTINE LE_Advance(xi,dtheta,poly_index,theta,step,nsteps,theta_target)
     146      REAL(KIND=QPREC):: xi, dtheta, poly_index, theta, h, step, theta_target
     147      REAL(KIND=QPREC) :: k1(2), k2(2),k3(2),k4(2)
     148      REAL(KIND=QPREC), Dimension(2) :: g
     149      REAL(KIND=QPREC) :: stepsize
     150      INTEGER :: i,nsteps
     151      IF (theta <= theta_target) THEN !use zero value for outside of polytrope
     152         xi=xi+step
     153         theta=theta_target
     154         RETURN
     155      END IF
     156      h=step/real(nsteps)
     157      DO i=1,nsteps
     158         g=(/theta, dtheta/)
     159         k1=h*f(poly_index,xi,g)
     160         k2=h*f(poly_index,xi+h/2.0,g+k1/2.0)
     161         k3=h*f(poly_index,xi+h/2.0,g+k2/2.0)
     162         k4=h*f(poly_index,xi+h,g+k3)
    146163
    147       g=(/theta, dtheta/)
    148       k1=h*f(poly_index,xi,g)
    149       k2=h*f(poly_index,xi+h/2.0,g+k1/2d0)
    150       k3=h*f(poly_index,xi+h/2.0,g+k2/2d0)
    151       k4=h*f(poly_index,xi+h,g+k3)
     164         g=g+(k1+2.0*k2+2.0*k3+k4)/6.0
    152165
    153       g=g+(k1+2.0*k2+2.0*k3+k4)/6.0
    154       theta=g(1)
    155       dtheta=g(2)
    156       xi=xi+h
     166         IF (g(1) <= theta_target .OR. ISNAN(g(1))) THEN !use previous values to estimate root
     167            xi = xi - ((theta-theta_target)/dtheta)
     168            theta=theta_target
     169            EXIT
     170         END IF
     171         theta=g(1)
     172         dtheta=g(2)
     173         xi=xi+h
     174      END DO
    157175   end SUBROUTINE LE_Advance
    158176
    … …  
    162180   function f(poly_index,xi,g)
    163181      implicit none
    164       real poly_index,xi,g(2), f(2)
     182      REAL(KIND=qPREC) poly_index,xi, f(2)
     183      REAL(KIND=qPREC), dimension(2) :: g
    165184      f(1)=g(2)
    166185      f(2)=-g(1)**poly_index-2.0*g(2)/xi
    167186   end function f
    168 
    169 
    170 
    171 
    172 
    173 
    174 
    175 
     187END MODULE LE_module }}}
     188