run star extras

run_star_extras.f — Fortran source code, 13 kB (14232 bytes)

File contents

! ***********************************************************************
!
!   Copyright (C) 2012  Bill Paxton
!
!   this file is part of mesa.
!
!   mesa is free software; you can redistribute it and/or modify
!   it under the terms of the gnu general library public license as published
!   by the free software foundation; either version 2 of the license, or
!   (at your option) any later version.
!
!   mesa is distributed in the hope that it will be useful, 
!   but without any warranty; without even the implied warranty of
!   merchantability or fitness for a particular purpose.  see the
!   gnu library general public license for more details.
!
!   you should have received a copy of the gnu library general public license
!   along with this software; if not, write to the free software
!   foundation, inc., 59 temple place, suite 330, boston, ma 02111-1307 usa
!
! ***********************************************************************
 
      module run_star_extras 

      use star_lib
      use star_def
      use const_def
      use const_def
      use chem_def
      use interp_1d_lib
      use interp_1d_def
      use binary_def
      
      implicit none
      
      integer :: time0, time1, clock_rate
      real(dp), parameter :: expected_runtime = 1d10 ! minutes

      integer, parameter :: restart_info_alloc = 1
      integer, parameter :: restart_info_get = 2
      integer, parameter :: restart_info_put = 3
      
      logical :: flash = .false.
      integer :: cur_model_number=1
      real(dp) :: dt_start

      real(dp) :: xfer_ratio
      
      contains

      subroutine extras_controls(s, ierr)
         type (star_info), pointer :: s
         integer, intent(out) :: ierr
         ierr = 0

         s% other_wind => rlo_other_wind

      end subroutine extras_controls
      
      
      integer function extras_startup(s, id, restart, ierr)
         type (star_info), pointer :: s
         integer, intent(in) :: id
         logical, intent(in) :: restart
         integer, intent(out) :: ierr
         integer :: restart_time, prev_time_used
         ierr = 0
         if (.not. restart) then
            call system_clock(time0,clock_rate)
            call alloc_restart_info(s)
         else
            call unpack_restart_info(s)
            call system_clock(restart_time,clock_rate)
            prev_time_used = time1 - time0
            time1 = restart_time
            time0 = time1 - prev_time_used
         end if
         extras_startup = keep_going

      end function extras_startup
      

      integer function extras_check_model(s, id, id_extra)
         type (star_info), pointer :: s
         type (binary_info), pointer :: b
         integer, intent(in) :: id, id_extra
         integer :: ierr = 0 
         real(dp) :: post_flash_mass

         call binary_ptr(b,ierr)

         s% envelope_mass_limit = -1d-1

         post_flash_mass = 0.8025d0 ! 1.0025 for 1.0 Msun WD

         if (.not. is_donor(s)) then
            if (flash) then
               flash = .false.
               cur_model_number = s% model_number
               call star_relax_mass(id, post_flash_mass, -5.1d0, ierr)
               b% angular_momentum_j = b% angular_momentum_j_old - &
               (b% m(2)-post_flash_mass*Msun)*(b% separation_old*b% m(1)/(b% m(1)+b% m(2)))**2*(2*pi/b% period)
               s% dt = 10*s% dt_next ! I don't know why this is needed, but it is
            endif
            
            if (s% L_nuc_burn_total > 1d1) then
               if (s% model_number - cur_model_number < 1000) then
                  flash = .false.
               else
               endif
            endif
         endif


         extras_check_model = keep_going
      end function extras_check_model


      integer function how_many_extra_history_columns(s, id, id_extra)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         how_many_extra_history_columns = 12
      end function how_many_extra_history_columns
      
      
      subroutine data_for_extra_history_columns(s, id, id_extra, n, names, vals, ierr)
