! $Id$
!
! This module takes care of everything related to internal particle number.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
!
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! MPVAR CONTRIBUTION 1
! CPARAM logical, parameter :: lparticles_number=.true.
!
!***************************************************************
module Particles_number
!
use Cparam
use Cdata
use General, only: keep_compiler_quiet
use General, only: itoa
use Messages
use Particles_cdata
use Particles_sub
!
implicit none
!
include 'particles_number.h'
!
real :: np_swarm0=1.0, rhop_swarm0=1.0, np_swarm0_luck=1.0
real :: vthresh_coagulation=0.0, deltavp22_floor=0.0
real :: tstart_fragmentation_par=0.0, cdtpf=0.2
real :: birthring_inner=0.0, birthring_outer=huge1
real :: depletion_rate=0.0
logical :: lfragmentation_par=.false.,lbirthring_depletion=.false.
character (len=labellen), dimension(ninit) :: initnpswarm='nothing'
character (len=intlen) :: sdust
!
integer :: idiag_npswarmm=0, idiag_dvp22mwnp=0, idiag_dvp22mwnp2=0, k
integer :: idiag_dtfragp=0, idiag_npsm=0
integer, parameter :: mmom=24
integer, dimension(0:mmom) :: idiag_admom=0
logical :: llog10_for_admom_above10=.true.
!
namelist /particles_number_init_pars/ &
initnpswarm, np_swarm0, rhop_swarm0, vthresh_coagulation, &
deltavp22_floor, lfragmentation_par, tstart_fragmentation_par, cdtpf, &
np_swarm0_luck
!
namelist /particles_number_run_pars/ &
initnpswarm, vthresh_coagulation, deltavp22_floor, &
lfragmentation_par, tstart_fragmentation_par, cdtpf, &
lbirthring_depletion, birthring_inner, birthring_outer, &
depletion_rate, llog10_for_admom_above10
!
contains
!***********************************************************************
subroutine register_particles_number()
!
! Set up indices for access to the fp and dfp arrays.
!
! 24-nov-05/anders: adapted
!
if (lroot) call svn_id( &
"$Id$")
!
! Index for particle internal number.
!
call append_npvar('inpswarm',inpswarm)
!
endsubroutine register_particles_number
!***********************************************************************
subroutine initialize_particles_number(f)
!
! Perform any post-parameter-read initialization i.e. calculate derived
! parameters.
!
! 24-nov-05/anders: adapted
!
real, dimension (mx,my,mz,mfarray) :: f
!
if (lfragmentation_par) then
if (.not. allocated(kneighbour)) allocate(kneighbour(mpar_loc))
lshepherd_neighbour=.true.
endif
!
if (lfragmentation_par .and..not.(lcartesian_coords.and.(all(lequidist)))) &
call fatal_error( 'initialize_particles_number', &
'fragmentation only '//'implemented for Cartesian equidistant grids.')
if (lbirthring_depletion .and..not.(lcylindrical_coords.or.lspherical_coords)) &
call fatal_error( 'initialize_particles_number', &
'birthring depletion only '//'implemented for cyl, sph polar grids.')
!
call keep_compiler_quiet(f)
!
endsubroutine initialize_particles_number
!***********************************************************************
subroutine init_particles_number(f,fp)
!
! Initial internal particle number.
!
! 24-nov-05/anders: adapted
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mpar_loc,mparray) :: fp
!
call set_particle_number(f,fp,1,npar_loc,init=.true.)
!
call keep_compiler_quiet(f)
!
endsubroutine init_particles_number
!***********************************************************************
subroutine set_particle_number(f,fp,npar_low,npar_high,init)
!
use General, only: loptest
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mpar_loc,mparray) :: fp
integer :: npar_low, npar_high
logical, optional :: init
logical :: initial
integer :: j
!
initial=lroot.and.loptest(init)
!
do j=1,ninit
select case (initnpswarm(j))
case ('nothing')
if (initial.and.j==1) print*, 'init_particles_number: nothing'
case ('constant')
if (initial) then
print*, 'init_particles_number: constant internal number'
print*, 'init_particles_number: np_swarm0=', np_swarm0
endif
fp(npar_low:npar_high,inpswarm)=np_swarm0
case ('constant_rhop_swarm')
if (initial) then
print*, 'init_particles_number: constant mass density per particle'
print*, 'init_particles_number: rhop_swarm0=', rhop_swarm0
endif
fp(npar_low:npar_high,inpswarm)=rhop_swarm0/ &
(four_pi_rhopmat_over_three*fp(npar_low:npar_high,iap)**3)
case ('constant_rhop')
if (initial) then
print*, 'init_particles_number: constant mass density'
print*, 'init_particles_number: rhop_swarm0=', rhop_swarm0
endif
fp(npar_low:npar_high,inpswarm)=rhop_swarm0/(float(npar)/nwgrid)/ &
(four_pi_rhopmat_over_three*fp(npar_low:npar_high,iap)**3)
case ('constant_luck')
if (initial) then
print*, 'init_particles_number: constant internal number, luck'
print*, 'init_particles_number: np_swarm0=', np_swarm0
print*, 'init_particles_number: np_swarm0_luck=', np_swarm0_luck
endif
do k = npar_low,npar_high
if (ipar(k) == 1) then
fp(k,inpswarm) = np_swarm0_luck
else
fp(k,inpswarm) = np_swarm0
endif
enddo
case ('constant_luck_proc')
if (initial) then
print*, 'init_particles_number: constant internal number, luck, proc'
print*, 'init_particles_number: np_swarm0=', np_swarm0
print*, 'init_particles_number: np_swarm0_luck=', np_swarm0_luck
endif
do k = npar_low,npar_high
fp(k,inpswarm) = np_swarm0
enddo
fp(npar_low,inpswarm) = np_swarm0_luck
endselect
enddo
!
call keep_compiler_quiet(f)
!
endsubroutine set_particle_number
!***********************************************************************
subroutine pencil_criteria_par_number()
!
! All pencils that the Particles_number module depends on are specified here.
!
! 21-nov-06/anders: coded
!
if (lfragmentation_par) then
lpenc_requested(i_rho1)=.true.
lpenc_requested(i_cc1)=.true.
endif
!
endsubroutine pencil_criteria_par_number
!***********************************************************************
subroutine dnpswarm_dt_pencil(f,df,fp,dfp,p,ineargrid)
!
! Evolution of internal particle number.
!
! 24-oct-05/anders: coded
!
use Diagnostics
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
real, dimension (mpar_loc,mparray) :: fp
real, dimension (mpar_loc,mpvar) :: dfp
type (pencil_case) :: p
integer, dimension (mpar_loc,3) :: ineargrid
!
real, dimension (nx) :: dt1_fragmentation
real :: deltavp, sigma_jk, cdot, deltavp_sum, npswarm_sum, np2swarm_sum
integer :: j, k, l, ncoll=-1
logical :: lheader, lfirstcall=.true.
!
intent (in) :: f, fp
intent (out) :: dfp
!
! Print out header information in first time step.
!
lheader=lfirstcall .and. lroot
!
! Identify module and boundary conditions.
!
if (lheader) print*,'dnpswarm_dt_pencil: Calculate dnpswarm_dt'
!
! Collisional fragmentation inside each superparticle.
!
if (lfragmentation_par .and. t>=tstart_fragmentation_par) then
!
if (lfirst.and.ldt) dt1_fragmentation=0.0
!
if (npar_imn(imn)/=0) then
do l=l1,l2
! Need to count collisions for diagnostics.
if (ldiagnos) then
deltavp_sum=0.0; npswarm_sum=0; np2swarm_sum=0; ncoll=0
endif
! Get index number of shepherd particle at grid point.
k=kshepherd(l-nghost)
! Continue only if there is actually a shepherd particle.
if (k>0) then
if (ip<=6.and.lroot) then
print*, 'dnpswarm_dt: l, m, n=', l, m, n
print*, 'dnpswarm_dt: kshepherd, np(l,m,n)=', k, f(l,m,n,inp)
print*, 'dnpswarm_dt: x(l), y(m), z(n) =', x(l), y(m), z(n)
endif
! Only continue of the shepherd particle has a neighbour.
do while (k/=0)
j=k
! Consider neighbours one at a time.
do while (kneighbour(j)/=0)
j=kneighbour(j)
if (ip<=6.and.lroot) then
print*, &
'dnpswarm_dt: collisions between particle ', k, 'and', j
endif
! Collision speed.
deltavp=sqrt( &
(fp(k,ivpx)-fp(j,ivpx))**2 + &
(fp(k,ivpy)-fp(j,ivpy))**2 + &
(fp(k,ivpz)-fp(j,ivpz))**2 )
if (deltavp22_floor/=0.0) &
deltavp=sqrt(deltavp**2+deltavp22_floor**2)
! Collision cross section.
sigma_jk=pi*(fp(j,iap)+fp(k,iap))**2
! Collision rate between two superparticles.
cdot=sigma_jk*fp(j,inpswarm)*fp(k,inpswarm)*deltavp
! Either coagulation... [warning: this coagulation scheme is a bit cheaty]
if (deltavp<=vthresh_coagulation) then
dfp(j,inpswarm) = dfp(j,inpswarm) - 0.5*cdot
dfp(k,inpswarm) = dfp(k,inpswarm) - 0.5*cdot
if (fp(k,inpswarm)/=0.0) &
dfp(k,iap) = dfp(k,iap) + &
1/3.*(0.5*cdot)*fp(k,iap)/fp(k,inpswarm)
if (fp(j,inpswarm)/=0.0) &
dfp(j,iap) = dfp(j,iap) + &
1/3.*(0.5*cdot)*fp(j,iap)/fp(j,inpswarm)
! ...or fragmentation.
else
dfp(j,inpswarm) = dfp(j,inpswarm) - cdot
dfp(k,inpswarm) = dfp(k,inpswarm) - cdot
if (lpscalar_nolog) then
df(l,m,n,icc) = df(l,m,n,icc) + &
p%rho1(l-nghost)*4/3.*pi*rhopmat* &
(fp(j,iap)**3+fp(k,iap)**3)*cdot
else
df(l,m,n,ilncc) = df(l,m,n,ilncc) + &
p%cc1(l-nghost)*p%rho1(l-nghost)*4/3.*pi*rhopmat* &
(fp(j,iap)**3+fp(k,iap)**3)*cdot
endif
endif ! fragmentation or coagulation
! Time-step contribution
! if (lfirst.and.ldt) then
! dt1_fragmentation(l-nghost) = dt1_fragmentation(l-nghost) +&
! p%cc1(l-nghost)*p%rho1(l-nghost)*4/3.*pi*rhopmat* &
! (fp(j,iap)**3+fp(k,iap)**3)*cdot
! endif
! Need to count collisions for diagnostics.
if (ldiagnos) then
deltavp_sum =deltavp_sum +deltavp
npswarm_sum =npswarm_sum +fp(j,inpswarm)+fp(k,inpswarm)
np2swarm_sum=np2swarm_sum+fp(j,inpswarm)*fp(k,inpswarm)
ncoll=ncoll+1
endif
enddo
! Subgrid model of collisions within a superparticle.
if (deltavp22_floor/=0.0) then
if (ip<=6.and.lroot) then
print*, 'dnpswarm_dt: collisions within particle ', k
endif
deltavp=deltavp22_floor
sigma_jk=pi*(fp(k,iap)+fp(k,iap))**2
cdot = sigma_jk*fp(k,inpswarm)*fp(k,inpswarm)*deltavp
dfp(k,inpswarm) = dfp(k,inpswarm) - cdot
if (lpscalar_nolog) then
df(l,m,n,icc) = df(l,m,n,icc) + &
p%rho1(l-nghost)*4/3.*pi*rhopmat*fp(k,iap)**3*cdot
else
df(l,m,n,ilncc) = df(l,m,n,ilncc) + &
p%cc1(l-nghost)*p%rho1(l-nghost)*4/3.*pi* &
rhopmat*fp(k,iap)**3*cdot
endif
! Need to count collisions for diagnostics.
if (ldiagnos) then
deltavp_sum =deltavp_sum +deltavp
npswarm_sum =npswarm_sum +fp(k,inpswarm)
np2swarm_sum=np2swarm_sum+fp(k,inpswarm)**2
ncoll=ncoll+1
endif
! Time-step contribution
! if (lfirst.and.ldt) then
! dt1_fragmentation(l-nghost) = dt1_fragmentation(l-nghost) +&
! p%cc1(l-nghost)*p%rho1(l-nghost)*4/3.*pi*rhopmat*fp(k,iap)**3*cdot
! endif
endif ! subgrid model
!
k=kneighbour(k)
enddo
! "if (k>0) then"
endif
!
if (ldiagnos) then
! Average collision speed per particle
if (idiag_dvp22mwnp/=0 .and. ncoll>=1) &
call sum_weighted_name((/ deltavp_sum/ncoll /), &
(/ npswarm_sum /),idiag_dvp22mwnp)
! Average collision speed per collision
if (idiag_dvp22mwnp2/=0 .and. ncoll>=1) &
call sum_weighted_name((/ deltavp_sum/ncoll /), &
(/ np2swarm_sum /),idiag_dvp22mwnp2)
endif
enddo ! l1,l2
!
if (lfirst.and.ldt) then
dt1_fragmentation=dt1_fragmentation/cdtpf
dt1_max=max(dt1_max,dt1_fragmentation)
if (ldiagnos.and.idiag_dtfragp/=0) &
call max_mn_name(dt1_fragmentation,idiag_dtfragp,l_dt=.true.)
endif
!
endif ! npar_imn/=0
endif ! lfragmentation_par
!
if (lbirthring_depletion) then
! deplete particles within a certain radial range
! (as in collisional grain destruction in debris disks)
where ((fp(1:npar_loc,ixp)>birthring_inner) .and. &
(fp(1:npar_loc,ixp)<birthring_outer)) &
dfp(1:npar_loc,inpswarm) = -depletion_rate
endif
!
lfirstcall=.false.
!
call keep_compiler_quiet(ineargrid)
!
endsubroutine dnpswarm_dt_pencil
!***********************************************************************
subroutine dnpswarm_dt(f,df,fp,dfp,ineargrid)
!
! Evolution of internal particle number.
!
! 21-nov-06/anders: coded
!
use Sub
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
real, dimension (mpar_loc,mparray) :: fp
real, dimension (mpar_loc,mpvar) :: dfp
integer, dimension (mpar_loc,3) :: ineargrid
!
! Diagnostic output
!
if (ldiagnos) then
if (idiag_npswarmm/=0) &
call sum_par_name(fp(1:npar_loc,inpswarm),idiag_npswarmm)
if (idiag_npsm/=0) &
call integrate_par_name(fp(1:npar_loc,inpswarm)/nwgrid,idiag_npsm)
endif
if (ldiagnos) then
do k=0,mmom
if(idiag_admom(k)/=0) then
if (llog10_for_admom_above10 .and. k>=12) then
call integrate_par_name(fp(1:npar_loc,inpswarm)*fp(1:npar_loc,iap)**k/nwgrid,idiag_admom(k),llog10=.true.)
! call integrate_par_name(fp(1:npar_loc,inpswarm)*fp(1:npar_loc,iap)**k/nwgrid,idiag_admom(k))
else
call integrate_par_name(fp(1:npar_loc,inpswarm)*fp(1:npar_loc,iap)**k/nwgrid,idiag_admom(k))
endif
endif
enddo
endif
!
call keep_compiler_quiet(f,df)
call keep_compiler_quiet(dfp)
call keep_compiler_quiet(ineargrid)
!
endsubroutine dnpswarm_dt
!***********************************************************************
subroutine read_particles_num_init_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=particles_number_init_pars, IOSTAT=iostat)
!
endsubroutine read_particles_num_init_pars
!***********************************************************************
subroutine write_particles_num_init_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=particles_number_init_pars)
!
endsubroutine write_particles_num_init_pars
!***********************************************************************
subroutine read_particles_num_run_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=particles_number_run_pars, IOSTAT=iostat)
!
endsubroutine read_particles_num_run_pars
!***********************************************************************
subroutine write_particles_num_run_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=particles_number_run_pars)
!
endsubroutine write_particles_num_run_pars
!***********************************************************************
subroutine rprint_particles_number(lreset,lwrite)
!
! Read and register print parameters relevant for internal particle number.
!
! 24-aug-05/anders: adapted
!
use Diagnostics
!
logical :: lreset
logical, optional :: lwrite
!
integer :: iname
!
! Reset everything in case of reset.
!
if (lreset) then
idiag_npswarmm=0; idiag_dvp22mwnp=0; idiag_dvp22mwnp2=0
idiag_dtfragp=0; idiag_npsm=0; idiag_admom=0
endif
!
! Run through all possible names that may be listed in print.in.
!
if (lroot.and.ip<14) &
print*, 'rprint_particles_number: run through parse list'
do iname=1,nname
call parse_name(iname,cname(iname),cform(iname), &
'npswarmm',idiag_npswarmm)
call parse_name(iname,cname(iname),cform(iname), &
'npsm',idiag_npsm)
do k=0,mmom
sdust=itoa(k)
call parse_name(iname,cname(iname),cform(iname),'admom',idiag_admom(k))
call parse_name(iname,cname(iname),cform(iname),'admom'//trim(sdust),idiag_admom(k))
enddo
call parse_name(iname,cname(iname),cform(iname), &
'dvp22mwnp',idiag_dvp22mwnp)
call parse_name(iname,cname(iname),cform(iname), &
'dvp22mwnp2',idiag_dvp22mwnp2)
call parse_name(iname,cname(iname),cform(iname),'dtfragp',idiag_dtfragp)
enddo
!
endsubroutine rprint_particles_number
!***********************************************************************
endmodule Particles_number