Wind test including radiation force + ieos = 17 + bug fixes

src/main/checksetup.f90

@@ -105,7 +105,7 @@ subroutine check_setup(nerror,nwarn,restart)
        nerror = nerror + 1
     endif
  else
-    if (polyk < tiny(0.) .and. ieos /= 2 .and. ieos /= 5) then
+    if (polyk < tiny(0.) .and. ieos /= 2 .and. ieos /= 5 .and. ieos /= 17) then
        print*,'WARNING! polyk = ',polyk,' in setup, speed of sound will be zero in equation of state'
        nwarn = nwarn + 1
     endif
@@ -239,7 +239,7 @@ subroutine check_setup(nerror,nwarn,restart)
        nerror = nerror + 1
     endif
  else
-    if (abs(gamma-1.) > tiny(gamma) .and. (ieos /= 2 .and. ieos /= 5 .and. ieos /=9)) then
+    if (abs(gamma-1.) > tiny(gamma) .and. (ieos /= 2 .and. ieos /= 5 .and. ieos /=9 .and. ieos /= 17)) then
        print*,'*** ERROR: using isothermal EOS, but gamma = ',gamma
        gamma = 1.
        print*,'*** Resetting gamma to 1, gamma = ',gamma

src/main/dust_formation.f90

@@ -372,7 +372,8 @@ end subroutine evol_K
 !----------------------------------------
 subroutine calc_muGamma(rho_cgs, T, mu, gamma, pH, pH_tot)
 ! all quantities are in cgs
- use io, only:fatal
+ use io,  only:fatal
+ use eos, only:ieos
 
  real, intent(in)    :: rho_cgs
  real, intent(inout) :: T, mu, gamma
@@ -388,8 +389,8 @@ subroutine calc_muGamma(rho_cgs, T, mu, gamma, pH, pH_tot)
  T_old = T
  if (T > 1.d4) then
     mu     = (1.+4.*eps(iHe))/(1.+eps(iHe))
-    gamma  = 5./3.
     pH     = pH_tot
+    if (ieos /= 17) gamma  = 5./3.
  elseif (T > 450.) then
 ! iterate to get consistently pH, T, mu and gamma
     tol       = 1.d-3
@@ -404,14 +405,15 @@ subroutine calc_muGamma(rho_cgs, T, mu, gamma, pH, pH_tot)
        pH        = solve_q(2.*KH2, 1., -pH_tot)
        pH2       = KH2*pH**2
        mu        = (1.+4.*eps(iHe))/(.5+eps(iHe)+0.5*pH/pH_tot)
+       if (ieos == 17) exit !only update mu, keep gamma constant
        x         = 2.*(1.+4.*eps(iHe))/mu
        gamma     = (3.*x+4.+4.*eps(iHe))/(x+4.+4.*eps(iHe))
        converged = (abs(T-T_old)/T_old) < tol
        if (i == 1) then
           mu_old = mu
           gamma_old = gamma
        else
-          T = 2.*T_old*mu/mu_old/(gamma_old-1.)*(x-eps(iHe))/(x+4.-eps(iHe))
+          T = T_old*mu/mu_old/(gamma_old-1.)*2.*x/(x+4.+4.*eps(iHe))
           if (i>=itermax .and. .not.converged) then
              if (isolve==0) then
                 isolve = isolve+1
@@ -431,7 +433,7 @@ subroutine calc_muGamma(rho_cgs, T, mu, gamma, pH, pH_tot)
     pH2    = pH_tot/2.
     pH     = 0.
     mu     = (1.+4.*eps(iHe))/(0.5+eps(iHe))
-    gamma  = (5.*eps(iHe)+3.5)/(3.*eps(iHe)+2.5)
+    if (ieos /= 17)  gamma  = (5.*eps(iHe)+3.5)/(3.*eps(iHe)+2.5)
  endif
 
 end subroutine calc_muGamma

src/main/energies.F90

