Merge pull request #591 from themikelau/CE-analysis

Update to analysis_common_envelope

scripts/test_analysis_ce.sh

@@ -29,6 +29,7 @@ SEP
 no
 2
 1.667
+0.6182
 0.6984
 0.0142
 BOUND
@@ -38,4 +39,5 @@ BOUND
 no
 2
 1.667
+0.6182
 ENERGIES

src/main/eos.f90

@@ -840,7 +840,7 @@ end subroutine calc_rec_ene
 !+
 !-----------------------------------------------------------------------
 subroutine calc_temp_and_ene(eos_type,rho,pres,ene,temp,ierr,guesseint,mu_local,X_local,Z_local)
- use physcon,          only:kb_on_mh
+ use physcon,          only:Rg
  use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_enfromtemp
  use eos_mesa,         only:get_eos_eT_from_rhop_mesa
  use eos_gasradrec,    only:calc_uT_from_rhoP_gasradrec
@@ -861,7 +861,7 @@ subroutine calc_temp_and_ene(eos_type,rho,pres,ene,temp,ierr,guesseint,mu_local,
  if (present(Z_local))  Z  = Z_local
  select case(eos_type)
  case(2,5,17) ! Ideal gas
-    temp = pres / (rho * kb_on_mh) * mu
+    temp = pres / (rho * Rg) * mu
     ene = pres / ( (gamma-1.) * rho)
  case(12) ! Ideal gas + radiation
     call get_idealgasplusrad_tempfrompres(pres,rho,mu,temp)
@@ -889,7 +889,7 @@ end subroutine calc_temp_and_ene
 !+
 !-----------------------------------------------------------------------
 subroutine calc_rho_from_PT(eos_type,pres,temp,rho,ierr,mu_local,X_local,Z_local)
- use physcon,          only:kb_on_mh
+ use physcon,          only:Rg
  use eos_idealplusrad, only:get_idealplusrad_rhofrompresT
  use eos_mesa,         only:get_eos_eT_from_rhop_mesa
  use eos_gasradrec,    only:calc_uT_from_rhoP_gasradrec
@@ -910,7 +910,7 @@ subroutine calc_rho_from_PT(eos_type,pres,temp,rho,ierr,mu_local,X_local,Z_local
  if (present(Z_local))  Z  = Z_local
  select case(eos_type)
  case(2) ! Ideal gas
-    rho = pres / (temp * kb_on_mh) * mu
+    rho = pres / (temp * Rg) * mu
  case(12) ! Ideal gas + radiation
     call get_idealplusrad_rhofrompresT(pres,temp,mu,rho)
  case default
@@ -925,36 +925,48 @@ end subroutine calc_rho_from_PT
 !  up to an additive integration constant, from density and pressure.
 !+
 !-----------------------------------------------------------------------
-function entropy(rho,pres,mu_in,ientropy,eint_in,ierr)
+function entropy(rho,pres,mu_in,ientropy,eint_in,ierr,T_in,Trad_in)
  use io,                only:fatal,warning
- use physcon,           only:radconst,kb_on_mh
+ use physcon,           only:radconst,kb_on_mh,Rg
  use eos_idealplusrad,  only:get_idealgasplusrad_tempfrompres
  use eos_mesa,          only:get_eos_eT_from_rhop_mesa
  use mesa_microphysics, only:getvalue_mesa
  real,    intent(in)            :: rho,pres,mu_in
- real,    intent(in),  optional :: eint_in
+ real,    intent(in),  optional :: eint_in,T_in,Trad_in
  integer, intent(in)            :: ientropy
  integer, intent(out), optional :: ierr
- real                           :: mu,entropy,logentropy,temp,eint
+ real                           :: mu,entropy,logentropy,temp,Trad,eint
 
  if (present(ierr)) ierr=0
 
  mu = mu_in
- select case(ientropy)
- case(1) ! Include only gas entropy (up to additive constants)
-    temp = pres * mu / (rho * kb_on_mh)
-    entropy = kb_on_mh / mu * log(temp**1.5/rho)
+ if (present(T_in)) then  ! is gas temperature specified?
+    temp = T_in
+ else
+    temp = pres * mu / (rho * Rg)  ! used as initial guess for case 2
+ endif
 
+ select case(ientropy)
+ case(1) ! Only include gas contribution
     ! check temp
     if (temp < tiny(0.)) call warning('entropy','temperature = 0 will give minus infinity with s entropy')
+    entropy = kb_on_mh / mu * log(temp**1.5/rho)
 
- case(2) ! Include both gas and radiation entropy (up to additive constants)
-    temp = pres * mu / (rho * kb_on_mh) ! Guess for temp
-    call get_idealgasplusrad_tempfrompres(pres,rho,mu,temp) ! First solve for temp from rho and pres
-    entropy = kb_on_mh / mu * log(temp**1.5/rho) + 4.*radconst*temp**3 / (3.*rho)
-
+ case(2) ! Include both gas and radiation contributions (up to additive constants)
+    if (present(T_in)) then
+       temp = T_in
+       if (present(Trad_in)) then
+          Trad = Trad_in
+       else
+          Trad = temp
+       endif
+    else
+       call get_idealgasplusrad_tempfrompres(pres,rho,mu,temp) ! First solve for temp from rho and pres
+       Trad = temp
+    endif
     ! check temp
     if (temp < tiny(0.)) call warning('entropy','temperature = 0 will give minus infinity with s entropy')
+    entropy = kb_on_mh / mu * log(temp**1.5/rho) + 4.*radconst*Trad**3 / (3.*rho)
 
  case(3) ! Get entropy from MESA tables if using MESA EoS
     if (ieos /= 10 .and. ieos /= 20) call fatal('eos','Using MESA tables to calculate S from rho and pres, but not using MESA EoS')
@@ -1039,7 +1051,7 @@ end subroutine get_rho_from_p_s
 !+
 !-----------------------------------------------------------------------
 subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp)
- use physcon, only:kb_on_mh,radconst,rg,mass_proton_cgs,kboltz
+ use physcon, only:radconst,Rg,mass_proton_cgs,kboltz
  use io,      only:fatal
  use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_pres
  use units,   only:unit_density,unit_pressure,unit_ergg
@@ -1060,7 +1072,7 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp)
  select case (ieos)
  case (2,5,17)
     temp = (cgsrho * exp(mu*cgss*mass_proton_cgs))**(2./3.)
