Skip to content
Open
Show file tree
Hide file tree
Changes from 6 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/main/dust_formation.f90
Original file line number Diff line number Diff line change
Expand Up @@ -25,11 +25,12 @@ module dust_formation
!

use part, only:idJstar,idK0,idK1,idK2,idK3,idmu,idgamma,idsat,idkappa
use physcon, only:kboltz,pi,atomic_mass_unit,mass_proton_cgs,patm
use physcon, only:kboltz,pi,atomic_mass_unit,mass_proton_cgs,patm,rho_Cdust,a0,mc
use dim, only:nElements

implicit none
integer, public :: idust_opacity = 0
! real, public :: rho_Cdust = 2.62, a0 = 1.28e-8, mc = 12.*atomic_mass_unit !cgs

public :: set_abundances,evolve_dust,evolve_chem,calc_kappa_dust,&
calc_kappa_bowen,chemical_equilibrium_light,psat_C,calc_nucleation,&
Expand Down Expand Up @@ -233,7 +234,6 @@ pure real function calc_kappa_dust(K3, Tdust, rho_cgs)
real, intent(in) :: K3, Tdust, rho_cgs

real :: kappa_cgs, fac
real, parameter :: rho_Cdust = 2.62, mc = 12.*atomic_mass_unit
real, parameter :: Qplanck_abs = 1.6846124267740528e+04
real, parameter :: Qross_ext = 9473.2722815583311

Expand Down
105 changes: 76 additions & 29 deletions src/main/growth.f90
Original file line number Diff line number Diff line change
Expand Up @@ -220,23 +220,27 @@ end subroutine print_growthinfo
!-----------------------------------------------------------------------
subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,filfac,dmdt,Vrel_disp)
use part, only:rhoh,idust,igas,iamtype,iphase,isdead_or_accreted,&
massoftype,Omega_k,dustfrac,tstop,deltav
massoftype,Omega_k,dustfrac,tstop,deltav,&
nucleation
use options, only:use_dustfrac,use_porosity
use physcon, only:fourpi
use eos, only:ieos,get_spsound
use timestep, only:dtmax
real, intent(inout) :: dustprop(:,:)
real, intent(in) :: Vrel_disp(:)
real, intent(inout) :: dustgasprop(:,:)
real, intent(in) :: xyzh(:,:)
real, intent(in) :: filfac(:)
real, intent(inout) :: VrelVf(:,:),vxyzu(:,:)
real, intent(out) :: dmdt(:)
integer, intent(in) :: npart
real :: rhog,rhod,vrel,rho,sdust
real :: massgrain,rhograin,filfaci
real :: dtarb,frac_masschange,att_factor
integer :: i,iam
use dim, only:do_nucleation
real, intent(in) :: dustprop(:,:)
real, intent(in) :: Vrel_disp(:)
real, intent(inout) :: dustgasprop(:,:)
real, intent(in) :: xyzh(:,:)
real, intent(in) :: filfac(:)
real, intent(inout) :: VrelVf(:,:),vxyzu(:,:)
real, intent(out) :: dmdt(:)
integer, intent(in) :: npart
!
real :: rhog,rhod,vrel,rho,sdust
real :: massgrain,rhograin,filfaci
real :: dtarb,frac_masschange,att_factor,tdyn
integer :: i,iam


vrel = 0.
rhod = 0.
Expand All @@ -249,23 +253,30 @@ subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,filfac,d
!$omp shared(npart,iphase,ieos,massoftype,use_dustfrac,dustfrac,use_porosity,grainsizemin) &
!$omp shared(ifrag,ieros,utime,umass,dsize,cohacc,dtmax) &
!$omp shared(xyzh,vxyzu,dustprop,dustgasprop,dmdt,filfac,VrelVf,tstop,deltav,Vrel_disp) &
!$omp private(i,iam,rho,rhog,rhod,vrel,sdust) &
!$omp shared(do_nucleation,nucleation) &
!$omp private(i,iam,rho,rhog,rhod,vrel,sdust,tdyn) &
!$omp private(massgrain,rhograin,dtarb,frac_masschange,att_factor) &
!$omp firstprivate(filfaci)
do i=1,npart
if (.not.isdead_or_accreted(xyzh(4,i))) then
iam = iamtype(iphase(i))

call get_dynamical_timescale(i,tdyn)

if (iam == idust .or. (iam == igas .and. use_dustfrac)) then

