Merge pull request #440 from lsiess/master

Bug Fix

docs/CE.rst

@@ -73,7 +73,7 @@ Use SETUP=star or SETUP=dustystar and if not specified, the default options.
   2.2 make setup
   2.3 ./phantomsetup star (option 5 MESA star, input profile = Jan_Star_Phantom_Profile.data, desired EOS = 10,
       use constant entropy profile, Relax star automatically = yes). The core radius is the softening radius (2-3Ro)
-      the core mass is the same as the one you have measured from MESA (0.46Mo in Jan_Star_Phantom_Profile).
+      the core mass is the same as the one you have measured from MESA (0.46Mo in Jan_Star_Phantom_Profile.data).
       This produces a file called star.setup - this file has all the options so you can edit it.
   2.4 vim  star.setup, (write_rho_to_file = T)
 Relaxation
@@ -124,3 +124,13 @@ if you come from 2.10, then use as initial model (hereafter initial_nnnnn) one o
         softening length for the primary core = 1., softening length for companion = 0.1)
   2.16 vim binary.in (optional, tmax=200.00, dtmax=0.100)
   2.17 ./phantom binary.in
+
+
+**D. Setup sink properties (luminosity)**
+
+::
+
+  2.18 ./phantommoddump binary_00000.tmp dusty_binary_00000.tmp 0.0
+       option 9, 12 lum  
+  2.19 vim dusty_binary.in (adapt isink_radiation, idust_opacity)
+  2.20 ./phantom dusty_binary.in

src/main/checksetup.F90

@@ -512,7 +512,8 @@ end function in_range
 !------------------------------------------------------------------
 subroutine check_setup_ptmass(nerror,nwarn,hmin)
  use dim,  only:maxptmass
- use part, only:nptmass,xyzmh_ptmass,ihacc,ihsoft,gr,iTeff,sinks_have_luminosity
+ use part, only:nptmass,xyzmh_ptmass,ihacc,ihsoft,gr,iTeff,sinks_have_luminosity,ilum
+ use ptmass_radiation, only:isink_radiation
  integer, intent(inout) :: nerror,nwarn
  real,    intent(in)    :: hmin
  integer :: i,j,n
@@ -589,6 +590,11 @@ subroutine check_setup_ptmass(nerror,nwarn,hmin)
  !
  !  check that radiation properties are sensible
  !
+ if (isink_radiation > 1 .and. xyzmh_ptmass(ilum,1) < 1e-10) then
+    nerror = nerror + 1
+    print*,'ERROR: isink_radiation > 1 and sink particle has no luminosity'
+    return
+ endif
  if (sinks_have_luminosity(nptmass,xyzmh_ptmass)) then
     if (any(xyzmh_ptmass(iTeff,1:nptmass) < 100.)) then
        print*,'WARNING: sink particle temperature less than 100K'

src/main/dust_formation.f90

@@ -32,7 +32,7 @@ module dust_formation
  integer, public :: idust_opacity = 0
 
  public :: set_abundances,evolve_dust,evolve_chem,calc_kappa_dust,&
-      calc_kappa_bowen,chemical_equilibrium_light,&
+      calc_kappa_bowen,chemical_equilibrium_light,psat_C,calc_nucleation,&
       read_options_dust_formation,write_options_dust_formation,&
       calc_Eddington_factor,calc_muGamma,init_muGamma,init_nucleation,&
       write_headeropts_dust_formation,read_headeropts_dust_formation
@@ -41,6 +41,7 @@ module dust_formation
 !
 
  real, public :: kappa_gas   = 2.d-4
+ real, public, parameter :: Scrit = 2. ! Critical saturation ratio
 
  private
 
@@ -82,7 +83,6 @@ module dust_formation
       -4.38897d+05, -1.58111d+05, 2.49224d+01, 1.08714d-03, -5.62504d-08, & !TiO
       -3.32351d+05, -3.04694d+05, 5.86984d+01, 1.17096d-03, -5.06729d-08, & !TiO2
        2.26786d+05, -1.43775d+05, 2.92429d+01, 1.69434d-04, -1.79867d-08], shape(coefs)) !C2
