! $Id$
!
! This modules solves for the gradient matrix of flow velocities.
! The gradient matrix Sigma obeys the equation:
! (d/dt) Sigma = (1/(taup))*(S - Sigma) - Sigma^2
! where A is the flow gradient matrix at the position of the
! particle; and taup is the Stokes time.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
!
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lparticles_caustics=.true.
!
! MAUX CONTRIBUTION 9
! COMMUNICATED AUXILIARIES 9
! MPVAR CONTRIBUTION 10
! MPAUX CONTRIBUTION 4
!
!***************************************************************
module Particles_caustics
!
use Cdata
use Messages
use Particles_cdata
use Particles_map
use Particles_mpicomm
use Particles_sub
!
implicit none
!
include 'particles_caustics.h'
!
integer :: dummy
real :: TrSigma_Cutoff=-1e10
!
namelist /particles_caustics_init_pars/ &
dummy
!
namelist /particles_caustics_run_pars/ &
TrSigma_Cutoff
!
! Diagnostic variables
!
integer :: idiag_TrSigmapm=0 ! DIAG_DOC: $\langle{\rm Tr}\left[\sigma\right]\rangle$
integer :: idiag_blowupm=0 ! DIAG_DOC: Mean no. of times $\sigma$ falls below cutoff
integer :: idiag_lnVpm=0 ! DIAG_DOC: Mean of (logarithm of) Volume around an inertial particle
integer, dimension(ndustrad) :: idiag_ncaustics=0
contains
!***********************************************************************
subroutine register_particles_caustics
!
! Set up indices for access to the fp and dfp arrays
!
! May-16/dhruba: coded
!
use FArrayManager, only: farray_register_auxiliary
!
if (lroot) call svn_id("particles_caustics")
!
! Indices for 9 elements of the sigmap matrix
!
call append_npvar('isigmap11',isigmap11)
call append_npvar('isigmap12',isigmap12)
call append_npvar('isigmap13',isigmap13)
call append_npvar('isigmap21',isigmap21)
call append_npvar('isigmap22',isigmap22)
call append_npvar('isigmap23',isigmap23)
call append_npvar('isigmap31',isigmap31)
call append_npvar('isigmap32',isigmap32)
call append_npvar('isigmap33',isigmap33)
!
! One extra variable that calculate the evolution
! of (logarithm of) volume around an inertial particle.
!
call append_npvar('lnVp',ilnVp)
!
! Now register the particle auxiliaries
!
call append_npaux('iPPp',iPPp)
call append_npaux('iQQp',iQQp)
call append_npaux('iRRp',iRRp)
call append_npaux('iblowup',iblowup)
!
! Set indices for velocity gradient matrix at grid points
!
call farray_register_auxiliary('guij',iguij,communicated=.true.,array=9)
igu11=iguij; igu12=iguij+1; igu13=iguij+2
igu21=iguij+3; igu22=iguij+4; igu23=iguij+5
igu31=iguij+6; igu32=iguij+7; igu33=iguij+8
!
endsubroutine register_particles_caustics
!***********************************************************************
subroutine initialize_particles_caustics(f)
!
! Perform any post-parameter-read initialization i.e. calculate derived
! parameters.
!
!
use General, only: keep_compiler_quiet
use FArrayManager
!
real, dimension (mx,my,mz,mfarray) :: f
!
! Stop if there is no velocity to calculate derivatives
!
if (.not. (lhydro.or.lhydro_kinematic)) &
call fatal_error('initialize_particles_caustics','you must select a hydro module')
!
call keep_compiler_quiet(f)
!
endsubroutine initialize_particles_caustics
!***********************************************************************
subroutine init_particles_caustics(f,fp,ineargrid)
!
use General, only: keep_compiler_quiet,random_number_wrapper
use Mpicomm, only: mpiallreduce_sum_int
use Hydro, only: calc_gradu
real, dimension (mx,my,mz,mfarray), intent (in) :: f
real, dimension (mpar_loc,mparray), intent (out) :: fp
integer, dimension (mpar_loc,3), intent (in) :: ineargrid
!
if (lroot) print*, 'init_particles_caustics: setting init. cond.'
call reinitialize_caustics(fp)
!
endsubroutine init_particles_caustics
!***********************************************************************
subroutine dcaustics_dt(f,df,fp,dfp,ineargrid)
!
use Particles_radius, only: get_stbin
use Diagnostics
use Particles_sub, only: sum_par_name
use Sub, only : linarray2matrix,Inv2_3X3mat,det3X3mat
!
real, dimension (mx,my,mz,mfarray), intent (in) :: f
real, dimension (mx,my,mz,mvar), intent (inout) :: df
real, dimension (mpar_loc,mparray), intent (inout) :: fp
real, dimension (mpar_loc,mpvar), intent (inout) :: dfp
integer, dimension (mpar_loc,3), intent (in) :: ineargrid
logical :: lheader, lfirstcall=.true.
real, dimension(3,3) :: Sigmap
real, dimension(9) :: Sigma_lin
real :: TrSigma,QSigma,detSigma
integer :: ip,iStbin,k
real :: prad
real,dimension(ndustrad) :: blowup_Stdep,np_Stdep
!
! Print out header information in first time step.
!
lheader=lfirstcall .and. lroot
if (lheader) then
print*,'dcaustics_dt: Calculate dcaustics_dt'
lfirstcall=.false.
endif
!
! Calculates the three invariants of the matrix sigma and stores them as auxiliary
! variable.
!
blowup_Stdep=0.
np_Stdep=0.
do ip=1,npar_loc
Sigma_lin=fp(ip,isigmap11:isigmap33)
call linarray2matrix(Sigma_lin,Sigmap)
TrSigma=Sigmap(1,1)+Sigmap(2,2)+Sigmap(3,3)
fp(ip,iPPp) = TrSigma
call Inv2_3X3mat(Sigmap,QSigma)
fp(ip,iQQp) = QSigma
call det3X3mat(Sigmap,detSigma)
fp(ip,iRRp) = detSigma
prad=fp(ip,iap)
call get_stbin(prad,iStbin)
blowup_Stdep(iStbin) = blowup_Stdep(iStbin)+fp(ip,iblowup)
np_Stdep(iStbin) = np_Stdep(iStbin)+1
enddo
!
if (ldiagnos) then
call sum_par_name(fp(1:npar_loc,iPPp),idiag_TrSigmapm)
call sum_par_name(fp(1:npar_loc,iblowup),idiag_blowupm)
call sum_par_name(fp(1:npar_loc,ilnVp),idiag_lnVpm)
do k=1,ndustrad
if (idiag_ncaustics(k)/=0) fname(idiag_ncaustics(k))=blowup_Stdep(k)/np_Stdep(k)
enddo
endif
!
if (lfirstcall) lfirstcall=.false.
!
endsubroutine dcaustics_dt
!***********************************************************************
subroutine dcaustics_dt_pencil(f,df,fp,dfp,p,ineargrid,k,taup1)
!
use Sub, only: linarray2matrix,matrix2linarray
!
real, dimension (mx,my,mz,mfarray),intent(in) :: f
real, dimension (mx,my,mz,mvar) :: df
type (pencil_case) :: p
real, dimension (mpar_loc,mparray), intent(in) :: fp
real, dimension (mpar_loc,mpvar) :: dfp
integer, dimension (mpar_loc,3) :: ineargrid
integer :: k
real :: taup1
intent (inout) :: df, dfp,ineargrid
intent (in) :: k,taup1
real, dimension(9) :: Sij_lin,Sigma_lin,dSigmap_lin
real, dimension(3,3) :: Sijp,Sigmap,dSigmap
!
! Identify module.
!
if (headtt) then
print*,'dcaustics_dt_pencil: Calculate dcaustics_dt_pencil'
print*,'called from dvvp_dt_pencil in the particles_dust module'
endif
!
! Loop over all particles in current pencil.
!
Sigma_lin = fp(k,isigmap11:isigmap33)
call linarray2matrix(Sigma_lin,Sigmap)
!
! interpolate the gradu matrix to particle positions
!
call interpolate_linear(f,igu11,igu33,fp(k,ixp:izp),Sij_lin,ineargrid(k,:),0,ipar(k))
call linarray2matrix(Sij_lin,Sijp)
!
! solve dynamical equation
!
dSigmap = taup1*(Sijp-Sigmap) - matmul(Sigmap,Sigmap)
call matrix2linarray(dSigmap,dSigmap_lin)
dfp(k,isigmap11:isigmap33) = dSigmap_lin
!
! update the (logarithm of) volume
!
dfp(k,ilnVp) = Sigmap(1,1)+Sigmap(2,2)+Sigmap(3,3)
!
endsubroutine dcaustics_dt_pencil
!***********************************************************************
subroutine read_pcaustics_init_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=particles_caustics_init_pars, IOSTAT=iostat)
!
endsubroutine read_pcaustics_init_pars
!***********************************************************************
subroutine write_pcaustics_init_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=particles_caustics_init_pars)
!
endsubroutine write_pcaustics_init_pars
!***********************************************************************
subroutine read_pcaustics_run_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=particles_caustics_run_pars, IOSTAT=iostat)
!
endsubroutine read_pcaustics_run_pars
!***********************************************************************
subroutine write_pcaustics_run_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=particles_caustics_run_pars)
!
endsubroutine write_pcaustics_run_pars
!***********************************************************************
subroutine rprint_particles_caustics(lreset,lwrite)
!
! Read and register print parameters relevant for particles.
!
! may-2016/dhruba+akshay: coded
!
use Diagnostics
use FArrayManager, only: farray_index_append
use General, only: itoa
!
logical :: lreset
logical, optional :: lwrite
!
integer :: iname
integer :: k
logical :: lwr
character (len=intlen) :: srad
!
! Write information to index.pro.
!
lwr = .false.
if (present(lwrite)) lwr=lwrite
!
if (lwr) then
call farray_index_append('isigmap11', isigmap11)
call farray_index_append('isigmap12', isigmap13)
call farray_index_append('isigmap13', isigmap13)
call farray_index_append('isigmap21', isigmap21)
call farray_index_append('isigmap22', isigmap22)
call farray_index_append('isigmap23', isigmap23)
call farray_index_append('isigmap31', isigmap31)
call farray_index_append('isigmap32', isigmap32)
call farray_index_append('isigmap33', isigmap33)
call farray_index_append('ilnVp', ilnVp)
call farray_index_append('iPPp', iPPp)
call farray_index_append('iQQp', iQQp)
call farray_index_append('iRRp', iRRp)
call farray_index_append('iblowup', iblowup)
endif
!
! Reset everything in case of reset.
!
if (lreset) then
idiag_TrSigmapm=0
idiag_blowupm=0
idiag_ncaustics=0
idiag_lnVpm=0
endif
!
! Run through all possible names that may be listed in print.in.
!
if (lroot .and. ip<14) print*,'rprint_particles: run through parse list'
do iname=1,nname
call parse_name(iname,cname(iname),cform(iname),'trsigmapm',idiag_TrSigmapm)
call parse_name(iname,cname(iname),cform(iname),'blowupm',idiag_blowupm)
call parse_name(iname,cname(iname),cform(iname),'lnVpm',idiag_lnVpm)
do k=1,ndustrad
srad=itoa(k)
call parse_name(iname,cname(iname),cform(iname), &
'ncaustics'//trim(srad),idiag_ncaustics(k))
enddo
enddo
!
endsubroutine rprint_particles_caustics
!***********************************************************************
subroutine reinitialize_caustics(fp)
!
real, dimension (mpar_loc,mparray), intent (out) :: fp
integer :: ip
if (lroot) then
print*, 'The sigmap matrix always starts from zero,'
print*, 'even when we restart from earlier runs.'
endif
!
! We set the gradient of the particle velocity field as zero.
!
do ip=1,npar_loc
fp(ip,isigmap11:isigmap33) = 0.
fp(ip,ilnVp) = 0.
fp(ip,iPPp) = 0.
fp(ip,iQQp) = 0.
fp(ip,iRRp) = 0.
fp(ip,iblowup) = 0.
enddo
endsubroutine reinitialize_caustics
!***********************************************************************
subroutine reset_caustics(fp)
!
real, dimension (mpar_loc,mparray), intent (out) :: fp
integer :: ip
real :: TrSigma, QSigma, RSigma
real, dimension(9) :: Sigma_lin
real, dimension(3,3) :: Sigmap
do ip=1,npar_loc
if (fp(ip,iPPp).lt.TrSigma_Cutoff) then
fp(ip,iblowup)=fp(ip,iblowup)+1.
fp(ip,isigmap11:isigmap33) = 0.
endif
enddo
!
endsubroutine reset_caustics
!***********************************************************************
endmodule Particles_caustics