(sink_radiation; nucleation) cleanup of unused code, use eos_vars to …

…get temperature instead of calling get_temperature; remove ifdefs in derivs (#55)

src/main/config.F90

@@ -283,18 +283,11 @@ module dim
  !number of elements considered in the nucleation chemical network
  integer, parameter :: nElements = 10
 #ifdef DUST_NUCLEATION
-#ifdef STAR
- logical :: star_radiation = .true.
-#else
- logical :: star_radiation = .false.
-#endif
  logical :: nucleation = .true.
- integer :: maxsp = maxp_hard
 #else
- logical :: star_radiation = .false.
  logical :: nucleation = .false.
- integer :: maxsp = 0
 #endif
+ integer :: maxp_nucleation = 0
 
 !--------------------
 ! MCFOST library
@@ -363,9 +356,8 @@ subroutine update_max_sizes(n,ntot)
  store_dust_temperature = .true.
 #endif
 
- if (store_dust_temperature) then
-    maxTdust = maxp
- endif
+ if (store_dust_temperature) maxTdust = maxp
+ if (do_nucleation) maxp_nucleation = maxp
 
 #ifdef NCELLSMAX
  ncellsmax       = NCELLSMAX

src/main/deriv.F90

@@ -38,14 +38,15 @@ module deriv
 subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
                   Bevol,dBevol,rad,drad,radprop,dustprop,ddustprop,&
                   dustevol,ddustevol,dustfrac,eos_vars,time,dt,dtnew,pxyzu,dens,metrics)
- use dim,            only:maxvxyzu,mhd,fast_divcurlB,gr
+ use dim,            only:maxvxyzu,maxp,mhd,fast_divcurlB,gr,periodic,&
+                          sink_radiation,use_dustgrowth
  use io,             only:iprint,fatal
  use linklist,       only:set_linklist
  use densityforce,   only:densityiterate
- use ptmass,         only:ipart_rhomax,ptmass_calc_enclosed_mass
+ use ptmass,         only:ipart_rhomax,ptmass_calc_enclosed_mass,ptmass_boundary_crossing
  use externalforces, only:externalforce
  use part,           only:dustgasprop,dvdx,Bxyz,set_boundaries_to_active,&
-                          nptmass,xyzmh_ptmass,sinks_have_heating
+                          nptmass,xyzmh_ptmass,sinks_have_heating,dust_temp,VrelVf
 #ifdef IND_TIMESTEPS
  use timestep_ind,   only:nbinmax
 #else
@@ -59,21 +60,11 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
  use photoevap,      only:find_ionfront,photo_ionize
  use part,           only:massoftype
 #endif
-#ifdef DUSTGROWTH
- use growth,         only:get_growth_rate
- use part,           only:VrelVf
-#endif
-#if defined(SINK_RADIATION) && !defined(ISOTHERMAL)
+ use growth,           only:get_growth_rate
  use ptmass_radiation, only:get_dust_temperature_from_ptmass
- use part,             only:dust_temp
-#endif
-#ifdef PERIODIC
- use ptmass,         only:ptmass_boundary_crossing
-#endif
- use part,           only:mhd,gradh,alphaind,igas
  use timing,         only:get_timings
  use forces,         only:force
- use part,           only:iradxi,ifluxx,ifluxy,ifluxz,ithick
+ use part,           only:mhd,gradh,alphaind,igas,iradxi,ifluxx,ifluxy,ifluxz,ithick
  use derivutils,     only:do_timing
  use cons2prim,      only:cons2primall,cons2prim_everything,prim2consall
  use metric_tools,   only:init_metric
@@ -131,9 +122,7 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
        call prim2consall(npart,xyzh,metrics,vxyzu,dens,pxyzu,use_dens=.false.)
     endif
 
-#ifdef PERIODIC
-    if (nptmass > 0) call ptmass_boundary_crossing(nptmass,xyzmh_ptmass)
-#endif
+    if (nptmass > 0 .and. periodic) call ptmass_boundary_crossing(nptmass,xyzmh_ptmass)
  endif
 
  call do_timing('link',tlast,tcpulast,start=.true.)
@@ -165,11 +154,11 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
     call do_timing('dens',tlast,tcpulast)
  endif
 
-#ifdef GR
- call cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars)
-#else
- call cons2prim_everything(npart,xyzh,vxyzu,dvdx,rad,eos_vars,radprop,Bevol,Bxyz,dustevol,dustfrac,alphaind)
-#endif
+ if (gr) then
+    call cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars)
+ else
+    call cons2prim_everything(npart,xyzh,vxyzu,dvdx,rad,eos_vars,radprop,Bevol,Bxyz,dustevol,dustfrac,alphaind)
+ endif
  call do_timing('cons2prim',tlast,tcpulast)
 !
 ! compute forces
