You are here: Home / Inlists of Published Results / AM Canum Venaticorum Progenitors with Low Mass Helium Star Donors and the resultant Explosions

AM Canum Venaticorum Progenitors with Low Mass Helium Star Donors and the resultant Explosions

Inlists, models, and run_star_extras.f for MESA simulations in Brooks, Bildsten, Marchant, & Paxton 2015, ApJ, 807, 74

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