!         use num_lib, only : safe_log10
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra, n
         character (len=maxlen_profile_column_name) :: names(n)
         real(dp) :: vals(n)
         integer, intent(out) :: ierr
         integer :: i
         real(dp) :: tb, pb, vc1, vc2, td, th, vc_div_cs, d1, d2, h, j, rhob
         ierr = 0
         tb = 0
         pb = 0
         vc_div_cs = 0
         td = 0
         th = 0
         d2 = 0
         h = 0
         vc2 = 0

         do i = 1800, s% nz, 1
            if (s% T(i) > tb) then
               tb = s% T(i)
               pb = s% P(i)
               rhob = s% Rho(i)
               th = s% Cp(i)*s% T(i)/(s% eps_nuc(i)+1d-10)
               td = s% scale_height(i)/s% csound(i)
            endif
            if (((s% conv_vel(i))/(s% csound(i))) > vc_div_cs .and. s% mixing_type(i) == 1) then
               vc_div_cs = ((s% conv_vel(i))/(s% csound(i)))
               d2 = s% D_mix(i)
               h = s% scale_height(i)
               vc2 = s% conv_vel(i)
               j = i
            endif
         end do

         vc1 = s% conv_vel(s% max_eps_nuc_k)
         d1 = s% D_mix(s% max_eps_nuc_k)

         names(1) = 'logTb'
         vals(1) = log10(tb)
         names(2) = 'logPb'
         vals(2) = log10(pb)
         names(3) = 'max_conv_vel_div_cs'
         vals(3) = vc_div_cs
         names(4) = 'heating_timescale'
         vals(4) = th
         names(5) = 'dyn_timescale'
         vals(5) = td
         names(6) = 'max_eps_nuc_vc'
         vals(6) = vc1
         names(7) = 'max_conv_vel'
         vals(7) = vc2
         names(8) = 'max_eps_nuc_D_mix'
         vals (8) = d1
         names(9) = 'max_conv_vel_D_mix'
         vals(9) = d2
         names(10) = 'max_conv_vel_H'
         vals(10) = h
         names(11) = 'max_conv_vel_zone'
         vals(11) = j
         names(12) = 'logRhob'
         vals(12) = log10(rhob)

      end subroutine data_for_extra_history_columns

      

      
      integer function how_many_extra_profile_columns(s, id, id_extra)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         how_many_extra_profile_columns = 0
      end function how_many_extra_profile_columns
      
      
      subroutine data_for_extra_profile_columns(s, id, id_extra, n, nz, names, vals, ierr)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra, n, nz
         character (len=maxlen_profile_column_name) :: names(n)
         real(dp) :: vals(nz,n)
         integer, intent(out) :: ierr
         integer :: k
         ierr = 0
      end subroutine data_for_extra_profile_columns
      

      integer function extras_finish_step(s, id, id_extra)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         extras_finish_step = keep_going
         call system_clock(time1,clock_rate)
         call store_restart_info(s)
      end function extras_finish_step
      
      
      subroutine extras_after_evolve(s, id, id_extra, ierr)
         type (star_info), pointer :: s
         integer, intent(in) :: id, id_extra
         integer, intent(out) :: ierr
         real(dp) :: dt
         ierr = 0
         dt = dble(time1 - time0) / clock_rate / 60
         if (dt > 10*expected_runtime) then
            write(*,'(/,a30,2f18.6,a,/)') '>>>>>>> EXCESSIVE runtime', &
               dt, expected_runtime, '   <<<<<<<<<  ERROR'
         else
            write(*,'(/,a50,2f18.6,99i10/)') 'runtime, retries, backups, steps', &
               dt, expected_runtime, s% num_retries, s% num_backups, s% model_number
         end if
      end subroutine extras_after_evolve
      
      
      ! routines for saving and restoring data so can do restarts

      
      subroutine alloc_restart_info(s)
         type (star_info), pointer :: s
         call move_restart_info(s,restart_info_alloc)
      end subroutine alloc_restart_info
      
      
      subroutine unpack_restart_info(s)
         type (star_info), pointer :: s
         call move_restart_info(s,restart_info_get)
      end subroutine unpack_restart_info
      
      
      subroutine store_restart_info(s)
         type (star_info), pointer :: s
         call move_restart_info(s,restart_info_put)
      end subroutine store_restart_info
      
      
      subroutine move_restart_info(s,op)
         type (star_info), pointer :: s
         integer, intent(in) :: op
         
         integer :: i, j, num_ints, num_dbls, ierr
         
         i = 0
         ! call move_int or move_flg 
         call move_int(time0)
         call move_int(time1)
         
         num_ints = i
         
         i = 0
         ! call move_dbl 
         
         num_dbls = i
         
         if (op /= restart_info_alloc) return
         if (num_ints == 0 .and. num_dbls == 0) return
         
         ierr = 0
         call star_alloc_extras(s% id, num_ints, num_dbls, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in star_alloc_extras'
            write(*,*) 'alloc_extras num_ints', num_ints
            write(*,*) 'alloc_extras num_dbls', num_dbls
            stop 1
         end if
         
         contains
         
         subroutine move_dbl(dbl)
            real(dp) :: dbl
            i = i+1
            select case (op)
            case (restart_info_get)
               dbl = s% extra_work(i)
            case (restart_info_put)
               s% extra_work(i) = dbl
            end select
         end subroutine move_dbl
         
         subroutine move_int(int)
            integer :: int
            include 'formats'
            i = i+1
            select case (op)
            case (restart_info_get)
               !write(*,3) 'restore int', i, s% extra_iwork(i)
               int = s% extra_iwork(i)
            case (restart_info_put)
               !write(*,3) 'save int', i, int
               s% extra_iwork(i) = int
            end select
         end subroutine move_int
         
         subroutine move_flg(flg)
            logical :: flg
            i = i+1
            select case (op)
            case (restart_info_get)
               flg = (s% extra_iwork(i) /= 0)
            case (restart_info_put)
               if (flg) then
                  s% extra_iwork(i) = 1
               else
                  s% extra_iwork(i) = 0
               end if
            end select
         end subroutine move_flg

      end subroutine move_restart_info
!'      
      subroutine rlo_other_wind(id, Lsurf, Msurf, Rsurf, Tsurf, w, ierr)
!         use binary_data,only: companion_mass, rl_rel_overlap_tolerance, &
!         r, rl, mass_transfer_rate, transfer_fraction, binary_separation
         use star_def
         integer, intent(in) :: id
         real(dp), intent(in) :: Lsurf, Msurf, Rsurf, Tsurf ! surface values (cgs)
         ! NOTE: surface is outermost cell. not necessarily at photosphere.
         ! NOTE: don't assume that vars are set at this point.
         ! so if you want values other than those given as args,
         ! you should use values from s% xh(:,:) and s% xa(:,:) only.
         ! rather than things like s% Teff or s% lnT(:) which have not been set yet.
         real(dp), intent(out) :: w ! wind in units of Msun/year (value is >= 0)
         integer, intent(out) :: ierr
         real(dp) :: roche_lobe_radius ! Rsun 
         real(dp) :: ratio, rho, p, grav, hp, v_th, h, rho_exponent, rho_rl, rho_rl0, mdot, xfer_ratio  
         real(dp) :: roche_lobe_xfer_full_on, roche_lobe_xfer_full_off, rlo_wind_base_mdot, rlo_wind_scale_height, rlo_wind_on
         real(dp) :: vx, rlobe2, donor_mass


         include 'formats'

         type (star_info), pointer :: s
         type (binary_info), pointer :: b

         call binary_ptr(b,ierr)

         call star_ptr(1, s, ierr)

         donor_mass = s% star_mass
         dt_start = s% time_step

         call star_ptr(2, s, ierr)

         vx = (s% star_mass/donor_mass)**one_third
         rlobe2 = b% separation*0.49d0*vx*vx/(0.6d0*vx*vx + log(1d0 + vx))
         ratio = (10**(s% log_surface_radius)*Rsun )/rlobe2
         roche_lobe_xfer_full_on = 0.05
         roche_lobe_xfer_full_off = 0.6
         xfer_ratio = 1d0

         if (ratio > roche_lobe_xfer_full_on .and. ratio < roche_lobe_xfer_full_off) then
         xfer_ratio = (roche_lobe_xfer_full_off - ratio) / &
            (roche_lobe_xfer_full_off - roche_lobe_xfer_full_on)
            xfer_ratio = 0.5d0*(1 - cos(pi*xfer_ratio))
         endif

         if (ratio > roche_lobe_xfer_full_off) then
            xfer_ratio = 0d0
         endif

         write(*,*) 'donor_mass', donor_mass, 'wd_mass', s% star_mass, 'ratio', ratio, 'xfer_ratio', xfer_ratio

         b% xfer_fraction = xfer_ratio

         w = 0
         ierr = 0

         rlo_wind_base_mdot = 1d-3  !s% x_ctrl(4)
         rlo_wind_scale_height = 1d-3  !s% x_ctrl(5)
         rlo_wind_on = 0.7  !s% x_ctrl(6)

         if (rlobe2 < 1) then
            ratio = 0d0
         endif

!         write(*,*) 'ratio', ratio, 'rlo_wind_on', rlo_wind_on

         if (ratio > rlo_wind_on) then

            w = rlo_wind_base_mdot* &
            exp((10**(s% log_surface_radius)*Rsun - rlobe2)/(Rsun*rlo_wind_scale_height))

!            if (w > 1d0) w = 1d0
!            if (s% surface_he4 > 0.5 ) w = 1d-2

!            write(*,1) 'log rlo mdot Msun/yr, R/R_L', log10(w), ratio
         endif
         if (.false. .and. (Rsurf/Rsun - roche_lobe_radius)/rlo_wind_scale_height > 100) then
            write(*,1) 'ratio', (Rsurf/Rsun - roche_lobe_radius)/rlo_wind_scale_height
            write(*,1) 'rlo_wind_scale_height', rlo_wind_scale_height
            write(*,1) 'R - roche_lobe_radius', Rsurf/Rsun - roche_lobe_radius
            write(*,1) 'R', Rsurf/Rsun
            write(*,1) 'exp(s% xh(s% i_lnR,1))/Rsun', &
               exp(s% xh(s% i_lnR,1))/Rsun
            write(*,1) 'roche_lobe_radius', roche_lobe_radius
            stop 'rlo_other_wind'
         end if


      end subroutine rlo_other_wind

      
      


      end module run_star_extras