@@ -186,15 +175,16 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
             ipart_rhomax,dt,stressmax,eos_vars,dens,metrics)
  call do_timing('force',tlast,tcpulast)
 
-#ifdef DUSTGROWTH
- ! compute growth rate of dust particles
- call get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,ddustprop(1,:))!--we only get ds/dt (i.e 1st dimension of ddustprop)
-#endif
+ if (use_dustgrowth) then ! compute growth rate of dust particles
+    call get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,ddustprop(1,:))!--we only get ds/dt (i.e 1st dimension of ddustprop)
+ endif
 
-#if defined(SINK_RADIATION) && !defined(ISOTHERMAL)
- !compute dust temperature
- call get_dust_temperature_from_ptmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,dust_temp)
-#endif
+ if (sink_radiation .and. maxvxyzu==maxp) then
+    !
+    ! compute dust temperature based on radiation from sink particles
+    !
+    call get_dust_temperature_from_ptmass(npart,xyzh,eos_vars,nptmass,xyzmh_ptmass,dust_temp)
+ endif
 !
 ! set new timestep from Courant/forces condition
 !
@@ -204,11 +194,8 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
  dtnew = min(dtforce,dtcourant,dtrad,dtmax)
 #endif
 
-
  call do_timing('total',t1,tcpu1,lunit=iprint)
 
- return
-
 end subroutine derivs
 
 !--------------------------------------

src/main/part.F90

@@ -30,7 +30,7 @@ module part
                maxgrav,ngradh,maxtypes,h2chemistry,gravity,maxp_dustfrac,&
                use_dust,use_dustgrowth,lightcurve,maxlum,nalpha,maxmhdni, &
                maxp_growth,maxdusttypes,maxdustsmall,maxdustlarge, &
-               maxphase,maxgradh,maxan,maxdustan,maxmhdan,maxneigh,maxprad,maxsp,&
+               maxphase,maxgradh,maxan,maxdustan,maxmhdan,maxneigh,maxprad,maxp_nucleation,&
                maxTdust,store_dust_temperature,use_krome,maxp_krome, &
                do_radiation,gr,maxgr,maxgran,n_nden_phantom,do_nucleation,inucleation
  use dtypekdtree, only:kdnode
@@ -487,7 +487,7 @@ subroutine allocate_part
  call allocate_array('ibelong', ibelong, maxp)
  call allocate_array('istsactive', istsactive, maxsts)
  call allocate_array('ibin_sts', ibin_sts, maxsts)
- call allocate_array('nucleation', nucleation, n_nucleation*inucleation, maxsp*inucleation)
+ call allocate_array('nucleation', nucleation, n_nucleation, maxp_nucleation*inucleation)
 #ifdef KROME
  call allocate_array('abundance', abundance, krome_nmols, maxp_krome)
 #else

src/main/ptmass.F90

@@ -52,9 +52,7 @@ module ptmass
  public :: update_ptmass
  public :: calculate_mdot
  public :: ptmass_calc_enclosed_mass
-#ifdef PERIODIC
  public :: ptmass_boundary_crossing
-#endif
 
  ! settings affecting routines in module (read from/written to input file)
  integer, public :: icreate_sinks = 0
@@ -423,10 +421,9 @@ end subroutine get_accel_sink_sink
 !  Update position of sink particles if they cross the periodic boundary
 !+
 !----------------------------------------------------------------
-#ifdef PERIODIC
 subroutine ptmass_boundary_crossing(nptmass,xyzmh_ptmass)
- use boundary, only:cross_boundary
- use mpidomain,only:isperiodic
+ use boundary,  only:cross_boundary
+ use mpidomain, only:isperiodic
  integer, intent(in)    :: nptmass
  real,    intent(inout) :: xyzmh_ptmass(:,:)
  integer                :: i,ncross
@@ -437,7 +434,7 @@ subroutine ptmass_boundary_crossing(nptmass,xyzmh_ptmass)
  enddo
 
 end subroutine ptmass_boundary_crossing
-#endif
+
 !----------------------------------------------------------------
 !+
 !  predictor step for the point masses

src/main/ptmass_radiation.F90

@@ -33,19 +33,13 @@ module ptmass_radiation
  real,    public  :: tdust_exp       = 0.5
  real,    public  :: alpha_rad       = 0.
 