if (use_dustfrac .and. iam == igas) then
!- no need for interpolations
rho = rhoh(xyzh(4,i),massoftype(igas))
rhog = rho*(1-dustfrac(1,i))
rhod = rho*dustfrac(1,i)
dustgasprop(1,i) = get_spsound(ieos,xyzh(:,i),rhog,vxyzu(:,i))
if (do_nucleation) then
dustgasprop(1,i) = get_spsound(ieos,xyzh(:,i),rhog,vxyzu(:,i),nucleation(7,i))
else
dustgasprop(1,i) = get_spsound(ieos,xyzh(:,i),rhog,vxyzu(:,i))
endif
dustgasprop(2,i) = rhog
dustgasprop(3,i) = tstop(1,i) * Omega_k(i)
dustgasprop(3,i) = tstop(1,i) / tdyn
dustgasprop(4,i) = sqrt(deltav(1,1,i)**2 + deltav(2,1,i)**2 + deltav(3,1,i)**2)
else
rhod = rhoh(xyzh(4,i),massoftype(idust))
Expand Down Expand Up @@ -632,31 +643,47 @@ end subroutine check_dustprop
!-----------------------------------------------------------------------
subroutine set_dustprop(npart,xyzh,sizedistrib,pwl_sizedistrib,R_ref,H_R_ref,q_index)
use dust, only:grainsizecgs,graindenscgs
use part, only:iamtype,iphase,idust,igas,dustprop,filfac,probastick
use part, only:iamtype,iphase,idust,igas,dustprop,filfac,probastick,nucleation
use physcon, only:fourpi
use options, only:use_dustfrac
use dim, only:do_nucleation
use physcon, only:rho_Cdust,a0
integer, intent(in) :: npart
real, intent(in) :: xyzh(:,:)
logical, intent(in), optional :: sizedistrib
real, intent(in), optional :: pwl_sizedistrib,R_ref,H_R_ref,q_index
integer :: i,iam
real :: r,h
real :: r,h,grainmassmin


do i=1,npart
iam = iamtype(iphase(i))
if (iam == idust .or. (iam == igas .and. use_dustfrac)) then
dustprop(2,i) = graindenscgs / unit_density
if (sizedistrib) then
r = sqrt(xyzh(1,i)**2 + xyzh(2,i)**2)
h = H_R_ref * R_ref * au / udist * (r * udist / au / R_ref)**(1.5-q_index)
dustprop(1,i) = grainsizecgs/udist * (r * udist / au / R_ref)**pwl_sizedistrib * exp(-0.5*xyzh(3,i)**2/h**2)
if (dustprop(1,i) < 2.e-5/udist) then
dustprop(1,i) = 2.e-5/udist
endif
dustprop(1,i) = fourpi/3. * dustprop(2,i) * (dustprop(1,i))**3

if (do_nucleation) then !used by moddump_nucleationtodustgrowth
dustprop(2,i) = rho_Cdust / unit_density
grainmassmin = fourpi/3 * (1.e-8/udist)**3 * dustprop(2,i) ! arbitrary threshold
if (nucleation(2,i) < 1e-30) then !arbitrary threshold in number of dust grain
dustprop(1,i) = grainmassmin !tiny(dustprop(1,i))
else
dustprop(1,i) = a0*nucleation(3,i)/nucleation(2,i) / udist !avg grain size = a0 K1/K0
dustprop(1,i) = fourpi/3. * dustprop(2,i) * (dustprop(1,i))**3 ! convert into dust mass
if (dustprop(1,i)<grainmassmin) dustprop(1,i) = grainmassmin
endif
else
dustprop(1,i) = fourpi/3. * dustprop(2,i) * (grainsizecgs / udist)**3
endif
dustprop(2,i) = graindenscgs / unit_density
if (sizedistrib) then
r = sqrt(xyzh(1,i)**2 + xyzh(2,i)**2)
h = H_R_ref * R_ref * au / udist * (r * udist / au / R_ref)**(1.5-q_index)
dustprop(1,i) = grainsizecgs/udist * (r * udist / au / R_ref)**pwl_sizedistrib * exp(-0.5*xyzh(3,i)**2/h**2)
if (dustprop(1,i) < 2.e-5/udist) then
dustprop(1,i) = 2.e-5/udist
endif
dustprop(1,i) = fourpi/3. * dustprop(2,i) * (dustprop(1,i))**3
else
dustprop(1,i) = fourpi/3. * dustprop(2,i) * (grainsizecgs / udist)**3
endif
endif
else
dustprop(:,i) = 0.
endif
Expand Down Expand Up @@ -1069,4 +1096,24 @@ real function get_size(mass,dens,filfac)

end function get_size


!-----------------------------------------------------------------------
!+
! Compute a characteristic dynamical timescale. Used in growth.f90 to get the St.
!+
!-----------------------------------------------------------------------
subroutine get_dynamical_timescale(i,tdyn)
use part, only:xyzmh_ptmass,ivwind,iReff,Omega_k
use dim, only:do_nucleation
integer, intent(in) :: i
real, intent(inout) :: tdyn