@@ -364,7 +364,7 @@ subroutine compute_energies(t)
              if (vxyzu(iu,i) < tiny(vxyzu(iu,i))) np_e_eq_0 = np_e_eq_0 + 1
              if (spsoundi < tiny(spsoundi) .and. vxyzu(iu,i) > 0. ) np_cs_eq_0 = np_cs_eq_0 + 1
           else
-             if ((ieos==2  .or. ieos == 5) .and. gammai > 1.001) then
+             if ((ieos==2 .or. ieos == 5 .or. ieos == 17) .and. gammai > 1.001) then
                 !--thermal energy using polytropic equation of state
                 etherm = etherm + pmassi*ponrhoi/(gammai-1.)*gasfrac
              elseif (ieos==9) then

src/main/eos.f90

@@ -12,7 +12,7 @@ module eos
 !     2 = adiabatic/polytropic eos
 !     3 = eos for a locally isothermal disc as in Lodato & Pringle (2007)
 !     4 = GR isothermal
-!     5 = polytropic EOS with vary mu and gamma depending on H2 formation
+!     5 = polytropic EOS with varying mu and gamma depending on H2 formation
 !     6 = eos for a locally isothermal disc as in Lodato & Pringle (2007),
 !         centered on a sink particle
 !     7 = z-dependent locally isothermal eos
@@ -25,6 +25,7 @@ module eos
 !    14 = locally isothermal prescription from Farris et al. (2014) for binary system
 !    15 = Helmholtz free energy eos
 !    16 = Shen eos
+!    17 = polytropic EOS with varying mu (depending on H2 formation)
 !    20 = Ideal gas + radiation + various forms of recombination energy from HORMONE (Hirai et al., 2020)
 !
 ! :References:
@@ -160,7 +161,7 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam
     spsoundi = sqrt(ponrhoi)
     tempi    = temperature_coef*mui*ponrhoi
 
- case(2,5)
+ case(2,5,17)
 !
 !--Adiabatic equation of state (code default)
 !
@@ -816,7 +817,7 @@ end subroutine calc_rec_ene
 !  pressure and density. Inputs and outputs are in cgs units.
 !
 !  Note on composition:
-!  For ieos=2, 5 and 12, mu_local is an input, X & Z are not used
+!  For ieos=2, 5, 12 and 17, mu_local is an input, X & Z are not used
 !  For ieos=10, mu_local is not used
 !  For ieos=20, mu_local is not used but available as an output
 !+
@@ -842,7 +843,7 @@ subroutine calc_temp_and_ene(eos_type,rho,pres,ene,temp,ierr,guesseint,mu_local,
  if (present(X_local))  X  = X_local
  if (present(Z_local))  Z  = Z_local
  select case(eos_type)
- case(2,5) ! Ideal gas
+ case(2,5,17) ! Ideal gas
     temp = pres / (rho * kb_on_mh) * mu
     ene = pres / ( (gamma-1.) * rho)
  case(12) ! Ideal gas + radiation
@@ -865,7 +866,7 @@ end subroutine calc_temp_and_ene
 !  are in cgs units.
 !
 !  Note on composition:
-!  For ieos=2 and 12, mu_local is an input, X & Z are not used
+!  For ieos=2, 5, 12 and 17, mu_local is an input, X & Z are not used
 !  For ieos=10, mu_local is not used
 !  For ieos=20, mu_local is not used but available as an output
 !+
@@ -1040,7 +1041,7 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp)
 
  niter = 0
  select case (ieos)
- case (2,5)
+ case (2,5,17)
     temp = (cgsrho * exp(mu*cgss*mass_proton_cgs))**(2./3.)
     cgsP = cgsrho*kb_on_mh*temp / mu
  case (12)
@@ -1145,7 +1146,7 @@ subroutine setpolyk(eos_type,iprint,utherm,xyzhi,npart)
        write(iprint,*) 'WARNING! different utherms but run is isothermal'
     endif
 
