Merge pull request #495 from danieljprice/star

Enhancements and cleanup related to star setup

src/main/eos.f90

@@ -59,7 +59,7 @@ module eos
  public  :: eos_is_non_ideal,eos_outputs_mu,eos_outputs_gasP
  public  :: get_local_u_internal,get_temperature_from_u
  public  :: calc_rec_ene,calc_temp_and_ene,entropy,get_rho_from_p_s,get_u_from_rhoT
- public  :: get_entropy,get_p_from_rho_s
+ public  :: calc_rho_from_PT,get_entropy,get_p_from_rho_s
  public  :: init_eos,finish_eos,write_options_eos,read_options_eos
  public  :: write_headeropts_eos, read_headeropts_eos
 
@@ -450,6 +450,7 @@ subroutine init_eos(eos_type,ierr)
  use dim,            only:maxvxyzu,do_radiation
  integer, intent(in)  :: eos_type
  integer, intent(out) :: ierr
+ integer              :: ierr_mesakapp
 
  ierr = 0
  !
@@ -526,7 +527,11 @@ subroutine init_eos(eos_type,ierr)
  end select
  done_init_eos = .true.
 
- if (do_radiation .and. iopacity_type==1) call init_eos_mesa(X_in,Z_in,ierr)
+ if (do_radiation .and. iopacity_type==1) then
+    write(*,'(1x,a,f7.5,a,f7.5)') 'Using radiation with MESA opacities. Initialising MESA EoS with X = ',X_in,', Z = ',Z_in
+    call init_eos_mesa(X_in,Z_in,ierr_mesakapp)
+    ierr = max(ierr,ierr_mesakapp)
+ endif
 
 end subroutine init_eos
 