if (do_nucleation) then
tdyn = 1/(xyzmh_ptmass(ivwind,1)/xyzmh_ptmass(iReff,1))
else
tdyn = 1/Omega_k(i)
endif

end subroutine get_dynamical_timescale
Comment on lines +1105 to +1117

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

high

In get_dynamical_timescale, the calculation 1/(xyzmh_ptmass(ivwind,1)/xyzmh_ptmass(iReff,1)) can be simplified to xyzmh_ptmass(iReff,1)/xyzmh_ptmass(ivwind,1) to avoid an extra division. Additionally, we should protect against out-of-bounds array access if nptmass < 1 (by importing and checking nptmass), and prevent division by zero if xyzmh_ptmass(ivwind,1) is zero. Finally, the intent of tdyn should be out rather than inout since it is only written to.

 subroutine get_dynamical_timescale(i,tdyn)
 use part,           only:xyzmh_ptmass,ivwind,iReff,Omega_k,nptmass
 use dim,            only:do_nucleation
 integer,   intent(in) :: i
 real,      intent(out) :: tdyn
 
 if (do_nucleation .and. nptmass >= 1) then
    if (xyzmh_ptmass(ivwind,1) > 0.0) then
       tdyn = xyzmh_ptmass(iReff,1)/xyzmh_ptmass(ivwind,1)
    else
       tdyn = 1/Omega_k(i)
    endif
 else
    tdyn = 1/Omega_k(i)
 endif

 end subroutine get_dynamical_timescale


end module growth
4 changes: 4 additions & 0 deletions src/main/physcon.f90
Original file line number Diff line number Diff line change
Expand Up @@ -101,4 +101,8 @@ module physcon
!
real(kind=8), parameter :: eVtoK = 1.1604519d4 !Degrees kelvin per eV K/eV

!
!--Dust nucleation
!
real(kind=8), parameter :: rho_Cdust = 2.62, a0 = 1.28e-8, mc = 12.*atomic_mass_unit
end module physcon
147 changes: 147 additions & 0 deletions src/utils/moddump_nucleationtogrowth.f90
Original file line number Diff line number Diff line change
@@ -0,0 +1,147 @@
!--------------------------------------------------------------------------!
! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. !
! Copyright (c) 2007-2026 The Authors (see AUTHORS) !
! See LICENCE file for usage and distribution conditions !
! http://phantomsph.github.io/ !
!--------------------------------------------------------------------------!
module moddump
!
! default moddump routine: does not make any modifications
!
! :References: None
!
! :Owner: Daniel Price
!
! :Runtime parameters: None
!
! :Dependencies: None
!
implicit none
character(len=*), parameter, public :: moddump_flags = ''

contains

subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu)
use dim, only:use_dust,maxdusttypes,maxdustlarge,maxdustsmall,use_dustgrowth,&
update_max_sizes,isothermal
use part, only:igas,idust,ndusttypes,ndustsmall,ndustlarge,&
grainsize,graindens,dustfrac,&
nucleation,dustprop
use dust_formation, only:mass_per_H
use set_dust, only:set_dustfrac,set_dustbinfrac
use options, only:use_dustfrac,use_porosity
use growth, only:set_dustprop,convert_to_twofluid,iporosity,init_growth,ifrag,gsizemincgs,&
ivrelkin
use prompting, only:prompt
use dust, only:grainsizecgs,graindenscgs
use table_utils, only:logspace
use units, only:umass,udist
use physcon, only:rho_Cdust,a0,mc
use cooling_solver, only:excitation_HI
integer, intent(inout) :: npart
integer, intent(inout) :: npartoftype(:)
real, intent(inout) :: massoftype(:)
real, intent(inout) :: xyzh(:,:),vxyzu(:,:)
integer :: i,ierr,j,itype,ipart,iloc,dust_method,np_ratio,np_gas,np_dust,maxdust
real :: dust_to_gas,smincgs,smaxcgs,sindex,dustbinfrac(maxdusttypes),udens

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

medium

Add a declaration for dust_to_gas_avg to support calculating the average dust-to-gas ratio across all gas particles.

 real    :: dust_to_gas,smincgs,smaxcgs,sindex,dustbinfrac(maxdusttypes),udens,dust_to_gas_avg

real :: inradius,outradius,pwl_sizedistrib,R_ref,H_R_ref,q_index
logical :: sizedistrib


! cooling values
excitation_HI = 1

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

high