- real, parameter :: Scrit = 2. ! Critical saturation ratio
  real, parameter :: vfactor = sqrt(kboltz/(2.*pi*atomic_mass_unit*12.01))
  !real, parameter :: vfactor = sqrt(kboltz/(8.*pi*atomic_mass_unit*12.01))
 
@@ -179,7 +179,7 @@ subroutine evolve_chem(dt, T, rho_cgs, JKmuS)
     S = pC/psat_C(T)
     if (S > Scrit) then
        !call nucleation(T, pC, pC2, 0., pC2H, pC2H2, S, JstarS, taustar, taugr)
-       call nucleation(T, pC, 0., 0., 0., pC2H2, S, JstarS, taustar, taugr)
+       call calc_nucleation(T, pC, 0., 0., 0., pC2H2, S, JstarS, taustar, taugr)
        JstarS = JstarS/ nH_tot
        call evol_K(JKmuS(idJstar), JKmuS(idK0:idK3), JstarS, taustar, taugr, dt, Jstar_new, K_new)
     else
@@ -272,7 +272,7 @@ end function calc_Eddington_factor
 !  Compute nucleation rate
 !
 !----------------------------
-subroutine nucleation(T, pC, pC2, pC3, pC2H, pC2H2, S, JstarS, taustar, taugr)
+subroutine calc_nucleation(T, pC, pC2, pC3, pC2H, pC2H2, S, JstarS, taustar, taugr)
 ! all quantities are in cgs
  real, intent(in)  :: T, pC, pC2, pC3, pC2H, pC2H2, S
  real, intent(out) :: JstarS, taustar, taugr
@@ -320,7 +320,7 @@ subroutine nucleation(T, pC, pC2, pC3, pC2H, pC2H2, S, JstarS, taustar, taugr)
     taustar = 1.d-30
  endif
  taugr = kboltz*T/(A0*v1*(alpha1*pC*(1.-1./S) + 2.*alpha2/sqrt(2.)*(pC2+pC2H+pC2H2)*(1.-1./S**2)))
-end subroutine nucleation
+end subroutine calc_nucleation
 
 !------------------------------------
 !

src/main/inject_wind.F90

@@ -28,7 +28,6 @@ module inject
 !   infile_utils, injectutils, io, options, part, partinject, physcon,
 !   ptmass_radiation, setbinary, timestep, units, wind, wind_equations
 !
- use physcon,           only:solarl
  use dim,               only:isothermal
 
  implicit none

src/main/readwrite_dumps_common.F90