- case(2,5)
+ case(2,5,17)
 !
 !--adiabatic/polytropic eos
 !  this routine is ONLY called if utherm is NOT stored, so polyk matters
@@ -1299,11 +1300,11 @@ subroutine eosinfo(eos_type,iprint)
     endif
  case(3)
     write(iprint,"(/,a,f10.6,a,f10.6)") ' Locally isothermal eq of state (R_sph): cs^2_0 = ',polyk,' qfac = ',qfacdisc
- case(5)
+ case(5,17)
     if (maxvxyzu >= 4) then
        write(iprint,"(' Adiabatic equation of state: P = (gamma-1)*rho*u, where gamma & mu depend on the formation of H2')")
     else
-       stop '[stop eos] eos = 5 cannot assume isothermal conditions'
+       stop '[stop eos] eos = 5,17 cannot assume isothermal conditions'
     endif
  case(6)
     write(iprint,"(/,a,i2,a,f10.6,a,f10.6)") ' Locally (on sink ',isink, &
@@ -1490,7 +1491,7 @@ subroutine read_options_eos(name,valstring,imatch,igotall,ierr)
     read(valstring,*,iostat=ierr) ieos
     ngot = ngot + 1
     if (ieos <= 0 .or. ieos > maxeos) call fatal(label,'equation of state choice out of range')
-    if (ieos == 5) then
+    if (ieos == 5 .or. ieos == 17) then
        store_dust_temperature = .true.
        update_muGamma = .true.
     endif

src/main/inject_wind.f90

@@ -78,7 +78,7 @@ subroutine init_inject(ierr)
  use options,           only:icooling,ieos
  use io,                only:fatal,iverbose
  use setbinary,         only:get_eccentricity_vector
- use timestep,          only:tmax,dtmax
+ use timestep,          only:tmax
  use wind_equations,    only:init_wind_equations
  use wind,              only:setup_wind,save_windprofile
  use physcon,           only:mass_proton_cgs, kboltz, Rg, days, km, au, years, solarm, pi, Gg
@@ -88,12 +88,11 @@ subroutine init_inject(ierr)
  use part,              only:xyzmh_ptmass,vxyz_ptmass,massoftype,igas,iboundary,imloss,ilum,iTeff,iReff,nptmass
  use injectutils,       only:get_sphere_resolution,get_parts_per_sphere,get_neighb_distance
  use cooling_molecular, only:do_molecular_cooling,fit_rho_power,fit_rho_inner,fit_vel,r_compOrb
- use ptmass_radiation,  only:alpha_rad
 
  integer, intent(out) :: ierr
- integer :: ires_min,nzones_per_sonic_point,new_nfill
+ integer :: nzones_per_sonic_point,new_nfill
  real :: mV_on_MdotR,initial_wind_velocity_cgs,dist_to_sonic_point,semimajoraxis_cgs
- real :: dr,dp,mass_of_particles1,tcross,tend,vesc,rsonic,tsonic,initial_Rinject,tboundary
+ real :: dr,dp,mass_of_particles1,tcross,tend,rsonic,tsonic,initial_Rinject,tboundary
  real :: separation_cgs,wind_mass_rate_cgs,wind_velocity_cgs,ecc(3),eccentricity,Tstar
 
  if (icooling > 0) nwrite = nwrite+1
@@ -232,6 +231,7 @@ subroutine init_inject(ierr)
  time_between_spheres  = mass_of_spheres / wind_mass_rate
  massoftype(iboundary) = mass_of_particles
  if (time_between_spheres > tmax)  then
+    call logging(initial_wind_velocity_cgs,rsonic,Tsonic,Tboundary)
     print *,'time_between_spheres = ',time_between_spheres,' < tmax = ',tmax
     call fatal(label,'no shell ejection : tmax < time_between_spheres')
  endif
@@ -268,8 +268,29 @@ subroutine init_inject(ierr)
     print*,'got dr/dp = ',dr/dp,' compared to desired dr on dp = ',wind_shell_spacing
  endif
 
