Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Dust growth: mass-weighted interpolations of dust-gas quantities #377

Merged
merged 1 commit into from
Mar 7, 2023
Merged

Dust growth: mass-weighted interpolations of dust-gas quantities #377

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
28 changes: 8 additions & 20 deletions src/main/force.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -883,7 +883,6 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g
use kernel, only:wkern_drag,cnormk_drag
use part, only:ndustsmall,grainsize,graindens
#ifdef DUSTGROWTH
use growth, only:wbymass
use kernel, only:wkern,cnormk
#endif
#endif
Expand Down Expand Up @@ -1831,18 +1830,12 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g
else
winter = wkern(q2j,qj)*hj21*hj1*cnormk
endif
!--following quantities are weighted by mass rather than mass/density
fsum(idensgasi) = fsum(idensgasi) + pmassj*winter
if (wbymass) then
fsum(idvix) = fsum(idvix) + pmassj*dvx*winter
fsum(idviy) = fsum(idviy) + pmassj*dvy*winter
fsum(idviz) = fsum(idviz) + pmassj*dvz*winter
fsum(icsi) = fsum(icsi) + pmassj*spsoundj*winter
else
fsum(idvix) = fsum(idvix) + pmassj/rhoj*dvx*winter
fsum(idviy) = fsum(idviy) + pmassj/rhoj*dvy*winter
fsum(idviz) = fsum(idviz) + pmassj/rhoj*dvz*winter
fsum(icsi) = fsum(icsi) + pmassj/rhoj*spsoundj*winter
endif
fsum(idvix) = fsum(idvix) + pmassj*dvx*winter
fsum(idviy) = fsum(idviy) + pmassj*dvy*winter
fsum(idviz) = fsum(idviz) + pmassj*dvz*winter
fsum(icsi) = fsum(icsi) + pmassj*spsoundj*winter
#endif
else
!--the following works for large grains only (not hybrid large and small grains)
Expand Down Expand Up @@ -2482,7 +2475,6 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv
use utils_gr, only:get_u0
use io, only:error
#ifdef DUSTGROWTH
use growth, only:wbymass
use dust, only:idrag,get_ts
use part, only:Omega_k
#endif
Expand Down Expand Up @@ -2954,13 +2946,9 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv
if (iamdusti) then
!- return interpolations to their respective arrays
dustgasprop(2,i) = fsum(idensgasi) !- rhogas
dustgasprop(4,i) = sqrt(fsum(idvix)**2 + fsum(idviy)**2 + fsum(idviz)**2) !- dv
dustgasprop(1,i) = fsum(icsi)
!- if interpolations are mass weigthed, divide result by rhog,i
if (wbymass) then
dustgasprop(1,i) = dustgasprop(1,i)/dustgasprop(2,i) !- sound speed
dustgasprop(4,i) = dustgasprop(4,i)/dustgasprop(2,i) !- |dv|
endif
!- interpolations are mass weigthed, divide result by rhog,i
dustgasprop(4,i) = sqrt(fsum(idvix)**2 + fsum(idviy)**2 + fsum(idviz)**2)/dustgasprop(2,i) !- |dv|
dustgasprop(1,i) = fsum(icsi)/dustgasprop(2,i) !- sound speed

!- get the Stokes number with get_ts using the interpolated quantities
rhoi = xpartveci(irhoi)
Expand Down
22 changes: 7 additions & 15 deletions src/main/growth.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,6 @@ module growth
! - vfragout : *inward fragmentation threshold in m/s*
! - cohacc : *strength of the cohesive acceleration in g/s^2*
! - dsize : *size of ejected grain during erosion in cm*
! - wbymass : *weight dustgasprops by mass rather than mass/density*
!
! :Dependencies: checkconserved, dim, dust, eos, infile_utils, io, options,
! part, physcon, table_utils, units, viscosity
Expand Down Expand Up @@ -63,9 +62,6 @@ module growth
real, public :: cohacc
real, public :: dsize


logical, public :: wbymass = .true.

#ifdef MCFOST
logical, public :: f_smax = .false.
real, public :: size_max = 0.2 !- cm
Expand Down Expand Up @@ -350,7 +346,6 @@ subroutine write_options_growth(iunit)
integer, intent(in) :: iunit

write(iunit,"(/,a)") '# options controlling growth'
call write_inopt(wbymass,'wbymass','weight dustgasprops by mass rather than mass/density',iunit)
if (nptmass > 1) call write_inopt(this_is_a_flyby,'flyby','use primary for keplerian freq. calculation',iunit)
call write_inopt(ifrag,'ifrag','dust fragmentation (0=off,1=on,2=Kobayashi)',iunit)
call write_inopt(ieros,'ieros','erosion of dust (0=off,1=on)',iunit)
Expand Down Expand Up @@ -434,9 +429,6 @@ subroutine read_options_growth(name,valstring,imatch,igotall,ierr)
read(valstring,*,iostat=ierr) this_is_a_flyby
ngot = ngot + 1
if (nptmass < 2) tmp = .true.
case('wbymass')
read(valstring,*,iostat=ierr) wbymass
ngot = ngot + 1
#ifdef MCFOST
case('force_smax')
read(valstring,*,iostat=ierr) f_smax
Expand All @@ -456,23 +448,23 @@ subroutine read_options_growth(name,valstring,imatch,igotall,ierr)
imcf = 3
#endif

if (ieros == 1) goteros = 2
if (ieros == 1) goteros = 3

if (nptmass > 1 .or. tmp) then
if ((ifrag <= 0) .and. ngot == 3+imcf+goteros) igotall = .true.
if ((ifrag <= 0) .and. ngot == 2+imcf+goteros) igotall = .true.
if (isnow == 0) then
if (ngot == 6+imcf+goteros) igotall = .true.
if (ngot == 5+imcf+goteros) igotall = .true.
elseif (isnow > 0) then
if (ngot == 8+imcf+goteros) igotall = .true.
if (ngot == 7+imcf+goteros) igotall = .true.
else
igotall = .false.
endif
else
if ((ifrag <= 0) .and. ngot == 2+imcf+goteros) igotall = .true.
if ((ifrag <= 0) .and. ngot == 1+imcf+goteros) igotall = .true.
if (isnow == 0) then
if (ngot == 5+imcf+goteros) igotall = .true.
if (ngot == 4+imcf+goteros) igotall = .true.
elseif (isnow > 0) then
if (ngot == 7+imcf+goteros) igotall = .true.
if (ngot == 6+imcf+goteros) igotall = .true.
else
igotall = .false.
endif
Expand Down
3 changes: 1 addition & 2 deletions src/tests/test_growth.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ subroutine test_farmingbox(ntests,npass,frag,onefluid)
use testutils, only:checkvalbuf,checkvalbuf_end
use eos, only:ieos,polyk,gamma,get_spsound
use dust, only:idrag,init_drag
use growth, only:ifrag,init_growth,isnow,vfrag,wbymass,gsizemincgs
use growth, only:ifrag,init_growth,isnow,vfrag,gsizemincgs
use options, only:alpha,alphamax,use_dustfrac
use unifdis, only:set_unifdis
use dim, only:periodic,mhd,use_dust,maxp,maxalpha
Expand Down Expand Up @@ -187,7 +187,6 @@ subroutine test_farmingbox(ntests,npass,frag,onefluid)
! initialise
!
this_is_a_test = .true.
wbymass = .false.

!
! setup for dustybox problem
Expand Down