Add a defensive check to ensure that the nucleation array is allocated before accessing it, preventing potential segmentation faults if this utility is run on a dump file without nucleation data.

 excitation_HI = 1

 if (.not. allocated(nucleation)) then
    print *, "Error: Nucleation array is not allocated. This dump does not contain nucleation data."
    stop
 endif


! dust default values
udens = umass/(udist**3)
dust_method = 1
np_ratio = 5
dust_to_gas = 0.01
ndusttypes = 1
smincgs = 1.e-5
smaxcgs = 1.
sindex = 3.5
dustbinfrac = 0.
pwl_sizedistrib = -2
R_ref = 100
H_R_ref = 0.0895
q_index = 0.25

sizedistrib = .false.

grainsizecgs = 1. !dummy
graindenscgs = 1.

gsizemincgs = 1.e-8
ivrelkin = 0 ! relative dust motions from 1-fluid dust velocities not yet available


! Interactive prompts
!call prompt('Which dust method do you want? (1=one fluid,2=two fluid)',dust_method,1,2)
dust_method = 1
if (dust_method==1) then
maxdust = maxdustsmall
elseif (dust_method==2) then
maxdust = maxdustlarge
call prompt('Enter ratio between number of gas particles and dust particles',np_ratio,1)
!--We do not care if modulo(npart,np_ratio) is stricly zero, since npart can
! be a weird value depdending on the simulation it comes from.
endif

! call prompt('Enter total dust to gas ratio',dust_to_gas,0.) computed from nuc.

!call prompt('How many grain sizes do you want?',ndusttypes,1,maxdust)
ndusttypes = 1
! if (ndusttypes > 1) then !uncomment if several dust populations
! !--grainsizes
! call prompt('Enter minimum grain size in cm',smincgs,0.)
! call prompt('Enter maximum grain size in cm',smaxcgs,0.)
! !--mass distribution
! call prompt('Enter power-law index, e.g. MRN',sindex)
! call set_dustbinfrac(smincgs/udist,smaxcgs/udist,sindex,dustbinfrac(1:ndusttypes),grainsize(1:ndusttypes))
! !--grain density
! call prompt('Enter grain density in g/cm^3',graindens(1),0.)
! graindens = graindens(1)/udens
! else
if (use_dustgrowth) then
ifrag = 1
iporosity = 0 !turn on use_porosity if you change that
grainsize(1) = grainsizecgs/udist
graindens(1) = graindenscgs/udens
endif

np_gas = npartoftype(igas)

if (dust_method == 1) then

use_dustfrac = .true.
ndustsmall = ndusttypes

do i=1,np_gas
if (ndusttypes > 1) then
dustfrac(1:ndusttypes,i) = dust_to_gas*dustbinfrac(1:ndusttypes)
else
if (nucleation(5,i)<tiny(nucleation(5,i))) then
dust_to_gas = tiny(dust_to_gas)
else
dust_to_gas = nucleation(5,i) * mc / mass_per_H
endif
call set_dustfrac(dust_to_gas,dustfrac(:,i))
endif
enddo

massoftype(igas) = massoftype(igas)*(1. + dust_to_gas)
Comment on lines +113 to +126

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

critical

In the loop over gas particles, dust_to_gas is a scalar variable that gets overwritten for each particle. Consequently, after the loop, massoftype(igas) is scaled using only the dust_to_gas ratio of the very last particle (np_gas), which is incorrect and can lead to wrong total gas mass scaling (especially if the last particle has zero dust). We should compute the average dust_to_gas ratio over all gas particles and use that to scale massoftype(igas) when ndusttypes == 1.

  dust_to_gas_avg = 0.
  do i=1,np_gas
     if (ndusttypes > 1) then
        dustfrac(1:ndusttypes,i) = dust_to_gas*dustbinfrac(1:ndusttypes)
     else
        if (nucleation(5,i)<tiny(nucleation(5,i))) then
            dust_to_gas = tiny(dust_to_gas)
        else
            dust_to_gas = nucleation(5,i) * mc / mass_per_H
        endif
        dust_to_gas_avg = dust_to_gas_avg + dust_to_gas
        call set_dustfrac(dust_to_gas,dustfrac(:,i))
     endif
  enddo

  if (ndusttypes > 1) then
     massoftype(igas) = massoftype(igas)*(1. + dust_to_gas)
  else
     dust_to_gas_avg = dust_to_gas_avg / np_gas
     massoftype(igas) = massoftype(igas)*(1. + dust_to_gas_avg)
  endif

npart = np_gas

if (use_dustgrowth) then
call set_dustprop(npart,xyzh,sizedistrib,pwl_sizedistrib,R_ref,H_R_ref,q_index)
endif
endif





end subroutine modify_dump

end module moddump

Loading