+ xyzmh_ptmass(imloss,wind_emitting_sink) = wind_mass_rate
 
 !logging
+ call logging(initial_wind_velocity_cgs,rsonic,Tsonic,Tboundary)
+
+end subroutine init_inject
+
+
+!-----------------------------------------------------------------------
+
+subroutine logging(initial_wind_velocity_cgs,rsonic,Tsonic,Tboundary)
+
+!-----------------------------------------------------------------------
+
+ use physcon,  only:pi,gg
+ use units,    only:utime,udist
+ use timestep, only:dtmax
+ use ptmass_radiation,  only:alpha_rad
+
+ real, intent(in) :: initial_wind_velocity_cgs,rsonic,Tsonic,Tboundary
+ integer :: ires_min
+ real :: vesc
+
  vesc = sqrt(2.*Gg*Mstar_cgs*(1.-alpha_rad)/Rstar_cgs)
  print*,'mass_of_particles          = ',mass_of_particles
  print*,'particles per sphere       = ',particles_per_sphere
@@ -309,9 +330,8 @@ subroutine init_inject(ierr)
     if (iwind_resolution < ires_min) print *,'WARNING! resolution too low to pass sonic point : iwind_resolution < ',ires_min
  endif
 
- xyzmh_ptmass(imloss,wind_emitting_sink) = wind_mass_rate
+end subroutine logging
 
-end subroutine init_inject
 
 !-----------------------------------------------------------------------
 !+
@@ -635,15 +655,15 @@ subroutine set_default_options_inject(flag)
     wind_mass_rate_Msun_yr = 8.2d-8
     wind_injection_radius_au = 0.
  else
-    !trans-sonic wind
     if (icase == 1) then
+       !trans-sonic wind
        sonic_type = 1
        wind_velocity_km_s = 0.
        wind_mass_rate_Msun_yr = 1.d-5
        wind_injection_radius_au = 2.
        wind_temperature = 50000.
-       !super sonic-wind
     else
+       !super sonic-wind
        sonic_type = 0
        wind_velocity_km_s = 20.
        wind_mass_rate_Msun_yr = 1.d-5

src/main/ptmass.F90

@@ -1212,7 +1212,7 @@ subroutine ptmass_create(nptmass,npart,itest,xyzh,vxyzu,fxyzu,fext,divcurlv,pote
           if (maxvxyzu >= 4) then
              etherm = etherm + pmassj*vxyzu(4,j)
           else
-             if (ieos==2 .and. gamma > 1.001) then
+             if ((ieos==2 .or. ieos==17) .and. gamma > 1.001) then
                 etherm = etherm + pmassj*(eos_vars(igasP,j)/rhoj)/(gamma - 1.)
              elseif (ieos==5 .and. gamma > 1.001) then
                 etherm = etherm + pmassj*(eos_vars(igasP,j)/rhoj)/(eos_vars(igamma,j) - 1.)

src/main/readwrite_infile.F90

@@ -215,7 +215,7 @@ subroutine write_infile(infile,logfile,evfile,dumpfile,iwritein,iprint)
  ! thermodynamics
  !
  call write_options_eos(iwritein)
- if (maxvxyzu >= 4 .and. (ieos==2 .or. ieos==5 .or. ieos==10 .or. ieos==15 .or. ieos==12 .or. ieos==16) ) then
+ if (maxvxyzu >= 4 .and. (ieos==2 .or. ieos==5 .or. ieos==10 .or. ieos==15 .or. ieos==12 .or. ieos==16 .or. ieos==17) ) then
     call write_inopt(ipdv_heating,'ipdv_heating','heating from PdV work (0=off, 1=on)',iwritein)
     call write_inopt(ishock_heating,'ishock_heating','shock heating (0=off, 1=on)',iwritein)
     if (mhd) then
@@ -685,14 +685,15 @@ subroutine read_infile(infile,logfile,evfile,dumpfile)
     if (beta > 4.)     call warn(label,'very high beta viscosity set')
 #ifndef MCFOST
     if (maxvxyzu >= 4 .and. (ieos /= 2 .and. ieos /= 5  .and. ieos /= 4  .and. ieos /= 10 .and. &
-             ieos /=11 .and. ieos /=12 .and. ieos /= 15 .and. ieos /= 16 .and. ieos /= 20)) &
+             ieos /=11 .and. ieos /=12 .and. ieos /= 15 .and. ieos /= 16 .and. ieos /= 17 .and. ieos /= 20)) &
        call fatal(label,'only ieos=2 makes sense if storing thermal energy')
 #endif
     if (irealvisc < 0 .or. irealvisc > 12)  call fatal(label,'invalid setting for physical viscosity')
     if (shearparam < 0.)                     call fatal(label,'stupid value for shear parameter (< 0)')
     if (irealvisc==2 .and. shearparam > 1) call error(label,'alpha > 1 for shakura-sunyaev viscosity')
     if (iverbose > 99 .or. iverbose < -9)   call fatal(label,'invalid verboseness setting (two digits only)')