-    cgsP = cgsrho*kb_on_mh*temp / mu
+    cgsP = cgsrho*Rg*temp / mu
  case (12)
     corr = huge(corr)
     do while (abs(corr) > eoserr .and. niter < nitermax)

src/main/eos_idealplusrad.f90

@@ -24,7 +24,7 @@ module eos_idealplusrad
 
  public :: get_idealplusrad_temp,get_idealplusrad_pres,get_idealplusrad_spsoundi,&
            get_idealgasplusrad_tempfrompres,get_idealplusrad_enfromtemp,&
-           get_idealplusrad_rhofrompresT
+           get_idealplusrad_rhofrompresT,egas_from_rhoT,erad_from_rhoT
 
  private
 
@@ -125,11 +125,37 @@ subroutine get_idealplusrad_enfromtemp(densi,tempi,mu,eni)
  real, intent(in)  :: densi,tempi,mu
  real, intent(out) :: eni
 
- eni = 1.5*Rg*tempi/mu + radconst*tempi**4/densi
+ eni = egas_from_rhoT(tempi,mu) + erad_from_rhoT(densi,tempi,mu)
 
 end subroutine get_idealplusrad_enfromtemp
 
 
+!----------------------------------------------------------------
+!+
+!  Calculates specific gas energy from density and temperature
+!+
+!----------------------------------------------------------------
+real function egas_from_rhoT(tempi,mu) result(egasi)
+ real, intent(in) :: tempi,mu
+
+ egasi = 1.5*Rg*tempi/mu
+
+end function egas_from_rhoT
+
+
+!----------------------------------------------------------------
+!+
+!  Calculates specific radiation energy from density and temperature
+!+
+!----------------------------------------------------------------
+real function erad_from_rhoT(densi,tempi,mu) result(eradi)
+ real, intent(in) :: densi,tempi,mu
+
+ eradi = radconst*tempi**4/densi
+
+end function erad_from_rhoT
+
+
 !----------------------------------------------------------------
 !+
 !  Calculates density from pressure and temperature

src/main/inject_rochelobe.f90

@@ -165,7 +165,6 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
     h = hfact*sw_chi/udist
     !add the particle
     call add_or_update_particle(part_type, xyzi, vxyz, h, u, i_part, npart, npartoftype, xyzh, vxyzu)