@@ -125,7 +125,7 @@ subroutine check_arrays(i1,i2,noffset,npartoftype,npartread,nptmass,nsinkpropert
  use dim,  only:maxp,maxvxyzu,maxalpha,maxBevol,mhd,h2chemistry,use_dustgrowth,gr,&
                 do_radiation,store_dust_temperature,do_nucleation
  use eos,  only:ieos,polyk,gamma,eos_is_non_ideal
- use part, only:maxphase,isetphase,set_particle_type,igas,ihacc,ihsoft,imacc,&
+ use part, only:maxphase,isetphase,set_particle_type,igas,ihacc,ihsoft,imacc,ilum,&
                 xyzmh_ptmass_label,vxyz_ptmass_label,get_pmass,rhoh,dustfrac,ndusttypes,norig
  use io,   only:warning,id,master
  use options,        only:alpha,use_dustfrac,use_var_comp
@@ -322,13 +322,13 @@ subroutine check_arrays(i1,i2,noffset,npartoftype,npartread,nptmass,nsinkpropert
     endif
     if (id==master .and. i1==1) then
        print "(2(a,i4),a)",' got ',nsinkproperties,' sink properties from ',nptmass,' sink particles'
-       if (nptmass > 0) print "(1x,47('-'),/,1x,a,'|',4(a9,1x,'|'),/,1x,47('-'))",&
-                              'ID',' Mass    ',' Racc    ',' Macc    ',' hsoft   '
+       if (nptmass > 0) print "(1x,58('-'),/,1x,a,'|',5(a9,1x,'|'),/,1x,58('-'))",&
+                              'ID',' Mass    ',' Racc    ',' Macc    ',' hsoft   ',' Lsink   '
        do i=1,min(nptmass,999)
-          if (xyzmh_ptmass(4,i) > 0.) print "(i3,'|',4(1pg9.2,1x,'|'))",i,xyzmh_ptmass(4,i), &
-                                            xyzmh_ptmass(ihacc,i),xyzmh_ptmass(imacc,i),xyzmh_ptmass(ihsoft,i)
+          if (xyzmh_ptmass(4,i) > 0.) print "(i3,'|',5(1pg9.2,1x,'|'))",i,xyzmh_ptmass(4,i),xyzmh_ptmass(ihacc,i),&
+                                            xyzmh_ptmass(imacc,i),xyzmh_ptmass(ihsoft,i),xyzmh_ptmass(ilum,i)
        enddo
-       if (nptmass > 0) print "(1x,47('-'))"
+       if (nptmass > 0) print "(1x,58('-'))"
     endif
  endif
  !

src/main/units.f90

@@ -35,7 +35,7 @@ module units
  public :: get_G_code, get_c_code, get_radconst_code, get_kbmh_code
  public :: c_is_unity, G_is_unity, in_geometric_units
  public :: is_time_unit, is_length_unit
- public :: in_solarr, in_solarm
+ public :: in_solarr, in_solarm, in_solarl
 
 contains
 
@@ -464,5 +464,17 @@ real(kind=8) function in_solarr(val) result(rval)
  rval = val*(udist/solarr)
 
 end function in_solarr
+!---------------------------------------------------------------------------
+!+
+!  function to convert a luminosity value from code units to solar luminosity
+!+
+!---------------------------------------------------------------------------
+real(kind=8) function in_solarl(val) result(rval)
+ use physcon, only:solarl
+ real, intent(in) :: val
+
+ rval = val*(unit_luminosity/solarl)
+
+end function in_solarl
 
 end module units

src/setup/set_star.f90

@@ -41,6 +41,7 @@ module setstar
     real :: initialtemp
     real :: rcore
     real :: mcore
+    real :: lcore
     real :: hsoft
     real :: hacc   ! accretion radius if star is a sink particle
     character(len=120) :: input_profile,dens_profile
@@ -80,6 +81,7 @@ subroutine set_defaults_star(star)
  star%hacc           = 1.
  star%rcore          = 0.
  star%mcore          = 0.
+ star%lcore          = 0.
  star%isofteningopt  = 1 ! By default, specify rcore
  star%np             = 1000
  star%input_profile  = 'P12_Phantom_Profile.data'
@@ -108,7 +110,7 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,&
                               write_kepler_comp
  use radiation_utils,    only:set_radiation_and_gas_temperature_equal
  use relaxstar,          only:relax_star
- use part,               only:ihsoft,igas,imu,set_particle_type
+ use part,               only:ihsoft,igas,imu,set_particle_type,ilum
  use extern_densprofile, only:write_rhotab
  use unifdis,            only:mask_prototype
  use physcon,            only:pi
@@ -185,10 +187,11 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,&
  !
  ! add sink particle stellar core
  !
- if (star%isinkcore) call set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,&
-                                           ihsoft,star%mcore,star%hsoft,ierr)
+ if (star%isinkcore) call set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,ihsoft,&
+                                           star%mcore,star%hsoft,ilum,star%lcore,ierr)
  if (ierr==1) call fatal('set_stellar_core','mcore <= 0')
  if (ierr==2) call fatal('set_stellar_core','hsoft <= 0')
+ if (ierr==3) call fatal('set_stellar_core','lcore < 0')
  !
  ! Write the desired profile to file (do this before relaxation)
  !