@@ -854,6 +859,48 @@ subroutine calc_temp_and_ene(eos_type,rho,pres,ene,temp,ierr,guesseint,mu_local,
 
 end subroutine calc_temp_and_ene
 
+!-----------------------------------------------------------------------
+!+
+!  Calculate density from pressure and temperature. Inputs and outputs
+!  are in cgs units.
+!
+!  Note on composition:
+!  For ieos=2 and 12, 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
+!+
+!-----------------------------------------------------------------------
+subroutine calc_rho_from_PT(eos_type,pres,temp,rho,ierr,mu_local,X_local,Z_local)
+ use physcon,          only:kb_on_mh
+ 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
+ integer, intent(in)              :: eos_type
+ real,    intent(in)              :: pres,temp
+ real,    intent(inout)           :: rho
+ real,    intent(in),    optional :: X_local,Z_local
+ real,    intent(inout), optional :: mu_local
+ integer, intent(out)             :: ierr
+ real                             :: mu,X,Z
+
+ ierr = 0
+ mu   = gmw
+ X    = X_in
+ Z    = Z_in
+ if (present(mu_local)) mu = mu_local
+ if (present(X_local))  X  = X_local
+ if (present(Z_local))  Z  = Z_local
+ select case(eos_type)
+ case(2) ! Ideal gas
+    rho = pres / (temp * kb_on_mh) * mu
+ case(12) ! Ideal gas + radiation
+    call get_idealplusrad_rhofrompresT(pres,temp,mu,rho)
+ case default
+    ierr = 1
+ end select
+
+end subroutine calc_rho_from_PT
+
 !-----------------------------------------------------------------------
 !+
 !  Calculates specific entropy (gas + radiation + recombination)

src/main/eos_idealplusrad.f90

@@ -23,7 +23,8 @@ module eos_idealplusrad
  real, parameter :: tolerance = 1e-15
 
  public :: get_idealplusrad_temp,get_idealplusrad_pres,get_idealplusrad_spsoundi,&
-           get_idealgasplusrad_tempfrompres,get_idealplusrad_enfromtemp
+           get_idealgasplusrad_tempfrompres,get_idealplusrad_enfromtemp,&
+           get_idealplusrad_rhofrompresT
 
  private
 
@@ -130,4 +131,18 @@ subroutine get_idealplusrad_enfromtemp(densi,tempi,mu,eni)
 
 end subroutine get_idealplusrad_enfromtemp
 
+
+!----------------------------------------------------------------
+!+
+!  Calculates density from pressure and temperature
+!+
+!----------------------------------------------------------------
+subroutine get_idealplusrad_rhofrompresT(presi,tempi,mu,densi)
+ real, intent(in)  :: presi,tempi,mu
+ real, intent(out) :: densi
+
+ densi = (presi - radconst*tempi**4 /3.) * mu / (Rg*tempi)
+
+end subroutine get_idealplusrad_rhofrompresT
+
 end module eos_idealplusrad

src/main/part.F90

@@ -759,6 +759,19 @@ logical function sinks_have_heating(nptmass,xyzmh_ptmass)
 
 end function sinks_have_heating
 
+!------------------------------------------------------------------------
+!+
+!  Query function to see if any sink particles have heating
+!+
+!------------------------------------------------------------------------
+logical function sink_has_heating(xyzmh_ptmassi)
+ real, intent(in) :: xyzmh_ptmassi(:)
+
+ sink_has_heating = xyzmh_ptmassi(iTeff) <= 0. .and. &
+                              xyzmh_ptmassi(ilum) > 0.
+
+end function sink_has_heating
+
 !----------------------------------------------------------------
 !+
 !  query function returning whether or not a particle is dead

src/main/ptmass.F90

@@ -1778,26 +1778,38 @@ end subroutine calculate_mdot
 
 !-----------------------------------------------------------------------
 !+
-!  calculate mass enclosed in sink softening radius
+!  calculate (weighted) sum of particle mass enclosed in sink softening radius
 !+
 !-----------------------------------------------------------------------
 subroutine ptmass_calc_enclosed_mass(nptmass,npart,xyzh)
- use part,   only:imassenc,ihsoft,massoftype,igas,xyzmh_ptmass
- use kernel, only:radkern2
+ use part, only:sink_has_heating,imassenc,ihsoft,massoftype,igas,xyzmh_ptmass,isdead_or_accreted
+ use ptmass_heating, only:isink_heating,heating_kernel
+ use kernel,         only:radkern2
  integer, intent(in) :: nptmass,npart
  real,    intent(in) :: xyzh(:,:)
- integer             :: i,j,ncount
- real                :: drj2
+ integer             :: i,j
+ real                :: wi,q2,x0,y0,z0,hsoft21
 
  do i = 1,nptmass
-    ncount = 0
+    if (.not. sink_has_heating(xyzmh_ptmass(:,i))) cycle
+    wi = 0.
+    x0 = xyzmh_ptmass(1,i)
+    y0 = xyzmh_ptmass(2,i)
+    z0 = xyzmh_ptmass(3,i)
+    hsoft21 = 1./xyzmh_ptmass(ihsoft,i)**2
+
+    !$omp parallel do default (none) &
+    !$omp reduction(+:wi) &
+    !$omp shared(npart,xyzh,x0,y0,z0,i,hsoft21,isink_heating) &
+    !$omp private(j,q2)
     do j = 1,npart
-       drj2 = (xyzh(1,j)-xyzmh_ptmass(1,i))**2 + (xyzh(2,j)-xyzmh_ptmass(2,i))**2 + (xyzh(3,j)-xyzmh_ptmass(3,i))**2
-       if (drj2 < radkern2*xyzmh_ptmass(ihsoft,i)**2) then
-          ncount = ncount + 1
+       if (.not. isdead_or_accreted(xyzh(4,j))) then
+          q2 = ((xyzh(1,j)-x0)**2 + (xyzh(2,j)-y0)**2 + (xyzh(3,j)-z0)**2)*hsoft21
+          if (q2 < radkern2) wi = wi + heating_kernel(q2,isink_heating)  ! wj = 1 for uniform heating
        endif
     enddo
-    xyzmh_ptmass(imassenc,i) = ncount * massoftype(igas)
+    !$omp end parallel do
+    xyzmh_ptmass(imassenc,i) = wi * massoftype(igas)
  enddo
 
 end subroutine ptmass_calc_enclosed_mass

src/main/ptmass_heating.f90

@@ -6,7 +6,7 @@
 !--------------------------------------------------------------------------!
 module ptmass_heating
 !
-! Heating of particles around softening radius of sink particle
+! Heating of particles around sink particles
 !
 ! :References: None
 !
@@ -18,14 +18,16 @@ module ptmass_heating
 !
 
  implicit none
- public :: energ_sinkheat
- real, public :: Lnuc
+ public :: energ_sinkheat,heating_kernel
+ real,    public :: Lnuc
+ integer, public :: isink_heating
+
  private
 
 contains
 !-----------------------------------------------------------------------
 !+
-!  Heat from point mass
+!  heating from point mass
 !+
 !-----------------------------------------------------------------------
 subroutine energ_sinkheat(nptmass,xyzmh_ptmass,xi,yi,zi,dudtheati)
@@ -35,17 +37,82 @@ subroutine energ_sinkheat(nptmass,xyzmh_ptmass,xi,yi,zi,dudtheati)
  real, intent(in)    :: xi,yi,zi,xyzmh_ptmass(:,:)
  real, intent(out)   :: dudtheati
  integer             :: i
- real                :: dri2
+ real                :: q2,dri2
 
  dudtheati = 0.
  do i = 1,nptmass
     dri2 = (xi-xyzmh_ptmass(1,i))**2 + (yi-xyzmh_ptmass(2,i))**2 + (zi-xyzmh_ptmass(3,i))**2
-    if (dri2 < radkern2*xyzmh_ptmass(ihsoft,i)**2) then
-       dudtheati = xyzmh_ptmass(iLum,i) / xyzmh_ptmass(imassenc,i)
+    q2 = dri2/xyzmh_ptmass(ihsoft,i)**2
+    if (q2 < radkern2) then
+       dudtheati = xyzmh_ptmass(iLum,i) / xyzmh_ptmass(imassenc,i) * heating_kernel(q2,isink_heating)
     endif
  enddo
 
 end subroutine energ_sinkheat
 
 
+!-----------------------------------------------------------------------
+!+
+!  heating weight function (note: arbitrary normalisation)
+!+
+!-----------------------------------------------------------------------
+real function heating_kernel(q2,kernel_type)
+ use kernel, only:wkern,cnormk
+ real, intent(in)    :: q2
+ integer, intent(in) :: kernel_type
+
+ select case(kernel_type)
+ case(1)
+    heating_kernel = cnormk*wkern(q2,sqrt(q2))
+ case default
+    heating_kernel = 1.
+ end select
+
+end function heating_kernel
+
+
+!-----------------------------------------------------------------------
+!+
+!  write options to input file (not used at the moment)
+!+
+!-----------------------------------------------------------------------
+! subroutine write_options_ptmass_heating(iunit)
+!  use infile_utils, only:write_inopt
+!  integer, intent(in) :: iunit
+
+!  call write_inopt(isink_heating,'isink_heating','sink heating distirbution (0=uniform,1=kernel)',iunit)
+
+! end subroutine write_options_ptmass_heating
+
+
+!-----------------------------------------------------------------------
+!+
+!  read options from input file (not used at the moment)
+!+
+!-----------------------------------------------------------------------
+! subroutine read_options_ptmass_heating(name,valstring,imatch,igotall,ierr)
+!  use io,  only:fatal
+!  character(len=*), intent(in)  :: name,valstring
+!  logical, intent(out)          :: imatch,igotall
+!  integer, intent(out)          :: ierr
+!  integer, save                 :: ngot = 0
+!  integer                       :: ni
+!  character(len=30), parameter  :: label = 'read_options_ptmass_heating'
+
+!  imatch  = .true.
+!  igotall = .false.
+!  select case(trim(name))
+!  case('isink_heating')
+!     read(valstring,*,iostat=ierr) isink_heating
+!     ngot = ngot + 1
+!     if (isink_heating < 0 .or. isink_heating > 1) call fatal(label,'invalid setting for isink_heating ([0,1])')
+!  case default
+!     imatch = .false.
+!  end select
+!  ni = 1
+!  igotall = (ngot >= ni)
+
+! end subroutine read_options_ptmass_heating
+
+
 end module ptmass_heating

src/main/utils_filenames.f90

@@ -591,11 +591,7 @@ subroutine get_column_labels(line,nlabels,labels,method,ndesired,csv)
     !
     istyle = 1
     i1 = max(index(line,'[')+1,i1)    ! strip leading square bracket
-    ! try with different number of spaces between brackets (if labels not found)
-    over_spaces1: do i=4,0,-1
-       call split(line(i1:),']'//spaces(1:i)//'[',labels,nlabels)
-       if (nlabels > 1) exit over_spaces1
-    enddo over_spaces1
+    call split(nospaces(line(i1:)),'][',labels,nlabels)
  elseif (index(line,',') > 1 .or. is_csv) then
     !
     ! format style 2: mylabel1,mylabel2,mylabel3

src/setup/readwrite_mesa.f90

@@ -198,39 +198,43 @@ subroutine write_mesa(outputpath,m,pres,temp,r,rho,ene,Xfrac,Yfrac,csound,mu)
  real, intent(in)                :: m(:),rho(:),pres(:),r(:),ene(:),temp(:)
  real, intent(in), optional      :: Xfrac(:),Yfrac(:),csound(:),mu(:)
  character(len=120), intent(in)  :: outputpath
- character(len=200)              :: headers
- integer                         :: i,noptionalcols,j,iu
+ character(len=200)              :: headers(100)
+ integer                         :: i,ncols,noptionalcols,j,iu
  real, allocatable               :: optionalcols(:,:)
- character(len=*), parameter     :: fmtstring = "(5(es13.6,2x),es13.6)"
+ character(len=*), parameter     :: fmtstring = "(5(es24.16e3,2x),es24.16e3)"
 
- headers = '[    Mass   ]  [  Pressure ]  [Temperature]  [   Radius  ]  [  Density  ]  [   E_int   ]'
+ ncols = 6
+ headers(1:ncols) = (/'       Mass','   Pressure','Temperature','     Radius','    Density','       Eint'/)
 
  ! Add optional columns
- allocate(optionalcols(size(r),10))
  noptionalcols = 0
+ allocate(optionalcols(size(r),10))
  if (present(Xfrac)) then
     noptionalcols = noptionalcols + 1
-    headers = trim(headers) // '  [   Xfrac   ]'
+    headers(noptionalcols+ncols) = '      Xfrac'
     optionalcols(:,noptionalcols) = Xfrac
  endif
  if (present(Yfrac)) then
     noptionalcols = noptionalcols + 1
-    headers = trim(headers) // '  [   Yfrac   ]'
+    headers(noptionalcols+ncols) = '      Yfrac'
     optionalcols(:,noptionalcols) = Yfrac
  endif
  if (present(mu)) then
     noptionalcols = noptionalcols + 1
-    headers = trim(headers) // '  [    mu     ]'
+    headers(noptionalcols+ncols) = '        mu'
     optionalcols(:,noptionalcols) = mu
  endif
  if (present(csound)) then
     noptionalcols = noptionalcols + 1
-    headers = trim(headers) // '  [Sound speed]'
+    headers(noptionalcols+ncols) = 'Sound speed'
     optionalcols(:,noptionalcols) = csound
  endif
 
  open(newunit=iu, file = outputpath, status = 'replace')
- write(iu,'(a)') headers
+ do i = 1,noptionalcols+ncols-1
+    write(iu,'(a24,2x)',advance="no") trim(headers(i))
+ enddo
+ write(iu,'(a24)') trim(headers(noptionalcols+ncols))
 
  do i=1,size(r)
     if (noptionalcols <= 0) then
@@ -239,9 +243,9 @@ subroutine write_mesa(outputpath,m,pres,temp,r,rho,ene,Xfrac,Yfrac,csound,mu)
        write(iu,fmtstring,advance="no") m(i),pres(i),temp(i),r(i),rho(i),ene(i)
        do j=1,noptionalcols
           if (j==noptionalcols) then
-             write(iu,'(2x,es13.6)') optionalcols(i,j)
+             write(iu,'(2x,es24.16e3)') optionalcols(i,j)
           else
-             write(iu,'(2x,es13.6)',advance="no") optionalcols(i,j)
+             write(iu,'(2x,es24.16e3)',advance="no") optionalcols(i,j)
           endif
        enddo
     endif