-    if (icooling > 0 .and. .not.(ieos == 2 .or. ieos == 5)) call fatal(label,'cooling requires adiabatic eos (ieos=2)')
+    if (icooling > 0 .and. .not.(ieos == 2 .or. ieos == 5 .or. ieos == 17)) &
+         call fatal(label,'cooling requires adiabatic eos (ieos=2)')
     if (icooling > 0 .and. (ipdv_heating <= 0 .or. ishock_heating <= 0)) &
          call fatal(label,'cooling requires shock and work contributions')
     if (((isink_radiation == 1 .or. isink_radiation == 3 ) .and. idust_opacity == 0 ) &

src/main/substepping.F90

@@ -1014,7 +1014,7 @@ subroutine cooling_abundances_update(i,pmassi,xyzh,vxyzu,eos_vars,abundance,nucl
  use cooling_ism,     only:nabn,dphotflag
  use options,         only:icooling
  use chem,            only:update_abundances,get_dphot
- use dust_formation,  only:evolve_dust
+ use dust_formation,  only:evolve_dust,calc_muGamma
  use cooling,         only:energ_cooling,cooling_in_step
  use part,            only:rhoh
 #ifdef KROME
@@ -1031,7 +1031,7 @@ subroutine cooling_abundances_update(i,pmassi,xyzh,vxyzu,eos_vars,abundance,nucl
  real,         intent(in)    :: dt,pmassi
  integer,      intent(in)    :: i
 
- real :: dudtcool,rhoi,dphot
+ real :: dudtcool,rhoi,dphot,pH,pH_tot
  real :: abundi(nabn)
 
  dudtcool = 0.
@@ -1059,6 +1059,8 @@ subroutine cooling_abundances_update(i,pmassi,xyzh,vxyzu,eos_vars,abundance,nucl
     call evolve_dust(dt, xyzh(:,i), vxyzu(4,i), nucleation(:,i), dust_temp(i), rhoi)
     eos_vars(imu,i)    = nucleation(idmu,i)
     eos_vars(igamma,i) = nucleation(idgamma,i)
+ elseif (update_muGamma) then
+    call calc_muGamma(rhoi, dust_temp(i),eos_vars(imu,i),eos_vars(igamma,i), pH, pH_tot)
  endif
  !
  ! COOLING

src/main/wind.F90

@@ -207,8 +207,7 @@ subroutine wind_step(state)
  use cooling_solver,   only:calc_cooling_rate
  use options,          only:icooling
  use units,            only:unit_ergg,unit_density
- use dim,              only:itau_alloc
- use eos,              only:ieos
+ use dim,              only:itau_alloc,update_muGamma
 
  type(wind_state), intent(inout) :: state
  real :: rvT(3), dt_next, v_old, dlnQ_dlnT, Q_code, pH, pH_tot
@@ -241,9 +240,9 @@ subroutine wind_step(state)
     state%gamma          = state%JKmuS(idgamma)
     state%kappa          = calc_kappa_dust(state%JKmuS(idK3), state%Tdust, state%rho)
     state%JKmuS(idalpha) = state%alpha_Edd+alpha_rad
- elseif (idust_opacity == 1) then
-    state%kappa = calc_kappa_bowen(state%Tdust)
-    if (ieos == 5) call calc_muGamma(state%rho, state%Tg,state%mu, state%gamma, pH, pH_tot)
+ else
+    if (idust_opacity == 1) state%kappa = calc_kappa_bowen(state%Tdust)
+    if (update_muGamma) call calc_muGamma(state%rho, state%Tg,state%mu, state%gamma, pH, pH_tot)
  endif
 
  if (itau_alloc == 1) then
@@ -345,13 +344,12 @@ subroutine wind_step(state)
  use ptmass_radiation, only:alpha_rad,iget_tdust,tdust_exp, isink_radiation
  use physcon,          only:pi,Rg
  use dust_formation,   only:evolve_chem,calc_kappa_dust,calc_kappa_bowen,&
-      calc_Eddington_factor,idust_opacity, calc_mugamma
+      calc_Eddington_factor,idust_opacity, calc_muGamma
  use part,             only:idK3,idmu,idgamma,idsat,idkappa
  use cooling_solver,   only:calc_cooling_rate
  use options,          only:icooling
  use units,            only:unit_ergg,unit_density
- use dim,              only:itau_alloc
- use eos,              only:ieos
+ use dim,              only:itau_alloc,update_muGamma
 
  type(wind_state), intent(inout) :: state
  real :: rvT(3), dt_next, v_old, dlnQ_dlnT, Q_code, pH,pH_tot
@@ -365,9 +363,9 @@ subroutine wind_step(state)
     state%mu        = state%JKmuS(idmu)
     state%gamma     = state%JKmuS(idgamma)
     state%kappa     = calc_kappa_dust(state%JKmuS(idK3), state%Tdust, state%rho)
- elseif (idust_opacity == 1) then
-    state%kappa     = calc_kappa_bowen(state%Tdust)
-    if (ieos == 5 ) call calc_muGamma(state%rho, state%Tg,state%mu, state%gamma, pH, pH_tot)
+ else
+    if (idust_opacity == 1) state%kappa     = calc_kappa_bowen(state%Tdust)
+    if (update_muGamma) call calc_muGamma(state%rho, state%Tg,state%mu, state%gamma, pH, pH_tot)
  endif
 
  if (itau_alloc == 1) then

src/main/wind_equations.f90

@@ -290,6 +290,7 @@ end subroutine RK4_step_dr
 subroutine calc_dvT_dr(r, v, T0, Rstar_cgs, Mdot_cgs, mu0, gamma0, alpha, dalpha_dr, Q, dQ_dr, dv_dr, dT_dr, numerator, denominator)
 !all quantities in cgs
  use physcon, only:Gg,Rg,pi
+ use dim,     only:update_muGamma
  use options, only:icooling,ieos
  use dust_formation,   only:calc_muGamma,idust_opacity
  real, intent(in) :: r, v, T0, mu0, gamma0, alpha, dalpha_dr, Q, dQ_dr, Rstar_cgs, Mdot_cgs
@@ -302,7 +303,7 @@ subroutine calc_dvT_dr(r, v, T0, Rstar_cgs, Mdot_cgs, mu0, gamma0, alpha, dalpha
  T = T0
  mu = mu0
  gamma = gamma0
- if (idust_opacity == 2) then
+ if (update_muGamma .or. idust_opacity == 2) then
     rho_cgs = Mdot_cgs/(4.*pi*r**2*v)
     call calc_muGamma(rho_cgs, T, mu, gamma, pH, pH_tot)
  endif