@@ -440,22 +443,23 @@ end subroutine set_defaults_given_profile
 subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk)
  use prompting,     only:prompt
  use setstar_utils, only:nprofile_opts,profile_opt,need_inputprofile,need_rstar
- use units,         only:in_solarm,in_solarr,udist,umass
- use physcon,       only:solarr,solarm
+ use units,         only:in_solarm,in_solarr,in_solarl,udist,umass,unit_luminosity
+ use physcon,       only:solarr,solarm,solarl
  type(star_t), intent(out)   :: star
  integer,      intent(in)    :: id,master
  logical,      intent(out)   :: use_var_comp
  integer,      intent(out)   :: need_iso
  integer,      intent(inout) :: ieos
  real,         intent(inout) :: polyk
  integer :: i
- real :: mstar_msun,rstar_rsun,rcore_rsun,mcore_msun,hsoft_rsun
+ real :: mstar_msun,rstar_rsun,rcore_rsun,mcore_msun,lcore_lsun,hsoft_rsun
 
  ! set defaults
  call set_defaults_star(star)
  mstar_msun = real(in_solarm(star%mstar))
  rstar_rsun = real(in_solarr(star%rstar))
  mcore_msun = real(in_solarm(star%mcore))
+ lcore_lsun = real(in_solarl(star%lcore))
  rcore_rsun = real(in_solarr(star%rcore))
  hsoft_rsun = real(in_solarr(star%hsoft))
 
@@ -506,8 +510,10 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk)
        if (star%isinkcore) then
           call prompt('Enter mass of the created sink particle core [Msun]',mcore_msun,0.)
           call prompt('Enter softening length of the sink particle core [Rsun]',hsoft_rsun,0.)
+          call prompt('Enter sink particle luminosity [Lsun]',lcore_lsun,0.)
           star%mcore = mcore_msun*real(solarm/umass)
           star%hsoft = hsoft_rsun*real(solarr/udist)
+          star%lcore = lcore_lsun*real(solarl/unit_luminosity)
        endif
     case(1)
        star%isinkcore = .true. ! Create sink particle core automatically
@@ -532,19 +538,20 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk)
           star%mcore = mcore_msun*real(solarm/umass)
           star%rcore = rcore_rsun*real(solarr/udist)
        end select
-
-       call prompt('Enter output file name of cored stellar profile:',star%outputfilename)
+       call prompt('Enter sink particle luminosity [Lsun]',lcore_lsun,0.)
+       star%lcore = lcore_lsun*real(solarl/unit_luminosity)
 
     case(2)
        star%isinkcore = .true. ! Create sink particle core automatically
        print*,'Specify core radius and initial guess for mass of sink particle core'
        call prompt('Enter core radius in Rsun : ',rcore_rsun,0.)
        call prompt('Enter guess for core mass in Msun : ',mcore_msun,0.)
+       call prompt('Enter sink particle luminosity [Lsun]',lcore_lsun,0.)
        call prompt('Enter output file name of cored stellar profile:',star%outputfilename)
        star%mcore = mcore_msun*real(solarm/umass)
        star%rcore = rcore_rsun*real(solarr/udist)
+       star%lcore = lcore_lsun*real(solarl/unit_luminosity)
     end select
-
  case(ievrard)
     call prompt('Enter the specific internal energy (units of GM/R) ',star%ui_coef,0.)
  case(:0)
@@ -561,7 +568,7 @@ end subroutine set_star_interactive
 subroutine write_options_star(star,iunit,label)
  use infile_utils,  only:write_inopt,get_optstring
  use setstar_utils, only:nprofile_opts,profile_opt,need_inputprofile,need_rstar
- use units,         only:in_solarm,in_solarr
+ use units,         only:in_solarm,in_solarr,in_solarl
  type(star_t),     intent(in) :: star
  integer,          intent(in) :: iunit
  character(len=*), intent(in), optional :: label