- public :: get_rad_accel_from_ptmass,read_options_ptmass_radiation,write_options_ptmass_radiation
+ public :: get_rad_accel_from_ptmass
+ public :: read_options_ptmass_radiation,write_options_ptmass_radiation
  public :: get_dust_temperature_from_ptmass
  public :: init_radiation_ptmass
 
- integer, parameter :: N = 1024
- real, parameter :: theta = 0., phi = 0.
- real, parameter :: u(3) = (/ sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta) /)
-
  private
 
- real  :: Lstar_lsun      = 5000.
- real  :: Mstar_msun      = 1.
-
 contains
 !-----------------------------------------------------------------------
 !+
@@ -66,39 +60,26 @@ end subroutine init_radiation_ptmass
 !-----------------------------------------------------------------------
 subroutine get_rad_accel_from_ptmass(nptmass,npart,xyzh,xyzmh_ptmass,fext)
  use part,    only:ilum
- use units,   only:umass,unit_energ,utime
- use dim,     only:star_radiation
- use physcon, only:solarl,solarm
+ use units,   only:umass,unit_luminosity
  integer,  intent(in)    :: nptmass,npart
  real,     intent(in)    :: xyzh(:,:)
  real,     intent(in)    :: xyzmh_ptmass(:,:)
  real,     intent(inout) :: fext(:,:)
  real                    :: xa,ya,za,Mstar_cgs,Lstar_cgs
  integer                 :: j
 
- if (star_radiation) then
-    Lstar_cgs  = Lstar_lsun*solarl
-    Mstar_cgs  = Mstar_msun*solarm
+ do j=1,nptmass
+    if (xyzmh_ptmass(4,j) < 0.) cycle
+    Mstar_cgs  = xyzmh_ptmass(4,j)*umass
+    Lstar_cgs  = xyzmh_ptmass(ilum,j)*unit_luminosity
+    !compute radiative acceleration if sink particle is assigned a non-zero luminosity
     if (Lstar_cgs > 0.d0) then
-       xa = xyzmh_ptmass(1,1)
-       ya = xyzmh_ptmass(2,1)
-       za = xyzmh_ptmass(3,1)
+       xa = xyzmh_ptmass(1,j)
+       ya = xyzmh_ptmass(2,j)
+       za = xyzmh_ptmass(3,j)
        call calc_rad_accel_from_ptmass(npart,xa,ya,za,Lstar_cgs,Mstar_cgs,xyzh,fext)
     endif
- else
-    do j=1,nptmass
-       if (xyzmh_ptmass(4,j) < 0.) cycle
-       Mstar_cgs  = xyzmh_ptmass(4,j)*umass
-       Lstar_cgs  = xyzmh_ptmass(ilum,j)*unit_energ/utime
-       !compute radiative acceleration if sink particle is assigned a non-zero luminosity
-       if (Lstar_cgs > 0.d0) then
-          xa = xyzmh_ptmass(1,j)
-          ya = xyzmh_ptmass(2,j)
-          za = xyzmh_ptmass(3,j)
-          call calc_rad_accel_from_ptmass(npart,xa,ya,za,Lstar_cgs,Mstar_cgs,xyzh,fext)
-       endif
-    enddo
- endif
+ enddo
 
 end subroutine get_rad_accel_from_ptmass
 
@@ -186,16 +167,14 @@ end subroutine get_radiative_acceleration_from_star
 !  through the gas, or by simpler approximations
 !+
 !-----------------------------------------------------------------------
-subroutine get_dust_temperature_from_ptmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,dust_temp)
- use part,    only:isdead_or_accreted,iLum,iTeff,iReff,rhoh,massoftype,igas,nucleation,idmu,idgamma
- use part,    only:eos_vars,itemp
- use eos,     only:ieos,get_temperature
- use dim,     only:do_nucleation
+subroutine get_dust_temperature_from_ptmass(npart,xyzh,eos_vars,nptmass,xyzmh_ptmass,dust_temp)
+ use part,    only:isdead_or_accreted,iLum,iTeff,iReff
+ use part,    only:itemp
  use io,      only:fatal
  integer,  intent(in)    :: nptmass,npart
- real,     intent(in)    :: xyzh(:,:),xyzmh_ptmass(:,:),vxyzu(:,:)
+ real,     intent(in)    :: xyzh(:,:),xyzmh_ptmass(:,:),eos_vars(:,:)
  real,     intent(out)   :: dust_temp(:)
