From f51413ccdc8f741312743a779f65a4874c78bf09 Mon Sep 17 00:00:00 2001 From: Stephane Michoulier Date: Tue, 7 Mar 2023 17:49:18 +1100 Subject: [PATCH] Dust growth: mass-weighted interpolations of dust-gas quantities are now default, option removed --- src/main/force.F90 | 28 ++++++++-------------------- src/main/growth.F90 | 22 +++++++--------------- src/tests/test_growth.F90 | 3 +-- 3 files changed, 16 insertions(+), 37 deletions(-) diff --git a/src/main/force.F90 b/src/main/force.F90 index 10345b224..7dcc4a344 100644 --- a/src/main/force.F90 +++ b/src/main/force.F90 @@ -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 @@ -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) @@ -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 @@ -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) diff --git a/src/main/growth.F90 b/src/main/growth.F90 index 4417089f4..e84c0b55e 100644 --- a/src/main/growth.F90 +++ b/src/main/growth.F90 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 diff --git a/src/tests/test_growth.F90 b/src/tests/test_growth.F90 index ee2813444..47623e125 100644 --- a/src/tests/test_growth.F90 +++ b/src/tests/test_growth.F90 @@ -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 @@ -187,7 +187,6 @@ subroutine test_farmingbox(ntests,npass,frag,onefluid) ! initialise ! this_is_a_test = .true. - wbymass = .false. ! ! setup for dustybox problem