@@ -614,6 +621,8 @@ subroutine write_options_star(star,iunit,label)
           call write_inopt(in_solarm(star%mcore),'mcore'//trim(c),&
                'Initial guess for mass of sink particle stellar core [Msun]',iunit)
        endif
+       call write_inopt(in_solarl(star%lcore),'lcore'//trim(c),&
+                              'Luminosity of point mass stellar core [Lsun]',iunit)
     else
        call write_inopt(star%isinkcore,'isinkcore'//trim(c),&
                'Add a sink particle stellar core',iunit)
@@ -623,6 +632,8 @@ subroutine write_options_star(star,iunit,label)
           call write_inopt(in_solarr(star%hsoft),'hsoft'//trim(c),&
                'Softening length of sink particle stellar core [Rsun]',iunit)
        endif
+       call write_inopt(in_solarl(star%lcore),'lcore'//trim(c),&
+               'Luminosity of sink core particle [Lsun]',iunit)
     endif
  case (ievrard)
     call write_inopt(star%ui_coef,'ui_coef'//trim(c),&
@@ -646,8 +657,8 @@ end subroutine write_options_star
 subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label)
  use infile_utils,  only:inopts,read_inopt
  use setstar_utils, only:need_inputprofile,need_rstar,nprofile_opts
- use units,         only:umass,udist
- use physcon,       only:solarm,solarr
+ use units,         only:umass,udist,unit_luminosity
+ use physcon,       only:solarm,solarr,solarl
  type(star_t),              intent(out)   :: star
  type(inopts), allocatable, intent(inout) :: db(:)
  integer,                   intent(out)   :: need_iso
@@ -656,7 +667,7 @@ subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label)
  integer,                   intent(inout) :: nerr
  character(len=*),          intent(in), optional :: label
  character(len=10) :: c
- real :: mcore_msun,rcore_rsun,mstar_msun,rstar_rsun,hsoft_rsun
+ real :: mcore_msun,rcore_rsun,lcore_lsun,mstar_msun,rstar_rsun,hsoft_rsun
 
  ! set defaults
  call set_defaults_star(star)
@@ -709,6 +720,8 @@ subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label)
           call read_inopt(mcore_msun,'mcore'//trim(c),db,errcount=nerr,min=0.)
           star%mcore = mcore_msun*real(solarm/umass)
        endif
+       call read_inopt(lcore_lsun,'lcore'//trim(c),db,errcount=nerr,min=0.)
+       star%lcore = lcore_lsun*real(solarl/unit_luminosity)
     endif
  case(ievrard)
     call read_inopt(star%ui_coef,'ui_coef'//trim(c),db,errcount=nerr,min=0.)

src/setup/set_star_utils.f90

@@ -310,11 +310,13 @@ end subroutine set_star_density
 !  Add a sink particle as a stellar core
 !+
 !-----------------------------------------------------------------------
-subroutine set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,ihsoft,mcore,hsoft,ierr)
+subroutine set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,ihsoft,mcore,&
+                            hsoft,ilum,lcore,ierr)
  integer, intent(out) :: nptmass,ierr
  real, intent(out)    :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:)
- real, intent(in)     :: mcore,hsoft
- integer              :: n,ihsoft
+ real, intent(in)     :: mcore,hsoft,lcore
+ integer, intent(in)  :: ihsoft,ilum
+ integer              :: n
 
  ierr = 0
  ! Check for sensible values
@@ -326,12 +328,17 @@ subroutine set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,ihsoft,mcore,hsoft,
     ierr = 2
     return
  endif
+ if (lcore < 0.) then
+    ierr = 3
+    return
+ endif
 
  nptmass                = 1
  n                      = nptmass
  xyzmh_ptmass(:,n)      = 0. ! zero all quantities by default
  xyzmh_ptmass(4,n)      = mcore
  xyzmh_ptmass(ihsoft,n) = hsoft
+ xyzmh_ptmass(ilum,n)   = lcore
  vxyz_ptmass(:,n)       = 0.
 
 end subroutine set_stellar_core