- real                    :: r,L_star,T_star,R_star,xa,ya,za,pmassi,vxyzui(4)
+ real                    :: r,L_star,T_star,R_star,xa,ya,za
  integer                 :: i,j
 
  !
@@ -233,21 +212,12 @@ subroutine get_dust_temperature_from_ptmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmas
     call fatal('ptmass_radiation','Lucy approximation not implemented')
  case default
     ! sets Tdust = Tgas
-    pmassi = massoftype(igas)
-    !$omp parallel  do default(none) &
-    !$omp shared(npart,ieos,xyzh,vxyzu,eos_vars,pmassi,dust_temp) &
-    !$omp shared(nucleation,do_nucleation) &
-    !$omp private(i,vxyzui)
+    !$omp parallel do default(none) &
+    !$omp shared(npart,xyzh,eos_vars,dust_temp) &
+    !$omp private(i)
     do i=1,npart
        if (.not.isdead_or_accreted(xyzh(4,i))) then
-          vxyzui= vxyzu(:,i)
-          if (do_nucleation) then
-             dust_temp(i) = get_temperature(ieos,xyzh(:,i),rhoh(xyzh(4,i),pmassi),vxyzui,&
-                  gammai=nucleation(idgamma,i),mui=nucleation(idmu,i))
-
-          else
-             dust_temp(i) = eos_vars(itemp,i)
-          endif
+          dust_temp(i) = eos_vars(itemp,i)
        endif
     enddo
     !$omp end parallel do
@@ -261,25 +231,19 @@ end subroutine get_dust_temperature_from_ptmass
 !+
 !-----------------------------------------------------------------------
 subroutine write_options_ptmass_radiation(iunit)
- use infile_utils, only: write_inopt
- use dim,          only: star_radiation!,store_dust_temperature
+ use infile_utils, only:write_inopt
  integer, intent(in) :: iunit
 
  call write_inopt(isink_radiation,'isink_radiation','sink radiation pressure method (0=off,1=alpha,2=dust,3=alpha+dust)',iunit)
  if (isink_radiation == 1 .or. isink_radiation == 3) then
     call write_inopt(alpha_rad,'alpha_rad','fraction of the gravitational acceleration imparted to the gas',iunit)
  endif
  if (isink_radiation >= 2) then
-    if (star_radiation) then
-       call write_inopt(Lstar_lsun,'Lstar','Stellar luminosity (for radiation pressure, in Lsun)',iunit)
-       call write_inopt(Mstar_msun,'Mstar','Stellar mass (in Msun)',iunit)
-    endif
     call write_inopt(iget_tdust,'iget_tdust','dust temperature (0:Tdust=Tgas 1:T(r) 2:Lucy (devel))',iunit)
  endif
- if (iget_tdust == 1 ) then
+ if (iget_tdust == 1) then
     call write_inopt(tdust_exp,'tdust_exp','exponent of the dust temperature profile',iunit)
  endif
- !if (iget_tdust > 0) store_dust_temperature = .true.
 
 end subroutine write_options_ptmass_radiation
 
@@ -289,8 +253,7 @@ end subroutine write_options_ptmass_radiation
 !+
 !-----------------------------------------------------------------------
 subroutine read_options_ptmass_radiation(name,valstring,imatch,igotall,ierr)
- use io,      only:fatal
- use dim,     only:star_radiation
+ use io, only:fatal
  character(len=*), intent(in)  :: name,valstring
  logical, intent(out) :: imatch,igotall
  integer,intent(out) :: ierr
@@ -302,14 +265,6 @@ subroutine read_options_ptmass_radiation(name,valstring,imatch,igotall,ierr)
  imatch  = .true.
  igotall = .false.
  select case(trim(name))
- case('Lstar')
-    read(valstring,*,iostat=ierr) Lstar_lsun
-    ngot = ngot + 1
-    if (Lstar_lsun < 0.) call fatal(label,'invalid setting for Lstar (must be >= 0)')
- case('Mstar')
-    read(valstring,*,iostat=ierr) Mstar_msun
-    ngot = ngot + 1
-    if (Mstar_msun < 0.) call fatal(label,'invalid setting for Mstar (must be >= 0)')
  case('alpha_rad')
     read(valstring,*,iostat=ierr) alpha_rad
     ngot = ngot + 1
@@ -331,7 +286,6 @@ subroutine read_options_ptmass_radiation(name,valstring,imatch,igotall,ierr)
  ni = 1
  if (isink_radiation > 0) then
     ni = ni+1
-    if (star_radiation)  ni = ni+2
  endif
  igotall = (ngot >= ni)