diff --git a/star/private/phase_separation.f90 b/star/private/phase_separation.f90
index 09b25374a..89a84b169 100644
--- a/star/private/phase_separation.f90
+++ b/star/private/phase_separation.f90
@@ -2,36 +2,40 @@
!
! Copyright (C) 2021 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.
+! MESA is free software; you can use it and/or modify
+! it under the combined terms and restrictions of the MESA MANIFESTO
+! and the GNU General Library Public License as published
+! by the Free Software Foundation; either version 2 of the License,
+! or (at your option) any later version.
!
-! This program is distributed in the hope that it will be useful,
+! You should have received a copy of the MESA MANIFESTO along with
+! this software; if not, it is available at the mesa website:
+! http://mesa.sourceforge.net/
+!
+! MESA is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-! See the GNU Lesser General Public License for more details.
+! See the GNU Library 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 .
+! You should have received a copy of the GNU Library General Public License
+! along with this software; if not, write to the Free Software
+! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!
! ***********************************************************************
-
module 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 +44,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 +58,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
+ integer, intent(out) :: ierr
+ 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
- end if
-
+ 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
@@ -82,19 +93,15 @@ subroutine do_2component_phase_separation(s, dt, components, ierr)
k_bound = k
exit
end if
- end do
-
+ 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 +111,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
@@ -112,138 +118,217 @@ subroutine do_2component_phase_separation(s, dt, components, ierr)
call update_model_(s,1,s%nz,.false.)
do k=1,s% nz
s% eps_phase_separation(k) = s% energy(k)
- end do
-
+ 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)
-
+ 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 .ne. '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
-
+ 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
-
- call update_model_(s,k-1,s%nz,.true.)
-
+ 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
-
+ 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 +337,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.)
-
+ 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)
-
+ 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 +369,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))
-
+ 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 +400,310 @@ 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 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
-
+ 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).
@@ -358,43 +713,36 @@ subroutine update_model_ (s, kc_t, kc_b, do_brunt)
write(*,*) 'phase_separation: error from call to set_eos_with_mask'
stop
end if
- s%fix_Pgas = .false.
-
+ 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)
+ ! 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
-
+ 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
-
- return
-
+ 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 100644
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 100644
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 100644
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 100644
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..796663911
--- /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
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..7e4dc99cb
--- /dev/null
+++ b/star/test_suite/wd_o_ne_3_phase/src/run_star_extras.f90
@@ -0,0 +1,450 @@
+! ***********************************************************************
+!
+! Copyright (C) 2010-2019 Bill Paxton & The MESA Team
+!
+! this file is part of mesa.
+!
+! mesa is free software; you can redistribute it and/or modify
+! it under the terms of the gnu general library public license as published
+! by the free software foundation; either version 2 of the license, or
+! (at your option) any later version.
+!
+! mesa is distributed in the hope that it will be useful,
+! but without any warranty; without even the implied warranty of
+! merchantability or fitness for a particular purpose. see the
+! gnu library general public license for more details.
+!
+! you should have received a copy of the gnu library general public license
+! along with this software; if not, write to the free software
+! foundation, inc., 59 temple place, suite 330, boston, ma 02111-1307 usa
+!
+! ***********************************************************************
+
+ module run_star_extras
+ use star_lib
+ use star_def
+ use const_def
+ use 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)) .eq. 'c12') then
+ iXC = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'o16') then
+ iXO = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'ne20') then
+ iXNe20 = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'ne22') then
+ iXNe22 = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'na23') then
+ iXNa = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'mg24') then
+ iXMg = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'he4') then
+ iXHe = k
+ end if
+ end do
+ if (iXC .eq. -1 .or. iXO .eq. -1 .or. iXNe20 .eq. -1 .or. iXNe22 .eq. -1 .or. iXNa .eq. -1 .or. iXMg .eq. -1 .or. iXHe .eq. -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) .ge. 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)) .eq. 'c12') then
+ iXC = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'o16') then
+ iXO = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'ne20') then
+ iXNe20 = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'ne22') then
+ iXNe22 = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'na23') then
+ iXNa = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'mg24') then
+ iXMg = k
+ end if
+ if (chem_isos% name(s% chem_id(k)) .eq. 'he4') then
+ iXHe = k
+ end if
+ end do
+ if (iXC .eq. -1 .or. iXO .eq. -1 .or. iXNe20 .eq. -1 .or. iXNe22 .eq. -1 .or. iXMg .eq. -1 .or. iXNa .eq. -1 .or. iXHe .eq. -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