-
  enddo
  !
  !-- no constraint on timestep

src/main/ionization.f90

@@ -339,28 +339,27 @@ end subroutine get_erec_components
 !  gas particle. Inputs and outputs in code units
 !+
 !----------------------------------------------------------------
-subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,ethi,radprop)
+subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,ethi,rad)
  use dim,              only:do_radiation
  use part,             only:rhoh,iradxi
  use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_enfromtemp
  use physcon,          only:radconst,Rg
  use units,            only:unit_density,unit_pressure,unit_ergg,unit_pressure
  integer, intent(in) :: ieos
  real, intent(in)    :: particlemass,presi,tempi,xyzh(4),vxyzu(4)
- real, intent(in), optional :: radprop(:)
+ real, intent(in), optional :: rad(:)
  real, intent(out)   :: ethi
- real                :: hi,densi_cgs,mui
+ real                :: densi_cgs,mui
 
  select case (ieos)
  case(10,20) ! calculate just gas + radiation thermal energy
-    hi = xyzh(4)
-    densi_cgs = rhoh(hi,particlemass)*unit_density
+    densi_cgs = rhoh(xyzh(4),particlemass)*unit_density
     mui = densi_cgs * Rg * tempi / (presi*unit_pressure - radconst * tempi**4 / 3.) ! Get mu from pres and temp
     call get_idealplusrad_enfromtemp(densi_cgs,tempi,mui,ethi)
     ethi = particlemass * ethi / unit_ergg
  case default ! assuming internal energy = thermal energy
     ethi = particlemass * vxyzu(4)
-    if (do_radiation) ethi  = ethi + particlemass*radprop(iradxi)
+    if (do_radiation) ethi  = ethi + particlemass*rad(iradxi)
  end select
 
 end subroutine calc_thermal_energy
@@ -420,9 +419,9 @@ subroutine ionisation_fraction(dens,temp,X,Y,xh0,xh1,xhe0,xhe1,xhe2)
     xhe2g = xhe2g + dx(3)
  enddo
 
- xh1 = xh1g * n / nh
- xhe1 = xhe1g * n / nhe
- xhe2 = xhe2g * n / nhe
+ xh1 = max(xh1g * n / nh,1.e-99)
+ xhe1 = max(xhe1g * n / nhe,1.e-99)
+ xhe2 = max(xhe2g * n / nhe,1.e-99)
  xh0 = ((nh/n) - xh1g) * n / nh
  xhe0 = ((nhe/n) - xhe1g - xhe2g) * n / nhe
 

src/main/radiation_utils.f90

@@ -23,8 +23,8 @@ module radiation_utils
  public :: get_rad_R
  public :: radiation_equation_of_state
  public :: T_from_Etot
- public :: radE_from_Trad
- public :: Trad_from_radE
+ public :: radxi_from_Trad
+ public :: Trad_from_radxi
  public :: ugas_from_Tgas
  public :: Tgas_from_ugas
  public :: get_opacity
@@ -136,29 +136,29 @@ end function T_from_Etot
 
 !---------------------------------------------------------
 !+
-!  get the radiation energy from the radiation temperature
+!  get specific radiation energy from radiation temperature
 !+
 !---------------------------------------------------------
-real function radE_from_Trad(Trad) result(radE)
+real function radxi_from_Trad(rho,Trad) result(radxi)
  use units, only:get_radconst_code
- real, intent(in) :: Trad
+ real, intent(in) :: rho,Trad
 
- radE = Trad**4*get_radconst_code()
+ radxi = Trad**4*get_radconst_code()/rho
 
-end function radE_from_Trad
+end function radxi_from_Trad
 
 !---------------------------------------------------------
 !+
-!  get the radiation temperature from the radiation energy per unit volume
+!  get radiation temperature from the specific radiation energy
 !+
 !---------------------------------------------------------
-real function Trad_from_radE(radE) result(Trad)
+real function Trad_from_radxi(rho,radxi) result(Trad)
  use units, only:get_radconst_code
- real, intent(in) :: radE
+ real, intent(in) :: rho,radxi
 
- Trad = (radE/get_radconst_code())**0.25
+ Trad = (rho*radxi/get_radconst_code())**0.25
 
-end function Trad_from_radE
+end function Trad_from_radxi
 
 !---------------------------------------------------------
 !+

src/setup/set_star_utils.f90

