diff --git a/star/private/phase_separation.f90 b/star/private/phase_separation.f90 index 09b25374a..05bb9ae99 100644 --- a/star/private/phase_separation.f90 +++ b/star/private/phase_separation.f90 @@ -20,18 +20,17 @@ module phase_separation use star_private_def - use const_def, only: dp + use const_def implicit none - private - public :: do_phase_separation - logical, parameter :: dbg = .false. ! offset to higher phase than 0.5 to avoid interference ! between phase separation mixing and latent heat for Skye. real(dp), parameter :: eos_phase_boundary = 0.9d0 + private + public :: do_phase_separation contains @@ -40,10 +39,13 @@ subroutine do_phase_separation(s, dt, ierr) real(dp), intent(in) :: dt integer, intent(out) :: ierr + ! 'CO' or 'ONe' will implement 2-species phase separation, for 'ONe' 22Ne is included if(s% phase_separation_option == 'CO') then call do_2component_phase_separation(s, dt, 'CO', ierr) else if(s% phase_separation_option == 'ONe') then call do_2component_phase_separation(s, dt, 'ONe', ierr) + else if(s% phase_separation_option == '3c') then + call do_2component_phase_separation(s, dt, '3c', ierr) else write(*,*) 'invalid phase_separation_option' stop @@ -51,30 +53,34 @@ subroutine do_phase_separation(s, dt, ierr) end subroutine do_phase_separation subroutine do_2component_phase_separation(s, dt, components, ierr) - use chem_def, only: ic12, io16, ine20 + use chem_def, only: chem_isos, ic12, io16, ine20, ine22, ina23, img24 use chem_lib, only: chem_get_iso_id type (star_info), pointer :: s real(dp), intent(in) :: dt character (len=*), intent(in) :: components integer, intent(out) :: ierr - - real(dp) :: XNe, XO, XC, pad - integer :: k, k_bound, kstart, net_ic12, net_io16, net_ine20 + real(dp) :: XNe20, XNe22, XO, XC, XNa, XMg , pad + integer :: k, k_bound, kstart, net_ic12, net_io16, net_ine20, net_ine22, net_ina23, net_img24 logical :: save_Skye_use_ion_offsets ! Set phase separation mixing mass negative at beginning of phase separation s% phase_sep_mixing_mass = -1d0 s% eps_phase_separation(1:s%nz) = 0d0 - - if(s% phase(s% nz) < eos_phase_boundary) then - s% crystal_core_boundary_mass = 0d0 - return + if(s% phase(s% nz) < eos_phase_boundary) then !!! prevent to move the core size inwards if the core is suddently "melted" leaving everything liquid under phi<0.9 + if (s% crystal_core_boundary_mass>0d0)then + s% crystal_core_boundary_mass=s% crystal_core_boundary_mass + return + else + s% crystal_core_boundary_mass = 0d0 + return + end if end if - net_ic12 = s% net_iso(ic12) net_io16 = s% net_iso(io16) net_ine20 = s% net_iso(ine20) - + net_ine22 = s% net_iso(ine22) + net_ina23 = s% net_iso(ina23) + net_img24 = s% net_iso(img24) ! Find zone of phase transition from liquid to solid k_bound = -1 do k = s%nz,1,-1 @@ -83,18 +89,14 @@ subroutine do_2component_phase_separation(s, dt, components, ierr) exit end if end do - XC = s% xa(net_ic12,k_bound) XO = s% xa(net_io16,k_bound) - XNe = s% xa(net_ine20,k_bound) - ! Check that we're still in C/O or O/Ne dominated material as appropriate, - ! otherwise skip phase separation - if(components == 'CO'.and. XO + XC < 0.9d0) return - if(components == 'ONe'.and. XNe + XO < 0.8d0) return ! O/Ne mixtures tend to have more byproducts of burning mixed in - + XNe20 = s% xa(net_ine20,k_bound) + XNe22 = s% xa(net_ine22,k_bound) + XNa = s% xa(net_ina23,k_bound) + XMg = s% xa(net_img24,k_bound) ! If there is a phase transition, reset the composition at the boundary if(k_bound > 0) then - ! core boundary needs to be padded by a minimal amount (less than a zone worth of mass) ! to account for loss of precision during remeshing. pad = s% min_dq * s% m(1) * 0.5d0 @@ -104,7 +106,6 @@ subroutine do_2component_phase_separation(s, dt, components, ierr) exit end if end do - ! calculate energy associated with phase separation, ignoring the ionization ! energy term that Skye sometimes calculates save_Skye_use_ion_offsets = s% eos_rq% Skye_use_ion_offsets @@ -113,137 +114,216 @@ subroutine do_2component_phase_separation(s, dt, components, ierr) do k=1,s% nz s% eps_phase_separation(k) = s% energy(k) end do - ! loop runs outward starting at previous crystallization boundary do k = kstart,1,-1 ! Start by checking if this material should be crystallizing if(s% phase(k) <= eos_phase_boundary) then - s% crystal_core_boundary_mass = s% m(k+1) - exit + if (s% crystal_core_boundary_mass>s% m(k+1)) then + s% crystal_core_boundary_mass=s% crystal_core_boundary_mass + exit + else + s% crystal_core_boundary_mass = s% m(k+1) + exit + end if end if - call move_one_zone(s,k,components) ! crystallized out to k now, liquid starts at k-1. ! now mix the liquid material outward until stably stratified call mix_outward(s, k-1, 0) - end do - call update_model_(s,1,s%nz,.false.) - - ! phase separation heating term for use by energy equation do k=1,s% nz s% eps_phase_separation(k) = (s% eps_phase_separation(k) - s% energy(k)) / dt end do s% eos_rq% Skye_use_ion_offsets = save_Skye_use_ion_offsets s% need_to_setvars = .true. end if - ierr = 0 end subroutine do_2component_phase_separation subroutine move_one_zone(s,k,components) - use chem_def, only: ic12, io16, ine20 + use chem_def, only: chem_isos, ic12, io16, ine20, ine22, ina23, img24 use chem_lib, only: chem_get_iso_id type(star_info), pointer :: s integer, intent(in) :: k + real(dp), dimension(2) :: dXNe + real(dp), dimension(4) :: Dd + real(dp) :: dx1_ character (len=*), intent(in) :: components - - real(dp) :: XC, XO, XNe, XC1, XO1, XNe1, dXO, dXNe, Xfac - integer :: net_ic12, net_io16, net_ine20 + real(dp) :: XC, XO, XNe20, XNe22, XNa, XMg, XC1, XO1, XNe120, XNe122, XNa1, XMg1, dXO, Xfac + integer :: net_ic12, net_io16, net_ine20, net_ine22, net_ina23, net_img24 net_ic12 = s% net_iso(ic12) net_io16 = s% net_iso(io16) net_ine20 = s% net_iso(ine20) - - if(components == 'CO') then - XO = s% xa(net_io16,k) - XC = s% xa(net_ic12,k) - + net_ine22 = s% net_iso(ine22) + net_ina23 = s% net_iso(ina23) + net_img24 = s% net_iso(img24) + XO = s% xa(net_io16,k) + XC = s% xa(net_ic12,k) + XNe20 = s% xa(net_ine20,k) + XNe22 = s% xa(net_ine22,k) + XNa = s% xa(net_ina23,k) + XMg = s% xa(net_img24,k) + if (components /= '3c') then + if(XO + XC > 0.7d0 .and. XC > XNe20 + XNe22) then ! Call Blouin phase diagram. ! Need to rescale temporarily because phase diagram assumes XO + XC = 1 Xfac = XO + XC XO = XO/Xfac XC = XC/Xfac - dXO = blouin_delta_xo(XO) - s% xa(net_io16,k) = Xfac*(XO + dXO) s% xa(net_ic12,k) = Xfac*(XC - dXO) - ! Redistribute change in C,O into zone k-1, ! conserving total mass of C,O XC1 = s% xa(net_ic12,k-1) XO1 = s% xa(net_io16,k-1) s% xa(net_ic12,k-1) = XC1 + Xfac*dXO * s% dq(k) / s% dq(k-1) s% xa(net_io16,k-1) = XO1 - Xfac*dXO * s% dq(k) / s% dq(k-1) - else if(components == 'ONe') then - XNe = s% xa(net_ine20,k) - XO = s% xa(net_io16,k) - + !write(*,*) 'phase CO',XO,XC,XNe20+XNe22 + else if(XO + XNe20 + XNe22> 0.7d0 .and. XNe20 + XNe22 > XC) then ! Call Blouin phase diagram. ! Need to rescale temporarily because phase diagram assumes XO + XNe = 1 - Xfac = XO + XNe + Xfac = XO + XNe20 + XNe22 XO = XO/Xfac - XNe = XNe/Xfac - - dXNe = blouin_delta_xne(XNe) - - s% xa(net_ine20,k) = Xfac*(XNe + dXNe) - s% xa(net_io16,k) = Xfac*(XO - dXNe) - + XNe20 = XNe20/Xfac + XNe22 = XNe22/Xfac + dXNe = blouin_delta_xne(XNe20,XNe22) + ! write(*,*) 'dXNe', dXNe + s% xa(net_ine20,k) = Xfac*(XNe20 + dXNe(1)) + s% xa(net_ine22,k) = Xfac*(XNe22 + dXNe(2)) + s% xa(net_io16,k) = Xfac*(XO - sum(dXNe)) ! Redistribute change in Ne,O into zone k-1, ! conserving total mass of Ne,O XO1 = s% xa(net_io16,k-1) - XNe1 = s% xa(net_ine20,k-1) - s% xa(net_io16,k-1) = XO1 + Xfac*dXNe * s% dq(k) / s% dq(k-1) - s% xa(net_ine20,k-1) = XNe1 - Xfac*dXNe * s% dq(k) / s% dq(k-1) - else - write(*,*) 'invalid components option in phase separation' - stop + XNe120 = s% xa(net_ine20,k-1) + XNe122 = s% xa(net_ine22,k-1) + s% xa(net_io16,k-1) = XO1 + Xfac*sum(dXNe) * s% dq(k) / s% dq(k-1) + s% xa(net_ine20,k-1) = XNe120 - Xfac*dXNe(1) * s% dq(k) / s% dq(k-1) + s% xa(net_ine22,k-1) = XNe122 - Xfac*dXNe(2) * s% dq(k) / s% dq(k-1) + !write(*,*) 'phase ONe',XO,XC,XNe20+XNe22 + end if + else if (components == '3c') then + ! check the abundances to decide which table use for interpolation + if (XO + XC + XNe20 + XNe22 > 0.7d0 .and. XC > XMg .and. XC > XNa) then + Xfac = XO + XC + XNe20 + XNe22 + XO = XO/Xfac + XC = XC/Xfac + XNe20 = XNe20/Xfac + XNe22 = XNe22/Xfac + ! call the deltas resulting from interpolation (in mass fraction) + call medin_cumming_3p_d_cone(XC,XO,XNe20,XNe22,Dd) + ! apply fractionation as given by the deltas from interpolation + s% xa(net_ic12,k) = Xfac*(XC + Dd(1)) + s% xa(net_io16,k) = Xfac*(XO + Dd(2)) + s% xa(net_ine20,k) = Xfac*(XNe20 + Dd(3)) + s% xa(net_ine22,k) = Xfac*(XNe22 + Dd(4)) + XC1 = s% xa(net_ic12,k-1) + XO1 = s% xa(net_io16,k-1) + XNe120 = s% xa(net_ine20,k-1) + XNe122 = s% xa(net_ine22,k-1) + s% xa(net_ic12,k-1) = XC1 - Xfac*Dd(1) * s% dq(k) / s% dq(k-1) + s% xa(net_io16,k-1) = XO1 - Xfac*Dd(2) * s% dq(k) / s% dq(k-1) + s% xa(net_ine20,k-1) = XNe120 - Xfac*(Dd(3)) * s% dq(k) / s% dq(k-1) + s% xa(net_ine22,k-1) = XNe122 - Xfac*(Dd(4)) * s% dq(k) / s% dq(k-1) + ! write(*,*) 'phase 3 CONe abundances',XC,XO,XNe20+XNe22 + else if (XO + XNe20 + XNe22 + XMg > 0.7d0 .and. XMg > XC .and. XMg > XNa) then + Xfac = XO + XNe20 + XNe22 + XMg + XMg = XMg/Xfac + XO = XO/Xfac + XNe20 = XNe20/Xfac + XNe22 = XNe22/Xfac + ! call the deltas resulting from interpolation (in mass fraction) + call medin_cumming_3p_d_neomg(XMg,XO,XNe20,XNe22,Dd) + ! apply fractionation as given by the deltas from interpolation + s% xa(net_img24,k) = Xfac*(XMg + Dd(1)) + s% xa(net_io16,k) = Xfac*(XO + Dd(2)) + s% xa(net_ine20,k) = Xfac*(XNe20 + Dd(3)) + s% xa(net_ine22,k) = Xfac*(XNe22 + Dd(4)) + XMg1 = s% xa(net_img24,k-1) + XO1 = s% xa(net_io16,k-1) + XNe120 = s% xa(net_ine20,k-1) + XNe122 = s% xa(net_ine22,k-1) + s% xa(net_img24,k-1) = XMg1 - Xfac*Dd(1) * s% dq(k) / s% dq(k-1) + s% xa(net_io16,k-1) = XO1 - Xfac*Dd(2) * s% dq(k) / s% dq(k-1) + s% xa(net_ine20,k-1) = XNe120 - Xfac*Dd(3) * s% dq(k) / s% dq(k-1) + s% xa(net_ine22,k-1) = XNe122 - Xfac*Dd(4) * s% dq(k) / s% dq(k-1) + ! write(*,*) 'phase 3 ONeMg abundances',XO,XNe20+XNe22,XMg + else if (XO + XNe20 + XNe22 + XNa > 0.7d0 .and. XNa > XC .and. XNa > XMg) then + Xfac = XO + XNe20 + XNe22 + XNa + XNa = XNa/Xfac + XO = XO/Xfac + XNe20 = XNe20/Xfac + XNe22 = XNe22/Xfac + ! call the deltas resulting from interpolation (in mass fraction) + call medin_cumming_3p_d_onena(XNa,XO,XNe20,XNe22,Dd) + ! apply fractionation as given by the deltas from interpolation + s% xa(net_ina23,k) = Xfac*(XNa + Dd(1)) + s% xa(net_io16,k) = Xfac*(XO + Dd(2)) + s% xa(net_ine20,k) = Xfac*(XNe20 + Dd(3)) + s% xa(net_ine22,k) = Xfac*(XNe22 + Dd(4)) + XNa1 = s% xa(net_ina23,k-1) + XO1 = s% xa(net_io16,k-1) + XNe120 = s% xa(net_ine20,k-1) + XNe122 = s% xa(net_ine22,k-1) + s% xa(net_ina23,k-1) = XNa1 - Xfac*Dd(1) * s% dq(k) / s% dq(k-1) + s% xa(net_io16,k-1) = XO1 - Xfac*Dd(2) * s% dq(k) / s% dq(k-1) + s% xa(net_ine20,k-1) = XNe120 - Xfac*Dd(3) * s% dq(k) / s% dq(k-1) + s% xa(net_ine22,k-1) = XNe122 - Xfac*Dd(4) * s% dq(k) / s% dq(k-1) + ! write(*,*) 'phase 3 ONeNa abundances',XO,XNe20+XNe22,XNa + else if (XC + XO + XMg > 0.7d0 .and. XMg > XNa .and. XMg > XNe20+XNe22) then + Xfac = XC + XO + XMg + XC = XC/Xfac + XO = XO/Xfac + XMg = XMg/Xfac + ! call the deltas resulting from interpolation (in mass fraction) + call medin_cumming_3p_d_comg(XC,XMg,XO,Dd) + ! apply fractionation as given by the deltas from interpolation + s% xa(net_ic12,k) = Xfac*(XC + Dd(1)) + s% xa(net_img24,k) = Xfac*(XMg + Dd(2)) + s% xa(net_io16,k) = Xfac*(XO - (Dd(1) + Dd(2))) + XC1 = s% xa(net_ic12,k-1) + XO1 = s% xa(net_io16,k-1) + XMg1 = s% xa(net_img24,k-1) + s% xa(net_ic12,k-1) = XC1 - Xfac*Dd(1) * s% dq(k) / s% dq(k-1) + s% xa(net_img24,k-1) = XMg1 - Xfac*Dd(2) * s% dq(k) / s% dq(k-1) + s% xa(net_io16,k-1) = XO1 + Xfac*(Dd(1)+Dd(2)) * s% dq(k) / s% dq(k-1) + ! write(*,*) 'phase 3 COMg abundances',XC,XO,XMg + end if end if - call update_model_(s,k-1,s%nz,.true.) - end subroutine move_one_zone ! mix composition outward until reaching stable composition profile subroutine mix_outward(s,kbot,min_mix_zones) type(star_info), pointer :: s integer, intent(in) :: kbot, min_mix_zones - real(dp) :: avg_xa(s%species) real(dp) :: mass, B_term, grada, gradr integer :: k, l, ktop logical :: use_brunt use_brunt = s% phase_separation_mixing_use_brunt - do k=kbot-min_mix_zones,1,-1 ktop = k - if (s% m(ktop) > s% phase_sep_mixing_mass) then s% phase_sep_mixing_mass = s% m(ktop) end if - mass = SUM(s%dm(ktop:kbot)) do l = 1, s%species avg_xa(l) = SUM(s%dm(ktop:kbot)*s%xa(l,ktop:kbot))/mass end do - ! some potential safeguards from conv_premix ! avg_xa = MAX(MIN(avg_xa, 1._dp), 0._dp) ! avg_xa = avg_xa/SUM(avg_xa) - do l = 1, s%species s%xa(l,ktop:kbot) = avg_xa(l) end do - ! updates, eos, opacities, mu, etc now that abundances have changed, ! but only in the cells near the boundary where we need to check here. ! Will call full update over mixed region after exiting loop. call update_model_(s, ktop-1, ktop+1, use_brunt) - if(use_brunt) then B_term = s% unsmoothed_brunt_B(ktop) grada = s% grada_face(ktop) @@ -252,36 +332,31 @@ subroutine mix_outward(s,kbot,min_mix_zones) ! stable against further mixing, so exit loop exit end if - else ! simpler calculation based on mu gradient + else ! simpler calculation based on mu gradient if(s% mu(ktop) >= s% mu(ktop-1)) then ! stable against further mixing, so exit loop exit end if end if - end do - ! Call a final update over all mixed cells now. call update_model_(s, ktop, kbot+1, .true.) - end subroutine mix_outward real(dp) function blouin_delta_xo(Xin) - real(dp), intent(in) :: Xin ! mass fraction - real(dp) :: Xnew ! mass fraction - real(dp) :: xo, dxo ! number fractions + real(dp), intent(in) :: Xin ! mass fraction + real(dp) :: Xnew ! mass fraction + real(dp) :: xo, dxo ! number fractions real(dp) :: a0, a1, a2, a3, a4, a5 ! Convert input mass fraction to number fraction, assuming C/O mixture xo = (Xin/16d0)/(Xin/16d0 + (1d0 - Xin)/12d0) - a0 = 0d0 a1 = -0.311540d0 a2 = 2.114743d0 a3 = -1.661095d0 a4 = -1.406005d0 a5 = 1.263897d0 - dxo = & a0 + & a1*xo + & @@ -289,31 +364,30 @@ real(dp) function blouin_delta_xo(Xin) a3*xo*xo*xo + & a4*xo*xo*xo*xo + & a5*xo*xo*xo*xo*xo - xo = xo + dxo - ! Convert back to mass fraction Xnew = 16d0*xo/(16d0*xo + 12d0*(1d0-xo)) - blouin_delta_xo = Xnew - Xin end function blouin_delta_xo - real(dp) function blouin_delta_xne(Xin) - real(dp), intent(in) :: Xin ! mass fraction - real(dp) :: Xnew ! mass fraction - real(dp) :: xne, dxne ! number fractions + function blouin_delta_xne(Xin20,Xin22) + real(dp), intent(in) :: Xin20, Xin22! mass fraction + real(dp) :: Xnew1, Xnew2 ! mass fraction + real(dp) :: xne, dxne, xne1, xne2 ! number fractions real(dp) :: a0, a1, a2, a3, a4, a5 + real(dp), dimension(2) :: blouin_delta_xne ! Convert input mass fraction to number fraction, assuming O/Ne mixture - xne = (Xin/20d0)/(Xin/20d0 + (1d0 - Xin)/16d0) - + xne1 =(Xin20/20d0)/(Xin20/20d0 + Xin22/22d0 + (1d0 - Xin20 - Xin22)/16d0) + xne2 =(Xin22/22d0)/(Xin20/20d0 + Xin22/22d0 + (1d0 - Xin20 - Xin22)/16d0) + ! isotope 22Ne is added to the Ne separation along with 20Ne + xne =((Xin22/22d0)+(Xin20/20d0))/(Xin20/20d0 + Xin22/22d0 + (1d0 - Xin20 - Xin22)/16d0) a0 = 0d0 a1 = -0.120299d0 a2 = 1.304399d0 a3 = -1.722625d0 a4 = 0.393996d0 a5 = 0.144529d0 - dxne = & a0 + & a1*xne + & @@ -321,34 +395,316 @@ real(dp) function blouin_delta_xne(Xin) a3*xne*xne*xne + & a4*xne*xne*xne*xne + & a5*xne*xne*xne*xne*xne - - xne = xne + dxne - + xne1 = xne1 + dxne*xne1/xne + xne2 = xne2 + dxne*xne2/xne + xne = xne1 + xne2 ! Convert back to mass fraction - Xnew = 20d0*xne/(20d0*xne + 16d0*(1d0-xne)) - - blouin_delta_xne = Xnew - Xin + Xnew1 = (20d0*xne1)/(20d0*xne1 + 22d0*xne2 + 16d0*(1d0-xne)) + Xnew2 = (22d0*xne2)/(20d0*xne1 + 22d0*xne2 + 16d0*(1d0-xne)) + blouin_delta_xne(1) = Xnew1 - (Xin20) + blouin_delta_xne(2) = Xnew2 - (Xin22) end function blouin_delta_xne - subroutine update_model_ (s, kc_t, kc_b, do_brunt) + subroutine tab_interp_medin_cumming_dx1(x1_,x2_,components,dx1_) + use interp_2D_lib_db, only: interp_mkbicub_db, interp_evbicub_db + use const_def, only: mesa_data_dir + use utils_lib, only: mesa_error, mkdir, is_bad + implicit none + integer, parameter :: num_x1 = 998, num_x2 = 998 + integer :: ilinx,iliny,ibcxmin,ibcxmax,ibcymin,ibcymax,iounit,ict(6),ierr,i,j,k + real(dp) :: bcxmin(num_x1), bcxmax(num_x1) + real(dp) :: bcymin(num_x2), bcymax(num_x2) + real(dp), pointer, dimension(:) :: x1_l, x2_l, deltax1_sob_f1, deltax1_sob_values + real(dp), pointer :: deltax1_sob_f(:,:,:) + real(dp) :: deltax1,x1l,x2l + real(dp), intent(in) :: x1_,x2_ ! target of this interpolation + character (len=*), intent(in) :: components + real(dp) :: fval(6) ! output data + real(dp), intent(out) :: dx1_ + integer :: ier + + ict = 0 + ict(1) = 1 + iounit=999 + ! setup interpolation table for x1 x2 dx1 + if (components=='CONe') then + open(unit=iounit, file='CONe_deltaC.dat', action='read',status='old') + else if (components=='NeOMg') then + open(unit=iounit, file='NeOMg_deltaMg.dat', action='read',status='old') + else if (components=='ONeNa') then + open(unit=iounit, file='ONeNa_deltaNa.dat', action='read',status='old') + else if (components=='COMg') then + open(unit=iounit, file='COMg_deltaC.dat', action='read',status='old') + end if + allocate(x1_l(num_x1), x2_l(num_x2), & + deltax1_sob_f1(4*num_x1*num_x2)) + deltax1_sob_f(1:4,1:num_x1,1:num_x2) => & + deltax1_sob_f1(1:4*num_x1*num_x2) + do j=1,num_x1 + do i=1,num_x2 + read(iounit,*) x1l, x2l, deltax1 + x1_l(j)=x1l + if (j == 1) then + x2_l(i) =x2l + end if + deltax1_sob_f(1,j,i) = deltax1 + end do + end do + close(iounit) + ! just use "not a knot" bc's at edges of tables + ibcxmin = 0; bcxmin(1:num_x1) = 0 + ibcxmax = 0; bcxmax(1:num_x1) = 0 + ibcymin = 0; bcymin(1:num_x2) = 0 + ibcymax = 0; bcymax(1:num_x2) = 0 + call interp_mkbicub_db( & + x1_l, num_x1, x2_l, num_x2, deltax1_sob_f1, num_x1, & + ibcxmin,bcxmin,ibcxmax,bcxmax, & + ibcymin,bcymin,ibcymax,bcymax, & + ilinx,iliny,ierr) + if (ierr /= 0) then + write(*,*) 'interp_mkbicub_db error' + ierr = -1 + call mesa_error(__FILE__,__LINE__) + end if + do j=1,num_x1 + do i=1,num_x2 + do k=1,4 + if (is_bad(deltax1_sob_f(k,j,i))) then + write(*,*) 'deltax1_sob_f', i, j, k, deltax1_sob_f(k,j,i) + end if + end do + end do + end do + call interp_evbicub_db( & + x1_, x2_, x1_l, num_x1, x2_l, num_x2, & + ilinx, iliny, deltax1_sob_f1, num_x1, ict, fval, ier) + dx1_=fval(1) ! delta_x1 from 2d interpolation + end subroutine tab_interp_medin_cumming_dx1 + + + subroutine tab_interp_medin_cumming_dx2(x1_,x2_,components,dx2_) + !use utils_lib + use interp_2D_lib_db, only: interp_mkbicub_db, interp_evbicub_db + use const_def, only: mesa_data_dir + use utils_lib, only: mesa_error, mkdir, is_bad + implicit none + integer, parameter :: num_x1 = 998, num_x2 = 998 + integer :: ilinx,iliny,ibcxmin,ibcxmax,ibcymin,ibcymax,iounit,ict(6),ierr,i,j,k + real(dp) :: bcxmin(num_x1), bcxmax(num_x1) + real(dp) :: bcymin(num_x2), bcymax(num_x2) + real(dp), pointer, dimension(:) :: x1_l, x2_l, deltax1_sob_f1, deltax1_sob_values + real(dp), pointer :: deltax1_sob_f(:,:,:) + real(dp) :: deltax1,x1l,x2l + real(dp), intent(in) :: x1_,x2_ ! target of this interpolation + character (len=*), intent(in) :: components + real(dp) :: fval(6) ! output data + real(dp), intent(out) :: dx2_ + integer :: ier + + ict = 0 + ict(1) = 1 + iounit=998 + ! setup interpolation table for tau sob eta + if (components=='CONe') then + open(unit=iounit, file='CONe_deltaO.dat', action='read',status='old') + else if (components=='NeOMg') then + open(unit=iounit, file='NeOMg_deltaO.dat', action='read',status='old') + else if (components=='ONeNa') then + open(unit=iounit, file='ONeNa_deltaO.dat', action='read',status='old') + else if (components=='COMg') then + open(unit=iounit, file='COMg_deltaMg.dat', action='read',status='old') + end if + allocate(x1_l(num_x1), x2_l(num_x2), & + deltax1_sob_f1(4*num_x1*num_x2)) + deltax1_sob_f(1:4,1:num_x1,1:num_x2) => & + deltax1_sob_f1(1:4*num_x1*num_x2) + do j=1,num_x1 + do i=1,num_x2 + read(iounit,*) x1l, x2l, deltax1 + x1_l(j)=x1l + if (j == 1) then + x2_l(i) =x2l + end if + deltax1_sob_f(1,j,i) = deltax1 + end do + end do + close(iounit) + ! just use "not a knot" bc's at edges of tables + ibcxmin = 0; bcxmin(1:num_x1) = 0 + ibcxmax = 0; bcxmax(1:num_x1) = 0 + ibcymin = 0; bcymin(1:num_x2) = 0 + ibcymax = 0; bcymax(1:num_x2) = 0 + call interp_mkbicub_db( & + x1_l, num_x1, x2_l, num_x2, deltax1_sob_f1, num_x1, & + ibcxmin,bcxmin,ibcxmax,bcxmax, & + ibcymin,bcymin,ibcymax,bcymax, & + ilinx,iliny,ierr) + if (ierr /= 0) then + write(*,*) 'interp_mkbicub_db error' + ierr = -1 + call mesa_error(__FILE__,__LINE__) + end if + do j=1,num_x1 + do i=1,num_x2 + do k=1,4 + if (is_bad(deltax1_sob_f(k,j,i))) then + write(*,*) 'deltax1_sob_f', i, j, k, deltax1_sob_f(k,j,i) + end if + end do + end do + end do + call interp_evbicub_db( & + x1_, x2_, x1_l, num_x1, x2_l, num_x2, & + ilinx, iliny, deltax1_sob_f1, num_x1, ict, fval, ier) + dx2_=fval(1) ! delta_x2 from 2d interpolation + end subroutine tab_interp_medin_cumming_dx2 + + + subroutine medin_cumming_3p_d_cone(X1,X2,X3_1,X3_2,Dd) + real(dp), intent(in) :: X1, X2, X3_1, X3_2 ! mass fraction + real(dp), dimension(4),intent(out) :: Dd + real(dp) :: Xnew1, Xnew2, Xnew3_1, Xnew3_2, Xfac ! mass fraction + real(dp) :: xc, dxc, xo, dxo, xne1, xne2 ! number fractions + real(dp) :: dx1_,dx2_ + integer :: i,j + + Xfac = X1 + X2 + X3_1 + X3_2 + xc = (X1/12)/(X1/12 + X2/16 + X3_1/20 + X3_2/22) + xo = (X2/16)/(X1/12 + X2/16 + X3_1/20 + X3_2/22) + xne1 = (X3_1/20)/(X1/12 + X2/16 + X3_1/20 + X3_2/22) + xne2 = (X3_2/22)/(X1/12 + X2/16 + X3_1/20 + X3_2/22) + call tab_interp_medin_cumming_dx1(xc,xo,'CONe',dx1_) + call tab_interp_medin_cumming_dx2(xc,xo,'CONe',dx2_) + dxc=dx1_ + dxo=dx2_ + !write(*,*) 'delta_xc: ',dxc,' delta_xo: ', dxo + xc = xc + dxc + xo = xo + dxo + ! convert deltas in number fraction to mass fraction + Xnew1 = 12*xc/(12*xc + 16*xo + 20*(1-xc-xo)*(xne1)/(xne1+xne2)+22*(1-xc-xo)*(xne2)/(xne1+xne2)) + Xnew2 = 16*xo/(12*xc + 16*xo + 20*(1-xc-xo)*(xne1)/(xne1+xne2)+22*(1-xc-xo)*(xne2)/(xne1+xne2)) + Xnew3_1 = (20*(1-xc-xo)*(xne1)/(xne1+xne2))/(12*xc + 16*xo + 20*(1-xc-xo)*(xne1) & + /(xne1+xne2)+22*(1-xc-xo)*(xne2)/(xne1+xne2)) + Xnew3_2 = (22*(1-xc-xo)*(xne2)/(xne1+xne2))/(12*xc + 16*xo + 20*(1-xc-xo)*(xne1) & + /(xne1+xne2)+22*(1-xc-xo)*(xne2)/(xne1+xne2)) + Dd=[0,0,0,0] + Dd(1)= Xnew1 - X1 + Dd(2)= Xnew2 - X2 + Dd(3)= Xnew3_1 - X3_1 + Dd(4)= Xnew3_2 - X3_2 + !write(*,*) 'delta_XC: ',Dd(1),' delta_XO: ', Dd(2), 'delta_XNe:', Dd(3)+Dd(4) + end subroutine medin_cumming_3p_d_cone + + subroutine medin_cumming_3p_d_neomg(X1,X2,X3_1,X3_2,Dd) + real(dp), intent(in) :: X1, X2, X3_1, X3_2 ! mass fraction + real(dp), dimension(4),intent(out) :: Dd + real(dp) :: Xnew1, Xnew2, Xnew3_1, Xnew3_2, Xfac ! mass fraction + real(dp) :: xmg, dxmg, xo, dxo, xne1, xne2 ! number fractions + real(dp) :: dx1_,dx2_ + integer :: i,j + + Xfac = X1 + X2 + X3_1 + X3_2 + xmg = (X1/24)/(X1/24 + X2/16 + X3_1/20 + X3_2/22) + xo = (X2/16)/(X1/24 + X2/16 + X3_1/20 + X3_2/22) + xne1 = (X3_1/20)/(X1/24 + X2/16 + X3_1/20 + X3_2/22) + xne2 = (X3_2/22)/(X1/24 + X2/16 + X3_1/20 + X3_2/22) + call tab_interp_medin_cumming_dx1(xmg,xo,'NeOMg',dx1_) + call tab_interp_medin_cumming_dx2(xmg,xo,'NeOMg',dx2_) + dxmg=dx1_ + dxo=dx2_ + xmg = xmg + dxmg + xo = xo + dxo + ! convert deltas in number fraction to mass fraction + Xnew1 = 24*xmg/(24*xmg + 16*xo + 20*(1-xmg-xo)*(xne1)/(xne1+xne2)+22*(1-xmg-xo)*(xne2)/(xne1+xne2)) + Xnew2 = 16*xo/(24*xmg + 16*xo + 20*(1-xmg-xo)*(xne1)/(xne1+xne2)+22*(1-xmg-xo)*(xne2)/(xne1+xne2)) + Xnew3_1 = (20*(1-xmg-xo)*(xne1)/(xne1+xne2))/(24*xmg + 16*xo + 20*(1-xmg-xo)*(xne1) & + /(xne1+xne2)+22*(1-xmg-xo)*(xne2)/(xne1+xne2)) + Xnew3_2 = (22*(1-xmg-xo)*(xne2)/(xne1+xne2))/(24*xmg + 16*xo + 20*(1-xmg-xo)*(xne1) & + /(xne1+xne2)+22*(1-xmg-xo)*(xne2)/(xne1+xne2)) + Dd=[0,0,0,0] + Dd(1)= Xnew1 - X1 + Dd(2)= Xnew2 - X2 + Dd(3)= Xnew3_1 - X3_1 + Dd(4)= Xnew3_2 - X3_2 + !write(*,*) 'delta_XMg: ',Dd(1),' delta_XO: ', Dd(2), 'delta_XNe:', Dd(3)+Dd(4) + end subroutine medin_cumming_3p_d_neomg + + subroutine medin_cumming_3p_d_onena(X1,X2,X3_1,X3_2,Dd) + real(dp), intent(in) :: X1, X2, X3_1, X3_2 ! mass fraction + real(dp), dimension(4),intent(out) :: Dd + real(dp) :: Xnew1, Xnew2, Xnew3_1, Xnew3_2, Xfac ! mass fraction + real(dp) :: xna, dxna, xo, dxo, xne1, xne2 ! number fractions + real(dp) :: dx1_,dx2_ + integer :: i,j + + Xfac = X1 + X2 + X3_1 + X3_2 + xna = (X1/23)/(X1/23 + X2/16 + X3_1/20 + X3_2/22) + xo = (X2/16)/(X1/23 + X2/16 + X3_1/20 + X3_2/22) + xne1 = (X3_1/20)/(X1/23 + X2/16 + X3_1/20 + X3_2/22) + xne2 = (X3_2/22)/(X1/23 + X2/16 + X3_1/20 + X3_2/22) + call tab_interp_medin_cumming_dx1(xna,xo,'ONeNa',dx1_) + call tab_interp_medin_cumming_dx2(xna,xo,'ONeNa',dx2_) + dxna=dx1_ + dxo=dx2_ + !write(*,*) xna,xo + !write(*,*) 'delta_xna: ',dxna,' delta_xo: ', dxo + xna = xna + dxna + xo = xo + dxo + ! convert deltas in number fraction to mass fraction + Xnew1 = 23*xna/(23*xna + 16*xo + 20*(1-xna-xo)*(xne1)/(xne1+xne2)+22*(1-xna-xo)*(xne2)/(xne1+xne2)) + Xnew2 = 16*xo/(23*xna + 16*xo + 20*(1-xna-xo)*(xne1)/(xne1+xne2)+22*(1-xna-xo)*(xne2)/(xne1+xne2)) + Xnew3_1 = (20*(1-xna-xo)*(xne1)/(xne1+xne2))/(23*xna + 16*xo + 20*(1-xna-xo)*(xne1) & + /(xne1+xne2)+22*(1-xna-xo)*(xne2)/(xne1+xne2)) + Xnew3_2 = (22*(1-xna-xo)*(xne2)/(xne1+xne2))/(23*xna + 16*xo + 20*(1-xna-xo)*(xne1) & + /(xne1+xne2)+22*(1-xna-xo)*(xne2)/(xne1+xne2)) + Dd=[0,0,0,0] + Dd(1)= Xnew1 - X1 + Dd(2)= Xnew2 - X2 + Dd(3)= Xnew3_1 - X3_1 + Dd(4)= Xnew3_2 - X3_2 + !write(*,*) 'delta_XNa: ',Dd(1),' delta_XO: ', Dd(2), 'delta_XNe:', Dd(3)+Dd(4) + end subroutine medin_cumming_3p_d_onena + + subroutine medin_cumming_3p_d_comg(X1,X2,X3,Dd) + real(dp), intent(in) :: X1, X2, X3 ! mass fraction + real(dp), dimension(4),intent(out) :: Dd + real(dp) :: Xnew1, Xnew2, Xfac ! mass fraction + real(dp) :: xc, dxc, xmg, dxmg, xo ! number fractions + real(dp) :: dx1_,dx2_ + integer :: i,j + + Xfac = X1 + X2 + X3 + xc = (X1/12)/(X1/12 + X2/24 + X3/16) + xmg = (X2/24)/(X1/12 + X2/24 + X3/16) + xo = (X3/16)/(X1/12 + X2/24 + X3/16) + call tab_interp_medin_cumming_dx1(xc,xmg,'COMg',dx1_) + call tab_interp_medin_cumming_dx2(xc,xmg,'COMg',dx2_) + dxc=dx1_ + dxmg=dx2_ + xc = xc + dxc + xmg = xmg + dxmg + ! convert deltas in number fraction to mass fraction + Xnew1 = 12*xc/(12*xc + 24*xmg + 16*(1-xc-xmg)) + Xnew2 = 24*xmg/(12*xc + 24*xmg + 16*(1-xc-xmg)) + Dd=[0,0,0,0] + Dd(1)= Xnew1 - X1 + Dd(2)= Xnew2 - X2 + end subroutine medin_cumming_3p_d_comg + subroutine update_model_ (s, kc_t, kc_b, do_brunt) use turb_info, only: set_mlt_vars use brunt, only: do_brunt_B use micro - type(star_info), pointer :: s integer, intent(in) :: kc_t integer, intent(in) :: kc_b logical, intent(in) :: do_brunt - integer :: ierr integer :: kf_t integer :: kf_b - logical :: mask(s%nz) mask(:) = .true. - ! Update the model to reflect changes in the abundances across ! cells kc_t:kc_b (the mask part of this call is unused, mask=true for all zones). ! Do updates at constant (P,T) rather than constant (rho,T). @@ -359,42 +715,35 @@ subroutine update_model_ (s, kc_t, kc_b, do_brunt) stop end if s%fix_Pgas = .false. - ! Update opacities across cells kc_t:kc_b (this also sets rho_face ! and related quantities on faces kc_t:kc_b) call set_micro_vars(s, kc_t, kc_b, & - skip_eos=.TRUE., skip_net=.TRUE., skip_neu=.TRUE., skip_kap=.FALSE., ierr=ierr) + skip_eos=.TRUE., skip_net=.TRUE., skip_neu=.TRUE., skip_kap=.TRUE., ierr=ierr) if (ierr /= 0) then write(*,*) 'phase_separation: error from call to set_micro_vars' stop end if - ! This is expensive, so only do it if we really need to. if(do_brunt) then ! Need to make sure we can set brunt for mix_outward calculation. if(.not. s% calculate_Brunt_B) then stop "phase separation requires s% calculate_Brunt_B = .true." end if - call do_brunt_B(s, kc_t, kc_b, ierr) ! for unsmoothed_brunt_B + call do_brunt_B(s, kc_t, kc_b, ierr) ! for unsmoothed_brunt_B if (ierr /= 0) then write(*,*) 'phase_separation: error from call to do_brunt_B' stop end if end if - ! Finally update MLT for interior faces - kf_t = kc_t kf_b = kc_b + 1 - call set_mlt_vars(s, kf_t+1, kf_b-1, ierr) if (ierr /= 0) then write(*,*) 'phase_separation: failed in call to set_mlt_vars during update_model_' stop - end if - + endif + ! Finish return - end subroutine update_model_ - - end module phase_separation +end module phase_separation diff --git a/star/test_suite/wd_o_ne_3_phase/.gitattributes b/star/test_suite/wd_o_ne_3_phase/.gitattributes new file mode 100644 index 000000000..f46bd6cb2 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/.gitattributes @@ -0,0 +1,2 @@ +*.dat filter=lfs diff=lfs merge=lfs -text +*.mod filter=lfs diff=lfs merge=lfs -text diff --git a/star/test_suite/wd_o_ne_3_phase/COMg_deltaC.dat b/star/test_suite/wd_o_ne_3_phase/COMg_deltaC.dat new file mode 100644 index 000000000..70f15a7b1 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/COMg_deltaC.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5c96bb17d819ca08bb9c1e3704849b78e892a04f6a38de93833bb17ebff15b66 +size 33341743 diff --git a/star/test_suite/wd_o_ne_3_phase/COMg_deltaMg.dat b/star/test_suite/wd_o_ne_3_phase/COMg_deltaMg.dat new file mode 100644 index 000000000..41e0538f1 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/COMg_deltaMg.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:417bd76545e6d43e2073e100b8db25914db0cb032999539b263c850ecea251d8 +size 33018225 diff --git a/star/test_suite/wd_o_ne_3_phase/CONe_deltaC.dat b/star/test_suite/wd_o_ne_3_phase/CONe_deltaC.dat new file mode 100644 index 000000000..8eb4add7c --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/CONe_deltaC.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:020c7a90ca01736aee534dbb3dc60a0dffa0fb5d800ef5c0b952d4f21a7b081a +size 33347923 diff --git a/star/test_suite/wd_o_ne_3_phase/CONe_deltaO.dat b/star/test_suite/wd_o_ne_3_phase/CONe_deltaO.dat new file mode 100644 index 000000000..a0cfb70c1 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/CONe_deltaO.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d2ff7d64b227caa2d021e09cfa43013b37d8c94e59edfc769ccc6454a1501acf +size 33106276 diff --git a/star/test_suite/wd_o_ne_3_phase/NeOMg_deltaMg.dat b/star/test_suite/wd_o_ne_3_phase/NeOMg_deltaMg.dat new file mode 100644 index 000000000..a760b0652 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/NeOMg_deltaMg.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ac05f770cbddf6089f515a6a75dbe90ba287ade216b8b6d953bcfbef545c1e0b +size 32912885 diff --git a/star/test_suite/wd_o_ne_3_phase/NeOMg_deltaO.dat b/star/test_suite/wd_o_ne_3_phase/NeOMg_deltaO.dat new file mode 100644 index 000000000..7b7b4b549 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/NeOMg_deltaO.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:30f8cda541ed77e92cd4bb7b9feaece2130cf6a7c58e6f210ec3fae771ead632 +size 33352468 diff --git a/star/test_suite/wd_o_ne_3_phase/ONeNa_deltaNa.dat b/star/test_suite/wd_o_ne_3_phase/ONeNa_deltaNa.dat new file mode 100644 index 000000000..14724dcd6 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/ONeNa_deltaNa.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ec03a91a0d63bd6e8041bcad27068e4ec7f6c8f2cd84b104fb4ca94c2ed4130b +size 32892398 diff --git a/star/test_suite/wd_o_ne_3_phase/ONeNa_deltaO.dat b/star/test_suite/wd_o_ne_3_phase/ONeNa_deltaO.dat new file mode 100644 index 000000000..3c03eb61e --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/ONeNa_deltaO.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:faefa3857cc651717ab03a0b8832edf0187eff99516c01dfb307c990e81db720 +size 33357817 diff --git a/star/test_suite/wd_o_ne_3_phase/clean b/star/test_suite/wd_o_ne_3_phase/clean new file mode 100755 index 000000000..95545a5c1 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/clean @@ -0,0 +1,4 @@ +#!/bin/bash + +cd make +make clean diff --git a/star/test_suite/wd_o_ne_3_phase/custom.net b/star/test_suite/wd_o_ne_3_phase/custom.net new file mode 100644 index 000000000..879f7013e --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/custom.net @@ -0,0 +1,8 @@ + + include 'basic.net' + + add_isos( + o18 + ne22 + na23 + ) diff --git a/star/test_suite/wd_o_ne_3_phase/history_columns.list b/star/test_suite/wd_o_ne_3_phase/history_columns.list new file mode 100644 index 000000000..d7dbb957c --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/history_columns.list @@ -0,0 +1,1053 @@ +! history_columns.list -- determines the contents of star history logs +! you can use a non-standard version by setting history_columns_file in your inlist + +! units are cgs unless otherwise noted. + +! reorder the following names as desired to reorder columns. +! comment out the name to omit a column (fewer columns => less IO => faster running). +! remove '!' to restore a column. + +! if you have a situation where you want a non-standard set of columns, +! make a copy of this file, edit as desired, and give the new filename in your inlist +! as history_columns_file. if you are just adding columns, you can 'include' this file, +! and just list the additions in your file. note: to include the standard default +! version, use include '' -- the 0 length string means include the default file. + +! blank lines and comments can be used freely. +! if a column name appears more than once in the list, only the first occurrence is used. + +! if you need to have something added to the list of options, let me know.... + + +! the first few lines of the log file contain a few items: + + ! version_number -- for the version of mesa being used + ! burn_min1 -- 1st limit for reported burning, in erg/g/s + ! burn_min2 -- 2nd limit for reported burning, in erg/g/s + + +!# other files + +! note: you can include another list by doing +! include 'filename' +! include '' means include the default standard list file + +! the following lines of the log file contain info about 1 model per row + +!---------------------------------------------------------------------------------------------- + +!# general info about the model + + model_number ! counting from the start of the run + num_zones ! number of zones in the model + + !## age + + star_age ! elapsed simulated time in years since the start of the run + !star_age_sec ! elapsed simulated time in seconds since the start of the run + !star_age_min ! elapsed simulated time in minutes since the start of the run + !star_age_hr ! elapsed simulated time in hours since the start of the run + !star_age_day ! elapsed simulated time in days since the start of the run + !day ! elapsed simulated time in days since the start of the run + + !log_star_age + !log_star_age_sec + + !## timestep + + !time_step ! timestep in years since previous model + time_step_sec ! timestep in seconds since previous model + !time_step_days + log_dt ! log10 time_step in years + !log_dt_sec ! log10 time_step in seconds + !log_dt_days ! log10 time_step in days + + !## mass + + star_mass ! in Msun units + !log_star_mass + + !star_gravitational_mass ! star_mass is baryonic mass + !star_mass_grav_div_mass + + !delta_mass ! star_mass - initial_mass in Msun units + log_xmstar ! log10 mass exterior to M_center (grams) + + !## mass change + + !star_mdot ! d(star_mass)/dt (in msolar per year) + log_abs_mdot ! log10(abs(star_mdot)) (in msolar per year) + + !## imposed surface conditions + !Tsurf_factor + !tau_factor + !tau_surface + + !## imposed center conditions + !m_center + !m_center_gm + !r_center + !r_center_cm + !r_center_km + !L_center + !log_L_center + !log_L_center_ergs_s + !v_center + !v_center_kms + + !logt_max + +!---------------------------------------------------------------------------------------------- + +!# mixing and convection + + !max_conv_vel_div_csound + !max_gradT_div_grada + !max_gradT_sub_grada + !min_log_mlt_Gamma + + + !## mixing regions + + mass_conv_core ! (Msun) mass coord of top of convective core. 0 if core is not convective + + ! mx1 refers to the largest (by mass) convective region. + ! mx2 is the 2nd largest. + + ! conv_mx1_top and conv_mx1_bot are the region where mixing_type == convective_mixing. + ! mx1_top and mx1_bot are the extent of all kinds of mixing, convective and other. + + ! values are m/Mstar + conv_mx1_top + conv_mx1_bot + conv_mx2_top + conv_mx2_bot + mx1_top + mx1_bot + mx2_top + mx2_bot + + ! radius -- values are radii in Rsun units + !conv_mx1_top_r + !conv_mx1_bot_r + !conv_mx2_top_r + !conv_mx2_bot_r + !mx1_top_r + !mx1_bot_r + !mx2_top_r + !mx2_bot_r + + ! you might want to get a more complete list of mixing regions by using the following + + !mixing_regions ! note: this includes regions where the mixing type is no_mixing. + + ! the is the number of regions to report + ! there will be 2* columns for this in the log file, 2 for each region. + ! the first column for a region gives the mixing type as defined in const/public/const_def.f90. + + ! the second column for a region gives the m/mstar location of the top of the region + ! entries for extra columns after the last region in the star will have an invalid mixing_type value of -1. + ! mstar is the total mass of the star, so these locations range from 0 to 1 + ! all regions are include starting from the center, so the bottom of one region + ! is the top of the previous one. since we start at the center, the bottom of the 1st region is 0. + + ! the columns in the log file will have names like 'mix_type_1' and 'mix_qtop_1' + + ! if the star has too many regions to report them all, + ! the smallest regions will be merged with neighbors for reporting purposes only. + + + !mix_relr_regions + ! same as above, but locations given as r/rstar instead of m/mstar. + ! the columns in the log file will have names like 'mix_relr_type_1' and 'mix_relr_top_1' + + + !## conditions at base of largest convection zone (by mass) + !cz_bot_mass ! mass coordinate of base (Msun) + !cz_mass ! mass coordinate of base (Msun) -- same as cz_bot_mass + !cz_log_xmass ! mass exterior to base (g) + !cz_log_xmsun ! mass exterior to base (Msun) + !cz_xm ! mass exterior to base (Msun) + !cz_logT + !cz_logRho + !cz_logP + !cz_bot_radius ! Rsun + !cz_log_column_depth + !cz_log_radial_depth + !cz_luminosity ! Lsun + !cz_opacity + !cz_log_tau + !cz_eta + !cz_log_eps_nuc ! log10(ergs/g/s) + !cz_t_heat ! Cp*T/eps_nuc (seconds) + + !cz_csound + !cz_scale_height + !cz_grav + + !cz_omega + !cz_omega_div_omega_crit + + !cz_zone + + ! mass fractions at base of largest convection zone (by mass) + !cz_log_xa h1 + !cz_log_xa he4 + + !## conditions at top of largest convection zone (by mass) + !cz_top_mass ! mass coordinate of top (Msun) + !cz_top_log_xmass ! mass exterior to top (g) + !cz_top_log_xmsun ! mass exterior to top (Msun) + !cz_top_xm ! mass exterior to top (Msun) + !cz_top_logT + !cz_top_logRho + !cz_top_logP + !cz_top_radius ! Rsun + !cz_top_log_column_depth + !cz_top_log_radial_depth + !cz_top_luminosity ! Lsun + !cz_top_opacity + !cz_top_log_tau + !cz_top_eta + !cz_top_log_eps_nuc ! log10(ergs/g/s) + !cz_top_t_heat ! Cp*T/eps_nuc (seconds) + + !cz_top_csound + !cz_top_scale_height + !cz_top_grav + + !cz_top_omega + !cz_top_omega_div_omega_crit + + !cz_top_zone + !cz_top_zone_logdq + + ! mass fractions at top of largest convection zone (by mass) + !cz_top_log_xa h1 + !cz_top_log_xa he4 + +!---------------------------------------------------------------------------------------------- + +!# nuclear reactions + + !## integrated quantities + + !power_h_burn ! total thermal power from PP and CNO, excluding neutrinos (in Lsun units) + !power_he_burn ! total thermal power from triple-alpha, excluding neutrinos (in Lsun units) + !power_photo + !power_z_burn + !log_power_nuc_burn ! total thermal power from all burning, excluding photodisintegrations + log_LH ! log10 power_h_burn + log_LHe ! log10 power_he_burn + log_LZ ! log10 total burning power including LC, but excluding LH and LHe and photodisintegrations + log_Lnuc ! log(LH + LHe + LZ) + !log_Lnuc_ergs_s + !log_Lnuc_sub_log_L + !lnuc_photo + + !extra_L ! integral of extra_heat in Lsun units + !log_extra_L ! log10 extra_L + + !## neutrino losses + !log_Lneu ! log10 power emitted in neutrinos, nuclear and thermal (in Lsun units) + !log_Lneu_nuc ! log10 power emitted in neutrinos, nuclear sources only (in Lsun units) + !log_Lneu_nonnuc ! log10 power emitted in neutrinos, thermal sources only (in Lsun units) + + !mass_loc_of_max_eps_nuc ! (in Msun units) + !mass_ext_to_max_eps_nuc ! (in Msun units) + !eps_grav_integral ! (in Lsun units) + !log_abs_Lgrav ! log10 abs(eps_grav_integral) (in Lsun units) + + !## information about reactions (by category) + + ! log10 total luminosity for reaction categories (Lsun units) + + pp + cno + !tri_alfa + !c_alpha + !n_alpha + !o_alpha + !ne_alpha + !na_alpha + !mg_alpha + !si_alpha + !s_alpha + !ar_alpha + !ca_alpha + !ti_alpha + !fe_co_ni + !c12_c12 + !c12_o16 + !o16_o16 + !photo + !pnhe4 + !other + + + + + !## nuclear reactions at center + + ! center log10 burn erg/g/s for reaction categories + + !c_log_eps_burn cno + !c_log_eps_burn tri_alfa + + ! center d_eps_nuc_dlnd for reaction categories + + !c_d_eps_dlnd cno + !c_d_eps_dlnd tri_alfa + + ! center d_eps_nuc_dlnT for reaction categories + + !c_d_eps_dlnT cno + !c_d_eps_dlnT tri_alfa + + !## regions of strong nuclear burning + + ! 2 zones where eps_nuc > burn_min1 erg/g/s + ! for each zone have 4 numbers: start1, start2, end2, end1 + ! start1 is mass of inner edge where first goes > burn_min1 (or -20 if none such) + ! start2 is mass of inner edge where first zone reaches burn_min2 erg/g/sec (or -20 if none such) + ! end2 is mass of outer edge where first zone drops back below burn_min2 erg/g/s + ! end1 is mass of outer edge where first zone ends (i.e. eps_nuc < burn_min1) + ! similar for the second zone + + epsnuc_M_1 ! start1 for 1st zone + epsnuc_M_2 ! start2 + epsnuc_M_3 ! end2 + epsnuc_M_4 ! end1 + + epsnuc_M_5 ! start1 for 2nd zone + epsnuc_M_6 ! start2 + epsnuc_M_7 ! end2 + epsnuc_M_8 ! end1 + + + ! you might want to get a more complete list of burning regions by using the following + + !burning_regions + ! the is the number of regions to report + ! there will be 2* columns for this in the log file, 2 for each region. + ! the first column for a region gives int(sign(val)*log10(max(1,abs(val)))) + ! where val = ergs/gm/sec nuclear energy minus all neutrino losses. + ! the second column for a region gives the q location of the top of the region + ! entries for extra columns after the last region in the star will have a value of -9999 + ! all regions are included starting from the center, so the bottom of one region + ! is the top of the previous one. + ! since we start at the center, the bottom of the 1st region is q=0 and top of last is q=1. + + ! the columns in the log file will have names like 'burn_type_1' and 'burn_qtop_1' + + !burn_relr_regions + ! same as above, but locations given as r/rstar instead of m/mstar. + ! the columns in the log file will have names like 'burn_relr_type_1' and 'burn_relr_top_1' + + + ! if the star has too many regions to report them all, + ! the smallest regions will be merged with neighbors for reporting purposes only. + +!---------------------------------------------------------------------------------------------- + +!# information about core and envelope + + !## helium core + he_core_mass + !he_core_radius + !he_core_lgT + !he_core_lgRho + !he_core_L + !he_core_v + !he_core_omega + !he_core_omega_div_omega_crit + !he_core_k + + !## CO core + co_core_mass + !CO_core + !co_core_radius + !co_core_lgT + !co_core_lgRho + !co_core_L + !co_core_v + !co_core_omega + !co_core_omega_div_omega_crit + !co_core_k + + !## ONe core + one_core_mass + !one_core_radius + !one_core_lgT + !one_core_lgRho + !one_core_L + !one_core_v + !one_core_omega + !one_core_omega_div_omega_crit + !one_core_k + + !## iron core + fe_core_mass + !fe_core_radius + !fe_core_lgT + !fe_core_lgRho + !fe_core_L + !fe_core_v + !fe_core_omega + !fe_core_omega_div_omega_crit + !fe_core_k + + !## neuton rich core + neutron_rich_core_mass + !neutron_rich_core_radius + !neutron_rich_core_lgT + !neutron_rich_core_lgRho + !neutron_rich_core_L + !neutron_rich_core_v + !neutron_rich_core_omega + !neutron_rich_core_omega_div_omega_crit + !neutron_rich_core_k + + !## envelope + + envelope_mass ! = star_mass - he_core_mass + !envelope_fraction_left ! = envelope_mass / (initial_mass - he_core_mass) + + !h_rich_layer_mass ! = star_mass - he_core_mass + !he_rich_layer_mass ! = he_core_mass - c_core_mass + !co_rich_layer_mass + +!---------------------------------------------------------------------------------------------- + +!# timescales + + !dynamic_timescale ! dynamic timescale (seconds) -- estimated by 2*pi*sqrt(r^3/(G*m)) + !kh_timescale ! kelvin-helmholtz timescale (years) + !mdot_timescale ! star_mass/abs(star_mdot) (years) + !kh_div_mdot_timescales ! kh_timescale/mdot_timescale + !nuc_timescale ! nuclear timescale (years) -- proportional to mass divided by luminosity + + !dt_cell_collapse ! min time for any cell to collapse at current velocities + !dt_div_dt_cell_collapse + + !dt_div_max_tau_conv ! dt/ maximum conv timescale + !dt_div_min_tau_conv ! dt/ minimum conv timescale + + + !min_dr_div_cs ! min over all cells of dr/csound (seconds) + !min_dr_div_cs_k ! location of min + !log_min_dr_div_cs ! log10 min dr_div_csound (seconds) + !min_dr_div_cs_yr ! min over all cells of dr/csound (years) + !log_min_dr_div_cs_yr ! log10 min dr_div_csound (years) + !dt_div_min_dr_div_cs + !log_dt_div_min_dr_div_cs + + !min_t_eddy ! minimum value of scale_height/conv_velocity + +!---------------------------------------------------------------------------------------------- + +!# conditions at or near the surface of the model + + !## conditions at the photosphere + !effective_T + !Teff + log_Teff ! log10 effective temperature + ! Teff is calculated using Stefan-Boltzmann relation L = 4 pi R^2 sigma Teff^4, + ! where L and R are evaluated at the photosphere (tau_factor < 1) + ! or surface of the model (tau_factor >= 1) when photosphere is not inside the model. + + !photosphere_black_body_T + !photosphere_cell_T ! temperature at model location closest to the photosphere, not necessarily Teff + !photosphere_cell_log_T + !photosphere_cell_density + !photosphere_cell_log_density + !photosphere_cell_opacity + !photosphere_cell_log_opacity + !photosphere_L ! Lsun units + !photosphere_log_L ! Lsun units + !photosphere_r ! Rsun units + !photosphere_log_r ! Rsun units + !photosphere_m ! Msun units + !photosphere_v_km_s + !photosphere_cell_k + !photosphere_column_density + !photosphere_csound + !photosphere_log_column_density + !photosphere_opacity + !photosphere_v_div_cs + !photosphere_xm + !photosphere_cell_free_e + !photosphere_cell_log_free_e + !photosphere_logg + !photosphere_T + + !## conditions at or near the surface of the model (outer edge of outer cell) + + !luminosity ! luminosity in Lsun units + luminosity_ergs_s ! luminosity in cgs units + log_L ! log10 luminosity in Lsun units + !log_L_ergs_s ! log10 luminosity in cgs units + !radius ! Rsun + log_R ! log10 radius in Rsun units + !radius_cm + !log_R_cm + + log_g ! log10 gravity + !gravity + !log_Ledd + log_L_div_Ledd ! log10(L/Leddington) + !lum_div_Ledd + !log_surf_optical_depth + !surface_optical_depth + + !log_surf_cell_opacity ! old name was log_surf_opacity + !log_surf_cell_P ! old name was log_surf_P + !log_surf_cell_pressure ! old name was log_surf_pressure + !log_surf_cell_density ! old name was log_surf_density + !log_surf_cell_temperature ! old name was log_surf_temperature + !surface_cell_temperature ! old name was surface_temperature + !log_surf_cell_z ! old name was log_surf_z + !surface_cell_entropy ! in units of kerg per baryon + ! old name was surface_entropy + + !v_surf ! (cm/s) + !v_surf_km_s ! (km/s) + v_div_csound_surf ! velocity divided by sound speed at outermost grid point + !v_div_csound_max ! max value of velocity divided by sound speed at face + !v_div_vesc + !v_phot_km_s + !v_surf_div_escape_v + + !v_surf_div_v_kh ! v_surf/(photosphere_r/kh_timescale) + + !surf_avg_j_rot + !surf_avg_omega + !surf_avg_omega_crit + !surf_avg_omega_div_omega_crit + !surf_avg_v_rot ! km/sec rotational velocity at equator + !surf_avg_v_crit ! critical rotational velocity at equator + !surf_avg_v_div_v_crit + !surf_avg_Lrad_div_Ledd + !surf_avg_logT + !surf_avg_logRho + !surf_avg_opacity + + ! Gravity Darkening, reports the surface averaged L/Lsun and Teff (K) caused by + ! gravity darkening in rotating stars. Based on the model of Espinosa Lara & Rieutord (2011) + ! 'polar' refers to the line of sight being directed along the rotation axis of the star + ! 'equatorial' refers to the line of sight coincident with the stellar equator + !grav_dark_L_polar !Lsun + !grav_dark_Teff_polar !K + !grav_dark_L_equatorial !Lsun + !grav_dark_Teff_equatorial !K + + !surf_escape_v ! cm/s + + !v_wind_Km_per_s ! Km/s + ! = 1d-5*s% opacity(1)*max(0d0,-s% mstar_dot)/ & + ! (4*pi*s% photosphere_r*Rsun*s% tau_base) + ! Lars says: + ! wind_mdot = 4*pi*R^2*rho*v_wind + ! tau = integral(opacity*rho*dr) from R to infinity + ! so tau = opacity*wind_mdot/(4*pi*R*v_wind) at photosphere + ! or v_wind = opacity*wind_mdot/(4*pi*R*tau) at photosphere + + !rotational_mdot_boost ! factor for increase in mass loss mdot due to rotation + !log_rotational_mdot_boost ! log factor for increase in mass loss mdot due to rotation + !surf_r_equatorial_div_r_polar + !surf_r_equatorial_div_r + !surf_r_polar_div_r + +!---------------------------------------------------------------------------------------------- + +!# conditions near center + + log_center_T ! temperature + log_center_Rho ! density + log_center_P ! pressure + + ! shorter names for above + log_cntr_P + log_cntr_Rho + log_cntr_T + + !center_T ! temperature + !center_Rho ! density + !center_P ! pressure + + center_degeneracy ! the electron chemical potential in units of k*T + center_gamma ! plasma interaction parameter + center_mu + center_ye + center_abar + !center_zbar + + !center_eps_grav + + !center_non_nuc_neu + !center_eps_nuc + !d_center_eps_nuc_dlnT + !d_center_eps_nuc_dlnd + !log_center_eps_nuc + + center_entropy ! in units of kerg per baryon + !max_entropy ! in units of kerg per baryon + !fe_core_infall + !non_fe_core_infall + !non_fe_core_rebound + !max_infall_speed + + !compactness_parameter ! (m/Msun)/(R(m)/1000km) for m = 2.5 Msun + !compactness + !m4 ! Mass co-ordinate where entropy=4 + !mu4 ! dM(Msun)/dr(1000km) where entropy=4 + + !center_omega + !center_omega_div_omega_crit + +!---------------------------------------------------------------------------------------------- + +!# abundances + + !species ! size of net + + !## mass fractions near center + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_center_abundances + !add_log_center_abundances + + ! individual central mass fractions (as many as desired) + center h1 + center he4 + center c12 + center o16 + + ! individual log10 central mass fractions (as many as desired) + !log_center h1 + !log_center he4 + ! etc. + + + !## mass fractions near surface + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_surface_abundances + !add_log_surface_abundances + + ! individual surface mass fractions (as many as desired) + !surface h1 + !surface he4 + surface c12 + surface o16 + ! etc. + + ! individual log10 surface mass fractions (as many as desired) + + !log_surface h1 + !log_surface he4 + + + !## mass fractions for entire star + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_average_abundances + !add_log_average_abundances + + ! individual average mass fractions (as many as desired) + !average h1 + !average he4 + ! etc. + + ! individual log10 average mass fractions (as many as desired) + !log_average h1 + !log_average he4 + ! etc. + + + !## mass totals for entire star (in Msun units) + + ! the following controls automatically add columns for all of the isos that are in the current net + !add_total_mass + !add_log_total_mass + + ! individual mass totals for entire star (as many as desired) + total_mass h1 + total_mass he4 + ! etc. + + ! individial log10 mass totals for entire star (in Msun units) + !log_total_mass h1 + !log_total_mass he4 + ! etc. + +!---------------------------------------------------------------------------------------------- + +!# info at specific locations + + !## info at location of max temperature + !max_T + !log_max_T + + +!---------------------------------------------------------------------------------------------- + +!# information about shocks + + !## info about outermost outward moving shock + ! excluding locations with q > max_q_for_outer_mach1_location + ! returns values at location of max velocity + !shock_mass ! baryonic (Msun) + !shock_mass_gm ! baryonic (grams) + !shock_q + !shock_radius ! (Rsun) + !shock_radius_cm ! (cm) + !shock_velocity + !shock_csound + !shock_v_div_cs + !shock_lgT + !shock_lgRho + !shock_lgP + !shock_gamma1 + !shock_entropy + !shock_tau + !shock_k + !shock_pre_lgRho + +!---------------------------------------------------------------------------------------------- + +!# asteroseismology + + !delta_nu ! large frequency separation for p-modes (microHz) + ! 1e6/(seconds for sound to cross diameter of star) + !delta_Pg ! g-mode period spacing for l=1 (seconds) + ! sqrt(2) pi^2/(integral of brunt_N/r dr) + !log_delta_Pg + !nu_max ! estimate from scaling relation (microHz) + ! nu_max = nu_max_sun * M/Msun / ((R/Rsun)^2 (Teff/Teff_sun)^0.5) + ! with nu_max_sun = 3100 microHz, Teff_sun = 5777 + !nu_max_3_4th_div_delta_nu ! nu_max^0.75/delta_nu + !acoustic_cutoff ! 0.5*g*sqrt(gamma1*rho/P) at surface + !acoustic_radius ! integral of dr/csound (seconds) + !ng_for_nu_max ! = 1 / (nu_max*delta_Pg) + ! period for g-mode with frequency nu_max = nu_max_ng*delta_Pg + !gs_per_delta_nu ! delta_nu / (nu_max**2*delta_Pg) + ! number of g-modes per delta_nu at nu_max + + !int_k_r_dr_nu_max_Sl1 ! integral of k_r*dr where nu < N < Sl for nu = nu_max, l=1 + !int_k_r_dr_2pt0_nu_max_Sl1 ! integral of k_r*dr where nu < N < Sl for nu = nu_max*2, l=1 + !int_k_r_dr_0pt5_nu_max_Sl1 ! integral of k_r*dr where nu < N < Sl for nu = nu_max/2, l=1 + !int_k_r_dr_nu_max_Sl2 ! integral of k_r*dr where nu < N < Sl for nu = nu_max, l=2 + !int_k_r_dr_2pt0_nu_max_Sl2 ! integral of k_r*dr where nu < N < Sl for nu = nu_max*2, l=2 + !int_k_r_dr_0pt5_nu_max_Sl2 ! integral of k_r*dr where nu < N < Sl for nu = nu_max/2, l=2 + !int_k_r_dr_nu_max_Sl3 ! integral of k_r*dr where nu < N < Sl for nu = nu_max, l=3 + !int_k_r_dr_2pt0_nu_max_Sl3 ! integral of k_r*dr where nu < N < Sl for nu = nu_max*2, l=3 + !int_k_r_dr_0pt5_nu_max_Sl3 ! integral of k_r*dr where nu < N < Sl for nu = nu_max/2, l=3 + +!---------------------------------------------------------------------------------------------- + +!# energy information + + !total_energy ! at end of step + !log_total_energy ! log(abs(total_energy)) + total_energy_after_adjust_mass ! after mass adjustments + + ! shorter versions of above + tot_E + log_tot_E + + + total_gravitational_energy + !log_total_gravitational_energy ! log(abs(total_gravitational_energy)) + !total_gravitational_energy_after_adjust_mass + + ! shorter versions of above + tot_PE + !log_tot_PE + + !total_internal_energy + !log_total_internal_energy + !total_internal_energy_after_adjust_mass + + ! shorter versions of above + tot_IE + !log_tot_IE + + !total_radial_kinetic_energy + !log_total_radial_kinetic_energy + !total_radial_kinetic_energy_after_adjust_mass + + ! shorter versions of above (does not include rot KE) + tot_KE + !log_tot_KE + + !total_turbulent_energy + !log_total_turbulent_energy + !total_turbulent_energy_after_adjust_mass + tot_Et + !log_tot_Et + + !total_energy_foe + + !tot_IE_div_IE_plus_KE + !total_IE_div_IE_plus_KE + + !total_entropy + total_eps_grav + + total_energy_sources_and_sinks ! for this step + !total_nuclear_heating + !total_non_nuc_neu_cooling + !total_irradiation_heating + total_extra_heating ! extra heat integrated over the model times dt (erg) + !total_WD_sedimentation_heating + + !rel_run_E_err + + !rel_E_err + !abs_rel_E_err + !log_rel_E_err + + !tot_e_equ_err + tot_e_err + + + error_in_energy_conservation ! for this step + ! = total_energy - (total_energy_start + total_energy_sources_and_sinks) + cumulative_energy_error ! = sum over all steps of error_in_energy_conservation + !rel_cumulative_energy_error ! = cumulative_energy_error/total_energy + !log_rel_cumulative_energy_error ! = log10 of rel_cumulative_energy_error + !log_rel_run_E_err ! shorter name for rel_cumulative_energy_error + + !rel_error_in_energy_conservation ! = error_in_energy_conservation/total_energy + !log_rel_error_in_energy_conservation + + !virial_thm_P_avg + !virial_thm_rel_err + !work_inward_at_center + !work_outward_at_surface + + +!---------------------------------------------------------------------------------------------- + + !# rotation + + !total_angular_momentum + !log_total_angular_momentum + !i_rot_total ! moment of inertia + + !total_rotational_kinetic_energy + !log_total_rotational_kinetic_energy + !total_rotational_kinetic_energy_after_adjust_mass + +!---------------------------------------------------------------------------------------------- + +!# velocities + + !avg_abs_v_div_cs + !log_avg_abs_v_div_cs + !max_abs_v_div_cs + !log_max_abs_v_div_cs + + !avg_abs_v + !log_avg_abs_v + !max_abs_v + !log_max_abs_v + + !u_surf + !u_surf_km_s + !u_div_csound_surf + !u_div_csound_max + + !infall_div_cs + +!---------------------------------------------------------------------------------------------- + +!# misc + + !e_thermal ! sum over all zones of Cp*T*dm + + !## eos + ! crystal_core_boundary_mass + !logQ_max ! logQ = logRho - 2*logT + 12 + !logQ_min + !gamma1_min + + !## core mixing + !mass_semiconv_core + + !## H-He boundary + + !diffusion_time_H_He_bdy + !temperature_H_He_bdy + + + !## optical depth and opacity + + !one_div_yphot + !log_one_div_yphot + + !log_min_opacity + !min_opacity + + !log_tau_center + + !log_max_tau_conv + !max_tau_conv + !log_min_tau_conv + !min_tau_conv + + !tau_qhse_yrs + + !## other + + !Lsurf_m + !dlnR_dlnM + !h1_czb_mass ! location (in Msun units) of base of 1st convection zone above he core + !kh_mdot_limit + !log_cntr_dr_cm + !min_Pgas_div_P + !surf_c12_minus_o16 ! this is useful for seeing effects of dredge up on AGB + !surf_num_c12_div_num_o16 + + !phase_of_evolution ! Integer mapping to the type of evolution see star_data/public/star_data_def.inc for definitions + + !## MLT++ + !gradT_excess_alpha + !gradT_excess_min_beta + !gradT_excess_max_lambda + + !max_L_rad_div_Ledd + !max_L_rad_div_Ledd_div_phi_Joss + + + !## RTI + !rti_regions + + !## Ni & Co + !total_ni_co_56 + + + !## internal structure constants + + ! this is evaluated assuming a spherical star and does not account for rotation + !apsidal_constant_k2 + + +!---------------------------------------------------------------------------------------------- + +!# accretion + + !k_below_const_q + !q_below_const_q + !logxq_below_const_q + + !k_const_mass + !q_const_mass + !logxq_const_mass + + !k_below_just_added + !q_below_just_added + !logxq_below_just_added + + !k_for_test_CpT_absMdot_div_L + !q_for_test_CpT_absMdot_div_L + !logxq_for_test_CpT_absMdot_div_L + +!---------------------------------------------------------------------------------------------- + +!# Color output + + ! Outputs the bolometric correction (bc) for the star in filter band ``filter'' (case sensitive) + !bc filter + + ! Outputs the absolute magnitude for the star in filter band ``filter'' (case sensitive) + !abs_mag filter + + ! Adds all the bc's to the output + !add_bc + + ! Adds all the absolute magnitudes to the output + !add_abs_mag + + ! Outputs luminosity in filter band ``filter'' (erg s^-1) (case sensitive) + ! lum_band filter + + ! Adds all the filter band luminosities to the output (erg s^-1) + ! add_lum_band + + ! Outputs log luminosity in filter band ``filter'' (log erg s^-1) (case sensitive) + ! log_lum_band filter + + ! Adds all the filter band luminosities to the output (log erg s^-1) + ! add_log_lum_band + +!---------------------------------------------------------------------------------------------- + +!# RSP + + !rsp_DeltaMag ! absolute magnitude difference between minimum and maximum light (mag) + !rsp_DeltaR ! R_max - R_min difference in the max and min radius (Rsun) + !rsp_GREKM ! fractional growth of the kinetic energy per pulsation period ("nonlinear growth rate") - see equation 5 in MESA5 + !rsp_num_periods ! Count of the number of pulsation cycles completed + !rsp_period_in_days ! Running period, ie., period between two consecutive values of R_max (days) + !rsp_phase ! Running pulsation phase for a cycle + +!---------------------------------------------------------------------------------------------- +!# debugging + + !## retries + num_retries ! total during the run + + !## solver iterations + + num_iters ! same as num_solver_iterations + !num_solver_iterations ! iterations at this step + !total_num_solver_iterations ! total iterations during the run + !avg_num_solver_iters + + !rotation_solver_steps + + !diffusion_solver_steps + !diffusion_solver_iters + + !avg_setvars_per_step + !avg_skipped_setvars_per_step + !avg_solver_setvars_per_step + + !burn_solver_maxsteps + + !total_num_solver_calls_converged + !total_num_solver_calls_failed + !total_num_solver_calls_made + !total_num_solver_relax_calls_converged + !total_num_solver_relax_calls_failed + !total_num_solver_relax_calls_made + !total_num_solver_relax_iterations + + !total_step_attempts + !total_step_redos + !total_step_retries + !total_steps_finished + !total_steps_taken + + !TDC_num_cells + + !## Relaxation steps + !total_relax_step_attempts + !total_relax_step_redos + !total_relax_step_retries + !total_relax_steps_finished + !total_relax_steps_taken + + !## conservation during mesh adjust + !log_mesh_adjust_IE_conservation + !log_mesh_adjust_KE_conservation + !log_mesh_adjust_PE_conservation + + !## amr + !num_hydro_merges + !num_hydro_splits + + !## timing + !elapsed_time ! time since start of run (seconds) diff --git a/star/test_suite/wd_o_ne_3_phase/inlist b/star/test_suite/wd_o_ne_3_phase/inlist new file mode 100644 index 000000000..460f300fb --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/inlist @@ -0,0 +1,228 @@ + +! inlist for cooling based on wd_cool_0.6M + + +&star_job + show_log_description_at_start = .false. + + load_saved_model = .true. + load_model_filename = 'wd1_10_mi8_3_ONeNa.mod'!'wd1_10_mi8_3_NeOMg.mod' + + set_initial_age = .true. + initial_age = 0d0 + set_initial_model_number = .true. + initial_model_number = 0 + + pgstar_flag = .true. + save_pgstar_files_when_terminate = .false. + + change_net = .true. + new_net_name = 'custom.net' + change_initial_net = .true. + + + + +/ ! end of star_job namelist + + +&eos + use_Skye = .true. + use_PC = .false. + + ! For more accurate C/O phase diagram + Skye_solid_mixing_rule = "Ogata" + +/ ! end of eos namelist + + + +&kap + Zbase = 0.02d0 + + kap_file_prefix = 'gs98' + use_Type2_opacities = .true. + include_electron_conduction = .true. + use_blouin_conductive_opacities = .true. + +/ ! end of kap name-list + +&controls + +! element diffusion + + do_element_diffusion = .true. + diffusion_use_full_net = .true. + + ! uncomment to add a He atmosphere + + ! accrete_same_as_surface = .true. + ! accrete_given_mass_fractions = .true. + ! num_accretion_species = 1 + ! accretion_species_xa(1) = 1.0 + ! accretion_species_id(1) = 'he4' + ! mass_change = 1d-8 + + !star_species_mass_max_limit = 1d-4 ! Msun + !star_species_mass_max_limit_iso = 'he4' + +! solver + + !to avoid problem with heating from crystallization at small timesteps + max_corr_jump_limit = 1d7!! + convergence_ignore_equL_residuals = .true.!! + convergence_ignore_alpha_RTI_residuals = .true.!! + make_gradr_sticky_in_solver_iters = .true. !! + hydro_mtx_min_allowed_logT = -10 + hydro_mtx_max_allowed_abs_dlogT = 99d0 + hydro_mtx_max_allowed_logT = 99d0 + tol_correction_high_T_limit = 1d6 + tol_correction_norm_high_T = 3d6 + tol_max_correction_high_T = 3d6 + + num_trace_history_values = 2 + trace_history_value_name(1) = 'rel_E_err' + trace_history_value_name(2) = 'log_eel_run_e_err' + + energy_eqn_option = 'eps_grav' + use_time_centered_eps_grav = .true. + + limit_for_rel_error_in_energy_conservation = 1d-4 + hard_limit_for_rel_error_in_energy_conservation = 1d-3 + + + +!phase separation + + do_phase_separation = .true. + do_phase_separation_heating = .true. + phase_separation_mixing_use_brunt = .true. + phase_separation_option = '3c'!'ONe' + + +!stop model + + !gamma_center_limit = 197.5d0! 190.5d0 ! Stop at the onset of crystallization using gamma center limit + !max_years_for_timestep = 1d6 ! Reduce timestep to 1d4 at the onset for resolving the core size better + + +! atm + + atm_option = 'table' + atm_table = 'DB_WD_tau_25'! 'photosphere' ! + atm_off_table_option = 'T_tau' + + + +! mesh + + varcontrol_target = 1d-4 + mesh_delta_coeff = 1.! 1.5 + +! mlt + + use_Ledoux_criterion = .true. + + thermohaline_coeff =3!1000 for a smoother profile before the onset of crystallization + thermohaline_option = 'Kippenhahn' + + MLT_option = 'ML2' + + mixing_length_alpha = 1.8 + + +! from wd_cool_0.6M test_suite model + + clip_D_limit = 0 ! zero mixing diffusion coeffs that are smaller than this + + prune_bad_cz_min_Hp_height = 0 ! lower limit on radial extent of cz + remove_mixing_glitches = .true. ! if true, then okay to remove gaps and singletons + + ! the following controls are for different kinds of "glitches" that can be removed + + + okay_to_remove_mixing_singleton = .true. + + min_convective_gap = 0 + ! close gap between convective regions if smaller than this (< 0 means skip this) + ! gap measured radially in units of pressure scale height + + min_thermohaline_gap = 0 + ! close gap between thermohaline mixing regions if smaller than this (< 0 means skip this) + ! gap measured radially in units of pressure scale height + + min_thermohaline_dropout = 0 + max_dropout_gradL_sub_grada = 1d-3 + ! if find radiative region embedded in thermohaline, + ! and max(gradL - grada) in region is everywhere < max_dropout_gradL_sub_grada + ! and region height is < min_thermohaline_dropout + ! then convert the region to thermohaline + + min_semiconvection_gap = 0 + ! close gap between semiconvective mixing regions if smaller than this (< 0 means skip this) + ! gap measured radially in units of pressure scale height + + remove_embedded_semiconvection = .false. + ! if have a semiconvection region bounded on each side by convection, + ! convert it to be convective too. + + + + photo_interval = 1 + profile_interval = 1 !10 + history_interval = 1 !10 + terminal_interval = 1 + write_header_frequency = 10 + log_directory = 'LOGS_110_Na_3c' + max_num_profile_models = 15000 + + +/ ! end of controls namelist + + +&pgstar + + file_white_on_black_flag = .true. ! white_on_black flags -- true means white foreground color on black background + file_device = 'png' ! png + file_extension = 'png' + + !file_device = 'vcps' ! postscript + !file_extension = 'ps' + + !Grid2_win_flag = .true. + Grid2_win_width = 20 + + !Grid2_file_flag = .true. + Grid2_file_width = 30 + !Grid2_file_interval = 10000 + + !Abundance_file_flag = .true. + + + TRho_Profile_win_flag = .true. + show_TRho_Profile_kap_regions = .false. + show_TRho_Profile_eos_regions = .true. + !show_TRho_Profile_degeneracy_line = .true. + !show_TRho_Profile_Pgas_Prad_line = .true. + !show_TRho_Profile_burn_lin.2es = .true. + show_TRho_Profile_burn_labels = .true. + show_TRho_Profile_text_info = .false. + + Abundance_win_flag = .true. + Abundance_xaxis_name = 'logxm' + Abundance_xaxis_reversed = .true. + Abundance_xmin = -7 + Abundance_xmax = 0.5 + Abundance_log_mass_frac_min = -3 + Abundance_log_mass_frac_max = 0 + + ! Mixing_win_flag = .true. + ! Mixing_win_width = 6 + ! Mixing_win_aspect_ratio = 0.75 + + + +/ ! end of pgstar namelist + + + diff --git a/star/test_suite/wd_o_ne_3_phase/inlist_create_wd b/star/test_suite/wd_o_ne_3_phase/inlist_create_wd new file mode 100644 index 000000000..43d5e3331 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/inlist_create_wd @@ -0,0 +1,155 @@ +! inlist to create a white dwarf based on Lauffer et al. (2018) + + +&star_job +show_log_description_at_start = .false. +save_model_when_terminate = .true. +save_model_filename = 'wd1_10_mi8_3_ONeNa.mod' + +change_initial_net = .true. +new_net_name = 'sagb_NeNa_MgAl.net' ! use co_burn_plus.net for Ne/O/Mg core +pgstar_flag = .true. +!steps_to_take_before_terminate = 1 + +/ ! end of star_job namelist +&eos + +!eos_file_prefix = 'mesa' +/ ! end of eos namelist + +&kap +Zbase = 0.02d0 +kap_file_prefix = 'gs98' +use_Type2_opacities = .true. +/ ! end of kap name-list + + +&controls + +initial_mass = 8.3! +initial_z = 0.02d0 + +! terminal output +photo_interval = 50 +profile_interval = 500 +history_interval = 100 +terminal_interval = 10 +write_header_frequency = 10 + +max_age = 10d9 +max_years_for_timestep = 1d9 +varcontrol_target = 1d-3 +dX_nuc_drop_limit = 1d-2 +mesh_delta_coeff = 2.5 +okay_to_reduce_gradT_excess = .true. + +! Stop when not more nuclear burning +power_nuc_burn_lower_limit = 1d-8 + +! convection +mixing_length_alpha = 2 +use_Ledoux_criterion = .true. +alpha_semiconvection = 0.01 +thermohaline_coeff = 3! + +! diffusion +do_element_diffusion = .true. +am_nu_visc_factor = 0 +am_D_mix_factor = 0.0333333333333333d0 +D_DSI_factor = 0 +D_SH_factor = 1 +D_SSI_factor = 1 +D_ES_factor = 1 +D_GSF_factor = 1 +D_ST_factor = 1 + +! mass Loss +cool_wind_full_on_T = 1d6 +hot_wind_full_on_T = 1d6 +cool_wind_RGB_scheme = 'Reimers' +cool_wind_AGB_scheme = 'Blocker' +RGB_to_AGB_wind_switch = 1d-4 +Reimers_scaling_factor = 0.5!0.1d0 !decrease if too much +Blocker_scaling_factor = 50d0!10d0 + +! stop if almost not He +star_species_mass_min_limit = 1d-20 ! Msun +star_species_mass_min_limit_iso = 'he4' + + +!overshoot + overshoot_scheme(1) = 'exponential' + overshoot_zone_type(1) = 'nonburn' + overshoot_zone_loc(1) = 'any' + overshoot_bdy_loc(1) = 'top'!'any'! + overshoot_f(1) = 0.01 + overshoot_f0(1) = 0.005 + + overshoot_scheme(2) = 'exponential' + overshoot_zone_type(2) = 'burn_H' + overshoot_zone_loc(2) = 'any' + overshoot_bdy_loc(2) = 'top'!'any' !'top' + overshoot_f(2) = 0.01 + overshoot_f0(2) = 0.005 + + overshoot_scheme(3) = 'exponential' + overshoot_zone_type(3) = 'burn_He' + overshoot_zone_loc(3) = 'any' + overshoot_bdy_loc(3) = 'top'!'any'!'top' + overshoot_f(3) = 0.01 + overshoot_f0(3) = 0.005 + + overshoot_scheme(4) = 'exponential' + overshoot_zone_type(4) = 'burn_Z' + overshoot_zone_loc(4) = 'shell' + overshoot_bdy_loc(4) = 'any' + overshoot_f(4) = 0.1 !0.5 increase if having problems to converge + overshoot_f0(4) = 0.01 !0.1 +/ ! end of controls namelist + +&pgstar + + ! show HR diagram + ! this plots the history of L,Teff over many timesteps + HR_win_flag = .true. + + ! set static plot bounds + HR_logT_min = 3.5 + HR_logT_max = 4.6 + HR_logL_min = 2.0 + HR_logL_max = 6.0 + + ! set window size (aspect_ratio = height/width) + HR_win_width = 6 + HR_win_aspect_ratio = 1.0 + + + ! show temperature/density profile + ! this plots the internal structure at single timestep + TRho_Profile_win_flag = .true. + + ! add legend explaining colors + show_TRho_Profile_legend = .true. + + ! display numerical info about the star + show_TRho_Profile_text_info = .true. + + ! set window size (aspect_ratio = height/width) + TRho_Profile_win_width = 8 + TRho_Profile_win_aspect_ratio = 0.75 + Abundance_win_flag = .true. + Abundance_xaxis_name = 'logxm' + Abundance_xaxis_reversed = .true. + Abundance_xmin = -4 + Abundance_xmax = 0.5 + Abundance_log_mass_frac_min = -4 + Abundance_log_mass_frac_max = -0.0 + + Mixing_win_flag = .true. + Mixing_win_width = 6 + Mixing_win_aspect_ratio = 0.75 + + + +/ ! end of pgstar namelist + diff --git a/star/test_suite/wd_o_ne_3_phase/make/makefile b/star/test_suite/wd_o_ne_3_phase/make/makefile new file mode 100644 index 000000000..89aa01f2c --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/make/makefile @@ -0,0 +1,9 @@ +ifeq ($(MESA_DIR),) +ifeq ($($MESA_DIR_INTENTIONALLY_EMPTY),) + $(error MESA_DIR enviroment variable is not set) +endif +endif + +STAR = star + +include $(MESA_DIR)/star/work_standard_makefile diff --git a/star/test_suite/wd_o_ne_3_phase/mk b/star/test_suite/wd_o_ne_3_phase/mk new file mode 100755 index 000000000..aec7a5195 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/mk @@ -0,0 +1,13 @@ +#!/bin/bash + +function check_okay { + if [ $? -ne 0 ] + then + echo + echo "FAILED" + echo + exit 1 + fi +} + +cd make; make; check_okay diff --git a/star/test_suite/wd_o_ne_3_phase/profile_columns.list b/star/test_suite/wd_o_ne_3_phase/profile_columns.list new file mode 100644 index 000000000..c05dc23a8 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/profile_columns.list @@ -0,0 +1,949 @@ +! profile_columns.list -- determines the contents of star model profiles +! you can use a non-standard version by setting profile_columns_file in your inlist + +! units are cgs unless otherwise noted. + +! reorder the following names as desired to reorder columns. +! comment out the name to omit a column (fewer columns => less IO => faster running). +! remove '!' to restore a column. + +! if you have a situation where you want a non-standard set of columns, +! make a copy of this file, edit as desired, and give the new filename in your inlist +! as profile_columns_file. if you are just adding columns, you can 'include' this file, +! and just list the additions in your file. note: to include the standard default +! version, use include '' -- the 0 length string means include the default file. + +! if you need to have something added to the list of options, let me know.... + +! the first few lines of the profile contain general info about the model. +! for completeness, those items are described at the end of this file. + + +! note: you can include another list by doing +! include 'filename' +! include '' means include the default standard list file + + +! the following lines of the profile contain info for 1 zone per row, surface to center. + +! minimal set of enabled columns: + + zone ! numbers start with 1 at the surface + mass ! m/Msun. mass coordinate of outer boundary of cell. + logR ! log10(radius/Rsun) at outer boundary of zone + logT ! log10(temperature) at center of zone + logRho ! log10(density) at center of zone + logP ! log10(pressure) at center of zone + x_mass_fraction_H + y_mass_fraction_He + z_mass_fraction_metals + + +! everything below this line is deactivated + + +!# Structure + !logM ! log10(m/Msun) + !log_mass + !dm ! cell mass (grams) + !dm_bar ! boundary mass (grams) average of adjacent dm's + logdq ! log10(dq) + !log_dq + dq_ratio ! dq(k-1)/dq(k) + q ! fraction of star mass interior to outer boundary of this zone + !log_q ! log10(q) + !xq + + grav ! gravitational acceleration (cm sec^2) + !log_g ! log10 gravitational acceleration (cm sec^2) + !g_div_r ! grav/radius (sec^2) + !r_div_g ! radius/grav (sec^-2) + !cgrav_factor ! = cgrav(k)/standard_cgrav + !vel_km_per_s ! velocity at outer boundary of zone (km/s) -- 0 if no velocity variable + + radius ! radius at outer boundary of zone (in Rsun units) + !radius_cm ! radius at outer boundary of zone (in centimeters) + !radius_km ! radius at outer boundary of zone (in kilometers) + !logR_cm ! log10 radius at outer boundary of zone (in centimeters) + !rmid ! radius at center by mass of zone (in Rsun units) + !r_div_R ! fraction of total radius + + velocity ! velocity at outer boundary of zone (cm/s) -- 0 if no velocity variable + !v_div_r ! velocity divided by radius + !v_times_t_div_r + !rho_times_r3 ! at face + !log_rho_times_r3 ! at face + !scale_height ! in Rsun units + pressure_scale_height ! in Rsun units + + !m_div_r ! gm/cm + !dmbar_m_div_r + !log_dmbar_m_div_r + !mass_grams ! mass coordinate of outer boundary of cell in grams + mmid ! mass at midpoint of cell (average of mass coords of the cell boundaries) Msun units. + + !m_grav ! total enclosed gravitational mass. Msun units. + !m_grav_div_m_baryonic ! mass_gravitational/mass at cell boundary + !mass_correction_factor ! dm_gravitational/dm (dm is baryonic mass of cell) + + !xm ! mass exterior to point (Msun units) + !dq ! mass of zone as a fraction of total star mass + logxq ! log10(1-q) + logxm ! log10(xm) + + !xr ! radial distance from point to surface (Rsun) + !xr_cm ! radial distance from point to surface (cm) + !xr_div_R ! radial distance from point to surface in units of star radius + !log_xr ! log10 radial distance from point to surface (Rsun) + !log_xr_cm ! log10 radial distance from point to surface (cm) + !log_xr_div_R ! log10 radial distance from point to surface in units of star radius + + !dr ! r(outer edge) - r(inner edge); radial extent of cell in cm. + !log_dr ! log10 cell width (cm) + !dv ! v(inner edge) - v(outer edge); rate at which delta_r is shrinking (cm/sec). + + !dt_dv_div_dr ! dt*dv/dr; need to have this << 1 for every cell + !dr_div_R ! cell width divided by star R + !log_dr_div_R ! log10 cell width divided by star R + !dr_div_rmid ! cell width divided by rmid + !log_dr_div_rmid ! log(dr_div_rmid) + + !dr_div_cs ! cell sound crossing time (sec) + !log_dr_div_cs ! log10 cell sound crossing time (sec) + !dr_div_cs_yr ! cell sound crossing time (years) + !log_dr_div_cs_yr ! log10 cell sound crossing time (years) + + !acoustic_radius ! sound time from center to outer cell boundary (sec) + !log_acoustic_radius ! log10(acoustic_radius) (sec) + !acoustic_depth ! sound time from surface to outer cell boundary (sec) + !log_acoustic_depth ! log10(acoustic_depth) (sec) + !acoustic_r_div_R_phot + + !cell_collapse_time ! only set if doing explicit hydro + ! time (seconds) for cell inner edge to catch cell outer edge at current velocities + ! 0 if distance between inner and outer is increasing + !log_cell_collapse_time ! log of cell_collapse_time + + !compression_gradient + + + +!# Thermodynamics + temperature ! temperature at center of zone + !logT_face ! log10(temperature) at outer boundary of zone + !logT_bb ! log10(black body temperature) at outer boundary of zone + !logT_face_div_logT_bb + + !energy ! internal energy (ergs/g) + !logE ! log10(specific internal energy) at center of zone + rho ! density + !density ! rho + + entropy ! specific entropy divided by (avo*kerg) + !logS ! log10(specific entropy) + !logS_per_baryon ! log10(specific entropy per baryon / kerg) + + pressure ! total pressure at center of zone (pgas + prad) + !prad ! radiation pressure at center of zone + !pgas ! gas pressure at center of zone (electrons and ions) + logPgas ! log10(pgas) + pgas_div_ptotal ! pgas/pressure + + eta ! electron degeneracy parameter (eta >> 1 for significant degeneracy) + mu ! mean molecular weight per gas particle (ions + free electrons) + + grada ! dlnT_dlnP at constant S + !dE_dRho ! at constant T + !cv ! specific heat at constant volume + cp ! specific heat at constant total pressure + + !log_CpT + !gamma1 ! dlnP_dlnRho at constant S + !gamma3 ! gamma3 - 1 = dlnT_dlnRho at constant S + gam ! plasma interaction parameter (> 160 or so means starting crystallization) + free_e ! free_e is mean number of free electrons per nucleon + !logfree_e ! log10(free_e), free_e is mean number of free electrons per nucleon + chiRho ! dlnP_dlnRho at constant T + chiT ! dlnP_dlnT at constant Rho + ! chiX + + csound ! sound speed + !log_csound + !csound_face ! sound speed (was previously called csound_at_face) + cs_at_cell_bdy ! sound speed at cell boundary (csound is at cell center) + !v_div_cs ! velocity divided by sound speed + v_div_csound ! velocity divided by sound speed + !div_v + + !thermal_time_to_surface ! in seconds + !log_thermal_time_to_surface + !t_rad + !log_t_rad + !log_t_sound + log_t_thermal + + eos_phase + !eos_frac_OPAL_SCVH + !eos_frac_HELM + eos_frac_Skye + !eos_frac_PC + !eos_frac_FreeEOS + !eos_frac_CMS + !eos_frac_ideal + + !pgas_div_p + !prad_div_pgas + !prad_div_pgas_div_L_div_Ledd + pressure_scale_height_cm + + eps_grav_composition_term + !eps_grav_plus_eps_mdot + + !chiRho_for_partials + !chiT_for_partials + !rel_diff_chiRho_for_partials + !rel_diff_chiT_for_partials + + latent_ddlnRho + latent_ddlnT + + !log_P_face + !log_Ptrb + !log_cp_T_div_t_sound + + !QQ + + +!# Mass accretion + eps_grav ! -T*ds/dt (negative for expansion) + !log_abs_eps_grav_dm_div_L + !log_abs_v ! log10(abs(velocity)) (cm/s) + !log_mdot_cs + !log_mdot_v + !eps_mdot + !env_eps_grav + !xm_div_delta_m + !log_xm_div_delta_m + + +!# Nuclear energy generation + signed_log_eps_grav ! sign(eps_grav)*log10(max(1,abs(eps_grav))) + !signed_log_eps_nuc + !net_nuclear_energy ! erg/gm/s from nuclear reactions minus all neutrino losses + ! The value plotted is net_nuclear_energy = sign(val)*log10(max(1,abs(val))) + ! where val = net nuclear energy minus all neutrino losses. + net_energy ! net_energy + eps_grav. + ! The value plotted is net_energy = sign(val)*log10(max(1,abs(val))) + ! where val = net nuclear energy plus eps_grav minus all neutrino losses. + !eps_nuc_plus_nuc_neu + !eps_nuc_minus_non_nuc_neu + !eps_nuc_start + + eps_nuc ! ergs/g/sec from nuclear reactions (including losses to reaction neutrinos) + !log_abs_eps_nuc + !d_lnepsnuc_dlnd + !d_epsnuc_dlnd + !deps_dlnd_face + ! (was previously called deps_dlnd_at_face) + !d_lnepsnuc_dlnT + !d_epsnuc_dlnT + !deps_dlnT_face + ! (was previously called deps_dlnT_at_face) + !eps_nuc_neu_total ! erg/gm/sec as neutrinos from nuclear reactions + + non_nuc_neu ! non-nuclear-reaction neutrino losses + !nonnucneu_plas ! plasmon neutrinos (for collective reactions like gamma_plasmon => nu_e + nubar_e) + !nonnucneu_brem ! bremsstrahlung (for reactions like e- + (z,a) => e- + (z,a) + nu + nubar) + !nonnucneu_phot ! photon neutrinos (for reactions like e- + gamma => e- + nu_e + nubar_e) + !nonnucneu_pair ! pair production (for reactions like e+ + e- => nu_e + nubar_e) + !nonnucneu_reco ! recombination neutrinos (for reactions like e- (continuum) => e- (bound) + nu_e + nubar_e) + + ! ergs/g/sec for reaction categories + add_reaction_categories ! this adds all the reaction categories + ! NOTE: you can list specific categories by giving their names (from chem_def) + pp + cno + !tri_alfa + !c_alpha + !n_alpha + !o_alpha + !ne_alpha + !na_alpha + !mg_alpha + !si_alpha + !s_alpha + !ar_alpha + !ca_alpha + !ti_alpha + !fe_co_ni + !c12_c12 + !c12_o16 + !o16_o16 + !photo + !pnhe4 + !other + + !burn_num_iters ! Number of split_burn iterations taken + !burn_avg_epsnuc + !log_burn_avg_epsnuc + +!# Composition + !x_mass_fraction_H + !y_mass_fraction_He + !z_mass_fraction_metals + abar ! average atomic weight (g/mole) + zbar ! average charge + !z2bar ! average charge^2 + ye ! average charge per baryon = proton fraction + + x ! hydrogen mass fraction + log_x + !y ! helium mass fraction + !log_y + z ! metallicity + !log_z ! metallicity + + add_abundances ! this adds all of the isos that are in the current net + ! NOTE: you can list specific isotopes by giving their names (from chem_def) + h1 + he3 + he4 + c12 + !n14 + o16 + + !add_log_abundances ! this adds log10 of all of the isos that are in the current net + ! NOTE: you can list specific isotopes by giving their names (from chem_def) + !log h1 + !log he3 + !log he4 + !log c12 + !log n14 + !log o16 + + ! log concentration of species + ! concentration = number density / number density of electrons + ! Ci = (Xi/Ai) / sum(Zi*Xi/Ai) [see Thoul et al, ApJ 421:828-842, 1994] + !log_concentration h1 + !log_concentration he4 + + + ! typical charge for given species + ! (used by diffusion) + !typical_charge he4 + typical_charge c12 + !typical_charge fe52 + typical_charge o16 + + ! ionization state for given species + ! (same as typical charge, except that it's unsmoothed) + !ionization he4 + !ionization c12 + !ionization fe52 + + !cno_div_z ! abundance of c12, n14, and o16 as a fraction of total z + + + + +!# Opacity + opacity ! opacity measured at center of zone + log_opacity ! log10(opacity) + !dkap_dlnrho_face ! partial derivative of opacity wrt. ln rho (at T=const) at outer edge of cell + ! (was previously called dkap_dlnrho_at_face) + !dkap_dlnT_face ! partial derivative of opacity wrt. ln T (at rho=const) at outer edge of cell + ! (was previously called dkap_dlnT_at_face) + !kap_frac_lowT ! fraction of opacity from lowT tables + !kap_frac_highT ! fraction of opacity from highT tables + !kap_frac_Type2 ! fraction of opacity from Type2 tables + !kap_frac_Compton ! fraction of opacity from Compton_Opacity + !kap_frac_op_mono ! fraction of opacity from OP mono + + log_kap + !log_kap_times_factor + + !log_c_div_tau + !xtau + !xlogtau + !logtau_sub_xlogtau + +!# Luminosity + luminosity ! luminosity at outer boundary of zone (in Lsun units) + !logL ! log10(max(1d-2,L/Lsun)) + !log_Lrad + !log_Ledd ! log10(Leddington/Lsun) -- local Ledd, 4 pi clight G m / kap + !log_L_div_Ledd ! log10(max(1d-12,L/Leddington)) + !log_Lrad_div_Ledd + !log_Lrad_div_L + signed_log_power ! sign(L)*log10(max(1,abs(L))) + + lum_adv + lum_conv + lum_conv_MLT + !lum_div_Ledd + lum_erg_s + !lum_plus_lum_adv + lum_rad + + !log_L_div_CpTMdot + !log_abs_lum_erg_s + + !L + !Lc + !Lc_div_L + !Lr + !Lr_div_L + !Lt + !Lt_div_L + +!# Energetics + total_energy ! specific total energy of cell (ergs/g). internal+potential+kinetic+rotation. + !cell_specific_IE + !cell_specific_KE + !cell_IE_div_IE_plus_KE + !cell_KE_div_IE_plus_KE + + !cell_ie_div_star_ie + !cell_internal_energy_fraction + !cell_internal_energy_fraction_start + !cell_specific_PE + + !log_cell_ie_div_star_ie + !log_cell_specific_IE + + !ergs_eps_grav_plus_eps_mdot + !ergs_error + !ergs_error_integral + !ergs_mdot + !ergs_rel_error_integral + !dm_eps_grav + + !dE + + !etrb + !log_etrb + !extra_grav + !log_rel_E_err + + !total_energy_sign + +!# Convection + !mlt_mixing_length ! mixing length for mlt (cm) + !mlt_mixing_type ! value returned by mlt + !mlt_Pturb + !alpha_mlt + + !conv_vel ! convection velocity (cm/sec) + !log_conv_vel ! log10 convection velocity (cm/sec) + + !conv_L_div_L + !log_conv_L_div_L + !lum_conv_div_lum_rad + !lum_rad_div_L_Edd + !lum_conv_div_lum_Edd + !lum_conv_div_L + !lum_rad_div_L + !lum_rad_div_L_Edd_sub_fourPrad_div_PchiT ! density increases outward if this is > 0 + ! see Joss, Salpeter, and Ostriker, "Critical Luminosity", ApJ 181:429-438, 1973. + + gradT ! mlt value for required temperature gradient dlnT/dlnP + + gradr ! dlnT/dlnP required for purely radiative transport + !grad_temperature ! smoothed dlnT/dlnP at cell boundary + !grad_density ! smoothed dlnRho/dlnP at cell boundary + + gradL ! gradient for Ledoux criterion for convection + !sch_stable ! 1 if grada > gradr, 0 otherwise + !ledoux_stable ! 1 if gradL > gradr, 0 otherwise + + !grada_sub_gradT + !gradT_sub_grada ! gradT-grada at cell boundary + !gradT_div_grada ! gradT/grada at cell boundary + + !gradr_sub_gradT + !gradT_sub_gradr ! gradT-gradr at cell boundary + !gradT_div_gradr ! gradT/gradr at cell boundary + + !log_gradT_div_gradr ! log10 gradT/gradr at cell boundary + !log_mlt_Gamma ! convective efficiency + ! conv_vel_div_csound ! convection velocity divided by sound speed + !conv_vel_div_L_vel ! L_vel is velocity needed to carry L by convection; L = 4*pi*r^2*rho*vel**3 + !log_mlt_D_mix ! log10 diffusion coefficient for mixing from mlt (cm^2/sec) + + !gradr_div_grada ! gradr/grada_face; > 1 => Schwarzschild unstable for convection + !gradr_sub_grada ! gradr - grada_face; > 0 => Schwarzschild unstable for convection + + !gradL_sub_gradr + !gradP_div_rho + !gradT_excess_effect + !gradT_rel_err + !gradT_sub_a + !grada_face + !grada_sub_gradr + !diff_grads + !log_diff_grads + + !mlt_D + !mlt_Gamma + !mlt_Y_face + !mlt_Zeta + !mlt_gradT + !mlt_log_abs_Y + !mlt_vc + !log_mlt_vc + + !superad_reduction_factor + !conv_vel_div_mlt_vc + + !log_Lconv + !log_Lconv_div_L + +!# Mixing + !mixing_type ! mixing types are defined in mesa/const/public/const_def + !log_D_mix ! log10 diffusion coefficient for mixing in units of cm^2/second (Eulerian) + !log_D_mix_non_rotation + !log_D_mix_rotation + + !log_D_conv ! D_mix for regions where mix_type = convective_mixing + !log_D_leftover ! D_mix for regions where mix_type = leftover_convective_mixing + !log_D_semi ! D_mix for regions where mix_type = semiconvective_mixing + !log_D_ovr ! D_mix for regions where mix_type = overshoot_mixing + !log_D_thrm ! D_mix for regions where mix_type = thermohaline_mixing + !log_D_minimum ! D_mix for regions where mix_type = minimum_mixing + !log_D_rayleigh_taylor ! D_mix for regions where mix_type = rayleigh_taylor_mixing + !log_D_anon ! D_mix for regions where mix_type = anonymous_mixing + !log_D_omega + + !log_sig_mix ! sig(k) is mixing flow across face k in (gm sec^1) + ! sig(k) = D_mix*(4*pi*r(k)**2*rho_face)**2/dmavg + + !dominant_isoA_for_thermohaline + !dominant_isoZ_for_thermohaline + !gradL_composition_term + + !mix_type + + + +!# Optical Depth + !tau ! optical depth + !log_column_depth ! log10 column depth, exterior mass / area (g cm^-2) + !log_radial_depth ! log10 radial distance to surface (cm) + logtau ! log10(optical depth) at center of zone + !tau_eff ! tau that gives the local P == P_atm if this location at surface + ! tau_eff = kap*(P/g - Pextra_factor*(L/M)/(6*pi*clight*cgrav)) + !tau_eff_div_tau + + + +!# Rotation + !omega ! angular velocity = j_rot/i_rot + !log_omega + !log_j_rot + !log_J_div_M53 ! J is j*1e-15 integrated from center; M53 is m^(5/3) + !log_J_inside ! J_inside is j_rot integrated from center + !shear ! -dlnomega/dlnR + !log_abs_shear ! log10(abs(dlnomega/dlnR)) + !richardson_number + !i_rot ! specific moment of inertia at cell boundary + !j_rot ! specific angular momentum at cell boundary + !v_rot ! rotation velocity at cell boundary (km/sec) + !w_div_w_crit_roche !ratio of rotational velocity to keplerian at the equator + !without the contribution from the Eddington factor + !fp_rot ! rotation factor for pressure + !ft_rot ! rotation factor for temperature + !ft_rot_div_fp_rot ! gradr factor + + !log_am_nu_non_rot ! log10(am_nu_non_rot) + !log_am_nu_rot ! log10(am_nu_rot) + !log_am_nu ! log10(am_nu_non_rot + am_nu_rot) + + !r_polar ! (Rsun) + !log_r_polar ! log10 (Rsun) + !r_equatorial ! (Rsun) + !log_r_equatorial ! log10 (Rsun) + !r_e_div_r_p ! equatorial/r_polar + !omega_crit ! breakup angular velocity = sqrt(G M / equatorial^3) + !omega_div_omega_crit + + !am_log_nu_omega ! for diffusion of omega + !am_log_nu_j ! for diffusion of angular momentum + + !am_log_nu_rot ! diffusion of angular momentum driven by rotation + !am_log_nu_non_rot ! diffusion driven by other sources, e.g. convection + + !am_log_sig_omega ! for diffusion of omega + !am_log_sig_j ! for diffusion of angular momentum + !am_log_sig ! == am_log_sig_omega + + !am_log_D_visc ! diffusion coeff for kinematic viscosity + !am_log_D_DSI ! diffusion coeff for dynamical shear instability + !am_log_D_SH ! diffusion coeff for Solberg-Hoiland instability + !am_log_D_SSI ! diffusion coeff for secular shear instability + !am_log_D_ES ! diffusion coeff for Eddington-Sweet circulation + !am_log_D_GSF ! diffusion coeff for Goldreich-Schubert-Fricke instability + !am_log_D_ST ! Spruit dynamo mixing diffusivity + !am_log_nu_ST ! Spruit dynamo effective viscosity + + !dynamo_log_B_r ! (Gauss) + !dynamo_log_B_phi ! (Gauss) + + !am_domega_dlnR + !log_abs_dlnR_domega + + !w_div_w_crit_roche2 + + +!# Diffusion + ! electric field from element diffusion calculation + !e_field + !log_e_field + + ! gravitational field from element diffusion calculation + !g_field_element_diffusion + !log_g_field_element_diffusion + + !eE_div_mg_element_diffusion + !log_eE_div_mg_element_diffusion + + ! element diffusion velocity for species + !edv h1 + !edv he4 + !edv o16 + + ! Energy generated by Ne22 sedimentation. + !eps_WD_sedimentation + !log_eps_WD_sedimentation + + !eps_diffusion + !log_eps_diffusion + + !diffusion_D h1 ! self diffusion coeff + !diffusion_dX h1 ! change in h1 mass fraction from diffusion + !diffusion_dX he4 ! change in he4 mass fraction from diffusion + !diffusion_dX n20 ! change in n20 mass fraction from diffusion + + !v_rad h1 ! velocity from radiative levitation + !v_rad he4 ! velocity from radiative levitation + !v_rad ne20 ! velocity from radiative levitation + + !log_g_rad h1 ! log10 acceleration from radiative levitation + !log_g_rad he4 ! log10 acceleration from radiative levitation + !log_g_rad ne20 ! log10 acceleration from radiative levitation + + + +!# Oscillations + brunt_N2 ! brunt-vaisala frequency squared + brunt_N2_structure_term + brunt_N2_composition_term + !log_brunt_N2_structure_term + !log_brunt_N2_composition_term + !brunt_A ! = N^2*r/g + !brunt_A_div_x2 ! x = r(k)/r(1) + !brunt_N2_dimensionless ! N2 in units of 3GM/R^3 + !brunt_N_dimensionless ! N in units of sqrt(3GM/R^3) + !brunt_frequency ! cycles per day + !brunt_N ! sqrt(abs(brunt_N2)) + !log_brunt_N ! log10(brunt_N) + !log_brunt_N2 ! log10(brunt_N2) + !log_brunt_N2_dimensionless ! log10(brunt_N2_dimensionless) + + brunt_B ! smoothed numerical difference + brunt_nonB ! = grada - gradT + !log_brunt_B ! smoothed numerical difference + !log_brunt_nonB ! = grada - gradT + + !sign_brunt_N2 ! sign of brunt_N2 (+1 for Ledoux stable; -1 for Ledoux unstable) + !brunt_nu ! brunt_frequency in microHz + !log_brunt_nu ! brunt_frequency in microHz + + !lamb_S ! lamb frequency for l=1: S = sqrt(2)*csound/r (rad/s) + !lamb_S2 ! squared lamb frequency for l=1: S2 = 2*(csound/r)^2 (rad^2/s^2) + + !lamb_Sl1 ! lamb frequency for l=1; = sqrt(2)*csound/r (microHz) + !lamb_Sl2 ! lamb frequency for l=2; = sqrt(6)*csound/r (microHz) + !lamb_Sl3 ! lamb frequency for l=3; = sqrt(12)*csound/r (microHz) + !lamb_Sl10 ! lamb frequency for l=10; = sqrt(110)*csound/r (microHz) + + !log_lamb_Sl1 ! log10(lamb_Sl1) + !log_lamb_Sl2 ! log10(lamb_Sl2) + !log_lamb_Sl3 ! log10(lamb_Sl3) + !log_lamb_Sl10 ! log10(lamb_Sl10) + + !brunt_N_div_r_integral ! integral from center of N*dr/r + !k_r_integral ! integral from center of k_r*dr + !brunt_N2_sub_omega2 + !sl2_sub_omega2 + + +!# RSP + + !rsp_Chi ! dlnP_dlnRho + !rsp_Et ! Specific turbulent energy + !rsp_logEt ! Log specific turbulent energy + !rsp_erad ! Specific internal (radiative) energy + !rsp_log_erad ! Log specific internal (radiative) energy + !rsp_Hp_face ! Pressure scale height at cell face + !rsp_Lc ! Convective luminosity + !rsp_Lc_div_L ! Convective luminosity div total luminosity + !rsp_Lr ! Radiative luminosity + !rsp_Lr_div_L ! Radiative luminosity div total luminosity + !rsp_Lt ! Turbulent luminosity + !rsp_Lt_div_L ! Turbulent luminosity div total luminosity + !rsp_Pt ! Turbulent pressure, p_t, see Table 1 in MESA5 + !rsp_Uq ! Viscous momentum transfer rate, U_q, see Table 1 in MESA5 + !rsp_Eq ! Viscous energy transfer rate, epsilon_q, see Table 1 in MESA5 + !rsp_Pvsc ! Artificial viscosity, p_av, see Table 1 in MESA5 + !rsp_gradT ! Temperature gradient + !rsp_Y_face ! Superadiabatic gradient at cell face, Y_sag, see Table 1 in MESA5 + !rsp_damp ! Turbulent dissipation, D, see Table 1 in MESA5 + !rsp_dampR ! Radiative cooling, D_r, see Table 1 in MESA5 + !rsp_sink ! Sum of turbulent dissipation and radiative cooling terms + !rsp_src ! Source function, S, see Table 1 in MESA5 + !rsp_src_snk ! Convective coupling, C, see Table 1 in MESA5 + !rsp_heat_exchange_timescale ! 1d0/(clight * opacity * density) + !rsp_log_heat_exchange_timescale + !rsp_log_dt_div_heat_exchange_timescale ! Ratio of time step to heat exchange timescale + !w + !log_w + + !COUPL + !DAMP + !DAMPR + !SOURCE + !Chi + !Eq + !Hp_face + !PII_face + !Ptrb + !Pvsc + !Uq + !Y_face + +!# RTI + + !RTI_du_diffusion_kick + !alpha_RTI + !boost_for_eta_RTI + !dedt_RTI + !dudt_RTI + !eta_RTI + !log_alpha_RTI + !log_boost_for_eta_RTI + !log_eta_RTI + !log_etamid_RTI + !log_lambda_RTI_div_Hrho + !log_sig_RTI + !log_sigmid_RTI + !log_source_RTI + !log_source_minus_alpha_RTI + !log_source_plus_alpha_RTI + !source_minus_alpha_RTI + !source_plus_alpha_RTI + !lambda_RTI + +!# Hydrodynamics + + + !v + !v_div_v_escape + !v_div_vesc + !v_kms + !log_v_escape + + !u + !u_face + + !P_face + + +!# Extras + !extra_heat + !extra_L ! extra_heat integrated from center (Lsun) + !log_extra_L ! log10 integrated from center (Lsun) + !log_irradiation_heat + + !extra_jdot ! set in other_torque routine + !extra_omegadot ! set in other_torque routine + + !extra_opacity_factor ! set in other_opacity_factor routine + + ! diffusion factor profile for species, set in other_diffusion_factor routine + !extra_diffusion_factor h1 + !extra_diffusion_factor he4 + !extra_diffusion_factor o16 + + + +!# Miscellaneous + + dlog_h1_dlogP ! (log(h1(k)) - log(h1(k-1)))/(log(P(k)) - log(P(k-1))) + !dlog_he3_dlogP + dlog_he4_dlogP + dlog_c12_dlogP + !dlog_c13_dlogP + !dlog_n14_dlogP + dlog_o16_dlogP + !dlog_ne20_dlogP + !dlog_mg24_dlogP + !dlog_si28_dlogP + + !dlog_pp_dlogP + !dlog_cno_dlogP + !dlog_3alf_dlogP + + !dlog_burn_c_dlogP + !dlog_burn_n_dlogP + !dlog_burn_o_dlogP + + !dlog_burn_ne_dlogP + !dlog_burn_na_dlogP + !dlog_burn_mg_dlogP + + !dlog_cc_dlogP + !dlog_co_dlogP + !dlog_oo_dlogP + + !dlog_burn_si_dlogP + !dlog_burn_s_dlogP + !dlog_burn_ar_dlogP + !dlog_burn_ca_dlogP + !dlog_burn_ti_dlogP + !dlog_burn_cr_dlogP + !dlog_burn_fe_dlogP + + !dlog_pnhe4_dlogP + !dlog_photo_dlogP + !dlog_other_dlogP + + !logR_kap ! logR = logRho - 3*logT + 18 ; used in kap tables + !logW ! logW = logPgas - 4*logT + !logQ ! logQ = logRho - 2*logT + 12 + !logV ! logV = logRho - 0.7*logE + 20 + + !log_CpT_absMdot_div_L ! log10(s% Cp(k)*s% T(k)*abs(s% mstar_dot)/s% L(k)) + + !delta_r ! r - r_start, change during step + !delta_L ! L - L_start, change during step + !delta_cell_vol ! cell_vol - cell_vol_start, change during step + !delta_entropy ! entropy - entropy_start, change during step (does not include effects of diffusion) + !delta_T ! T - T_start, change during step + !delta_rho ! rho - rho_start, change during step + !delta_eps_nuc ! eps_nuc - eps_nuc_start, change during step + !delta_mu ! mu - mu_start, change during step + + !zFe ! mass fraction of "Fe" = Fe+Co+Ni + !log_zFe + !dPdr_dRhodr_info + !log_sig_raw_mix + + !d_u_div_rmid + !d_u_div_rmid_start + !d_v_div_r_dm + !d_v_div_r_dr + + !dlnP_dlnR + !dlnRho_dlnR + !dlnRho_dr + !dlnX_dr + !dlnY_dr + !dlogR + !dPdr_div_grav + !dPdr_info + !dRhodr_info + !dRstar_div_dr + !dr_ratio + !dm_eps_grav + !dr_ratio + !dt_cs_div_dr + !dt_div_tau_conv + !dt_times_conv_vel_div_mixing_length + !log_dt_cs_div_dr + !log_dt_div_tau_conv + !log_dt_times_conv_vel_div_mixing_length + !log_du_kick_div_du + !du + !dvdt_dPdm + !dvdt_grav + + !tau_conv + !tau_cool + !tau_epsnuc + !tau_qhse + + !max_abs_xa_corr + + !tdc_num_iters + + !k + + +! the first few lines of the profile contain general info about the model. +! for completeness, those items are described here. + + ! initial mass and Z + ! initial_mass + ! initial_z + ! general properties of the current state + ! model_number + ! num_zones + ! star_age + ! time_step + ! properties at the photosphere + ! Teff + ! photosphere_L + ! photosphere_r + ! properties at the outermost zone of the model + ! log_surface_L + ! log_surface_radius + ! log_surface_temp + ! properties near the center of the model + ! log_center_temp + ! log_center_density + ! log_center_P + ! center_eta + ! abundances near the center + ! center_h1 + ! center_he3 + ! center_he4 + ! center_c12 + ! center_n14 + ! center_o16 + ! center_ne20 + ! information about total mass + ! star_mass + ! star_mdot + ! star_mass_h1 + ! star_mass_he3 + ! star_mass_he4 + ! star_mass_c12 + ! star_mass_n14 + ! star_mass_o16 + ! star_mass_ne20 + ! locations of abundance transitions + ! he_core_mass + ! c_core_mass + ! o_core_mass + ! si_core_mass + ! fe_core_mass + ! location of optical depths 10 and 100 + ! tau10_mass + ! tau10_radius + ! tau100_mass + ! tau100_radius + ! time scales + ! dynamic_time + ! kh_timescale + ! nuc_timescale + ! various kinds of total power + ! power_nuc_burn + ! power_h_burn + ! power_he_burn + ! power_neu + ! a few control parameter values + ! h1_boundary_limit + ! he4_boundary_limit + ! c12_boundary_limit + ! burn_min1 + ! burn_min2 diff --git a/star/test_suite/wd_o_ne_3_phase/re b/star/test_suite/wd_o_ne_3_phase/re new file mode 100755 index 000000000..c9ef26f96 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/re @@ -0,0 +1,33 @@ +#!/bin/bash + +shopt -u expand_aliases + +photo_directory=photos + +function most_recent_photo { + ls -tp "$photo_directory" | grep -v / | head -1 +} + +if [ $# -eq 0 ] +then + photo=$(most_recent_photo) +else + photo=$1 +fi + +if [ -z "$photo" ] || ! [ -f "$photo_directory/$photo" ] +then + echo "specified photo ($photo) does not exist" + exit 1 +fi + +echo "restart from $photo" +if ! cp "$photo_directory/$photo" restart_photo +then + echo "failed to copy photo ($photo)" + exit 1 +fi + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" +./star +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" diff --git a/star/test_suite/wd_o_ne_3_phase/rn b/star/test_suite/wd_o_ne_3_phase/rn new file mode 100755 index 000000000..25590040a --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/rn @@ -0,0 +1,7 @@ +#!/bin/bash + +rm -f restart_photo + +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" +./star +date "+DATE: %Y-%m-%d%nTIME: %H:%M:%S" diff --git a/star/test_suite/wd_o_ne_3_phase/src/run.f90 b/star/test_suite/wd_o_ne_3_phase/src/run.f90 new file mode 100644 index 000000000..2016c4921 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/src/run.f90 @@ -0,0 +1,13 @@ + program run + use run_star_support, only: do_read_star_job + use run_star, only: do_run_star + implicit none + integer :: ierr + character (len=32) :: inlist_fname + + ierr = 0 + inlist_fname = 'inlist' + call do_read_star_job(inlist_fname, ierr) + if (ierr /= 0) stop 1 + call do_run_star(inlist_fname) + end program run diff --git a/star/test_suite/wd_o_ne_3_phase/src/run_star_extras.f90 b/star/test_suite/wd_o_ne_3_phase/src/run_star_extras.f90 new file mode 100644 index 000000000..735b536d5 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/src/run_star_extras.f90 @@ -0,0 +1,447 @@ +! *********************************************************************** +! +! Copyright (C) 2010-2019 Bill Paxton & The MESA Team +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU Lesser General Public License +! as published by the Free Software Foundation, +! either version 3 of the License, or (at your option) any later version. +! +! This program 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 Lesser General Public License for more details. +! +! You should have received a copy of the GNU Lesser General Public License +! along with this program. If not, see . +! +! *********************************************************************** + + module run_star_extras + use star_lib + use star_def + use const_def + use math_lib + implicit none + ! these routines are called by the standard run_star check_model + contains + + subroutine extras_controls(id, ierr) + integer, intent(in) :: id + integer, intent(out) :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + ! this is the place to set any procedure pointers you want to change + ! e.g., other_wind, other_mixing, other_energy (see star_data.inc) + ! the extras functions in this file will not be called + ! unless you set their function pointers as done below. + ! otherwise we use a null_ version which does nothing (except warn). + s% extras_startup => extras_startup + s% extras_start_step => extras_start_step + s% extras_check_model => extras_check_model + s% extras_finish_step => extras_finish_step + s% extras_after_evolve => extras_after_evolve + s% how_many_extra_history_columns => how_many_extra_history_columns + s% data_for_extra_history_columns => data_for_extra_history_columns + s% how_many_extra_profile_columns => how_many_extra_profile_columns + s% data_for_extra_profile_columns => data_for_extra_profile_columns + s% how_many_extra_history_header_items => how_many_extra_history_header_items + s% data_for_extra_history_header_items => data_for_extra_history_header_items + s% how_many_extra_profile_header_items => how_many_extra_profile_header_items + s% data_for_extra_profile_header_items => data_for_extra_profile_header_items + end subroutine extras_controls + + subroutine extras_startup(id, restart, ierr) + integer, intent(in) :: id + logical, intent(in) :: restart + integer, intent(out) :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + end subroutine extras_startup + + integer function extras_start_step(id) + integer, intent(in) :: id + integer :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + extras_start_step = 0 + end function extras_start_step + + ! returns either keep_going, retry, or terminate. + integer function extras_check_model(id) + use chem_def + use eos_def + integer, intent(in) :: id + integer :: ierr, k, i_accr, iXC, iXO, iXne20, iXNe22, iXNa, iXMg, iXHe + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + extras_check_model = keep_going + if (.false. .and. s% star_mass_h1 < 0.35d0) then + ! stop when star hydrogen mass drops to specified level + extras_check_model = terminate + write(*, *) 'have reached desired hydrogen mass' + return + end if + ! if you want to check multiple conditions, it can be useful + ! to set a different termination code depending on which + ! condition was triggered. MESA provides 9 customizeable + ! termination codes, named t_xtra1 .. t_xtra9. You can + ! customize the messages that will be printed upon exit by + ! setting the corresponding termination_code_str value. + ! termination_code_str(t_xtra1) = 'my termination condition' + ! by default, indicate where (in the code) MESA terminated + if (extras_check_model == terminate) s% termination_code = t_extras_check_model + ! Set iX to the index of o16 and iY to ne20 + iXC = -1 + iXO = -1 + iXNe20 = -1 + iXNe22 = -1 + iXNa = -1 + iXMg = -1 + iXHe = -1 + do k = 1,s% species + if (chem_isos% name(s% chem_id(k)) == 'c12') then + iXC = k + end if + if (chem_isos% name(s% chem_id(k)) == 'o16') then + iXO = k + end if + if (chem_isos% name(s% chem_id(k)) == 'ne20') then + iXNe20 = k + end if + if (chem_isos% name(s% chem_id(k)) == 'ne22') then + iXNe22 = k + end if + if (chem_isos% name(s% chem_id(k)) == 'na23') then + iXNa = k + end if + if (chem_isos% name(s% chem_id(k)) == 'mg24') then + iXMg = k + end if + if (chem_isos% name(s% chem_id(k)) == 'he4') then + iXHe = k + end if + end do + if (iXC == -1 .or. iXO == -1 .or. iXNe20 == -1 .or. iXNe22 == -1 .or. iXNa == -1 .or. iXMg == -1 .or. iXHe == -1) then + write (*,*) 'Could not find elements specified!!' + end if + ! Find the base of the oxygen layer + !do k = 1, s% nz + ! if (s% xa(iX, k) .le. 0.1) then + !write(*,*) log10(s% rho(k)), s% xa(iX, k), s% gam(k), s% chiRho(k), s% chiT(k), s% d_eos_dxa(i_lnPgas,iX,k), s% d_eos_dxa(i_lnPgas,iY,k) + ! exit + !end if + !end do + ! terminate the model when the base density reaches a particular value + if (.false.) then + write(*,*) 'base density=', log10(s% rho(s% nz)) + if (s% rho(s% nz) > 1d10) then + extras_check_model = terminate + termination_code_str(t_xtra1) = 'base density' + s% termination_code = t_xtra1 + end if + end if + end function extras_check_model + + integer function how_many_extra_history_columns(id) + integer, intent(in) :: id + integer :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + how_many_extra_history_columns = 3 + end function how_many_extra_history_columns + + subroutine data_for_extra_history_columns(id, n, names, vals, ierr) + use chem_def, only: i3alf + integer, intent(in) :: id, n + character (len=maxlen_history_column_name) :: names(n) + real(dp) :: vals(n), rcore + integer, intent(out) :: ierr + type (star_info), pointer :: s + integer :: i_max,kc,k + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + ! note: do NOT add the extras names to history_columns.list + ! the history_columns.list is only for the built-in history column options. + ! it must not include the new column names you are adding here. + ! find the maximum in the helium burning rate + i_max = maxloc(s% eps_nuc_categories(i3alf,:), dim=1) + ! to find the maximum in the total epsilon: + !i_max = maxloc(s% eps_nuc, dim=1) + names(1) = "max_eps_he_lgT" + vals(1) = log10(s% T(i_max)) + names(2) = "mass_core_cryst" + vals(2) = s% crystal_core_boundary_mass + do k = s%nz,1,-1 + if(s% m(k) >= s% crystal_core_boundary_mass) then + kc = k + exit + end if + end do + names(3) = "r_core_cryst" + vals(3) = exp(s% lnR(kc)) + end subroutine data_for_extra_history_columns + + integer function how_many_extra_profile_columns(id) + integer, intent(in) :: id + integer :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + how_many_extra_profile_columns = 14 + end function how_many_extra_profile_columns + + subroutine data_for_extra_profile_columns(id, n, nz, names, vals, ierr) + use chem_def + use eos_def + use eos_lib + integer :: iXC, iXO, iXne20, iXNe22, iXNa, iXMg, iXHe + integer, intent(in) :: id, n, nz + character (len=maxlen_profile_column_name) :: names(n) + real(dp) :: vals(nz,n) + real(dp), allocatable :: xa1_c12(:), xa1_o16(:), xa1_ne20(:), xa1_ne22(:), xa1_na23(:), xa1_mg24(:), xa1_he4(:) + real(dp), dimension(num_eos_basic_results) :: res, d_dlnd, d_dlnT + integer, intent(out) :: ierr + real(dp), allocatable, dimension(:,:) :: d_dxa + integer :: k + real(dp) :: P1, S1, my_chiT, eps, chiX_C12, chiX_O16, chiX_Ne20, chiX_Ne22, chiX_Na23, chiX_Mg24,chiX_He4, bs_C12, bs_O16, bs_Ne20, bs_Ne22, bs_Na23, bs_Mg24, bs_He4, mu1 !!! + !real(dp) :: chimu_c12,chimu_o16,chimu_ne20,chimu_mg24 + real(dp) :: plnxc_plnye, plnxo_plnye, plnxne20_plnye, plnxne22_plnye, plnxmg_plnye, ln_ye + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + allocate(d_dxa(num_eos_d_dxa_results,s% species)) + allocate(xa1_c12(s% species)) + allocate(xa1_o16(s% species)) + allocate(xa1_ne20(s% species)) + allocate(xa1_ne22(s% species)) + allocate(xa1_na23(s% species)) + allocate(xa1_mg24(s% species)) + allocate(xa1_he4(s% species)) + names(1) = 'chiX_C12' + names(2) = 'chiX_O16' + names(3) = 'chiX_Ne20' + names(4) = 'chiX_Ne22' + names(5) = 'chiX_Na23' + names(6) = 'chiX_Mg24' + names(7) = 'chiX_He4' + names(8) = 'bs_C12' + names(9) = 'bs_O16' + names(10) = 'bs_Ne20' + names(11) = 'bs_Ne22' + names(12) = 'bs_N23' + names(13) = 'bs_Mg24' + names(14) = 'bs_He4' + iXC = -1 + iXO = -1 + iXNe20 = -1 + iXNe22 = -1 + iXNa = -1 + iXMg = -1 + iXHe = -1 + do k = 1,s% species + if (chem_isos% name(s% chem_id(k)) == 'c12') then + iXC = k + end if + if (chem_isos% name(s% chem_id(k)) == 'o16') then + iXO = k + end if + if (chem_isos% name(s% chem_id(k)) == 'ne20') then + iXNe20 = k + end if + if (chem_isos% name(s% chem_id(k)) == 'ne22') then + iXNe22 = k + end if + if (chem_isos% name(s% chem_id(k)) == 'na23') then + iXNa = k + end if + if (chem_isos% name(s% chem_id(k)) == 'mg24') then + iXMg = k + end if + if (chem_isos% name(s% chem_id(k)) == 'he4') then + iXHe = k + end if + end do + if (iXC == -1 .or. iXO == -1 .or. iXNe20 == -1 .or. iXNe22 == -1 .or. iXMg == -1 .or. iXNa == -1 .or. iXHe == -1) then + write (*,*) 'Could not find elements specified!!' + end if + do k = 1, nz + eps = 1d-4 + call eosDT_get( & + s% eos_handle, s% species, s% chem_id, s% net_iso, s% xa(:, k), & + s% rho(k), log10(s% rho(k)), s% T(k), log10(s% T(k)), & + res, d_dlnd, d_dlnT, d_dxa, ierr) + P1 = res(i_lnPgas) + S1 = exp(res(i_lnS)) + mu1= res(i_mu)!s% mu(k) + ln_ye=res(i_lnfree_e) + xa1_c12 = s% xa(:, k) + xa1_c12(iXC) = xa1_c12(iXC)*(1d0+eps) + xa1_o16 = s% xa(:, k) + xa1_o16(iXO) = xa1_o16(iXO)*(1d0+eps) + xa1_ne20 = s% xa(:, k) + xa1_ne20(iXNe20) = xa1_ne20(iXNe20)*(1d0+eps) + xa1_ne22 = s% xa(:, k) + xa1_ne22(iXNe22) = xa1_ne22(iXNe22)*(1d0+eps) + xa1_na23 = s% xa(:, k) + xa1_na23(iXNa) = xa1_na23(iXNa)*(1d0+eps) + xa1_mg24 = s% xa(:, k) + xa1_mg24(iXMg) = xa1_mg24(iXMg)*(1d0+eps) + xa1_he4 = s% xa(:, k) + xa1_he4(iXHe) = xa1_he4(iXHe)*(1d0+eps) + call eosDT_get( & + s% eos_handle, s% species, s% chem_id, s% net_iso, xa1_c12, & + s% rho(k), log10(s% rho(k)), s% T(k), log10(s% T(k)), & + res, d_dlnd, d_dlnT, d_dxa, ierr) + chiX_C12 = (res(i_lnPgas)-P1)/eps + bs_C12 = - s% xa(iXC,k)*(exp(res(i_lnS))-S1)/(eps* s% xa(iXC,k)) !!-X ds/dX as in Medin & Cumming (2015) + plnxc_plnye= (log(s% xa(iXC, k))-log(xa1_c12(iXC)))/(ln_ye-res(i_lnfree_e)) + call eosDT_get( & + s% eos_handle, s% species, s% chem_id, s% net_iso, xa1_o16, & + s% rho(k), log10(s% rho(k)), s% T(k), log10(s% T(k)), & + res, d_dlnd, d_dlnT, d_dxa, ierr) + chiX_O16 = (res(i_lnPgas)-P1)/eps + bs_O16 = - s% xa(iXO,k)*(exp(res(i_lnS))-S1)/(eps* s% xa(iXO,k)) + plnxo_plnye= (log(s% xa(iXO, k))-log(xa1_o16(iXO)))/(ln_ye-res(i_lnfree_e)) + call eosDT_get( & + s% eos_handle, s% species, s% chem_id, s% net_iso, xa1_ne20, & + s% rho(k), log10(s% rho(k)), s% T(k), log10(s% T(k)), & + res, d_dlnd, d_dlnT, d_dxa, ierr) + chiX_Ne20 = (res(i_lnPgas)-P1)/eps + bs_Ne20 = - s% xa(iXNe20,k)*(exp(res(i_lnS))-S1)/(eps* s% xa(iXNe20,k)) + plnxne20_plnye= (log(s% xa(iXNe20, k))-log(xa1_ne20(iXNe20)))/(ln_ye-res(i_lnfree_e)) + call eosDT_get( & + s% eos_handle, s% species, s% chem_id, s% net_iso, xa1_ne22, & + s% rho(k), log10(s% rho(k)), s% T(k), log10(s% T(k)), & + res, d_dlnd, d_dlnT, d_dxa, ierr) + chiX_Ne22 = (res(i_lnPgas)-P1)/eps + bs_Ne22 = - s% xa(iXNe22,k)*(exp(res(i_lnS))-S1)/(eps* s% xa(iXNe22,k)) + plnxne22_plnye= (log(s% xa(iXNe22, k))-log(xa1_ne20(iXNe22)))/(ln_ye-res(i_lnfree_e)) + call eosDT_get( & + s% eos_handle, s% species, s% chem_id, s% net_iso, xa1_na23, & + s% rho(k), log10(s% rho(k)), s% T(k), log10(s% T(k)), & + res, d_dlnd, d_dlnT, d_dxa, ierr) + chiX_Na23 = (res(i_lnPgas)-P1)/eps + bs_Na23 = - s% xa(iXNa,k)*(exp(res(i_lnS))-S1)/(eps* s% xa(iXNa,k)) + call eosDT_get( & + s% eos_handle, s% species, s% chem_id, s% net_iso, xa1_mg24, & + s% rho(k), log10(s% rho(k)), s% T(k), log10(s% T(k)), & + res, d_dlnd, d_dlnT, d_dxa, ierr) + chiX_Mg24 = (res(i_lnPgas)-P1)/eps + bs_Mg24 = - s% xa(iXMg,k)*(exp(res(i_lnS))-S1)/(eps* s% xa(iXMg,k)) + plnxmg_plnye= (log(s% xa(iXMg, k))-log(xa1_mg24(iXMg)))/(ln_ye-res(i_lnfree_e)) + call eosDT_get( & + s% eos_handle, s% species, s% chem_id, s% net_iso, xa1_he4, & + s% rho(k), log10(s% rho(k)), s% T(k), log10(s% T(k)), & + res, d_dlnd, d_dlnT, d_dxa, ierr) + chiX_He4 = (res(i_lnPgas)-P1)/eps + bs_He4 = - s% xa(iXHe,k)*(exp(res(i_lnS))-S1)/(eps* s% xa(iXHe,k)) + vals(k,1) = chiX_C12 + vals(k,2) = chiX_O16 + vals(k,3) = chiX_Ne20 + vals(k,4) = chiX_Ne22 + vals(k,5) = chiX_Na23 + vals(k,6) = chiX_Mg24 + vals(k,7) = chiX_He4 + vals(k,8) = bs_C12 + vals(k,9) = bs_O16 + vals(k,10) = bs_Ne20 + vals(k,11) = bs_Ne22 + vals(k,12) = bs_Na23 + vals(k,13) = bs_Mg24 + vals(k,14) = bs_He4 + end do + end subroutine data_for_extra_profile_columns + + integer function how_many_extra_history_header_items(id) + integer, intent(in) :: id + integer :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + how_many_extra_history_header_items = 0 + end function how_many_extra_history_header_items + + subroutine data_for_extra_history_header_items(id, n, names, vals, ierr) + integer, intent(in) :: id, n + character (len=maxlen_history_column_name) :: names(n) + real(dp) :: vals(n) + type(star_info), pointer :: s + integer, intent(out) :: ierr + ierr = 0 + call star_ptr(id,s,ierr) + if(ierr/=0) return + ! here is an example for adding an extra history header item + ! also set how_many_extra_history_header_items + ! names(1) = 'mixing_length_alpha' + ! vals(1) = s% mixing_length_alpha + end subroutine data_for_extra_history_header_items + + integer function how_many_extra_profile_header_items(id) + integer, intent(in) :: id + integer :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + how_many_extra_profile_header_items = 0 + end function how_many_extra_profile_header_items + + subroutine data_for_extra_profile_header_items(id, n, names, vals, ierr) + integer, intent(in) :: id, n + character (len=maxlen_profile_column_name) :: names(n) + real(dp) :: vals(n) + type(star_info), pointer :: s + integer, intent(out) :: ierr + ierr = 0 + call star_ptr(id,s,ierr) + if(ierr/=0) return + ! here is an example for adding an extra profile header item + ! also set how_many_extra_profile_header_items + ! names(1) = 'mixing_length_alpha' + ! vals(1) = s% mixing_length_alpha + end subroutine data_for_extra_profile_header_items + + ! returns either keep_going or terminate. + ! note: cannot request retry; extras_check_model can do that. + integer function extras_finish_step(id) + integer, intent(in) :: id + integer :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + extras_finish_step = keep_going + ! to save a profile, + ! s% need_to_save_profiles_now = .true. + ! to update the star log, + ! s% need_to_update_history_now = .true. + ! see extras_check_model for information about custom termination codes + ! by default, indicate where (in the code) MESA terminated + if (extras_finish_step == terminate) s% termination_code = t_extras_finish_step + end function extras_finish_step + + subroutine extras_after_evolve(id, ierr) + integer, intent(in) :: id + integer, intent(out) :: ierr + type (star_info), pointer :: s + ierr = 0 + call star_ptr(id, s, ierr) + if (ierr /= 0) return + end subroutine extras_after_evolve + end module run_star_extras + diff --git a/star/test_suite/wd_o_ne_3_phase/wd1_10_mi8_3_NeOMg.mod b/star/test_suite/wd_o_ne_3_phase/wd1_10_mi8_3_NeOMg.mod new file mode 100644 index 000000000..a25c75bbb --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/wd1_10_mi8_3_NeOMg.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:123e893958689c6803f62b68fc657fa247a53e1b4430747776ab12bd64797cdc +size 565419 diff --git a/star/test_suite/wd_o_ne_3_phase/wd1_10_mi8_3_ONeNa.mod b/star/test_suite/wd_o_ne_3_phase/wd1_10_mi8_3_ONeNa.mod new file mode 100644 index 000000000..acb31441d --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/wd1_10_mi8_3_ONeNa.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b3d4ae9a436c801bc38b733fe1d978caa677a4bc3099e3603426ac144718f48c +size 509615 diff --git a/star/test_suite/wd_o_ne_3_phase/wd1_30_mi9_4_NeOMg.mod b/star/test_suite/wd_o_ne_3_phase/wd1_30_mi9_4_NeOMg.mod new file mode 100644 index 000000000..5fd907111 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/wd1_30_mi9_4_NeOMg.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:207b4cf56aec6edb46af9ecfc81bb0ed416b1ab695f12fd5e3e0e5e91534748c +size 794400 diff --git a/star/test_suite/wd_o_ne_3_phase/wd1_30_mi9_4_ONeNa.mod b/star/test_suite/wd_o_ne_3_phase/wd1_30_mi9_4_ONeNa.mod new file mode 100644 index 000000000..c8d22fe41 --- /dev/null +++ b/star/test_suite/wd_o_ne_3_phase/wd1_30_mi9_4_ONeNa.mod @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7039a22c97e99fa3cdbcb3a6c99322c7db96416ce85b772c9b71d889ca33d2ec +size 508217