@@ -439,7 +439,7 @@ subroutine set_star_thermalenergy(ieos,den,pres,r,npts,npart,xyzh,vxyzu,rad,eos_
                                   relaxed,use_var_comp,initialtemp,npin)
  use part,            only:do_radiation,rhoh,massoftype,igas,itemp,igasP,iX,iZ,imu,iradxi
  use eos,             only:equationofstate,calc_temp_and_ene,gamma,gmw
- use radiation_utils, only:ugas_from_Tgas,radE_from_Trad
+ use radiation_utils, only:ugas_from_Tgas,radxi_from_Trad
  use table_utils,     only:yinterp
  use units,           only:unit_density,unit_ergg,unit_pressure
  integer, intent(in)    :: ieos,npart,npts
@@ -497,7 +497,7 @@ subroutine set_star_thermalenergy(ieos,den,pres,r,npts,npart,xyzh,vxyzu,rad,eos_
        endif
        if (do_radiation) then
           vxyzu(4,i) = ugas_from_Tgas(tempi,gamma,gmw)
-          rad(iradxi,i) = radE_from_Trad(tempi)/densi
+          rad(iradxi,i) = radxi_from_Trad(densi,tempi)
        else
           vxyzu(4,i) = eni / unit_ergg
        endif
@@ -514,7 +514,7 @@ end subroutine set_star_thermalenergy
 !-----------------------------------------------------------------------
 subroutine solve_uT_profiles(eos_type,r,den,pres,Xfrac,Yfrac,regrid_core,temp,en,mu)
  use eos,     only:get_mean_molecular_weight,calc_temp_and_ene
- use physcon, only:radconst,kb_on_mh
+ use physcon, only:radconst,Rg
  integer, intent(in) :: eos_type
  real, intent(in)    :: r(:),den(:),pres(:),Xfrac(:),Yfrac(:)
  logical, intent(in) :: regrid_core
@@ -531,7 +531,7 @@ subroutine solve_uT_profiles(eos_type,r,den,pres,Xfrac,Yfrac,regrid_core,temp,en
     mu(i) = get_mean_molecular_weight(Xfrac(i),1.-Xfrac(i)-Yfrac(i))  ! only used in u, T calculation if ieos==2,12
     if (i==1) then
        guessene = 1.5*pres(i)/den(i)  ! initial guess
-       tempi = min((3.*pres(i)/radconst)**0.25, pres(i)*mu(i)/(den(i)*kb_on_mh)) ! guess for temperature
+       tempi = min((3.*pres(i)/radconst)**0.25, pres(i)*mu(i)/(den(i)*Rg)) ! guess for temperature
     else
        guessene = en(i-1)
        tempi = temp(i-1)

src/setup/setup_star.f90

@@ -33,7 +33,7 @@ module setup
 !
  use io,             only:fatal,error,warning,master
  use part,           only:gravity,gr
- use physcon,        only:solarm,solarr,km,pi,c,kb_on_mh,radconst
+ use physcon,        only:solarm,solarr,km,pi,c,radconst
  use options,        only:nfulldump,iexternalforce,calc_erot,use_var_comp
  use timestep,       only:tmax,dtmax
  use eos,                only:ieos

src/tests/test_gr.f90

@@ -468,7 +468,7 @@ subroutine test_cons2prim_i(x,v,dens,u,p,ncheck,nfail,errmax,tol)
  use part,            only:ien_entropy,ien_etotal,ien_entropy_s
  use metric_tools,    only:pack_metric,unpack_metric
  use eos,             only:ieos,equationofstate,calc_temp_and_ene
- use physcon,         only:radconst,kb_on_mh
+ use physcon,         only:radconst
 
  real, intent(in) :: x(1:3),v(1:3),dens,p,tol
  real,    intent(inout) :: u

src/tests/test_sedov.F90

@@ -54,7 +54,7 @@ subroutine test_sedov(ntests,npass)
  use mpidomain, only:i_belong
  use checkconserved,  only:etot_in,angtot_in,totmom_in,mdust_in
  use radiation_utils, only:set_radiation_and_gas_temperature_equal,&
-                           T_from_Etot,Tgas_from_ugas,ugas_from_Tgas,radE_from_Trad,Trad_from_radE
+                           T_from_Etot,Tgas_from_ugas,ugas_from_Tgas
  use readwrite_dumps, only:write_fulldump
  use step_lf_global,  only:init_step
  integer, intent(inout) :: ntests,npass