diff --git a/AUTHORS b/AUTHORS index dc6272205..4b972f1d7 100644 --- a/AUTHORS +++ b/AUTHORS @@ -22,8 +22,8 @@ Elisabeth Borchert Ward Homan Christophe Pinte Terrence Tricco -Simone Ceppi Stephane Michoulier +Simone Ceppi Spencer Magnall Caitlyn Hardiman Enrico Ragusa @@ -31,8 +31,8 @@ Sergei Biriukov Cristiano Longarini Giovanni Dipierro Roberto Iaconi -Amena Faruqi Hauke Worpel +Amena Faruqi Alison Young Stephen Neilson <36410751+s-neilson@users.noreply.github.com> Martina Toscani @@ -49,21 +49,23 @@ Phantom benchmark bot Kieran Hirsh Nicole Rodrigues David Trevascus -Farzana Meru Nicolás Cuello +Farzana Meru +Mike Lau Chris Nixon Miguel Gonzalez-Bolivar -Benoit Commercon -Giulia Ballabio -Joe Fisher -Maxime Lombart -Mike Lau Orsola De Marco +Maxime Lombart +Joe Fisher Zachary Pellow +Benoit Commercon +Giulia Ballabio s-neilson <36410751+s-neilson@users.noreply.github.com> -Cox, Samuel +MICHOULIER Stephane +Steven Rieder Jeremy Smallwood +Cox, Samuel Jorge Cuadra -Steven Rieder Stéven Toupin Taj Jankovič +Chunliang Mu diff --git a/build/Makefile b/build/Makefile index 219f2eb0b..32282cb91 100644 --- a/build/Makefile +++ b/build/Makefile @@ -96,7 +96,7 @@ ifndef SRCNIMHD endif ifndef SRCDUST - SRCDUST = dust.f90 growth.F90 + SRCDUST = dust.f90 growth.f90 porosity.f90 endif ifdef SMOL @@ -111,10 +111,6 @@ else SRCINJECT=utils_binary.f90 set_binary.f90 inject_wind.f90 endif -#ifndef SRCGROWTH -# SRCGROWTH = growth.F90 -#endif - #--- live feedback from mcfost ifeq ($(MCFOST), yes) ANALYSIS= analysis_mcfost.f90 @@ -469,7 +465,7 @@ SRCPOTS= extern_gr.F90 \ extern_spiral.f90 \ extern_lensethirring.f90 \ extern_gnewton.f90 \ - lumin_nsdisc.F90 extern_prdrag.f90 \ + lumin_nsdisc.f90 extern_prdrag.f90 \ extern_Bfield.f90 \ extern_densprofile.f90 \ extern_staticsine.f90 \ @@ -512,9 +508,9 @@ SRCMESA= eos_mesa_microphysics.f90 eos_mesa.f90 SRCEOS = eos_barotropic.f90 eos_stratified.f90 eos_piecewise.f90 ${SRCMESA} eos_shen.f90 eos_helmholtz.f90 eos_idealplusrad.f90 ionization.F90 eos_gasradrec.f90 eos.f90 ifeq ($(HDF5), yes) - SRCREADWRITE_DUMPS= utils_hdf5.f90 utils_dumpfiles_hdf5.f90 readwrite_dumps_common.F90 readwrite_dumps_fortran.F90 readwrite_dumps_hdf5.F90 readwrite_dumps.F90 + SRCREADWRITE_DUMPS= utils_hdf5.f90 utils_dumpfiles_hdf5.f90 readwrite_dumps_common.f90 readwrite_dumps_fortran.F90 readwrite_dumps_hdf5.F90 readwrite_dumps.F90 else - SRCREADWRITE_DUMPS= readwrite_dumps_common.F90 readwrite_dumps_fortran.F90 readwrite_dumps.F90 + SRCREADWRITE_DUMPS= readwrite_dumps_common.f90 readwrite_dumps_fortran.F90 readwrite_dumps.F90 endif ifeq ($(KROME), krome) @@ -538,12 +534,12 @@ SOURCES= physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 units.f90 \ partinject.F90 utils_inject.f90 dust_formation.f90 ptmass_radiation.f90 ptmass_heating.f90 \ utils_deriv.f90 utils_implicit.f90 radiation_implicit.f90 ${SRCTURB} \ ${SRCINJECT_DEPS} wind_equations.f90 wind.F90 ${SRCINJECT} \ - ${SRCKROME} memory.F90 ${SRCREADWRITE_DUMPS} \ + ${SRCKROME} memory.f90 ${SRCREADWRITE_DUMPS} \ quitdump.f90 ptmass.F90 \ readwrite_infile.F90 dens.F90 force.F90 deriv.F90 energies.F90 sort_particles.f90 \ utils_shuffleparticles.F90 evwrite.f90 step_leapfrog.F90 writeheader.F90 ${SRCAN} step_supertimestep.F90 \ mf_write.f90 evolve.F90 utils_orbits.f90 utils_linalg.f90 \ - checksetup.F90 initial.F90 + checksetup.f90 initial.F90 # Needed as einsteintk_wrapper depends on initial ifeq ($(GR),yes) @@ -610,20 +606,19 @@ edit: checksetup #---------------------------------------------------- # these are the sources for anything which uses the readwrite_dumps module # -SRCDUMP= physcon.f90 ${CONFIG} ${SRCKERNEL} utils_omp.F90 io.F90 units.f90 \ - boundary.f90 boundary_dynamic.f90 mpi_utils.F90 \ +SRCDUMP= physcon.f90 ${CONFIG} ${SRCKERNEL} utils_omp.F90 io.F90 units.f90 mpi_utils.F90 \ utils_timing.f90 utils_infiles.f90 dtype_kdtree.f90 utils_allocate.f90 part.F90 \ - ${DOMAIN} mpi_dens.F90 mpi_force.F90 \ + ${DOMAIN} mpi_dens.F90 mpi_force.F90 boundary.f90 boundary_dynamic.f90 \ mpi_balance.F90 mpi_memory.f90 mpi_derivs.F90 mpi_tree.F90 kdtree.F90 linklist_kdtree.F90 \ utils_dumpfiles.f90 utils_vectors.f90 utils_mathfunc.f90 \ utils_datafiles.f90 utils_filenames.f90 utils_system.f90 utils_tables.f90 datafiles.f90 gitinfo.f90 \ centreofmass.f90 \ timestep.f90 ${SRCEOS} cullendehnen.f90 dust_formation.f90 \ ${SRCGR} ${SRCPOT} \ - memory.F90 \ + memory.f90 \ utils_sphNG.f90 \ setup_params.f90 ${SRCFASTMATH} checkoptions.F90 \ - viscosity.f90 damping.f90 options.f90 checkconserved.f90 prompting.f90 ${SRCDUST} \ + viscosity.f90 damping.f90 options.f90 checkconserved.f90 prompting.f90 dust.f90 \ ${SRCREADWRITE_DUMPS} \ utils_sort.f90 sort_particles.f90 OBJDUMP1= $(SRCDUMP:.f90=.o) @@ -685,7 +680,7 @@ else endif SRCTESTS=utils_testsuite.f90 ${TEST_FASTMATH} test_kernel.f90 \ - test_dust.F90 test_growth.f90 test_smol.F90 \ + test_dust.f90 test_growth.f90 test_smol.F90 \ test_nonidealmhd.F90 directsum.f90 test_gravity.f90 \ test_derivs.F90 test_cooling.f90 test_eos_stratified.f90 \ test_eos.f90 test_externf.f90 test_rwdump.f90 \ diff --git a/build/Makefile_setups b/build/Makefile_setups index ae562595d..85ccc0e0e 100644 --- a/build/Makefile_setups +++ b/build/Makefile_setups @@ -420,7 +420,7 @@ endif ifeq ($(SETUP), shock) # shock tube tests PERIODIC=yes - SETUPFILE= setup_shock.F90 + SETUPFILE= setup_shock.f90 KERNEL=quintic KNOWN_SETUP=yes endif @@ -428,7 +428,7 @@ endif ifeq ($(SETUP), dustyshock) # shock tube tests with dust PERIODIC=yes - SETUPFILE= setup_shock.F90 + SETUPFILE= setup_shock.f90 DUST=yes KERNEL=quintic KNOWN_SETUP=yes @@ -437,7 +437,7 @@ endif ifeq ($(SETUP), mhdshock) # Ryu & Brio-Wu shock tube tests PERIODIC=yes - SETUPFILE= setup_shock.F90 + SETUPFILE= setup_shock.f90 MHD=yes KERNEL=quintic KNOWN_SETUP=yes @@ -446,7 +446,7 @@ endif ifeq ($(SETUP), nimhdshock) # non-ideal mhd standing and C shock tests PERIODIC=yes - SETUPFILE= setup_shock.F90 + SETUPFILE= setup_shock.f90 MHD=yes STS_TIMESTEPS=no NONIDEALMHD=yes @@ -460,7 +460,7 @@ ifeq ($(SETUP), radshock) # shock tube in radiation hydrodynamics PERIODIC=yes RADIATION=yes - SETUPFILE= setup_shock.F90 + SETUPFILE= setup_shock.f90 KERNEL=quintic KNOWN_SETUP=yes endif @@ -468,7 +468,7 @@ endif ifeq ($(SETUP), srshock) # special relativistic sod shock tube test PERIODIC=yes - SETUPFILE= setup_shock.F90 + SETUPFILE= setup_shock.f90 KERNEL=quintic GR=yes METRIC=minkowski @@ -479,7 +479,7 @@ endif ifeq ($(SETUP), testparticles) # test particles - SETUPFILE= setup_testparticles.F90 + SETUPFILE= setup_testparticles.f90 KNOWN_SETUP=yes MAXP=500000 ANALYSIS= analysis_1particle.f90 @@ -1037,7 +1037,7 @@ ifeq ($(SETUP), testgr) endif ifeq ($(SETUP), flrw) -# constant density FLRW cosmology with perturbations +# constant density FLRW cosmology with perturbations GR=yes KNOWN_SETUP=yes IND_TIMESTEPS=no diff --git a/src/main/checksetup.F90 b/src/main/checksetup.f90 similarity index 100% rename from src/main/checksetup.F90 rename to src/main/checksetup.f90 diff --git a/src/main/config.F90 b/src/main/config.F90 index 92faa467c..069005dd4 100644 --- a/src/main/config.F90 +++ b/src/main/config.F90 @@ -137,7 +137,7 @@ module dim ! xpartveci integer, parameter :: maxxpartvecidens = 14 + radenxpartvetden - integer, parameter :: maxxpartvecvars = 61 ! Number of scalars in xpartvec + integer, parameter :: maxxpartvecvars = 62 ! Number of scalars in xpartvec integer, parameter :: maxxpartvecarrs = 2 ! Number of arrays in xpartvec integer, parameter :: maxxpartvecGR = 33 ! Number of GR values in xpartvec (1 for dens, 16 for gcov, 16 for gcon) integer, parameter :: maxxpartveciforce = maxxpartvecvars + & ! Total number of values diff --git a/src/main/deriv.F90 b/src/main/deriv.F90 index f86a8ba63..fd6b06ceb 100644 --- a/src/main/deriv.F90 +++ b/src/main/deriv.F90 @@ -14,10 +14,10 @@ module deriv ! ! :Runtime parameters: None ! -! :Dependencies: cons2prim, densityforce, derivutils, dim, externalforces, -! forces, forcing, growth, io, linklist, metric_tools, options, part, -! ptmass, ptmass_radiation, radiation_implicit, timestep, timestep_ind, -! timing +! :Dependencies: cons2prim, densityforce, derivutils, dim, dust_formation, +! externalforces, forces, forcing, growth, io, linklist, metric_tools, +! options, part, porosity, ptmass, ptmass_radiation, radiation_implicit, +! timestep, timestep_ind, timing ! implicit none @@ -36,7 +36,7 @@ module deriv !------------------------------------------------------------- subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& Bevol,dBevol,rad,drad,radprop,dustprop,ddustprop,& - dustevol,ddustevol,dustfrac,eos_vars,time,dt,dtnew,pxyzu,dens,metrics) + dustevol,ddustevol,filfac,dustfrac,eos_vars,time,dt,dtnew,pxyzu,dens,metrics) use dim, only:maxvxyzu,mhd,fast_divcurlB,gr,periodic,do_radiation,& sink_radiation,use_dustgrowth,ind_timesteps use io, only:iprint,fatal,error @@ -51,7 +51,8 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& #ifdef DRIVING use forcing, only:forceit #endif - use growth, only:get_growth_rate + use growth, only:get_growth_rate + use porosity, only:get_disruption,get_probastick use ptmass_radiation, only:get_dust_temperature use timing, only:get_timings use forces, only:force @@ -60,7 +61,7 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& use cons2prim, only:cons2primall,cons2prim_everything,prim2consall use metric_tools, only:init_metric use radiation_implicit, only:do_radiation_implicit,ierr_failed_to_converge - use options, only:implicit_radiation,implicit_radiation_store_drad + use options, only:implicit_radiation,implicit_radiation_store_drad,use_porosity integer, intent(in) :: icall integer, intent(inout) :: npart integer, intent(in) :: nactive @@ -80,6 +81,7 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& real, intent(inout) :: dustprop(:,:) real, intent(out) :: dustfrac(:,:) real, intent(out) :: ddustevol(:,:),ddustprop(:,:) + real, intent(inout) :: filfac(:) real, intent(in) :: time,dt real, intent(out) :: dtnew real, intent(inout) :: pxyzu(:,:), dens(:) @@ -121,6 +123,10 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& call do_timing('link',tlast,tcpulast,start=.true.) + ! + ! compute disruption of dust particles + ! + if (use_dustgrowth .and. use_porosity) call get_disruption(npart,xyzh,filfac,dustprop,dustgasprop) ! ! calculate density by direct summation ! @@ -170,7 +176,9 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& call do_timing('force',tlast,tcpulast) if (use_dustgrowth) then ! compute growth rate of dust particles - call get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,ddustprop(1,:))!--we only get ds/dt (i.e 1st dimension of ddustprop) + call get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,filfac,ddustprop(1,:))!--we only get dm/dt (i.e 1st dimension of ddustprop) + ! compute growth rate and probability of sticking/bouncing of porous dust + if (use_porosity) call get_probastick(npart,xyzh,ddustprop(1,:),dustprop,dustgasprop,filfac) endif ! @@ -214,7 +222,7 @@ end subroutine derivs !-------------------------------------- subroutine get_derivs_global(tused,dt_new,dt) use part, only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& - Bevol,dBevol,rad,drad,radprop,dustprop,ddustprop,& + Bevol,dBevol,rad,drad,radprop,dustprop,ddustprop,filfac,& dustfrac,ddustevol,eos_vars,pxyzu,dens,metrics,dustevol use timing, only:printused,getused use io, only:id,master @@ -230,8 +238,8 @@ subroutine get_derivs_global(tused,dt_new,dt) if (present(dt)) dti = dt call getused(t1) call derivs(1,npart,npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,& - rad,drad,radprop,dustprop,ddustprop,dustevol,ddustevol,dustfrac,eos_vars,& - time,dti,dtnew,pxyzu,dens,metrics) + rad,drad,radprop,dustprop,ddustprop,dustevol,ddustevol,filfac,dustfrac,& + eos_vars,time,dti,dtnew,pxyzu,dens,metrics) call getused(t2) if (id==master .and. present(tused)) call printused(t1) if (present(tused)) tused = t2 - t1 diff --git a/src/main/dust.f90 b/src/main/dust.f90 index fda11c216..8270c5ea9 100644 --- a/src/main/dust.f90 +++ b/src/main/dust.f90 @@ -47,6 +47,7 @@ module dust public :: print_dustinfo public :: write_options_dust public :: read_options_dust + public :: get_viscmol_nu real, private :: cste_mu,coeff_gei_1,seff private @@ -413,4 +414,11 @@ subroutine read_options_dust(name,valstring,imatch,igotall,ierr) end subroutine read_options_dust +real function get_viscmol_nu(spsoundgas,rhogas) + real,intent(in) :: spsoundgas,rhogas + + get_viscmol_nu = cste_mu*seff*spsoundgas/rhogas + +end function get_viscmol_nu + end module dust diff --git a/src/main/extern_gr.F90 b/src/main/extern_gr.f90 similarity index 78% rename from src/main/extern_gr.F90 rename to src/main/extern_gr.f90 index 8696ffd10..3d3aacdb2 100644 --- a/src/main/extern_gr.F90 +++ b/src/main/extern_gr.f90 @@ -6,9 +6,12 @@ !--------------------------------------------------------------------------! module extern_gr ! -! None +! Compute terms related to derivatives of the metric which appear +! on the right hand side of the momentum equation ! -! :References: None +! :References: +! Liptai & Price (2019), MNRAS 485, 819 +! Magnall, Price, Lasky & Macpherson (2023), Phys. Rev D. 108, 103534 ! ! :Owner: Spencer Magnall ! @@ -49,6 +52,12 @@ subroutine get_grforce(xyzhi,metrici,metricderivsi,veli,densi,ui,pi,fexti,dtf) end subroutine get_grforce +!--------------------------------------------------------------- +!+ +! Wrapper of the above, computing accelerations due to metric +! gradients on all particles +!+ +!--------------------------------------------------------------- subroutine get_grforce_all(npart,xyzh,metrics,metricderivs,vxyzu,dens,fext,dtexternal) use timestep, only:C_force use eos, only:ieos,get_pressure @@ -77,8 +86,12 @@ subroutine get_grforce_all(npart,xyzh,metrics,metricderivs,vxyzu,dens,fext,dtext end subroutine get_grforce_all -!--- Subroutine to calculate the timestep constraint from the 'external force' -! this is multiplied by the safety factor C_force elsewhere +!--------------------------------------------------------------------------- +!+ +! Subroutine to calculate the timestep constraint from the 'external force' +! this is multiplied by the safety factor C_force elsewhere +!+ +!--------------------------------------------------------------------------- subroutine dt_grforce(xyzh,fext,dtf) use physcon, only:pi use metric_tools, only:imetric,imet_schwarzschild,imet_kerr @@ -107,7 +120,6 @@ subroutine dt_grforce(xyzh,fext,dtf) end subroutine dt_grforce - !---------------------------------------------------------------- !+ ! Compute the source terms required on the right hand side of @@ -139,20 +151,19 @@ pure subroutine forcegr(x,metrici,metricderivsi,v,dens,u,p,fterm,ierr) ! energy-momentum tensor times sqrtg on 2rho* do i=0,3 - term(0:3,i) = 0.5*(enth*uzero*v4(0:3)*v4(i) + P*gcon(0:3,i)/(dens*uzero)) + term(0:3,i) = 0.5*(enth*uzero*v4(0:3)*v4(i) + P*gcon(0:3,i)/(dens*uzero)) enddo ! source term fterm = 0. do i=0,3 - fterm(1) = fterm(1) + dot_product(term(:,i),metricderivsi(:,i,1)) - fterm(2) = fterm(2) + dot_product(term(:,i),metricderivsi(:,i,2)) - fterm(3) = fterm(3) + dot_product(term(:,i),metricderivsi(:,i,3)) + fterm(1) = fterm(1) + dot_product(term(:,i),metricderivsi(:,i,1)) + fterm(2) = fterm(2) + dot_product(term(:,i),metricderivsi(:,i,2)) + fterm(3) = fterm(3) + dot_product(term(:,i),metricderivsi(:,i,3)) enddo end subroutine forcegr - !-------- I don't think this is actually being used at the moment.... subroutine update_grforce_leapfrog(vhalfx,vhalfy,vhalfz,fxi,fyi,fzi,fexti,dt,xi,yi,zi,densi,ui,pi) use io, only:fatal @@ -220,82 +231,45 @@ subroutine update_grforce_leapfrog(vhalfx,vhalfy,vhalfz,fxi,fyi,fzi,fexti,dt,xi, end subroutine update_grforce_leapfrog +!---------------------------------------------------------------- +!+ +! compute stress energy tensor on all particles +!+ +!---------------------------------------------------------------- subroutine get_tmunu_all(npart,xyzh,metrics,vxyzu,metricderivs,dens,tmunus) - use eos, only:ieos,get_pressure - use part, only:isdead_or_accreted + use eos, only:ieos,get_pressure + use part, only:isdead_or_accreted integer, intent(in) :: npart real, intent(in) :: xyzh(:,:), metrics(:,:,:,:), metricderivs(:,:,:,:), dens(:) real, intent(inout) :: vxyzu(:,:),tmunus(:,:,:) real :: pi integer :: i - logical :: verbose - verbose = .false. - ! TODO write openmp parallel code !$omp parallel do default(none) & !$omp shared(npart,xyzh,metrics,vxyzu,dens,ieos,tmunus) & - !$omp private(i,pi,verbose) - do i=1, npart - !print*, "i: ", i - if (i==1) then - verbose = .true. - else - verbose = .false. - endif + !$omp private(i,pi) + do i=1,npart if (.not.isdead_or_accreted(xyzh(4,i))) then pi = get_pressure(ieos,xyzh(:,i),dens(i),vxyzu(:,i)) - call get_tmunu(xyzh(:,i),metrics(:,:,:,i),& - vxyzu(1:3,i),dens(i),vxyzu(4,i),pi,tmunus(:,:,i),verbose) + call get_tmunu(xyzh(:,i),metrics(:,:,:,i),vxyzu(1:3,i),& + dens(i),vxyzu(4,i),pi,tmunus(:,:,i)) endif enddo !$omp end parallel do - !print*, "tmunu calc val is: ", tmunus(0,0,5) -end subroutine get_tmunu_all - -subroutine get_tmunu_all_exact(npart,xyzh,metrics,vxyzu,metricderivs,dens,tmunus) - use eos, only:ieos,get_pressure - use part, only:isdead_or_accreted - integer, intent(in) :: npart - real, intent(in) :: xyzh(:,:), metrics(:,:,:,:), metricderivs(:,:,:,:), dens(:) - real, intent(inout) :: vxyzu(:,:),tmunus(:,:,:) - real :: pi - integer :: i - logical :: firstpart - real :: tmunu(4,4) - !print*, "entered get tmunu_all_exact" - tmunu = 0. - firstpart = .true. - ! TODO write openmp parallel code - do i=1, npart - if (.not.isdead_or_accreted(xyzh(4,i)) .and. firstpart) then - pi = get_pressure(ieos,xyzh(:,i),dens(i),vxyzu(:,i)) - call get_tmunu_exact(xyzh(:,i),metrics(:,:,:,i), metricderivs(:,:,:,i), & - vxyzu(1:3,i),dens(i),vxyzu(4,i),pi,tmunus(:,:,i)) - !print*, "finished get_tmunu call!" - firstpart = .false. - !print*, "tmunu: ", tmunu - !print*, "tmunus: ", tmunus(:,:,i) - tmunu(:,:) = tmunus(:,:,i) - !print*, "Got tmunu val: ", tmunu - !stop - else - !print*, "setting tmunu for part: ", i - tmunus(:,:,i) = tmunu(:,:) - endif - - enddo - !print*, "tmunu calc val is: ", tmunus(0,0,5) -end subroutine get_tmunu_all_exact +end subroutine get_tmunu_all -! Subroutine to calculate the covariant form of the stress energy tensor -! For a particle at position p -subroutine get_tmunu(x,metrici,v,dens,u,p,tmunu,verbose) - use metric_tools, only:unpack_metric +!------------------------------------------------------------------------- +!+ +! calculate the covariant form of the stress energy tensor +! for a particle at position x +!+ +!------------------------------------------------------------------------- +subroutine get_tmunu(x,metrici,v,dens,u,p,tmunu) + use metric_tools, only:unpack_metric use utils_gr, only:get_u0 real, intent(in) :: x(3),metrici(:,:,:),v(3),dens,u,p real, intent(out) :: tmunu(0:3,0:3) - logical, optional, intent(in) :: verbose real :: w,v4(0:3),uzero,u_upper(0:3),u_lower(0:3) real :: gcov(0:3,0:3), gcon(0:3,0:3) real :: gammaijdown(1:3,1:3),betadown(3),alpha @@ -318,42 +292,7 @@ subroutine get_tmunu(x,metrici,v,dens,u,p,tmunu,verbose) ! Get cov and con versions of the metric + spatial metric and lapse and shift ! Not entirely convinced that the lapse and shift calculations are acccurate for the general case!! - !print*, "Before unpack metric " call unpack_metric(metrici,gcov=gcov,gcon=gcon,gammaijdown=gammaijdown,alpha=alpha,betadown=betadown) - !print*, "After unpack metric" - -! if (present(verbose) .and. verbose) then -! ! Do we get sensible values -! print*, "Unpacked metric quantities..." -! print*, "gcov: ", gcov -! print*, "gcon: ", gcon -! print*, "gammaijdown: ", gammaijdown -! print* , "alpha: ", alpha -! print*, "betadown: ", betadown -! print*, "v4: ", v4 -! endif - - ! ! Need to change Betadown to betaup - ! ! Won't matter at this point as it is allways zero - ! ! get big V - ! bigV(:) = (v(:) + betadown)/alpha - - ! ! We need the covariant version of the 3 velocity - ! ! gamma_ij v^j = v_i where gamma_ij is the spatial metric - ! do i=1, 3 - ! vcov(i) = gammaijdown(i,1)*bigv(1) + gammaijdown(i,2)*bigv(2) + gammaijdown(i,3)*bigv(3) - ! enddo - - - ! ! Calculate the lorentz factor - ! lorentz = (1. - (vcov(1)*bigv(1) + vcov(2)*bigv(2) + vcov(3)*bigv(3)))**(-0.5) - - ! ! Calculate the 4-velocity - ! velshiftterm = vcov(1)*betadown(1) + vcov(2)*betadown(2) + vcov(3)*betadown(3) - ! v4(0) = lorentz*(-alpha + velshiftterm) - ! ! This should be vcov not v - ! v4(1:3) = lorentz*vcov(1:3) - ! We are going to use the same Tmunu calc as force GR ! And then lower it using the metric @@ -363,7 +302,6 @@ subroutine get_tmunu(x,metrici,v,dens,u,p,tmunu,verbose) v4(0) = 1. v4(1:3) = v(:) - ! first component of the upper-case 4-velocity (contravariant) call get_u0(gcov,v,uzero,ierr) @@ -380,34 +318,14 @@ subroutine get_tmunu(x,metrici,v,dens,u,p,tmunu,verbose) enddo enddo - -! if (present(verbose) .and. verbose) then -! ! Do we get sensible values -! print*, "Unpacked metric quantities..." -! print*, "gcov: ", gcov -! print*, "gcon: ", gcon -! print*, "gammaijdown: ", gammaijdown -! print* , "alpha: ", alpha -! print*, "betadown: ", betadown -! print*, "v4: ", v4 -! endif - -! if (verbose) then -! print*, "tmunu part: ", tmunu -! print*, "dens: ", dens -! print*, "w: ", w -! print*, "p: ", p -! print*, "gcov: ", gcov -! endif - - ! print*, "tmunu part: ", tmunu - ! print*, "dens: ", dens - ! print*, "w: ", w - ! print*, "p: ", p - ! print*, "gcov: ", gcov - ! stop end subroutine get_tmunu +!------------------------------------------------------------------------- +!+ +! the following two routines are for testing purposes +! and could be deleted at some stage (as used in Magnall et al. 2023) +!+ +!------------------------------------------------------------------------- subroutine get_tmunu_exact(x,metrici,metricderivsi,v,dens,u,p,tmunu) use metric_tools, only:unpack_metric use utils_gr, only:get_sqrtg @@ -443,10 +361,10 @@ subroutine get_tmunu_exact(x,metrici,metricderivsi,v,dens,u,p,tmunu) rhostar = 13.294563008157013D0 call get_sqrtg(gcov,negsqrtg) + ! Set/Calculate primitive density using rhostar exactly rhoprim = rhostar/(negsqrtg/alpha) - ! Stress energy tensor do j=0,3 do i=0,3 @@ -454,8 +372,40 @@ subroutine get_tmunu_exact(x,metrici,metricderivsi,v,dens,u,p,tmunu) enddo enddo +end subroutine get_tmunu_exact +!------------------------------------------------------------------------- +!+ +! see above: for testing purposes and could be deleted at some stage +! (as used in Magnall et al. 2023) +!+ +!------------------------------------------------------------------------- +subroutine get_tmunu_all_exact(npart,xyzh,metrics,vxyzu,metricderivs,dens,tmunus) + use eos, only:ieos,get_pressure + use part, only:isdead_or_accreted + integer, intent(in) :: npart + real, intent(in) :: xyzh(:,:), metrics(:,:,:,:), metricderivs(:,:,:,:), dens(:) + real, intent(inout) :: vxyzu(:,:),tmunus(:,:,:) + real :: pi + integer :: i + logical :: firstpart + real :: tmunu(4,4) -end subroutine get_tmunu_exact + tmunu = 0. + firstpart = .true. + ! TODO write openmp parallel code + do i=1, npart + if (.not.isdead_or_accreted(xyzh(4,i)) .and. firstpart) then + pi = get_pressure(ieos,xyzh(:,i),dens(i),vxyzu(:,i)) + call get_tmunu_exact(xyzh(:,i),metrics(:,:,:,i), metricderivs(:,:,:,i), & + vxyzu(1:3,i),dens(i),vxyzu(4,i),pi,tmunus(:,:,i)) + firstpart = .false. + tmunu(:,:) = tmunus(:,:,i) + else + tmunus(:,:,i) = tmunu(:,:) + endif + enddo + +end subroutine get_tmunu_all_exact end module extern_gr diff --git a/src/main/force.F90 b/src/main/force.F90 index 7a9cb72f2..7586fc82a 100644 --- a/src/main/force.F90 +++ b/src/main/force.F90 @@ -39,11 +39,11 @@ module forces ! ! :Runtime parameters: None ! -! :Dependencies: boundary, cooling, dim, dust, eos, eos_shen, fastmath, io, -! io_summary, kdtree, kernel, linklist, metric_tools, mpiderivs, -! mpiforce, mpimemory, mpiutils, nicil, omputils, options, part, physcon, -! ptmass, ptmass_heating, radiation_utils, timestep, timestep_ind, -! timestep_sts, timing, units, utils_gr, viscosity +! :Dependencies: boundary, cooling, dim, dust, eos, eos_shen, fastmath, +! growth, io, io_summary, kdtree, kernel, linklist, metric_tools, +! mpiderivs, mpiforce, mpimemory, mpiutils, nicil, omputils, options, +! part, physcon, ptmass, ptmass_heating, radiation_utils, timestep, +! timestep_ind, timestep_sts, timing, units, utils_gr, viscosity ! use dim, only:maxfsum,maxxpartveciforce,maxp,ndivcurlB,ndivcurlv,& maxdusttypes,maxdustsmall,do_radiation @@ -104,22 +104,23 @@ module forces icurlBxi = 42, & icurlByi = 43, & icurlBzi = 44, & - igrainsizei = 45, & + igrainmassi = 45, & igraindensi = 46, & - ifxi_drag = 47, & - ifyi_drag = 48, & - ifzi_drag = 49, & - idti = 50, & - idvxdxi = 51, & - idvzdzi = 59, & + ifilfaci = 47, & + ifxi_drag = 48, & + ifyi_drag = 49, & + ifzi_drag = 50, & + idti = 51, & + idvxdxi = 52, & + idvzdzi = 60, & !--dust arrays initial index - idustfraci = 60, & + idustfraci = 61, & !--dust arrays final index - idustfraciend = 60 + (maxdusttypes - 1), & - itstop = 61 + (maxdusttypes - 1), & - itstopend = 61 + 2*(maxdusttypes - 1), & + idustfraciend = 61 + (maxdusttypes - 1), & + itstop = 62 + (maxdusttypes - 1), & + itstopend = 62 + 2*(maxdusttypes - 1), & !--final dust index - lastxpvdust = 61 + 2*(maxdusttypes - 1), & + lastxpvdust = 62 + 2*(maxdusttypes - 1), & iradxii = lastxpvdust + 1, & iradfxi = lastxpvdust + 2, & iradfyi = lastxpvdust + 3, & @@ -184,13 +185,13 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,& rad,drad,radprop,dustprop,dustgasprop,dustfrac,ddustevol,fext,fxyz_drag,& ipart_rhomax,dt,stressmax,eos_vars,dens,metrics) - use dim, only:maxvxyzu,maxneigh,mhd,mhd_nonideal,lightcurve,mpi + use dim, only:maxvxyzu,maxneigh,mhd,mhd_nonideal,lightcurve,mpi,use_dust use io, only:iprint,fatal,iverbose,id,master,real4,warning,error,nprocs use linklist, only:ncells,get_neighbour_list,get_hmaxcell,get_cell_location,listneigh use options, only:iresistive_heating use part, only:rhoh,dhdrho,rhoanddhdrho,alphaind,iactive,gradh,& - hrho,iphase,igas,maxgradh,dvdx,eta_nimhd,deltav,poten,iamtype,use_dust,& - fxyz_dragold + hrho,iphase,igas,maxgradh,dvdx,eta_nimhd,deltav,poten,iamtype,& + dragreg,filfac,fxyz_dragold use timestep, only:dtcourant,dtforce,dtrad,bignumber,dtdiff use io_summary, only:summary_variable, & iosumdtf,iosumdtd,iosumdtv,iosumdtc,iosumdto,iosumdth,iosumdta, & @@ -397,6 +398,8 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,& !$omp shared(ncells,ifirstincell) & !$omp shared(xyzh) & !$omp shared(dustprop) & +!$omp shared(dragreg) & +!$omp shared(filfac) & !$omp shared(dustgasprop) & !$omp shared(fxyz_drag) & !$omp shared(fxyz_dragold) & @@ -530,8 +533,8 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,& call write_cell(stack_waiting,cell) else call finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dvdx,& - divBsymm,divcurlv,dBevol,ddustevol,deltav,dustgasprop,fxyz_drag,fext, & - dtcourant,dtforce,dtvisc,dtohm,dthall,dtambi,dtdiff,dtmini,dtmaxi, & + divBsymm,divcurlv,dBevol,ddustevol,deltav,dustgasprop,fxyz_drag,fext,dragreg,& + filfac,dtcourant,dtforce,dtvisc,dtohm,dthall,dtambi,dtdiff,dtmini,dtmaxi, & #ifdef IND_TIMESTEPS nbinmaxnew,nbinmaxstsnew,ncheckbin, & ndtforce,ndtforceng,ndtcool,ndtdrag,ndtdragd, & @@ -621,8 +624,8 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,& cell = get_cell(stack_waiting,i) call finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dvdx, & - divBsymm,divcurlv,dBevol,ddustevol,deltav,dustgasprop,fxyz_drag,fext, & - dtcourant,dtforce,dtvisc,dtohm,dthall,dtambi,dtdiff,dtmini,dtmaxi, & + divBsymm,divcurlv,dBevol,ddustevol,deltav,dustgasprop,fxyz_drag,fext,dragreg, & + filfac,dtcourant,dtforce,dtvisc,dtohm,dthall,dtambi,dtdiff,dtmini,dtmaxi, & #ifdef IND_TIMESTEPS nbinmaxnew,nbinmaxstsnew,ncheckbin, & ndtforce,ndtforceng,ndtcool,ndtdrag,ndtdragd, & @@ -910,7 +913,9 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g use dim, only:use_dust,use_dustgrowth,ind_timesteps use dust, only:get_ts,idrag,icut_backreaction,ilimitdustflux,irecon,drag_implicit use kernel, only:wkern_drag,cnormk_drag,wkern,cnormk - use part, only:ndustsmall,grainsize,graindens + use part, only:ndustsmall,grainsize,graindens,ndustsmall,grainsize,graindens,filfac + use options, only:use_porosity + use growth, only:get_size use kernel, only:wkern,cnormk #ifdef IND_TIMESTEPS use part, only:ibin_old,iamboundary @@ -1001,7 +1006,7 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g real :: pri,pro2i real :: etaohmi,etahalli,etaambii real :: jcbcbi(3),jcbi(3) - real :: alphai,grainsizei,graindensi + real :: alphai,grainmassi,graindensi,filfaci logical :: usej integer :: iamtypei real :: radFi(3),radFj(3),radRj,radDFWi,radDFWj,c_code,radkappai,radkappaj,& @@ -1071,8 +1076,9 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g endif if (use_dustgrowth) then - grainsizei = xpartveci(igrainsizei) + grainmassi = xpartveci(igrainmassi) graindensi = xpartveci(igraindensi) + filfaci = xpartveci(ifilfaci) endif fxi_drag = xpartveci(ifxi_drag) fyi_drag = xpartveci(ifyi_drag) @@ -1741,7 +1747,13 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g do l=1,ndustsmall ! get stopping time - for one fluid dust we do not know deltav, but it is small by definition if (use_dustgrowth) then !- only work for ndustsmall=1 though - call get_ts(idrag,l,grainsizei,graindensi,rhogasj,rhoj*dustfracjsum,spsoundj,0.,tsj(l),iregime) + if (use_porosity) then + call get_ts(idrag,l,get_size(grainmassi,graindensi,filfaci),& + graindensi*filfaci,rhogasj,rhoj*dustfracjsum,spsoundj,0.,tsj(l),iregime) + else + call get_ts(idrag,l,get_size(grainmassi,graindensi),& + graindensi,rhogasj,rhoj*dustfracjsum,spsoundj,0.,tsj(l),iregime) + endif else call get_ts(idrag,l,grainsize(l),graindens(l),rhogasj,rhoj*dustfracjsum,spsoundj,0.,tsj(l),iregime) endif @@ -1807,7 +1819,13 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g wdrag = wkern_drag(q2j,qj)*hj21*hj1*cnormk_drag endif if (use_dustgrowth) then - call get_ts(idrag,1,dustprop(1,j),dustprop(2,j),rhoi,rhoj,spsoundi,dv2,tsijtmp,iregime) + if (use_porosity) then + call get_ts(idrag,1,get_size(dustprop(1,j),dustprop(2,j),filfac(j)),& + dustprop(2,j)*filfac(j),rhoi,rhoj,spsoundi,dv2,tsijtmp,iregime) + else + call get_ts(idrag,1,get_size(dustprop(1,j),dustprop(2,j)),& + dustprop(2,j),rhoi,rhoj,spsoundi,dv2,tsijtmp,iregime) + endif else !--the following works for large grains only (not hybrid large and small grains) idusttype = iamtypej - idust + 1 @@ -1852,7 +1870,12 @@ subroutine compute_forces(i,iamgasi,iamdusti,xpartveci,hi,hi1,hi21,hi41,gradhi,g wdrag = wkern_drag(q2j,qj)*hj21*hj1*cnormk_drag endif if (use_dustgrowth) then - call get_ts(idrag,1,grainsizei,graindensi,rhoj,rhoi,spsoundj,dv2,tsijtmp,iregime) + if (use_porosity) then + call get_ts(idrag,1,get_size(grainmassi,graindensi,filfaci),& + graindensi*filfaci,rhoj,rhoi,spsoundj,dv2,tsijtmp,iregime) + else + call get_ts(idrag,1,get_size(grainmassi,graindensi),graindensi,rhoj,rhoi,spsoundj,dv2,tsijtmp,iregime) + endif if (q2i < q2j) then winter = wkern(q2i,qi)*hi21*hi1*cnormk else @@ -2065,8 +2088,10 @@ subroutine start_cell(cell,iphase,xyzh,vxyzu,gradh,divcurlv,divcurlB,dvdx,Bevol, iohm,ihall,iambi,ndustsmall,iradP,igasP,ics,itemp use viscosity, only:irealvisc,bulkvisc use dust, only:get_ts,idrag - use part, only:grainsize,graindens - use part, only:ibin_old + use options, only:use_porosity + use part, only:grainsize,graindens,filfac + use growth, only:get_size + use part, only:ibin_old use timestep_ind, only:get_dt use nicil, only:nimhd_get_jcbcb use radiation_utils, only:get_rad_R @@ -2211,7 +2236,13 @@ subroutine start_cell(cell,iphase,xyzh,vxyzu,gradh,divcurlv,divcurlB,dvdx,Bevol, if (use_dust .and. use_dustfrac .and. iamgasi) then do j=1,ndustsmall if (use_dustgrowth) then - call get_ts(idrag,j,dustprop(1,i),dustprop(2,i),rhogasi,rhoi*dustfracisum,spsoundi,0.,tstopi(j),iregime) + if (use_porosity) then + call get_ts(idrag,j,get_size(dustprop(1,i),dustprop(2,i),filfac(j)),& + dustprop(2,i)*filfac(j),rhogasi,rhoi*dustfracisum,spsoundi,0.,tstopi(j),iregime) + else + call get_ts(idrag,j,get_size(dustprop(1,i),dustprop(2,i)),& + dustprop(2,i),rhogasi,rhoi*dustfracisum,spsoundi,0.,tstopi(j),iregime) + endif else call get_ts(idrag,j,grainsize(j),graindens(j),rhogasi,rhoi*dustfracisum,spsoundi,0.,tstopi(j),iregime) endif @@ -2343,8 +2374,9 @@ subroutine start_cell(cell,iphase,xyzh,vxyzu,gradh,divcurlv,divcurlB,dvdx,Bevol, endif if (use_dustgrowth) then - cell%xpartvec(igrainsizei,cell%npcell) = dustprop(1,i) + cell%xpartvec(igrainmassi,cell%npcell) = dustprop(1,i) cell%xpartvec(igraindensi,cell%npcell) = dustprop(2,i) + cell%xpartvec(ifilfaci,cell%npcell) = filfac(i) endif if (use_dust) then cell%xpartvec(ifxi_drag,cell%npcell) = fxyz_drag(1,i) @@ -2469,8 +2501,8 @@ subroutine compute_cell(cell,listneigh,nneigh,Bevol,xyzh,vxyzu,fxyzu, & end subroutine compute_cell subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dvdx,& - divBsymm,divcurlv,dBevol,ddustevol,deltav,dustgasprop,fxyz_drag,fext, & - dtcourant,dtforce,dtvisc,dtohm,dthall,dtambi,dtdiff,dtmini,dtmaxi, & + divBsymm,divcurlv,dBevol,ddustevol,deltav,dustgasprop,fxyz_drag,fext,dragreg, & + filfac,dtcourant,dtforce,dtvisc,dtohm,dthall,dtambi,dtdiff,dtmini,dtmaxi, & #ifdef IND_TIMESTEPS nbinmaxnew,nbinmaxstsnew,ncheckbin, & ndtforce,ndtforceng,ndtcool,ndtdrag,ndtdragd, & @@ -2516,7 +2548,10 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv use metric_tools, only:unpack_metric use utils_gr, only:get_u0 use io, only:error + use growth, only:get_size use dust, only:idrag,get_ts + use physcon, only:fourpi + use options, only:use_porosity use part, only:Omega_k use io, only:warning use physcon, only:c,kboltz @@ -2534,6 +2569,8 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv real, intent(out) :: ddustevol(:,:) real, intent(out) :: deltav(:,:,:) real, intent(out) :: dustgasprop(:,:) + real, intent(in) :: filfac(:) + integer, intent(out) :: dragreg(:) real, intent(inout) :: fxyz_drag(:,:) real, intent(in) :: fext(:,:) real, intent(inout) :: dtcourant,dtforce,dtvisc @@ -2581,7 +2618,7 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv logical :: allow_decrease,dtcheck character(len=16) :: dtchar #endif - real :: tstopint,gsizei,gdensi + real :: tstopint,gmassi,gdensi integer :: ireg integer :: ip,i real :: densi, vxi,vyi,vzi,u0i,dudtcool,dudtheat @@ -3004,9 +3041,15 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv !- get the Stokes number with get_ts using the interpolated quantities rhoi = xpartveci(irhoi) gdensi = xpartveci(igraindensi) - gsizei = xpartveci(igrainsizei) - call get_ts(idrag,1,gsizei,gdensi,dustgasprop(2,i),rhoi,dustgasprop(1,i),& - dustgasprop(4,i)**2,tstopint,ireg) + gmassi = xpartveci(igrainmassi) + if (use_porosity) then + call get_ts(idrag,1,get_size(gmassi,gdensi,filfac(i)),gdensi*filfac(i),& + dustgasprop(2,i),rhoi,dustgasprop(1,i),dustgasprop(4,i)**2,tstopint,ireg) + dragreg(i) = ireg + else + call get_ts(idrag,1,get_size(gmassi,gdensi),gdensi,& + dustgasprop(2,i),rhoi,dustgasprop(1,i),dustgasprop(4,i)**2,tstopint,ireg) + endif dustgasprop(3,i) = tstopint * Omega_k(i) !- Stokes number endif diff --git a/src/main/growth.F90 b/src/main/growth.f90 similarity index 76% rename from src/main/growth.F90 rename to src/main/growth.f90 index a8ad0f4a0..a2d35aec7 100644 --- a/src/main/growth.F90 +++ b/src/main/growth.f90 @@ -28,6 +28,7 @@ module growth ! - isnow : *snow line (0=off,1=position based,2=temperature based)* ! - rsnow : *snow line position in AU* ! - size_max_user : *(mcfost) maximum size for binning in cm* +! - tsmincgs : *minimum allowed stopping time* ! - vfrag : *uniform fragmentation threshold in m/s* ! - vfragin : *inward fragmentation threshold in m/s* ! - vfragout : *inward fragmentation threshold in m/s* @@ -47,6 +48,7 @@ module growth integer, public :: ieros = 0 real, public :: gsizemincgs = 5.e-3 + real, public :: tsmincgs = 1.e5 real, public :: rsnow = 100. real, public :: Tsnow = 150. real, public :: vfragSI = 15. @@ -60,6 +62,7 @@ module growth real, public :: vfragin real, public :: vfragout real, public :: grainsizemin + real, public :: tsmin real, public :: cohacc real, public :: dsize @@ -70,7 +73,7 @@ module growth public :: get_growth_rate,get_vrelonvfrag,check_dustprop public :: write_options_growth,read_options_growth,print_growthinfo,init_growth public :: vrelative,read_growth_setup_options,write_growth_setup_options - public :: comp_snow_line,bin_to_multi,convert_to_twofluid + public :: comp_snow_line,bin_to_multi,convert_to_twofluid,get_size contains @@ -81,7 +84,10 @@ module growth !------------------------------------------------ subroutine init_growth(ierr) use io, only:error + use physcon, only:fourpi use viscosity, only:irealvisc,shearparam + use dust, only:grainsizecgs + use options, only:use_porosity integer, intent(out) :: ierr ierr = 0 @@ -95,12 +101,17 @@ subroutine init_growth(ierr) grainsizemin = gsizemincgs / udist cohacc = cohacccgs * utime * utime / umass dsize = dsizecgs / udist + tsmin = tsmincgs / utime if (ifrag > 0) then if (grainsizemin < 0.) then call error('init_growth','grainsizemin < 0',var='grainsizemin',val=grainsizemin) ierr = 1 endif + if (gsizemincgs > grainsizecgs .and. .not. use_porosity) then + call error('init_growth','grainsizemin > grainsize',var='grainsizemin',val=grainsizemin) + ierr = 1 + endif select case(isnow) case(0) !-- uniform vfrag if (vfrag <= 0.) then @@ -163,8 +174,8 @@ subroutine print_growthinfo(iprint) integer, intent(in) :: iprint - if (ifrag == 0) write(iprint,"(a)") ' Using pure growth model where ds = + vrel*rhod/graindens*dt ' - if (ifrag == 1) write(iprint,"(a)") ' Using growth/frag where ds = (+ or -) vrel*rhod/graindens*dt ' + if (ifrag == 0) write(iprint,"(a)") ' Using pure growth model where dm/dt = + 4pi*rhod*s**2*vrel*dt ' + if (ifrag == 1) write(iprint,"(a)") ' Using growth/frag where dm/dt = (+ or -) 4pi*rhod*s**2*vrel*dt ' if (ifrag == 2) write(iprint,"(a)") ' Using growth/frag with Kobayashi fragmentation model ' if (ifrag > -1) write(iprint,"((a,1pg10.3))")' Computing Vrel with alphaSS = ',shearparam if (ifrag > 0) then @@ -185,44 +196,47 @@ subroutine print_growthinfo(iprint) endif endif if (ieros == 1) then - write(iprint,"(a)") ' Using aeolian-erosion model where ds = -rhog*(deltav**3)*(d**2)/(3*cohacc*s)*dt ' + write(iprint,"(a)") ' Using aeolian-erosion model where ds = -fourpi*rhos*rhog*s*(deltav**3)*(dsize**2)/(3*cohacc)*dt ' write(iprint,"(2(a,1pg10.3),a)")' dsize = ',dsizecgs,' cm = ',dsize,' (code units)' endif + end subroutine print_growthinfo !----------------------------------------------------------------------- !+ -! Main routine that returns ds/dt and calculate Vrel/Vfrag. +! Main routine that returns dm/dt and calculate Vrel/Vfrag. ! This growth model is currently only available for the ! two-fluid dust method. !+ !----------------------------------------------------------------------- -subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,dsdt) +subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,filfac,dmdt) use part, only:rhoh,idust,igas,iamtype,iphase,isdead_or_accreted,& massoftype,Omega_k,dustfrac,tstop,deltav - use options, only:use_dustfrac + use options, only:use_dustfrac,use_porosity + use physcon, only:fourpi use eos, only:ieos,get_spsound real, intent(in) :: dustprop(:,:) real, intent(inout) :: dustgasprop(:,:) real, intent(in) :: xyzh(:,:) + real, intent(in) :: filfac(:) real, intent(inout) :: VrelVf(:),vxyzu(:,:) - real, intent(out) :: dsdt(:) + real, intent(out) :: dmdt(:) integer, intent(in) :: npart ! - real :: rhog,rhod,vrel,rho + real :: rhog,rhod,vrel,rho,sdust integer :: i,iam vrel = 0. rhod = 0. rho = 0. - !--get ds/dt over all particles + !--get dm/dt over all particles !$omp parallel do default(none) & - !$omp shared(npart,iphase,ieos,massoftype,use_dustfrac,dustfrac) & + !$omp shared(npart,iphase,ieos,massoftype,use_dustfrac,dustfrac,use_porosity) & !$omp shared(ifrag,ieros,utime,umass,dsize,cohacc) & - !$omp shared(xyzh,vxyzu,dustprop,dustgasprop,dsdt,VrelVf,tstop,deltav) & - !$omp private(i,iam,rho,rhog,rhod,vrel) + !$omp shared(xyzh,vxyzu,dustprop,dustgasprop,dmdt,filfac,VrelVf,tstop,deltav) & + !$omp private(i,iam,rho,rhog,rhod,vrel,sdust) do i=1,npart if (.not.isdead_or_accreted(xyzh(4,i))) then iam = iamtype(iphase(i)) @@ -242,30 +256,40 @@ subroutine get_growth_rate(npart,xyzh,vxyzu,dustgasprop,VrelVf,dustprop,dsdt) rhod = rhoh(xyzh(4,i),massoftype(idust)) endif + !--dust size from mass and filling factor + if (use_porosity) then + sdust = get_size(dustprop(1,i),dustprop(2,i),filfac(i)) + else + sdust = get_size(dustprop(1,i),dustprop(2,i)) + endif + call get_vrelonvfrag(xyzh(:,i),vxyzu(:,i),vrel,VrelVf(i),dustgasprop(:,i)) ! - !--dustprop(1)= size, dustprop(2) = intrinsic density, + !--dustprop(1) = mass, dustprop(2) = intrinsic density, ! - !--if statements to compute ds/dt + !--if statements to compute dm/dt ! - if (ifrag == -1) dsdt(i) = 0. + if (ifrag == -1) dmdt(i) = 0. if ((VrelVf(i) < 1. .or. ifrag == 0) .and. ifrag /= -1) then ! vrel/vfrag < 1 or pure growth --> growth - dsdt(i) = rhod/dustprop(2,i)*vrel + dmdt(i) = fourpi*sdust**2*rhod*vrel elseif (VrelVf(i) >= 1. .and. ifrag > 0) then ! vrel/vfrag > 1 --> fragmentation select case(ifrag) case(1) - dsdt(i) = -rhod/dustprop(2,i)*vrel ! Symmetrical of Stepinski & Valageas + dmdt(i) = -fourpi*sdust**2*rhod*vrel ! Symmetrical of Stepinski & Valageas case(2) - dsdt(i) = -rhod/dustprop(2,i)*vrel*(VrelVf(i)**2)/(1+VrelVf(i)**2) ! Kobayashi model + dmdt(i) = -fourpi*sdust**2*rhod*vrel*(VrelVf(i)**2)/(1+VrelVf(i)**2) ! Kobayashi model end select - endif !sqrt(0.0123)=0.110905 !1.65 -> surface energy in cgs - if (ieros == 1 .and. (dustgasprop(4,i) >= 0.110905*utime*sqrt(1.65/umass/dustgasprop(2,i)/dsize))) then - dsdt(i) = dsdt(i) - dustgasprop(2,i)*(dustgasprop(4,i)**3)*(dsize**2)/(3.*cohacc*dustprop(1,i)) ! Erosion model + endif + if (ieros == 1) then !sqrt(0.0123)=0.110905 !1.65 -> surface energy in cgs + ! Erosion model of Rozner, Grishin & Perets (2020) + if (dustgasprop(4,i) >= 0.110905*sqrt(1.65*utime*utime/umass/dustprop(2,i)/dsize)) then + dmdt(i) = dmdt(i) - fourpi*sdust*dustprop(2,i)*dustgasprop(2,i)*(dustgasprop(4,i)**3)*(dsize**2)/(3.*cohacc) + endif endif endif else - dsdt(i) = 0. + dmdt(i) = 0. endif enddo !$omp end parallel do @@ -350,6 +374,7 @@ end subroutine comp_snow_line !----------------------------------------------------------------------- subroutine write_options_growth(iunit) use infile_utils, only:write_inopt + use options, only:use_porosity integer, intent(in) :: iunit write(iunit,"(/,a)") '# options controlling growth' @@ -357,7 +382,11 @@ subroutine write_options_growth(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) if (ifrag /= 0) then - call write_inopt(gsizemincgs,'grainsizemin','minimum grain size in cm',iunit) + if (use_porosity) then + call write_inopt(tsmincgs,'tsmincgs','minimum allowed stopping time',iunit) + else + call write_inopt(gsizemincgs,'grainsizemin','minimum grain size in cm',iunit) + endif call write_inopt(isnow,'isnow','snow line (0=off,1=position based,2=temperature based)',iunit) if (isnow == 1) call write_inopt(rsnow,'rsnow','position of the snow line in AU',iunit) if (isnow == 2) call write_inopt(Tsnow,'Tsnow','snow line condensation temperature in K',iunit) @@ -408,6 +437,9 @@ subroutine read_options_growth(name,valstring,imatch,igotall,ierr) case('grainsizemin') read(valstring,*,iostat=ierr) gsizemincgs ngot = ngot + 1 + case('tsmincgs') + read(valstring,*,iostat=ierr) tsmincgs + ngot = ngot + 1 case('isnow') read(valstring,*,iostat=ierr) isnow ngot = ngot + 1 @@ -482,6 +514,7 @@ end subroutine read_options_growth !----------------------------------------------------------------------- subroutine write_growth_setup_options(iunit) use infile_utils, only:write_inopt + use options, only:use_porosity integer, intent(in) :: iunit write(iunit,"(/,a)") '# options for growth and fragmentation of dust' @@ -494,7 +527,11 @@ subroutine write_growth_setup_options(iunit) call write_inopt(vfragSI,'vfrag','uniform fragmentation threshold in m/s',iunit) call write_inopt(vfraginSI,'vfragin','inward fragmentation threshold in m/s',iunit) call write_inopt(vfragoutSI,'vfragout','inward fragmentation threshold in m/s',iunit) - call write_inopt(gsizemincgs,'grainsizemin','minimum allowed grain size in cm',iunit) + if (use_porosity) then + call write_inopt(tsmincgs,'tsmincgs','minimum allowed stopping time',iunit) + else + call write_inopt(gsizemincgs,'grainsizemin','minimum allowed grain size in cm',iunit) + endif end subroutine write_growth_setup_options @@ -505,6 +542,7 @@ end subroutine write_growth_setup_options !----------------------------------------------------------------------- subroutine read_growth_setup_options(db,nerr) use infile_utils, only:read_inopt,inopts + use options, only:use_porosity type(inopts), allocatable, intent(inout) :: db(:) integer, intent(inout) :: nerr @@ -512,7 +550,11 @@ subroutine read_growth_setup_options(db,nerr) call read_inopt(ieros,'ieros',db,min=0,max=1,errcount=nerr) if (ifrag > 0) then call read_inopt(isnow,'isnow',db,min=0,max=2,errcount=nerr) - call read_inopt(gsizemincgs,'grainsizemin',db,min=1.e-5,errcount=nerr) + if (use_porosity) then + call read_inopt(tsmincgs,'tsmincgs',db,min=1.e-5,errcount=nerr) + else + call read_inopt(gsizemincgs,'grainsizemin',db,min=1.e-5,errcount=nerr) + endif select case(isnow) case(0) call read_inopt(vfragSI,'vfrag',db,min=0.,errcount=nerr) @@ -531,22 +573,43 @@ end subroutine read_growth_setup_options !----------------------------------------------------------------------- !+ -! In case of fragmentation, limit sizes to a minimum value +! In case of fragmentation, limit masses to a minimum value !+ !----------------------------------------------------------------------- -subroutine check_dustprop(npart,size) - use part, only:iamtype,iphase,idust,igas - use options, only:use_dustfrac - real,intent(inout) :: size(:) +subroutine check_dustprop(npart,dustprop,filfac,mprev,filfacprev) + use part, only:iamtype,iphase,idust,igas,dustgasprop,Omega_k + use options, only:use_dustfrac,use_porosity + real,intent(inout) :: dustprop(:,:) integer,intent(in) :: npart + real, intent(in) :: filfac(:),mprev(:),filfacprev(:) integer :: i,iam + real :: tsnew,sdustprev,sdustmin,sdust + !$omp parallel do default(none) & + !$omp shared(iphase,dustgasprop,use_dustfrac,use_porosity) & + !$omp shared(npart,ifrag,dustprop,filfac,mprev,filfacprev) & + !$omp shared(tsmin,grainsizemin) & + !$omp private(i,iam,tsnew,sdustprev,sdustmin,sdust) do i=1,npart iam = iamtype(iphase(i)) - if (iam==idust .or. (use_dustfrac .and. iam==igas)) then - if (ifrag > 0 .and. size(i) < grainsizemin) size(i) = grainsizemin + if ((iam == idust .or. (use_dustfrac .and. iam == igas)) .and. ifrag > 0 .and. dustprop(1,i) <= mprev(i)) then + if (use_porosity) then + sdustprev = get_size(mprev(i),dustprop(2,i),filfacprev(i)) + sdust = get_size(dustprop(1,i),dustprop(2,i),filfac(i)) + tsnew = dustgasprop(3,i)*sdustprev*filfacprev(i)/sdust/filfac(i)/Omega_k(i) + if (tsnew < tsmin) then + sdustmin = tsmin*sdustprev*filfacprev(i)*Omega_k(i)/filfac(i)/dustgasprop(3,i) + dustprop(1,i) = dustprop(1,i) * (sdustmin/sdust)**3. + endif + else + sdust = get_size(dustprop(1,i),dustprop(2,i)) + if (sdust < grainsizemin) then + dustprop(1,i) = dustprop(1,i) * (grainsizemin/sdust)**3. ! fragmentation at constant density and filling factor + endif + endif endif enddo + !$omp end parallel do end subroutine check_dustprop @@ -555,15 +618,38 @@ end subroutine check_dustprop ! Set dustprop (used by moddump) !+ !----------------------------------------------------------------------- -subroutine set_dustprop(npart) - use dust, only:grainsizecgs,graindenscgs - use part, only:dustprop - integer,intent(in) :: npart - integer :: i +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 physcon, only:fourpi + use options, only:use_dustfrac + integer, intent(in) :: npart + real, intent(in) :: xyzh(:,:) + integer :: i,iam + real :: r,h + logical, optional, intent(in) :: sizedistrib + real, optional, intent(in) :: pwl_sizedistrib,R_ref,H_R_ref,q_index do i=1,npart - dustprop(1,i) = grainsizecgs / udist - dustprop(2,i) = graindenscgs / unit_density + 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 + else + dustprop(1,i) = fourpi/3. * dustprop(2,i) * (grainsizecgs / udist)**3 + endif + else + dustprop(:,i) = 0. + endif + filfac(i) = 0. + probastick(i) = 1. enddo end subroutine set_dustprop @@ -575,8 +661,9 @@ end subroutine set_dustprop !----------------------------------------------------------------------- subroutine bin_to_multi(bins_per_dex,force_smax,smax_user,verbose) use part, only:npart,npartoftype,massoftype,ndusttypes,& - ndustlarge,grainsize,dustprop,graindens,& - iamtype,iphase,set_particle_type,idust + ndustlarge,grainsize,dustprop,& + iamtype,iphase,set_particle_type,idust,filfac + use options, only:use_porosity use units, only:udist use table_utils, only:logspace use io, only:fatal @@ -585,38 +672,51 @@ subroutine bin_to_multi(bins_per_dex,force_smax,smax_user,verbose) real, intent(in) :: smax_user logical, intent(inout) :: force_smax logical, intent(in) :: verbose - real :: smaxtmp,smintmp,smax,smin,tolm,& - mdustold,mdustnew,code_to_mum + real :: smaxtmp,smintmp,smax,smin,tolm + real :: mdustold,mdustnew,code_to_mum logical :: init - integer :: nbins,nbinmax,i,j,itype,ndustold,ndustnew,npartmin,imerge,iu - real, allocatable, dimension(:) :: grid - character(len=20) :: outfile = "bin_distrib.dat" + integer :: nbinsize,nbinsizemax,i,j,itype,ndustold,ndustnew,npartmin,imerge,iu + integer :: nbinfilfacmax,ndustsizetypes + real, allocatable, dimension(: ) :: grid + real, allocatable, dimension(:,:) :: dustpropmcfost !dustpropmcfost(1=size,2=filfac) + character(len=20) :: outfile = "bin_distrib.dat" !- initialise - code_to_mum = udist*1.e4 - tolm = 1.e-5 - smaxtmp = 0. - smintmp = 1.e26 - ndustold = 0 - ndustnew = 0 - mdustold = 0. - mdustnew = 0. - nbinmax = 25 - npartmin = 50 !- limit to find neighbours - init = .false. - graindens = maxval(dustprop(2,:)) + code_to_mum = udist*1.e4 + tolm = 1.e-5 + smaxtmp = 0. + smintmp = 1.e26 + ndustold = 0 + ndustnew = 0 + mdustold = 0. + mdustnew = 0. + nbinsizemax = 25 + nbinfilfacmax = 10 + npartmin = 50 !- limit to find neighbours + init = .false. + allocate(dustpropmcfost(2,npart)) + !graindens = maxval(dustprop(2,:)) !- loop over particles, find min and max on non-accreted dust particles do i = 1,npart itype = iamtype(iphase(i)) if (itype==idust) then - if (dustprop(1,i) < smintmp) smintmp = dustprop(1,i) - if (dustprop(1,i) > smaxtmp) smaxtmp = dustprop(1,i) + if (use_porosity) then + dustpropmcfost(1,i) = get_size(dustprop(1,i),dustprop(2,i),filfac(i)) + dustpropmcfost(2,i) = filfac(i) + else + dustpropmcfost(1,i) = get_size(dustprop(1,i),dustprop(2,i)) + dustpropmcfost(2,i) = 1 + endif + if (dustpropmcfost(1,i) < smintmp) smintmp = dustpropmcfost(1,i) + if (dustpropmcfost(1,i) > smaxtmp) smaxtmp = dustpropmcfost(1,i) + !if (dustpropmcfost(2,i) < fmintmp) fmintmp = dustpropmcfost(2,i) + !if (dustpropmcfost(2,i) > fmaxtmp) fmaxtmp = dustpropmcfost(2,i) endif enddo !- overrule force_smax if particles are small, avoid empty bins - if ((maxval(dustprop(1,:))*udist < smax_user) .and. force_smax) then + if ((maxval(dustpropmcfost(1,1:npart))*udist < smax_user) .and. force_smax) then force_smax = .false. write(*,*) "Overruled force_smax from T to F" endif @@ -628,28 +728,31 @@ subroutine bin_to_multi(bins_per_dex,force_smax,smax_user,verbose) smax = smaxtmp else init = .true. - write(*,*) "Detected initial condition, restraining nbins = 1" + write(*,*) "Detected initial condition, restraining nbinsize = 1" endif if (.not. init) then smin = smintmp !- set ndusttypes based on desired log size spacing - nbins = int((log10(smax)-log10(smin))*bins_per_dex + 1.) - ndusttypes = min(nbins, nbinmax) !- prevent memory allocation errors - ndustlarge = ndusttypes !- this is written to the header + nbinsize = int((log10(smax)-log10(smin))*bins_per_dex + 1.) + !nbinfilfac = int((log10(smax)-log10(smin))*bins_per_dex + 1.) + ndustsizetypes = min(nbinsize, nbinsizemax) !- prevent memory allocation errors + !ndustfilfactypes = min(nbinfilfac,nbinfilfacmax) + ndustlarge = ndustsizetypes!*ndustfilfactypes !- this is written to the header !- allocate memory for a grid of ndusttypes+1 elements - allocate(grid(ndusttypes+1)) + allocate(grid(ndustsizetypes+1))!,ndustfilfactypes)) !- bin sizes in ndusttypes bins write(*,"(a,f10.1,a,f10.1,a,i3,a)") "Binning sizes between ",smin*code_to_mum, " (µm) and ",& smax*code_to_mum," (µm) in ",ndusttypes, " bins" - call logspace(grid(1:ndusttypes+1),smin,smax) !- bad for live mcfost, need to compile it before growth.F90 + call logspace(grid(1:ndustsizetypes+1),smin,smax) !- bad for live mcfost, need to compile it before growth.F90 + !call logspace(grid(2,1:ndustfilfactypes),fmin,fmax) !- find representative size for each bin - do i = 1,ndusttypes + do i = 1,ndustsizetypes grainsize(i) = sqrt(grid(i)*grid(i+1)) enddo @@ -661,7 +764,7 @@ subroutine bin_to_multi(bins_per_dex,force_smax,smax_user,verbose) if (itype==idust) then !- figure out which bin do j=1,ndusttypes - if ((dustprop(1,i) >= grid(j)) .and. (dustprop(1,i) < grid(j+1))) then + if ((dustpropmcfost(1,i) >= grid(j)) .and. (dustpropmcfost(1,i) < grid(j+1))) then if (j > 1) then npartoftype(idust+j-1) = npartoftype(idust+j-1) + 1 npartoftype(idust) = npartoftype(idust) - 1 @@ -669,7 +772,7 @@ subroutine bin_to_multi(bins_per_dex,force_smax,smax_user,verbose) endif endif !- if smax has been forced, put larger grains inside last bin - if ((j==ndusttypes) .and. force_smax .and. (dustprop(1,i) >= grid(j+1))) then + if ((j==ndusttypes) .and. force_smax .and. (dustpropmcfost(1,i) >= grid(j+1))) then npartoftype(idust+j-1) = npartoftype(idust+j-1) + 1 npartoftype(idust) = npartoftype(idust) - 1 call set_particle_type(i,idust+j-1) @@ -720,6 +823,10 @@ subroutine bin_to_multi(bins_per_dex,force_smax,smax_user,verbose) grainsize(ndusttypes) = smaxtmp !- only 1 bin, all particles have same size endif + ! clean up + if (allocated(dustpropmcfost)) deallocate(dustpropmcfost) + if (allocated(grid)) deallocate(grid) + end subroutine bin_to_multi !----------------------------------------------------------------------- @@ -821,8 +928,8 @@ end subroutine merge_bins !----------------------------------------------------------------------- subroutine convert_to_twofluid(npart,xyzh,vxyzu,massoftype,npartoftype,np_ratio,dust_to_gas) use part, only: dustprop,dustgasprop,ndustlarge,ndustsmall,igas,idust,VrelVf,& - dustfrac,iamtype,iphase,deltav,set_particle_type - use options, only: use_dustfrac + dustfrac,iamtype,iphase,deltav,set_particle_type,filfac + use options, only: use_dustfrac,use_porosity use dim, only: update_max_sizes integer, intent(inout) :: npart,npartoftype(:) real, intent(inout) :: xyzh(:,:),vxyzu(:,:),massoftype(:) @@ -861,6 +968,9 @@ subroutine convert_to_twofluid(npart,xyzh,vxyzu,massoftype,npartoftype,np_ratio, dustgasprop(3,ipart) = dustgasprop(3,iloc) dustgasprop(4,ipart) = dustgasprop(4,iloc) VrelVf(ipart) = VrelVf(iloc) + if (use_porosity) then + filfac(ipart) = filfac(iloc) + endif call set_particle_type(ipart,idust) enddo @@ -875,6 +985,19 @@ subroutine convert_to_twofluid(npart,xyzh,vxyzu,massoftype,npartoftype,np_ratio, vxyzu(1,j) = vxyzu(1,j) - dustfrac(1,j) * deltav(1,1,j) vxyzu(2,j) = vxyzu(2,j) - dustfrac(1,j) * deltav(2,1,j) vxyzu(3,j) = vxyzu(3,j) - dustfrac(1,j) * deltav(3,1,j) + + !- reset dust quantities of the mixture + dustprop(1,j) = 0. + dustprop(2,j) = 0. + dustgasprop(1,j) = 0. + dustgasprop(2,j) = 0. + dustgasprop(3,j) = 0. + dustgasprop(4,j) = 0. + VrelVf(j) = 0. + if (use_porosity) then + filfac(j) = 0. + endif + endif enddo @@ -909,4 +1032,25 @@ real function vrelative(dustgasprop,Vt) end function vrelative +!--Compute size from mass and filling factor +real function get_size(mass,dens,filfac) + use physcon, only:fourpi + real, intent(in) :: mass,dens + real, optional, intent(in) :: filfac + real :: f + + if (present(filfac)) then + f = filfac + else + f = 1.0 + endif + + if (dens > 0. .and. f > 0.) then + get_size = ( 3.*mass / (fourpi*dens*f) )**(1./3.) + else + get_size = 0. + endif + +end function get_size + end module growth diff --git a/src/main/initial.F90 b/src/main/initial.F90 index bbcd5a3d7..314b72644 100644 --- a/src/main/initial.F90 +++ b/src/main/initial.F90 @@ -21,9 +21,9 @@ module initial ! fastmath, fileutils, forcing, growth, inject, io, io_summary, ! krome_interface, linklist, metric_tools, mf_write, mpibalance, ! mpidomain, mpimemory, mpitree, mpiutils, nicil, nicil_sup, omputils, -! options, part, partinject, ptmass, radiation_utils, readwrite_dumps, -! readwrite_infile, timestep, timestep_ind, timestep_sts, timing, -! tmunu2grid, units, writeheader +! options, part, partinject, porosity, ptmass, radiation_utils, +! readwrite_dumps, readwrite_infile, timestep, timestep_ind, +! timestep_sts, timing, tmunu2grid, units, writeheader ! implicit none @@ -130,7 +130,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,& epot_sinksink,get_ntypes,isdead_or_accreted,dustfrac,ddustevol,& nden_nimhd,dustevol,rhoh,gradh, & - Bevol,Bxyz,dustprop,ddustprop,ndustsmall,iboundary,eos_vars,dvdx + Bevol,Bxyz,dustprop,filfac,ddustprop,ndustsmall,iboundary,eos_vars,dvdx use part, only:pxyzu,dens,metrics,rad,radprop,drad,ithick use densityforce, only:densityiterate use linklist, only:set_linklist @@ -167,6 +167,8 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) #endif use dust, only:init_drag use growth, only:init_growth + use porosity, only:init_porosity,init_filfac + use options, only:use_porosity #ifdef MFLOW use mf_write, only:mflow_write,mflow_init use io, only:imflow @@ -312,7 +314,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) ! !--get total number of particles (on all processors) ! - ntot = reduceall_mpi('+',npart) + ntot = reduceall_mpi('+',npart) call update_npartoftypetot if (id==master) write(iprint,"(a,i12)") ' npart total = ',ntot if (npart > 0) then @@ -333,6 +335,11 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) if (use_dustgrowth) then call init_growth(ierr) if (ierr /= 0) call fatal('initial','error initialising growth variables') + if (use_porosity) then + call init_porosity(ierr) + if (ierr /= 0) call fatal('initial','error initialising porosity variables') + call init_filfac(npart,xyzh,vxyzu) + endif endif endif ! @@ -594,14 +601,14 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) do j=1,nderivinit if (ntot > 0) call derivs(1,npart,npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,& - rad,drad,radprop,dustprop,ddustprop,dustevol,ddustevol,dustfrac,& - eos_vars,time,0.,dtnew_first,pxyzu,dens,metrics) + rad,drad,radprop,dustprop,ddustprop,dustevol,ddustevol,filfac,& + dustfrac,eos_vars,time,0.,dtnew_first,pxyzu,dens,metrics) #ifdef LIVE_ANALYSIS call do_analysis(dumpfile,numfromfile(dumpfile),xyzh,vxyzu, & massoftype(igas),npart,time,ianalysis) call derivs(1,npart,npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& Bevol,dBevol,rad,drad,radprop,dustprop,ddustprop,dustevol,& - ddustevol,dustfrac,eos_vars,time,0.,dtnew_first,pxyzu,dens,metrics) + ddustevol,filfac,dustfrac,eos_vars,time,0.,dtnew_first,pxyzu,dens,metrics) if (do_radiation) call set_radiation_and_gas_temperature_equal(npart,xyzh,vxyzu,massoftype,rad) #endif diff --git a/src/main/interp_metric.F90 b/src/main/interp_metric.f90 similarity index 91% rename from src/main/interp_metric.F90 rename to src/main/interp_metric.f90 index 362eb129f..a6037037d 100644 --- a/src/main/interp_metric.F90 +++ b/src/main/interp_metric.f90 @@ -6,9 +6,10 @@ !--------------------------------------------------------------------------! module metric_interp ! -! metric_interp +! Interpolate a tabulated metric onto the particle positions ! -! :References: None +! :References: +! Magnall, Price, Lasky & Macpherson (2023), Phys. Rev D. 108, 103534 ! ! :Owner: Spencer Magnall ! @@ -16,13 +17,16 @@ module metric_interp ! ! :Dependencies: einsteintk_utils ! + implicit none interface trilinear_interp module procedure interp_g, interp_sqrtg, interp_gderiv end interface trilinear_interp + contains subroutine interp_g() + end subroutine interp_g subroutine interp_sqrtg() @@ -54,7 +58,6 @@ pure subroutine get_grid_neighbours(position,dx,xlower,ylower,zlower) ylower = ylower + 1 zlower = zlower + 1 - end subroutine get_grid_neighbours end module metric_interp diff --git a/src/main/lumin_nsdisc.F90 b/src/main/lumin_nsdisc.f90 similarity index 100% rename from src/main/lumin_nsdisc.F90 rename to src/main/lumin_nsdisc.f90 diff --git a/src/main/memory.F90 b/src/main/memory.f90 similarity index 99% rename from src/main/memory.F90 rename to src/main/memory.f90 index 5275c132a..b20dae9f4 100644 --- a/src/main/memory.F90 +++ b/src/main/memory.f90 @@ -6,7 +6,7 @@ !--------------------------------------------------------------------------! module memory ! -! None +! Wrapper routines for memory allocation ! ! :References: None ! diff --git a/src/main/options.f90 b/src/main/options.f90 index 85887742a..36bc7e5eb 100644 --- a/src/main/options.f90 +++ b/src/main/options.f90 @@ -45,7 +45,7 @@ module options real, public :: rhofinal_cgs,rhofinal1 ! dust method - logical, public :: use_dustfrac, use_hybrid + logical, public :: use_dustfrac, use_hybrid, use_porosity ! mcfost logical, public :: use_mcfost, use_Voronoi_limits_file, use_mcfost_stellar_parameters, mcfost_computes_Lacc diff --git a/src/main/part.F90 b/src/main/part.F90 index 27c6391f7..652935668 100644 --- a/src/main/part.F90 +++ b/src/main/part.F90 @@ -70,12 +70,20 @@ module part ! !--storage of dust growth properties ! - real, allocatable :: dustprop(:,:) !- size and intrinsic density + real, allocatable :: dustprop(:,:) !- mass and intrinsic density real, allocatable :: dustgasprop(:,:) !- gas related quantites interpolated on dust particles (see Force.F90) real, allocatable :: VrelVf(:) - character(len=*), parameter :: dustprop_label(2) = (/'grainsize','graindens'/) + character(len=*), parameter :: dustprop_label(2) = (/'grainmass','graindens'/) character(len=*), parameter :: dustgasprop_label(4) = (/'csound','rhogas','St ','dv '/) character(len=*), parameter :: VrelVf_label = 'Vrel/Vfrag' + + !- porosity + integer, allocatable :: dragreg(:) !- drag regime + real, allocatable :: mprev(:) !- previous mass + real, allocatable :: filfac(:) !- filling factor + real, allocatable :: filfacprev(:) !- previous filling factor needed for minimum St condition + real, allocatable :: probastick(:) !-probabily of sticking, when bounce is on + character(len=*), parameter :: filfac_label = 'filfac' !- options logical, public :: this_is_a_test = .false. logical, public :: this_is_a_flyby = .false. @@ -290,9 +298,9 @@ module part real(kind=4), allocatable :: divBsymm(:) real, allocatable :: fext(:,:) real, allocatable :: ddustevol(:,:) - real, allocatable :: ddustprop(:,:) !--grainsize is the only prop that evolves for now + real, allocatable :: ddustprop(:,:) !--grainmass is the only prop that evolves for now real, allocatable :: drad(:,:) - character(len=*), parameter :: ddustprop_label(2) = (/' ds/dt ','drho/dt'/) + character(len=*), parameter :: ddustprop_label(2) = (/' dm/dt ','drho/dt'/) ! !--storage associated with/dependent on timestepping ! @@ -301,6 +309,7 @@ module part real, allocatable :: dustpred(:,:) real, allocatable :: Bpred(:,:) real, allocatable :: dustproppred(:,:) + real, allocatable :: filfacpred(:) real, allocatable :: radpred(:,:) integer(kind=1), allocatable :: ibin(:) integer(kind=1), allocatable :: ibin_old(:) @@ -334,7 +343,7 @@ module part ! information between MPI threads ! integer, parameter, private :: usedivcurlv = min(ndivcurlv,1) - integer, parameter :: ipartbufsize = 128 + integer, parameter :: ipartbufsize = 129 real :: hfact,Bextx,Bexty,Bextz integer :: npart @@ -406,6 +415,11 @@ subroutine allocate_part call allocate_array('dustevol', dustevol, maxdustsmall, maxp_dustfrac) call allocate_array('ddustevol', ddustevol, maxdustsmall, maxdustan) call allocate_array('ddustprop', ddustprop, 2, maxp_growth) + call allocate_array('dragreg', dragreg, maxp_growth) + call allocate_array('filfac', filfac, maxp_growth) + call allocate_array('mprev', mprev, maxp_growth) + call allocate_array('probastick', probastick, maxp_growth) + call allocate_array('filfacprev', filfacprev, maxp_growth) call allocate_array('deltav', deltav, 3, maxdustsmall, maxp_dustfrac) call allocate_array('pxyzu', pxyzu, maxvxyzu, maxgr) call allocate_array('dens', dens, maxgr) @@ -434,6 +448,7 @@ subroutine allocate_part call allocate_array('dustpred', dustpred, maxdustsmall, maxdustan) call allocate_array('Bpred', Bpred, maxBevol, maxmhdan) call allocate_array('dustproppred', dustproppred, 2, maxp_growth) + call allocate_array('filfacpred', filfacpred, maxp_growth) call allocate_array('dust_temp',dust_temp,maxTdust) call allocate_array('rad', rad, maxirad, maxprad) call allocate_array('radpred', radpred, maxirad, maxprad) @@ -486,6 +501,11 @@ subroutine deallocate_part if (allocated(dustevol)) deallocate(dustevol) if (allocated(ddustevol)) deallocate(ddustevol) if (allocated(ddustprop)) deallocate(ddustprop) + if (allocated(dragreg)) deallocate(dragreg) + if (allocated(filfac)) deallocate(filfac) + if (allocated(mprev)) deallocate(mprev) + if (allocated(filfacprev)) deallocate(filfacprev) + if (allocated(probastick)) deallocate(probastick) if (allocated(deltav)) deallocate(deltav) if (allocated(pxyzu)) deallocate(pxyzu) if (allocated(dens)) deallocate(dens) @@ -514,6 +534,7 @@ subroutine deallocate_part if (allocated(dustpred)) deallocate(dustpred) if (allocated(Bpred)) deallocate(Bpred) if (allocated(dustproppred)) deallocate(dustproppred) + if (allocated(filfacpred)) deallocate(filfacpred) if (allocated(ibin)) deallocate(ibin) if (allocated(ibin_old)) deallocate(ibin_old) if (allocated(ibin_wake)) deallocate(ibin_wake) @@ -1256,6 +1277,7 @@ subroutine copy_particle_all(src,dst,new_part) dustgasprop(:,dst) = dustgasprop(:,src) VrelVf(dst) = VrelVf(src) dustproppred(:,dst) = dustproppred(:,src) + filfacpred(dst) = filfacpred(src) endif fxyz_drag(:,dst) = fxyz_drag(:,src) fxyz_dragold(:,dst) = fxyz_dragold(:,src) @@ -1710,18 +1732,44 @@ end subroutine delete_particles_outside_box ! Delete particles outside (or inside) of a defined sphere !+ !---------------------------------------------------------------- -subroutine delete_particles_outside_sphere(center,radius,np) +subroutine delete_particles_outside_sphere(center,radius,np,revert,mytype) use io, only:fatal real, intent(in) :: center(3), radius integer, intent(inout) :: np + logical, intent(in), optional :: revert + integer, intent(in), optional :: mytype + integer :: i - real :: r(3), radius_squared + real :: r(3), radius_squared + logical :: use_revert + + if (present(revert)) then + use_revert = revert + else + use_revert = .false. + endif radius_squared = radius**2 - do i=1,np - r = xyzh(1:3,i) - center - if (dot_product(r,r) > radius_squared) call kill_particle(i,npartoftype) - enddo + + if (present(mytype)) then + do i=1,np + r = xyzh(1:3,i) - center + if (use_revert) then + if (dot_product(r,r) < radius_squared .and. iamtype(iphase(i)) == mytype) call kill_particle(i,npartoftype) + else + if (dot_product(r,r) > radius_squared .and. iamtype(iphase(i)) == mytype) call kill_particle(i,npartoftype) + endif + enddo + else + do i=1,np + r = xyzh(1:3,i) - center + if (use_revert) then + if (dot_product(r,r) < radius_squared) call kill_particle(i,npartoftype) + else + if (dot_product(r,r) > radius_squared) call kill_particle(i,npartoftype) + endif + enddo + endif call shuffle_part(np) if (np /= sum(npartoftype)) call fatal('del_part_outside_sphere','particles not conserved') diff --git a/src/main/porosity.f90 b/src/main/porosity.f90 new file mode 100755 index 000000000..547584f34 --- /dev/null +++ b/src/main/porosity.f90 @@ -0,0 +1,806 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2024 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.github.io/ ! +!--------------------------------------------------------------------------! +module porosity +! +! Contains routine for porosity evolution (growth, bouncing, +! fragmentation, compaction, disruption) +! +! :References: +! Okuzumi et al. (1997), ApJ 752, 106 +! Garcia, Gonzalez (2020), MNRAS 493, 1788 +! Tatsuuma et Kataoka (2021), ApJ 913, 132 +! Michoulier & Gonzalez (2022), MNRAS 517, 3064 +! +! :Owner: Stephane Michoulier +! +! :Runtime parameters: +! - gammaft : *Force to torque efficient of gas flow on dust* +! - ibounce : *bouncing (0=Off,1=On)* +! - icompact : *Compaction during fragmentation (ifrag > 0) (0=off,1=on)* +! - idisrupt : *disruption (0=Off,1=On)* +! - iporosity : *porosity (0=Off,1=On)* +! - smonocgs : *Monomer size in cm (smaller or equal to 1.e-4 cm)* +! - surfenergSI : *Monomer surface energy in J/m**2* +! - youngmodSI : *Monomer young modulus in Pa* +! +! :Dependencies: dim, dust, eos, growth, infile_utils, io, options, part, +! physcon, random, units, viscosity +! + use units, only:umass,udist,unit_energ,unit_pressure,unit_density + use physcon, only:Ro,pi,fourpi,roottwo + implicit none + + !--Default values + + integer, public :: iporosity = 0 !--0=Off 1=On (-1=On for checkup, filfac is initialized but does not evolve) + integer, public :: icompact = 1 !--0=off 1=on (Compaction of dust grain during fragmentation) + integer, public :: ibounce = 0 !--0=off 1=on (Allow dust grains to bounce) + integer, public :: idisrupt = 0 !--0=off 1=on (Rotational disruption by gas flow: Tatsuuma et al. 2021) + real, public :: smonocgs = 1e-4 !--monomer size in cm + real, public :: surfenergSI = 0.20 !--surface energy of monomers in SI: J/m**2 (here for Si: Kimura et al. 2020) + real, public :: youngmodSI = 72e9 !--young modulus of monomers in SI: Pa (here for Si: Yamamoto et al. 2014) + real, public :: gammaft = 0.1 !--force-to-torque efficiency (Tatsuuma et al. 2021) + + real, parameter :: cratio = -0.5801454844 !--common ratio for a power + real, parameter :: b_oku = 0.15 !--parameter b (Okuzumi et al. 2012) + real, parameter :: maxpacking = 0.74048 !--max sphere packing for hexagonal close packing + + real, public :: smono !--monomer size + real, public :: mmono !--monomer mass + real, public :: surfenerg + real, public :: youngmod + real :: eroll !--rolling + real :: grainmassminlog + real :: Yd0 !test for compaction + real :: Ydpow !test for compaction + + public :: get_filfac,init_filfac,get_disruption,get_probastick + public :: init_porosity,print_porosity_info,write_options_porosity,read_options_porosity + public :: write_porosity_setup_options,read_porosity_setup_options + + private + +contains + +!------------------------------------------------ +!+ +! Initialise variables for computing porosity +!+ +!------------------------------------------------ +subroutine init_porosity(ierr) + use io, only:error + use dust, only:idrag,grainsizecgs,graindenscgs + integer, intent(out) :: ierr + + ierr = 0 + + !--initialise variables in code units + smono = smonocgs / udist + mmono = fourpi/3 * (graindenscgs / unit_density) * smono**3 + surfenerg = surfenergSI * udist * udist * 1000 / unit_energ + youngmod = youngmodSI * 10 / unit_pressure + eroll = 302.455974078*(surfenerg**5 * smono**4 / youngmod**2)**(1./3.) + + Yd0 = 9.5e6 *10/unit_pressure ! for water+silicate; 9.8e6 for water only + Ydpow = 6.4 !for silicate+water, 4 for water only + + grainmassminlog = log10(50.*mmono) + + if (smono <= 0.) then + call error('init_porosity','smonocgs <= 0',var='smonocgs',val=smonocgs) + ierr = 1 + endif + + if (grainsizecgs < smonocgs) then + call error('init_porosity','grainsizecgs < smonocgs',var='smonocgs',val=smonocgs) + ierr = 1 + endif + + if (surfenerg <= 0.) then + call error('init_porosity','surfenerg <= 0',var='surfenerg',val=surfenerg) + ierr = 2 + endif + + if (youngmod <= 0.) then + call error('init_porosity','youngmod <= 0',var='youngmod',val=youngmod) + ierr = 3 + endif + + if (idrag /= 1) then + call error('init_porosity','idrag = 1 should be used for porosity',var='idrag',val=real(idrag)) + ierr = 4 + endif + +end subroutine init_porosity + +!----------------------------------------------------------------------- +!+ +! Compute the initial filling factor +!+ +!----------------------------------------------------------------------- +subroutine init_filfac(npart,xyzh,vxyzu) + use options, only:use_dustfrac + use viscosity, only:shearparam + use part, only:idust,igas,iamtype,iphase,massoftype,& + rhoh,dustfrac,dustprop,filfac,Omega_k + use dust, only:get_viscmol_nu!,grainsizecgs + use eos, only:gamma,get_spsound + integer, intent(in) :: npart + real, intent(in) :: xyzh(:,:) + real, intent(inout) :: vxyzu(:,:) + + integer :: i,iam + real :: rho,rhogas,cs,cparam,coeff_gei,nu + real :: sfrac,s1,s2,s3,filfacmax +! real :: mfrac,m1,m2,m3 + + + select case (iporosity) ! add other case for other models here + case (1) + + !--initialize filling factor (Garcia & Gonzalez 2020, Suyama et al. 2008, Okuzumi et al. 2012) + + if (all(filfac(:) == 0.)) then ! check if filfac(i) was already initialize by init_filfac or not + coeff_gei = sqrt(8./(pi*gamma)) + do i=1,npart + iam = iamtype(iphase(i)) + if (iam == idust .or. (iam == igas .and. use_dustfrac)) then + sfrac = (dustprop(1,i)/mmono)**(1./3.) + if (sfrac > 1.) then ! if grainsize > monomer size, compute filling factor + !- compute rho, rhogas and cs + if (iam == igas .and. use_dustfrac) then + rho = rhoh(xyzh(4,i),massoftype(igas)) + rhogas = rho*(1-dustfrac(1,i)) + cs = get_spsound(3,xyzh(:,i),rhogas,vxyzu(:,i)) + else + rhogas = rhoh(xyzh(4,i),massoftype(igas)) + cs = get_spsound(3,xyzh(:,i),rhogas,vxyzu(:,i)) + rho = rhogas + rhoh(xyzh(4,i),massoftype(idust)) + endif + + !- molecular viscosity + nu = get_viscmol_nu(cs,rhogas) + + !- shared parameter for the following filling factors + cparam = (243.*pi*roottwo/15625.)*(Ro*shearparam*smono**4*dustprop(2,i)*dustprop(2,i)*cs & + *Omega_k(i))/(rho*b_oku*eroll) + + !--transition masses m1/mmono and m2/mmono between hit&stick and Epstein/Stokes regimes with St < 1 + s1 = (cparam/(2.*(2.**0.075 - 1.)*coeff_gei))**((1.-cratio)/(1.+8.*cratio)) + s2 = (cparam*cs*smono/(9.*nu*(2.**0.2 - 1.)))**((1.-cratio)/(9.*cratio)) + + !--we assume St < 1 here (grainsizecgs < 100-1000 cm) + if (s1 < s2) then + if (sfrac < s1) then ! filling factor: hit&stick regime + filfac(i) = sfrac**(3.*cratio/(1.-cratio)) + else + !- transition masses m3/mmono between Epstein and Stokes regimes with St < 1 + s3 = s1**((1.+8.*cratio)/(1.-cratio)) / s2**(9.*cratio/(1.-cratio)) + + if (sfrac < s3) then ! filling factor: Epstein regime - St<1 + filfac(i) = s1**((1.+8.*cratio)/(3.-3.*cratio))/sfrac**(1./3.) + else ! filling factor: Stokes regime - St<1 + filfac(i) = s2**(3.*cratio/(1.-cratio)) + endif + endif + else + if (sfrac < s2) then ! filling factor: hit&stick regime + filfac(i) = sfrac**(3.*cratio/(1.-cratio)) + else ! filling factor: Stokes regime - St<1 + filfac(i) = s2**(3.*cratio/(1.-cratio)) + endif + endif + + !- max value of filfac is maxpacking == max compaction + filfacmax = 0.5*maxpacking *(1+ sqrt(1 + 4*(1.-maxpacking)/maxpacking/maxpacking*sfrac**(-3))) + if (filfac(i) > filfacmax) filfac(i) = filfacmax + + !- Compute grain mass of the grain using grain size and filfac + dustprop(1,i) = filfac(i) * dustprop(1,i) + else + filfac(i) = 1. + dustprop(1,i) = mmono + endif + endif + enddo + endif + case (-1) + !--initialize filling factor for compact grains + if (all(filfac(:) == 0.)) then ! check if filfac(i) was already initialize by init_filfac or not + do i=1,npart + iam = iamtype(iphase(i)) + if (iam == idust .or. (iam == igas .and. use_dustfrac)) then + sfrac = (dustprop(1,i)/mmono)**(1./3.) + if (sfrac > 1.) then ! if grainsize > monomer size, compute filling factor + filfac(i) = 1. + else + filfac(i) = 1. + dustprop(1,i) = mmono + endif + endif + enddo + endif + end select + +end subroutine init_filfac + +!---------------------------------------------------------- +!+ +! print information about porosity +!+ +!---------------------------------------------------------- +subroutine print_porosity_info(iprint) + integer, intent(in) :: iprint + + if (iporosity == 1) then + write(iprint,"(a)") ' Using porosity ' + if (icompact == 1) then + write(iprint,"(a)") ' Using compaction during fragmentation ' + endif + write(iprint,"(2(a,1pg10.3),a)")' Monomer size = ',smonocgs,' cm = ',smono,' (code units)' + write(iprint,"(2(a,1pg10.3),a)")' Surface energy = ',surfenergSI,' J/m**2 = ',surfenerg,' (code units)' + write(iprint,"(2(a,1pg10.3),a)")' Young modulus = ',youngmodSI,' Pa = ',youngmod,' (code units)' + endif + +end subroutine print_porosity_info + +!----------------------------------------------------------------------- +!+ +! Compute the final filling factor +!+ +!----------------------------------------------------------------------- +subroutine get_filfac(npart,xyzh,mprev,filfac,dustprop,dt) + use dim, only:use_dustgrowth + use options, only:use_dustfrac + use part, only:rhoh,idust,igas,iamtype,iphase,isdead_or_accreted,& + massoftype,dustfrac,dustgasprop,VrelVf,probastick + integer, intent(in) :: npart + real, intent(in) :: dt + real, intent(inout) :: filfac(:),dustprop(:,:) + real, intent(in) :: xyzh(:,:),mprev(:) + integer :: i,iam + real :: filfacevol,filfacmin,filfacmax + real :: rho,rhod + + select case (iporosity) ! add other cases for other models here + case (1) + !$omp parallel do default(none) & + !$omp shared(xyzh,npart,iphase,massoftype,use_dustfrac,dustfrac,icompact) & + !$omp shared(mprev,filfac,dustprop,dustgasprop,VrelVf,probastick,mmono,dt,ibounce) & + !$omp private(i,iam,rho,rhod,filfacevol,filfacmin,filfacmax) + do i=1,npart + if (.not.isdead_or_accreted(xyzh(4,i))) then + iam = iamtype(iphase(i)) + + if (iam == idust .or. (iam == igas .and. use_dustfrac)) then + if (dustprop(1,i) > mmono) then + !- compute rho = rho_gas + rho_dust + + if (use_dustfrac .and. iam == igas) then + rho = rhoh(xyzh(4,i),massoftype(igas)) + rhod = rho*dustfrac(1,i) + else + rhod = rhoh(xyzh(4,i),massoftype(idust)) + rho = dustgasprop(2,i) + rhod + endif + + call get_filfac_min(i,rho,dustprop(1,i)/mmono,dustprop(2,i),dustgasprop(:,i),filfacmin) + !--if new mass > previous mass, compute the new filling factor due to growth + if (dustprop(1,i) > mprev(i)) then + call get_filfac_growth(mprev(i),dustprop(1,i),filfac(i),dustgasprop(:,i),filfacevol) + if (ibounce == 1) call get_filfac_bounce(mprev(i),dustprop(2,i),filfac(i),& + dustgasprop(:,i),probastick(i),rhod,dt,filfacevol,filfacmin) + !--if new mass < previous mass, compute the new filling factor due to fragmentation + else + call get_filfac_frag(mprev(i),dustprop(:,i),filfac(i),dustgasprop(:,i),rhod,VrelVf(i),dt,filfacevol) + endif + filfac(i) = filfacevol + + !--check if the filling factor is smaller than the minimum filling factor + filfac(i) = max(filfac(i),filfacmin) + !-- max value of filfac is maxpacking == max compaction + filfacmax = maxpacking + (1.-maxpacking)*mmono/dustprop(1,i) + filfac(i) = min(filfac(i),filfacmax) + else + filfac(i) = 1. + dustprop(1,i) = mmono + endif + endif + else + filfac(i) = 0. + endif + enddo + !$omp end parallel do + end select + +end subroutine get_filfac + +!----------------------------------------------------------------------- +!+ +! Compute the filling factor during growth +!+ +!----------------------------------------------------------------------- +subroutine get_filfac_growth(mprev,mass,filfac,dustgasprop,filfacgrowth) + use viscosity, only:shearparam + use growth, only:vrelative + real, intent(in) :: mprev,mass,filfac + real, intent(in) :: dustgasprop(:) + real, intent(out) :: filfacgrowth + real :: ekincdt,vrel,vt + real :: j ! Power of the filling factor dependency in mass + + vt = sqrt(roottwo*Ro*shearparam)*dustgasprop(1) + vrel = vrelative(dustgasprop,vt) + + !- kinetic energy condition Ekin/(3*b_oku/eroll) + ekincdt = mprev*vrel*vrel/(12.*b_oku*eroll) + + !-choose power according to the value of ekincdt + if (ekincdt <= 1.) then + j = cratio + else + j = -0.2 + endif + + !- filling factor due to growth + filfacgrowth = filfac*(mass/mprev)**j + +end subroutine get_filfac_growth + +!----------------------------------------------------------------------- +!+ +! Compute the filling factor during bounce +!+ +!----------------------------------------------------------------------- +subroutine get_filfac_bounce(mprev,graindens,filfac,dustgasprop,probastick,rhod,dt,filfacevol,filfacmin) + use viscosity, only:shearparam + use growth, only:vrelative,get_size + use physcon, only:fourpi + real, intent(in) :: mprev,graindens,filfac,probastick,rhod,dt + real, intent(in) :: dustgasprop(:),filfacmin + real, intent(inout) :: filfacevol + real :: sdust,vrel,ncoll,vol,deltavol + real :: ekin,pdyn,coeffrest,filfacbnc + real :: vstick,vyield,vend,vt + + if (probastick < 1.) then + vt = sqrt(roottwo*Ro*shearparam)*dustgasprop(1) + vrel = vrelative(dustgasprop,vt) + sdust = get_size(mprev,graindens,filfac) + vstick = compute_vstick(mprev,sdust) !-compute vstick, i.e. max velocity before bouncing appears + + if (vrel >= vstick) then !-if vrel>=vstick -> bouncing + vyield = compute_vyield(vstick) !-compute vyield, i.e. max velocity before inelastic collisions appear + vend = compute_vend(vstick) !-compute vend, i.e. max velocity before there is only bouncing => no growth + + if (vrel < vyield) then !-elastic collision, no compaction + filfacbnc = filfac + else !-inelastic collision, compaction + vol = fourpi/3. * sdust**3 + ncoll = fourpi*sdust**2*rhod*vrel*dt/mprev !-number of collision in dt + ekin = mprev*vrel*vrel/4. + coeffrest = get_coeffrest(vstick/vrel,vyield/vrel) !-coefficient of restitution + !pdyn = eroll * (filfac/(maxpacking - filfac)/smono)**3 + pdyn = eroll /((1./filfac - 1./maxpacking)*smono)**3 + deltavol = (1.-coeffrest*coeffrest)*ekin/pdyn + if (deltavol > vol) deltavol = vol + + filfacbnc = filfac *(1./(1.-0.5*(deltavol/vol)))**ncoll + if (filfacbnc > maxpacking) filfacbnc = maxpacking + endif + + if (vrel < vend) then !-final filfac is a combination of filfac due to growth + bouncing + if (filfacevol < filfacmin) filfacevol = filfacmin + filfacevol = filfacevol*probastick + (1-probastick)*filfacbnc + else + filfacevol = filfacbnc + endif + endif + endif + +end subroutine get_filfac_bounce + +!----------------------------------------------------------------------- +!+ +! Compute the filling factor during fragmentation +!+ +!----------------------------------------------------------------------- +subroutine get_filfac_frag(mprev,dustprop,filfac,dustgasprop,rhod,VrelVf,dt,filfacfrag) + use viscosity, only:shearparam + use growth, only:vrelative,get_size + use physcon, only:fourpi + real, intent(in) :: mprev,filfac,rhod,VrelVf,dt + real, intent(in) :: dustprop(:),dustgasprop(:) + real, intent(out) :: filfacfrag + real :: sdust,vrel,ncoll,vol,deltavol!,compfactor + real :: ekin,pdyn,vt + + select case (icompact) + case (1) + ! model Garcia + Kataoka mod + sdust = get_size(mprev,dustprop(2),filfac) + vol = fourpi/3. * sdust**3 + vt = sqrt(roottwo*Ro*shearparam)*dustgasprop(1) + vrel = vrelative(dustgasprop,vt) + ncoll = fourpi*sdust**2*rhod*vrel*dt/mprev !number of collisions in dt + + ekin = mprev*vrel*vrel/4. - (2.*mprev - dustprop(1))*0.85697283*eroll/mmono !0.856973 = 3* 1.8 * 48/302.46 + pdyn = eroll /((1./filfac - 1./maxpacking)*smono)**3 + deltavol = ekin/pdyn !-ekin is kinetic energy - all energy needed to break monomers + + if (deltavol < 0) deltavol = 0. + if (deltavol > vol) deltavol = vol + + filfacfrag = filfac *(1./(1.-0.5*exp(1-VrelVf**2.)*deltavol/vol))**ncoll + case default ! (0) + ! Fragmentation at constant filling factor + filfacfrag = filfac + end select + +end subroutine get_filfac_frag + +!----------------------------------------------------------------------- +!+ +! Compute the filling factor in the collisional compression regime +!+ +!----------------------------------------------------------------------- +subroutine get_filfac_col(i,rho,mfrac,graindens,dustgasprop,filfaccol) + use part, only:Omega_k,dragreg + use viscosity, only:shearparam + use dust, only:get_viscmol_nu + use eos, only:gamma + integer, intent(in) :: i + real, intent(in) :: rho,mfrac,graindens + real, intent(in) :: dustgasprop(:) + real, intent(out) :: filfaccol + real :: cparam,coeff_gei,nu,kwok + real :: m1,m2,m3,m4,m5 + + !--compute filling factor due to collisions (Garcia & Gonzalez 2020, Suyama et al. 2008, Okuzumi et al. 2012) + + !- shared parameter for the following filling factors + cparam = (243.*pi*roottwo/15625.)*(Ro*shearparam*smono**4*graindens*graindens*dustgasprop(1) & + *Omega_k(i))/(rho*b_oku*eroll) + + coeff_gei = sqrt(8./(pi*gamma)) + + !- molecular viscosity + nu = get_viscmol_nu(dustgasprop(1),dustgasprop(2)) + + !- Kwok (1975) correction for supersonic drag is important + if (dragreg(i) == 2) then + kwok = sqrt(1.+9.*pi/128.*dustgasprop(4)*dustgasprop(4)/(dustgasprop(1)*dustgasprop(1))) + else + kwok = 1. + endif + + !--transition sizes m1/mmono and m2/mmono between hit&stick and Epstein/Stokes regimes with St < 1 + m1 = (cparam/(2.*(2.**0.075 - 1.)*coeff_gei*kwok))**(0.375/(cratio+0.125)) + m2 = (cparam*dustgasprop(1)*smono/(9.*nu*(2.**0.2 - 1.)))**(1./(3.*cratio)) + + if (dustgasprop(3) <= 1) then !- Stokes < 1 + if (m1 < m2) then + if (mfrac < m1) then !- filling factor: hit&stick regime + filfaccol = mfrac**cratio + else + !- transition masses m3/mmono between Epstein and Stokes regimes with St < 1 + m3 = m1**(8.*cratio+1.) / m2**(8*cratio) + if (mfrac < m3) then !- filling factor: Epstein regime - St<1 + filfaccol = m1**(cratio+0.125)/mfrac**(0.125) + else !- filling factor: Stokes regime - St<1 + filfaccol = m2**cratio + endif + endif + else + if (mfrac < m2) then !- filling factor: hit&stick regime + filfaccol = mfrac**cratio + else !- filling factor: Stokes regime - St<1 + filfaccol = m2**cratio + endif + endif + else !- Stokes > 1 + !--transition masses m4/mmono and m5/mmono between hit&stick and Epstein/Stokes regimes with St > 1 + m4 = (rho*coeff_gei*kwok*dustgasprop(1)/(graindens*smono*Omega_k(i)))**4 / m1**((cratio+0.125)/0.375) + m5 = (9.*nu*rho/(2.*graindens*smono**2*Omega_k(i)))**1.5 / m2**(0.5*cratio) + + if (m4 < m5) then !- filling factor: Epstein regime - St>1 + filfaccol = m1**(cratio+0.125) * m4**0.075 / mfrac**0.2 + else !- filling factor: Stokes regime - St>1 + filfaccol = m2**cratio * (m5/mfrac)**0.2 + endif + endif + +end subroutine get_filfac_col + +!----------------------------------------------------------------------- +!+ +! Compute the minimum filling factor +!+ +!----------------------------------------------------------------------- +subroutine get_filfac_min(i,rho,mfrac,graindens,dustgasprop,filfacmin) + use part, only:Omega_k + integer, intent(in) :: i + real, intent(in) :: rho,mfrac,graindens + real, intent(in) :: dustgasprop(:) + real, intent(out) :: filfacmin + real :: filfaccol,filfacgas,filfacgrav + + call get_filfac_col(i,rho,mfrac,graindens,dustgasprop,filfaccol) + + !--compute filling factor due to gas drag compression (Garcia & Gonzalez 2020, Kataoka et al. 2013a) + filfacgas = ((mmono*smono*dustgasprop(4)*Omega_k(i))/(pi*eroll*dustgasprop(3)))**(3./7.) * mfrac**(1./7.) + + !--compute filling factor due to self-gravity (Garcia & Gonzalez 2020, Kataoka et al. 2013b) + filfacgrav = (mmono*mmono/(pi*smono*eroll))**0.6 * mfrac**0.4 + + !--return the maximum filling factor between filfaccol, filfacgas and filfacgrav + filfacmin = max(filfaccol,filfacgas,filfacgrav) + +end subroutine get_filfac_min + + +subroutine get_disruption(npart,xyzh,filfac,dustprop,dustgasprop) + use options, only:use_dustfrac + use part, only:idust,igas,iamtype,iphase,massoftype,isdead_or_accreted,rhoh + use growth, only:check_dustprop,get_size + use random, only:ran2 + integer, intent(in) :: npart + real, intent(in) :: xyzh(:,:),dustgasprop(:,:) + real, intent(inout) :: dustprop(:,:),filfac(:) + integer :: i,iam,seed + real :: stress,strength,filfacmin,rho + real :: grainmasscurlog,grainmassmaxlog,randmass + + select case (idisrupt) + case(1) + !$omp parallel do default(none) & + !$omp shared(xyzh,npart,massoftype,iphase,use_dustfrac) & + !$omp shared(filfac,dustprop,dustgasprop,mmono,smono,grainmassminlog,surfenerg,gammaft) & + !$omp private(grainmasscurlog,grainmassmaxlog,randmass,seed) & + !$omp private(i,iam,rho,filfacmin,stress,strength) + do i=1, npart + if (.not.isdead_or_accreted(xyzh(4,i))) then + iam = iamtype(iphase(i)) + if (iam == idust .or. (iam == igas .and. use_dustfrac)) then + + stress = 25./36. * dustprop(2,i) * filfac(i) * gammaft**2 * dustgasprop(4,i)**2 + strength = 0.6*filfac(i)**(1.8)*surfenerg/smono + seed = int(stress) + + if (stress >= strength) then !-grain is rotationnaly disrupted + !-compute rho to compute filfacmin + if (use_dustfrac .and. iam == igas) then + rho = rhoh(xyzh(4,i),massoftype(igas)) + else + rho = dustgasprop(2,i) + rhoh(xyzh(4,i),massoftype(idust)) + endif + + !-compute current, current/2 and min mass in log10 + grainmasscurlog = log10(dustprop(1,i)) + grainmassmaxlog = log10(dustprop(1,i)/(2.)) + + !--call random number between 2 float values to assign a random mass to dustprop(1) + if (grainmassmaxlog > grainmassminlog) then + randmass = (grainmassmaxlog - grainmassminlog) * ran2(seed) + grainmassminlog + else + if (grainmasscurlog > grainmassminlog) then + randmass = grainmassminlog + else + randmass = grainmasscurlog + endif + endif + + dustprop(1,i) = 10.**randmass + + !-compute filfacmin and compare it to filfac(i) + call get_filfac_min(i,rho,dustprop(1,i)/mmono,dustprop(2,i),dustgasprop(:,i),filfacmin) + filfac(i) = max(filfac(i),filfacmin) + endif + endif + endif + enddo + !$omp end parallel do + end select + +end subroutine get_disruption + +!----------------------------------------------------------------------- +!+ +! Compute the probability of bounce and associated growth rate +!+ +!----------------------------------------------------------------------- + +subroutine get_probastick(npart,xyzh,dmdt,dustprop,dustgasprop,filfac) + use options, only:use_dustfrac + use part, only:idust,igas,iamtype,iphase,isdead_or_accreted,rhoh,probastick + use viscosity, only:shearparam + use growth, only:vrelative,get_size + integer, intent(in) :: npart + real, intent(in) :: filfac(:) + real, intent(in) :: xyzh(:,:),dustprop(:,:),dustgasprop(:,:) + real, intent(inout) :: dmdt(:) + integer :: i,iam + real :: vrel,vstick,vend,sdust,vt + + if (ibounce == 1) then + !$omp parallel do default(none) & + !$omp shared(xyzh,npart,iphase,use_dustfrac) & + !$omp shared(filfac,dmdt,dustprop,dustgasprop,probastick,shearparam) & + !$omp private(i,iam,vrel,vstick,vend,sdust,vt) + do i=1, npart + if (.not.isdead_or_accreted(xyzh(4,i))) then + iam = iamtype(iphase(i)) + if ((iam == idust .or. (iam == igas .and. use_dustfrac))) then + if (filfac(i) >= 0.3 .and. dmdt(i) >= 0.) then + vt = sqrt(roottwo*Ro*shearparam)*dustgasprop(1,i) + vrel = vrelative(dustgasprop(:,i),vt) + sdust = get_size(dustprop(1,i),dustprop(2,i),filfac(i)) + vstick = compute_vstick(dustprop(1,i),sdust) + vend = compute_vend(vstick) + + !compute the probability of bounce depending on the velocity + if (vrel >= vstick) then + if (vrel < vend) then + probastick(i) = (log(vrel)-log(vend))/(log(vstick)-log(vend)) + else + probastick(i) = 0. !full bounce -> no growth + endif + else + probastick(i) = 1. + endif + else + probastick(i) = 1. + endif + !compute new growth rate + dmdt(i) = dmdt(i)*probastick(i) + endif + endif + enddo + !$omp end parallel do + endif + +end subroutine get_probastick + +!----------------------------------------------------------------------- +!+ +! Write porosity options in the input file +!+ +!----------------------------------------------------------------------- +subroutine write_options_porosity(iunit) + use infile_utils, only:write_inopt + integer, intent(in) :: iunit + + write(iunit,"(/,a)") '# options controlling porosity (require idrag=1)' + call write_inopt(iporosity,'iporosity','porosity (0=off,1=on) ',iunit) + if (iporosity == 1 .or. iporosity == -1) then + call write_inopt(icompact, 'icompact', 'Compaction during fragmentation (ifrag > 0) (0=off,1=on)', iunit) + call write_inopt(ibounce, 'ibounce', 'Dust bouncing (0=off,1=on)', iunit) + call write_inopt(idisrupt, 'idisrupt', 'Rotational disruption (0=off,1=on)', iunit) + call write_inopt(smonocgs,'smonocgs','Monomer size in cm (smaller or equal to 1.e-4 cm)',iunit) + call write_inopt(surfenergSI,'surfenergSI','Monomer surface energy in J/m**2',iunit) + call write_inopt(youngmodSI,'youngmodSI','Monomer young modulus in Pa',iunit) + call write_inopt(gammaft,'gammaft','Force to torque efficient of gas flow on dust',iunit) + endif + +end subroutine write_options_porosity + +!----------------------------------------------------------------------- +!+ +! Read porosity options from the input file +!+ +!----------------------------------------------------------------------- +subroutine read_options_porosity(name,valstring,imatch,igotall,ierr) + use options, only: use_porosity + character(len=*), intent(in) :: name,valstring + logical, intent(out) :: imatch,igotall + integer, intent(out) :: ierr + + integer, save :: ngot = 0 + + imatch = .true. + igotall = .false. + + select case(trim(name)) + case('iporosity') + read(valstring,*,iostat=ierr) iporosity + ngot = ngot + 1 + if (iporosity == 1 .or. iporosity == -1) use_porosity = .true. + case('icompact') + read(valstring,*,iostat=ierr) icompact + ngot = ngot + 1 + case('ibounce') + read(valstring,*,iostat=ierr) ibounce + ngot = ngot + 1 + case('idisrupt') + read(valstring,*,iostat=ierr) idisrupt + ngot = ngot + 1 + case('smonocgs') + read(valstring,*,iostat=ierr) smonocgs + ngot = ngot + 1 + case('surfenergSI') + read(valstring,*,iostat=ierr) surfenergSI + ngot = ngot + 1 + case('youngmodSI') + read(valstring,*,iostat=ierr) youngmodSI + ngot = ngot + 1 + case('gammaft') + read(valstring,*,iostat=ierr) gammaft + ngot = ngot + 1 + case default + imatch = .false. + end select + + if ((iporosity == 0) .and. ngot == 1) igotall = .true. + if ((iporosity /= 0) .and. ngot == 8) igotall = .true. + +end subroutine read_options_porosity + +!----------------------------------------------------------------------- +!+ +! Write porosity options to the .setup file +!+ +!----------------------------------------------------------------------- +subroutine write_porosity_setup_options(iunit) + use infile_utils, only:write_inopt + integer, intent(in) :: iunit + + write(iunit,"(/,a)") '# options for porosity' + call write_inopt(iporosity,'iporosity','porosity (0=Off,1=On)',iunit) + call write_inopt(ibounce,'ibounce','bouncing (0=Off,1=On)',iunit) + call write_inopt(idisrupt,'idisrupt','disruption (0=Off,1=On)',iunit) + +end subroutine write_porosity_setup_options + +!----------------------------------------------------------------------- +!+ +! Read growth options from the .setup file +!+ +!----------------------------------------------------------------------- +subroutine read_porosity_setup_options(db,nerr) + use options, only:use_porosity + use infile_utils, only:read_inopt,inopts + type(inopts), allocatable, intent(inout) :: db(:) + integer, intent(inout) :: nerr + + call read_inopt(iporosity,'iporosity',db,min=-1,max=1,errcount=nerr) + if (iporosity == 1 .or. iporosity == -1) use_porosity = .true. + call read_inopt(ibounce,'ibounce',db,min=0,max=1,errcount=nerr) + call read_inopt(idisrupt,'idisrupt',db,min=0,max=1,errcount=nerr) + +end subroutine read_porosity_setup_options + +real function get_coeffrest(vstickvrel,vyieldvrel) + real, intent(in) :: vstickvrel,vyieldvrel + + if (vyieldvrel >= 1.) then + get_coeffrest = sqrt(1.-vstickvrel*vstickvrel) + else + get_coeffrest = sqrt(1.2*sqrt(3.)*(1.-(vyieldvrel*vyieldvrel/6.))*& + sqrt(1./(1.+2.*sqrt((1.2/(vyieldvrel*vyieldvrel))-0.2)))-(vstickvrel*vstickvrel)) + endif + +end function get_coeffrest + +!--velocity limit between full sticking regime and partial sticking + bouncing regime +real function compute_vstick(mass,size) + real, intent(in) ::mass,size + compute_vstick = 8.76*((surfenerg**5 * size**4)/(mass**3*youngmod**2))**(1./6.) +end function + +!--velocity limit between elastic and inelastic bouncing regime +real function compute_vyield(vstick) + real, intent(in) ::vstick + compute_vyield = 10.*vstick +end function + +!--velocity limit between partial sticking + bouncing regime and full bouncing regime +real function compute_vend(vstick) + real, intent(in) ::vstick + compute_vend = 24343220.*vstick +end function + +end module porosity diff --git a/src/main/readwrite_dumps_common.F90 b/src/main/readwrite_dumps_common.f90 similarity index 99% rename from src/main/readwrite_dumps_common.F90 rename to src/main/readwrite_dumps_common.f90 index 998e45e61..5cf289f29 100644 --- a/src/main/readwrite_dumps_common.F90 +++ b/src/main/readwrite_dumps_common.f90 @@ -285,7 +285,7 @@ subroutine check_arrays(i1,i2,noffset,npartoftype,npartread,nptmass,nsinkpropert dustfrac = 0. endif if (use_dustgrowth .and. .not.got_dustprop(1)) then - if (id==master) write(*,*) 'ERROR! using dustgrowth, but no grain size found in dump file' + if (id==master) write(*,*) 'ERROR! using dustgrowth, but no grain mass found in dump file' ierr = ierr + 1 endif if (use_dustgrowth .and. .not.got_dustprop(2)) then diff --git a/src/main/readwrite_dumps_fortran.F90 b/src/main/readwrite_dumps_fortran.F90 index 065057bca..2d80c89b4 100644 --- a/src/main/readwrite_dumps_fortran.F90 +++ b/src/main/readwrite_dumps_fortran.F90 @@ -214,12 +214,12 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG) maxptmass,get_pmass,nabundances,abundance,abundance_label,mhd,& divcurlv,divcurlv_label,divcurlB,divcurlB_label,poten,dustfrac,deltav,deltav_label,tstop,& dustfrac_label,tstop_label,dustprop,dustprop_label,eos_vars,eos_vars_label,ndusttypes,ndustsmall,VrelVf,& - VrelVf_label,dustgasprop,dustgasprop_label,dust_temp,pxyzu,pxyzu_label,dens,& !,dvdx,dvdx_label + VrelVf_label,dustgasprop,dustgasprop_label,filfac,filfac_label,dust_temp,pxyzu,pxyzu_label,dens,& !,dvdx,dvdx_label rad,rad_label,radprop,radprop_label,do_radiation,maxirad,maxradprop,itemp,igasP,igamma,& iorig,iX,iZ,imu,nucleation,nucleation_label,n_nucleation,tau,itau_alloc,tau_lucy,itauL_alloc,& luminosity,eta_nimhd,eta_nimhd_label use part, only:metrics,metricderivs,tmunus - use options, only:use_dustfrac,use_var_comp,icooling + use options, only:use_dustfrac,use_porosity,use_var_comp,icooling use dump_utils, only:tag,open_dumpfile_w,allocate_header,& free_header,write_header,write_array,write_block_header use mpiutils, only:reduce_mpi,reduceall_mpi @@ -354,6 +354,7 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG) call write_array(1,dustprop,dustprop_label,2,npart,k,ipass,idump,nums,ierrs(3)) call write_array(1,VrelVf,VrelVf_label,npart,k,ipass,idump,nums,ierrs(3)) call write_array(1,dustgasprop,dustgasprop_label,4,npart,k,ipass,idump,nums,ierrs(3)) + if (use_porosity) call write_array(1,filfac,filfac_label,npart,k,ipass,idump,nums,ierrs(3)) endif if (h2chemistry) call write_array(1,abundance,abundance_label,nabundances,npart,k,ipass,idump,nums,ierrs(5)) if (use_dust) call write_array(1,dustfrac,dustfrac_label,ndusttypes,npart,k,ipass,idump,nums,ierrs(7)) @@ -517,13 +518,15 @@ end subroutine write_fulldump_fortran subroutine write_smalldump_fortran(t,dumpfile) use dim, only:maxp,maxtypes,use_dust,lightcurve,use_dustgrowth,h2chemistry + use options, only:use_porosity use io, only:idump,iprint,real4,id,master,error,warning,nprocs use part, only:xyzh,xyzh_label,npart,Bxyz,Bxyz_label,& npartoftypetot,update_npartoftypetot,& maxphase,iphase,nabundances,& nptmass,nsinkproperties,xyzmh_ptmass,xyzmh_ptmass_label,& abundance,abundance_label,mhd,dustfrac,iamtype_int11,& - dustprop,dustprop_label,dustfrac_label,ndusttypes,& + dustprop,dustprop_label,dustfrac_label,& + filfac,filfac_label,ndusttypes,& rad,rad_label,do_radiation,maxirad,luminosity use dump_utils, only:open_dumpfile_w,dump_h,allocate_header,free_header,& write_header,write_array,write_block_header @@ -596,6 +599,7 @@ subroutine write_smalldump_fortran(t,dumpfile) call write_array(1,xyzh,xyzh_label,3,npart,k,ipass,idump,nums,ierr,singleprec=.true.) if (use_dustgrowth) then call write_array(1,dustprop,dustprop_label,2,npart,k,ipass,idump,nums,ierr,singleprec=.true.) + if (use_porosity) call write_array(1,filfac,filfac_label,npart,k,ipass,idump,nums,ierr,singleprec=.true.) endif if (h2chemistry .and. nabundances >= 1) & call write_array(1,abundance,abundance_label,1,npart,k,ipass,idump,nums,ierr,singleprec=.true.) @@ -1133,11 +1137,12 @@ subroutine read_phantom_arrays(i1,i2,noffset,narraylengths,nums,npartread,nparto alphaind,poten,xyzmh_ptmass,xyzmh_ptmass_label,vxyz_ptmass,vxyz_ptmass_label, & Bevol,Bxyz,Bxyz_label,nabundances,iphase,idust, & eos_vars,eos_vars_label,maxeosvars,dustprop,dustprop_label,divcurlv,divcurlv_label,iX,iZ,imu, & - VrelVf,VrelVf_label,dustgasprop,dustgasprop_label,pxyzu,pxyzu_label,dust_temp, & + VrelVf,VrelVf_label,dustgasprop,dustgasprop_label,filfac,filfac_label,pxyzu,pxyzu_label,dust_temp, & rad,rad_label,radprop,radprop_label,do_radiation,maxirad,maxradprop,ifluxx,ifluxy,ifluxz, & nucleation,nucleation_label,n_nucleation,ikappa,tau,itau_alloc,tau_lucy,itauL_alloc,& ithick,ilambda,iorig,dt_in,krome_nmols,T_gas_cool use sphNGutils, only:mass_sphng,got_mass,set_gas_particle_mass + use options, only:use_porosity integer, intent(in) :: i1,i2,noffset,narraylengths,nums(:,:),npartread,npartoftype(:),idisk1,iprint real, intent(in) :: massoftype(:) integer, intent(in) :: nptmass,nsinkproperties @@ -1151,7 +1156,7 @@ subroutine read_phantom_arrays(i1,i2,noffset,narraylengths,nums,npartread,nparto logical :: got_krome_mols(krome_nmols),got_krome_T,got_krome_gamma,got_krome_mu logical :: got_eosvars(maxeosvars),got_nucleation(n_nucleation),got_ray_tracer logical :: got_psi,got_Tdust,got_dustprop(2),got_VrelVf,got_dustgasprop(4) - logical :: got_divcurlv(4),got_rad(maxirad),got_radprop(maxradprop),got_pxyzu(4),got_iorig + logical :: got_filfac,got_divcurlv(4),got_rad(maxirad),got_radprop(maxradprop),got_pxyzu(4),got_iorig character(len=lentag) :: tag,tagarr(64) integer :: k,i,iarr,ik,ndustfraci @@ -1172,6 +1177,7 @@ subroutine read_phantom_arrays(i1,i2,noffset,narraylengths,nums,npartread,nparto got_eosvars = .false. got_dustprop = .false. got_VrelVf = .false. + got_filfac = .false. got_dustgasprop = .false. got_divcurlv = .false. got_Tdust = .false. @@ -1215,6 +1221,7 @@ subroutine read_phantom_arrays(i1,i2,noffset,narraylengths,nums,npartread,nparto call read_array(dustprop,dustprop_label,got_dustprop,ik,i1,i2,noffset,idisk1,tag,match,ierr) call read_array(VrelVf,VrelVf_label,got_VrelVf,ik,i1,i2,noffset,idisk1,tag,match,ierr) call read_array(dustgasprop,dustgasprop_label,got_dustgasprop,ik,i1,i2,noffset,idisk1,tag,match,ierr) + if (use_porosity) call read_array(filfac,filfac_label,got_filfac,ik,i1,i2,noffset,idisk1,tag,match,ierr) endif if (use_dust) then if (any(tag == dustfrac_label)) then diff --git a/src/main/readwrite_infile.F90 b/src/main/readwrite_infile.F90 index 83278feaf..c7aec5707 100644 --- a/src/main/readwrite_infile.F90 +++ b/src/main/readwrite_infile.F90 @@ -16,6 +16,8 @@ module readwrite_infile ! :Runtime parameters: ! - C_cour : *Courant number* ! - C_force : *dt_force number* +! - X : *hydrogen mass fraction for MESA opacity table* +! - Z : *metallicity for MESA opacity table* ! - alpha : *shock viscosity parameter* ! - alphaB : *shock resistivity parameter* ! - alphamax : *MAXIMUM shock viscosity parameter* @@ -66,9 +68,9 @@ module readwrite_infile ! ! :Dependencies: boundary_dyn, cooling, damping, dim, dust, dust_formation, ! eos, externalforces, forcing, gravwaveutils, growth, infile_utils, -! inject, io, linklist, metric, nicil_sup, options, part, ptmass, -! ptmass_radiation, radiation_implicit, radiation_utils, timestep, -! viscosity +! inject, io, linklist, metric, nicil_sup, options, part, porosity, +! ptmass, ptmass_radiation, radiation_implicit, radiation_utils, +! timestep, viscosity ! use timestep, only:dtmax_dratio,dtmax_max,dtmax_min use options, only:nfulldump,nmaxdumps,twallmax,iexternalforce,tolh, & @@ -105,6 +107,7 @@ subroutine write_infile(infile,logfile,evfile,dumpfile,iwritein,iprint) use linklist, only:write_inopts_link use dust, only:write_options_dust use growth, only:write_options_growth + use porosity, only:write_options_porosity #ifdef INJECT_PARTICLES use inject, only:write_options_inject #endif @@ -260,7 +263,10 @@ subroutine write_infile(infile,logfile,evfile,dumpfile,iwritein,iprint) #endif if (use_dust) call write_options_dust(iwritein) - if (use_dustgrowth) call write_options_growth(iwritein) + if (use_dustgrowth) then + call write_options_growth(iwritein) + call write_options_porosity(iwritein) + endif write(iwritein,"(/,a)") '# options for injecting/removing particles' #ifdef INJECT_PARTICLES @@ -323,6 +329,8 @@ subroutine read_infile(infile,logfile,evfile,dumpfile) use linklist, only:read_inopts_link use dust, only:read_options_dust use growth, only:read_options_growth + use options, only:use_porosity + use porosity, only:read_options_porosity use metric, only:read_options_metric #ifdef INJECT_PARTICLES use inject, only:read_options_inject @@ -349,7 +357,7 @@ subroutine read_infile(infile,logfile,evfile,dumpfile) integer :: ierr,ireaderr,line,idot,ngot,nlinesread real :: ratio logical :: imatch,igotallrequired,igotallturb,igotalllink,igotloops - logical :: igotallbowen,igotallcooling,igotalldust,igotallextern,igotallinject,igotallgrowth + logical :: igotallbowen,igotallcooling,igotalldust,igotallextern,igotallinject,igotallgrowth,igotallporosity logical :: igotallionise,igotallnonideal,igotalleos,igotallptmass,igotalldamping logical :: igotallprad,igotalldustform,igotallgw,igotallgr,igotallbdy integer, parameter :: nrequired = 1 @@ -366,6 +374,7 @@ subroutine read_infile(infile,logfile,evfile,dumpfile) igotallturb = .true. igotalldust = .true. igotallgrowth = .true. + igotallporosity = .true. igotalllink = .true. igotallextern = .true. igotallinject = .true. @@ -540,6 +549,7 @@ subroutine read_infile(infile,logfile,evfile,dumpfile) !--Extract if one-fluid dust is used from the fileid if (.not.imatch .and. use_dust) call read_options_dust(name,valstring,imatch,igotalldust,ierr) if (.not.imatch .and. use_dustgrowth) call read_options_growth(name,valstring,imatch,igotallgrowth,ierr) + if (.not.imatch .and. use_porosity) call read_options_porosity(name,valstring,imatch,igotallporosity,ierr) if (.not.imatch .and. gr) call read_options_metric(name,valstring,imatch,igotallgr,ierr) #ifdef INJECT_PARTICLES if (.not.imatch) call read_options_inject(name,valstring,imatch,igotallinject,ierr) @@ -575,7 +585,7 @@ subroutine read_infile(infile,logfile,evfile,dumpfile) igotallrequired = (ngot >= nrequired) .and. igotalllink .and. igotallbowen .and. igotalldust & .and. igotalleos .and. igotallcooling .and. igotallextern .and. igotallturb & .and. igotallptmass .and. igotallinject .and. igotallionise .and. igotallnonideal & - .and. igotallgrowth .and. igotalldamping .and. igotallprad & + .and. igotallgrowth .and. igotallporosity .and. igotalldamping .and. igotallprad & .and. igotalldustform .and. igotallgw .and. igotallgr .and. igotallbdy if (ierr /= 0 .or. ireaderr > 0 .or. .not.igotallrequired) then @@ -595,6 +605,8 @@ subroutine read_infile(infile,logfile,evfile,dumpfile) if (.not.igotalldust) write(*,*) 'missing dust options' if (.not.igotallgr) write(*,*) 'missing metric parameters (eg, spin, mass)' if (.not.igotallgrowth) write(*,*) 'missing growth options' + if (.not.igotallporosity) write(*,*) 'missing porosity options' + if (.not.igotallextern) then if (gr) then write(*,*) 'missing GR quantities (eg: accretion radius)' diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 70e4376dd..a211f5dd9 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -25,7 +25,7 @@ module step_lf_global ! :Dependencies: boundary_dyn, chem, cons2prim, cons2primsolver, cooling, ! cooling_ism, damping, deriv, dim, dust_formation, eos, extern_gr, ! externalforces, growth, io, io_summary, krome_interface, metric_tools, -! mpiutils, options, part, ptmass, ptmass_radiation, timestep, +! mpiutils, options, part, porosity, ptmass, ptmass_radiation, timestep, ! timestep_ind, timestep_sts, timing, units ! use dim, only:maxp,maxvxyzu,do_radiation,ind_timesteps @@ -98,7 +98,8 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) iphase,iamtype,massoftype,maxphase,igas,idust,mhd,& iamboundary,get_ntypes,npartoftypetot,& dustfrac,dustevol,ddustevol,eos_vars,alphaind,nptmass,& - dustprop,ddustprop,dustproppred,pxyzu,dens,metrics,ics + dustprop,ddustprop,dustproppred,pxyzu,dens,metrics,ics,& + filfac,filfacpred,mprev,filfacprev use options, only:avdecayconst,alpha,ieos,alphamax use deriv, only:derivs use timestep, only:dterr,bignumber,tolv @@ -118,6 +119,8 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) use cooling, only:ufloor,cooling_in_step use timing, only:increment_timer,get_timings,itimer_extf use growth, only:check_dustprop + use options, only:use_porosity + use porosity, only:get_filfac use damping, only:idamp use cons2primsolver, only:conservative2primitive,primitive2conservative use eos, only:equationofstate @@ -167,6 +170,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) !$omp shared(rad,drad,pxyzu) & !$omp shared(Bevol,dBevol,dustevol,ddustevol,use_dustfrac) & !$omp shared(dustprop,ddustprop,dustproppred,ufloor) & + !$omp shared(mprev,filfacprev,filfac,use_porosity) & !$omp shared(ibin,ibin_old,twas,timei) & !$omp firstprivate(itype) & !$omp private(i,hdti) & @@ -200,7 +204,10 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) nvfloorp = nvfloorp + 1 endif endif - + if (use_porosity) then + mprev(i) = dustprop(1,i) + filfacprev(i) = filfac(i) + endif if (itype==idust .and. use_dustgrowth) then dustprop(:,i) = dustprop(:,i) + hdti*ddustprop(:,i) endif @@ -215,7 +222,12 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) endif enddo predictor !omp end parallel do - if (use_dustgrowth) call check_dustprop(npart,dustprop(1,:)) + if (use_dustgrowth) then + if (use_porosity) then + call get_filfac(npart,xyzh,mprev,filfac,dustprop,hdti) + endif + call check_dustprop(npart,dustprop,filfac,mprev,filfacprev) + endif !---------------------------------------------------------------------- @@ -254,6 +266,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) !$omp shared(pxyzu,ppred) & !$omp shared(Bevol,dBevol,Bpred,dtsph,massoftype,iphase) & !$omp shared(dustevol,ddustprop,dustprop,dustproppred,dustfrac,ddustevol,dustpred,use_dustfrac) & +!$omp shared(filfac,filfacpred,use_porosity) & !$omp shared(alphaind,ieos,alphamax,ialphaloc) & !$omp shared(eos_vars,ufloor) & !$omp shared(twas,timei) & @@ -310,7 +323,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) nvfloorps = nvfloorps + 1 endif endif - + if (use_porosity) filfacpred(i) = filfac(i) if (use_dustgrowth .and. itype==idust) then dustproppred(:,i) = dustprop(:,i) + hdti*ddustprop(:,i) endif @@ -350,9 +363,12 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) endif enddo predict_sph !$omp end parallel do - - if (use_dustgrowth) call check_dustprop(npart,dustproppred(1,:)) - + if (use_dustgrowth) then + if (use_porosity) then + call get_filfac(npart,xyzh,dustprop(1,:),filfacpred,dustproppred,hdti) + endif + call check_dustprop(npart,dustproppred(:,:),filfacpred,dustprop(1,:),filfac) + endif ! ! recalculate all SPH forces, and new timestep ! @@ -365,7 +381,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) call derivs(1,npart,nactive,xyzh,vpred,fxyzu,fext,divcurlv,& divcurlB,Bpred,dBevol,radpred,drad,radprop,dustproppred,ddustprop,& - dustpred,ddustevol,dustfrac,eos_vars,timei,dtsph,dtnew,& + dustpred,ddustevol,filfacpred,dustfrac,eos_vars,timei,dtsph,dtnew,& ppred,dens,metrics) if (do_radiation .and. implicit_radiation) then @@ -569,7 +585,12 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) enddo corrector !$omp enddo !$omp end parallel - if (use_dustgrowth) call check_dustprop(npart,dustprop(1,:)) + if (use_dustgrowth) then + if (use_porosity) then + call get_filfac(npart,xyzh,mprev,filfac,dustprop,dtsph) + endif + call check_dustprop(npart,dustprop,filfac,mprev,filfacprev) + endif if (gr) then call check_velocity_error(errmax,p2mean,np,its,tolv,dtsph,timei,idamp,dterr,errmaxmean,converged) @@ -584,6 +605,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) !$omp shared(store_itype,vxyzu,fxyzu,vpred,iphase) & !$omp shared(Bevol,dBevol,Bpred,pxyzu,ppred) & !$omp shared(dustprop,ddustprop,dustproppred,use_dustfrac,dustevol,dustpred,ddustevol) & +!$omp shared(filfac,filfacpred,use_porosity) & !$omp shared(rad,drad,radpred) & !$omp firstprivate(itype) & !$omp schedule(static) @@ -600,6 +622,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) vpred(:,i) = vxyzu(:,i) endif if (use_dustgrowth) dustproppred(:,i) = dustprop(:,i) + if (use_porosity) filfacpred(i) = filfac(i) if (mhd) Bpred(:,i) = Bevol(:,i) if (use_dustfrac) dustpred(:,i) = dustevol(:,i) if (do_radiation) radpred(:,i) = rad(:,i) @@ -611,6 +634,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) vpred(:,i) = vxyzu(:,i) endif if (use_dustgrowth) dustproppred(:,i) = dustprop(:,i) + if (use_porosity) filfacpred(i) = filfac(i) if (mhd) Bpred(:,i) = Bevol(:,i) if (use_dustfrac) dustpred(:,i) = dustevol(:,i) if (do_radiation) radpred(:,i) = rad(:,i) @@ -635,15 +659,19 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) enddo until_converged !$omp end parallel do - if (use_dustgrowth) call check_dustprop(npart,dustprop(1,:)) - + if (use_dustgrowth) then + if (use_porosity) then + call get_filfac(npart,xyzh,mprev,filfac,dustprop,dtsph) + endif + call check_dustprop(npart,dustprop,filfac,mprev,filfacprev) + endif ! ! get new force using updated velocity: no need to recalculate density etc. ! if (gr) vpred = vxyzu ! Need primitive utherm as a guess in cons2prim call derivs(2,npart,nactive,xyzh,vpred,fxyzu,fext,divcurlv,divcurlB, & - Bpred,dBevol,radpred,drad,radprop,dustproppred,ddustprop,dustpred,ddustevol,dustfrac,& - eos_vars,timei,dtsph,dtnew,ppred,dens,metrics) + Bpred,dBevol,radpred,drad,radprop,dustproppred,ddustprop,dustpred,ddustevol,filfacpred,& + dustfrac,eos_vars,timei,dtsph,dtnew,ppred,dens,metrics) if (gr) vxyzu = vpred ! May need primitive variables elsewhere? if (do_radiation .and. implicit_radiation) then rad = radpred diff --git a/src/main/utils_dumpfiles_hdf5.f90 b/src/main/utils_dumpfiles_hdf5.f90 index cffcc3b32..581d06cd8 100644 --- a/src/main/utils_dumpfiles_hdf5.f90 +++ b/src/main/utils_dumpfiles_hdf5.f90 @@ -468,7 +468,7 @@ subroutine write_hdf5_arrays( & ! Dust growth if (array_options%use_dustgrowth) then - call write_to_hdf5(dustprop(1,1:npart), 'grainsize', group_id, error) + call write_to_hdf5(dustprop(1,1:npart), 'grainmass', group_id, error) call write_to_hdf5(dustprop(2,1:npart), 'graindens', group_id, error) call write_to_hdf5(VrelVf(1:npart), 'vrel_on_vfrag', group_id, error) call write_to_hdf5(dustgasprop(3,1:npart), 'St', group_id, error) @@ -614,7 +614,7 @@ subroutine write_hdf5_arrays_small( & ! Dustgrowth if (array_options%use_dustgrowth) then - call write_to_hdf5(real(dustprop(1,1:npart), kind=4), 'grainsize', group_id, error) + call write_to_hdf5(real(dustprop(1,1:npart), kind=4), 'grainmass', group_id, error) call write_to_hdf5(real(dustprop(2,1:npart), kind=4), 'graindens', group_id, error) endif @@ -933,7 +933,7 @@ subroutine read_hdf5_arrays( & ! Dust growth if (array_options%use_dustgrowth) then - call read_from_hdf5(dustprop(1,:), 'grainsize', group_id, got_arrays%got_dustprop(1), error) + call read_from_hdf5(dustprop(1,:), 'grainmass', group_id, got_arrays%got_dustprop(1), error) call read_from_hdf5(dustprop(2,:), 'graindens', group_id, got_arrays%got_dustprop(2), error) call read_from_hdf5(VrelVf(:), 'vrel_on_vfrag', group_id, got_arrays%got_VrelVf, error) call read_from_hdf5(dustgasprop(3,:), 'St', group_id, got_arrays%got_St, error) diff --git a/src/main/writeheader.F90 b/src/main/writeheader.F90 index 256b0dbbe..0d32e639b 100644 --- a/src/main/writeheader.F90 +++ b/src/main/writeheader.F90 @@ -80,7 +80,7 @@ subroutine write_header(icall,infile,evfile,logfile,dumpfile,ntot) use io, only:iprint use boundary, only:xmin,xmax,ymin,ymax,zmin,zmax use boundary_dyn, only:dynamic_bdy,rho_thresh_bdy,width_bkg - use options, only:tolh,alpha,alphau,alphaB,ieos,alphamax,use_dustfrac + use options, only:tolh,alpha,alphau,alphaB,ieos,alphamax,use_dustfrac,use_porosity use part, only:hfact,massoftype,mhd,gravity,periodic,massoftype,npartoftypetot,& labeltype,maxtypes use mpiutils, only:reduceall_mpi @@ -186,6 +186,7 @@ subroutine write_header(icall,infile,evfile,logfile,dumpfile,ntot) endif endif if (use_dustgrowth) write(iprint,"(1x,a)") 'Dust growth is ON' + if (use_porosity) write(iprint,"(1x,a)") 'Dust porosity is ON' if (cooling_in_step) then write(iprint,"(1x,a)") 'Cooling is calculated in step' else diff --git a/src/setup/set_dust_options.f90 b/src/setup/set_dust_options.f90 index e3d548a6b..c97bbf0ca 100644 --- a/src/setup/set_dust_options.f90 +++ b/src/setup/set_dust_options.f90 @@ -27,7 +27,7 @@ module set_dust_options ! - ndusttypesinp : *number of grain sizes* ! ! :Dependencies: dim, dust, eos, fileutils, growth, infile_utils, io, -! options, part, prompting +! options, part, porosity, prompting ! use dim, only:maxdusttypes,maxdustsmall,maxdustlarge,use_dustgrowth use prompting, only:prompt @@ -236,6 +236,7 @@ end subroutine set_log_dist_options !-------------------------------------------------------------------------- subroutine read_dust_setup_options(db,nerr) use growth, only:read_growth_setup_options + use porosity, only:read_porosity_setup_options use infile_utils, only:inopts,read_inopt use io, only:error use fileutils, only:make_tags_unique @@ -432,8 +433,10 @@ subroutine read_dust_setup_options(db,nerr) call read_inopt(graindensinp(1),'graindensinp',db,min=0.,errcount=nerr) endif - if (use_dustgrowth) call read_growth_setup_options(db,nerr) - + if (use_dustgrowth) then + call read_growth_setup_options(db,nerr) + call read_porosity_setup_options(db,nerr) + endif end subroutine read_dust_setup_options !-------------------------------------------------------------------------- @@ -517,6 +520,7 @@ end subroutine read_log_dist_options !-------------------------------------------------------------------------- subroutine write_dust_setup_options(iunit) use growth, only:write_growth_setup_options + use porosity, only:write_porosity_setup_options use infile_utils, only:write_inopt use fileutils, only:make_tags_unique @@ -714,7 +718,10 @@ subroutine write_dust_setup_options(iunit) call write_inopt(isetdust,'isetdust', & 'how to set dust density profile (0=equal to gas,1=custom,2=equal to gas with cutoffs)',iunit) - if (use_dustgrowth) call write_growth_setup_options(iunit) + if (use_dustgrowth) then + call write_growth_setup_options(iunit) + call write_porosity_setup_options(iunit) + endif end subroutine write_dust_setup_options diff --git a/src/setup/setup_disc.f90 b/src/setup/setup_disc.f90 index b1df4b584..590f8a2b4 100644 --- a/src/setup/setup_disc.f90 +++ b/src/setup/setup_disc.f90 @@ -90,9 +90,9 @@ module setup ! :Dependencies: centreofmass, dim, dust, eos, extern_binary, ! extern_corotate, extern_lensethirring, externalforces, fileutils, ! growth, infile_utils, io, kernel, memory, options, part, physcon, -! prompting, radiation_utils, set_dust, set_dust_options, setbinary, -! setdisc, setflyby, sethierarchical, spherical, timestep, units, -! vectorutils +! porosity, prompting, radiation_utils, set_dust, set_dust_options, +! setbinary, setdisc, setflyby, sethierarchical, spherical, timestep, +! units, vectorutils ! use dim, only:use_dust,maxalpha,use_dustgrowth,maxdusttypes,& maxdustlarge,maxdustsmall,compiled_with_mcfost @@ -102,15 +102,16 @@ module setup use extern_binary, only:mass2,accradius1,accradius2,ramp,surface_force,eps_soft1 use fileutils, only:make_tags_unique use growth, only:ifrag,isnow,rsnow,Tsnow,vfragSI,vfraginSI,vfragoutSI,gsizemincgs + use porosity, only:iporosity use io, only:master,warning,error,fatal use kernel, only:hfact_default - use options, only:use_dustfrac,iexternalforce,use_hybrid + use options, only:use_dustfrac,iexternalforce,use_hybrid,use_porosity use options, only:use_mcfost,use_mcfost_stellar_parameters use part, only:xyzmh_ptmass,maxvxyzu,vxyz_ptmass,ihacc,ihsoft,& iJ2,ispinx,ispinz,iReff,igas,& idust,iphase,dustprop,dustfrac,ndusttypes,ndustsmall,& ndustlarge,grainsize,graindens,nptmass,iamtype,dustgasprop,& - VrelVf,rad,radprop,ikappa,iradxi + VrelVf,filfac,probastick,rad,radprop,ikappa,iradxi use physcon, only:au,solarm,jupiterm,earthm,pi,twopi,years,hours,deg_to_rad use setdisc, only:scaled_sigma,get_disc_mass,maxbins use set_dust_options, only:set_dust_default_options,dust_method,dust_to_gas,& @@ -1607,6 +1608,7 @@ end subroutine set_planet_atm ! !-------------------------------------------------------------------------- subroutine initialise_dustprop(npart) + use physcon, only:fourpi integer, intent(in) :: npart integer :: i,iam @@ -1615,11 +1617,13 @@ subroutine initialise_dustprop(npart) do i=1,npart iam = iamtype(iphase(i)) if (iam==idust .or. (use_dustfrac .and. iam==igas)) then - dustprop(1,i) = grainsize(1) + dustprop(1,i) = fourpi/3.*graindens(1)*grainsize(1)**3 dustprop(2,i) = graindens(1) else dustprop(:,i) = 0. endif + filfac(i) = 0. + probastick(i) = 1. dustgasprop(:,i) = 0. VrelVf(i) = 0. enddo @@ -2372,6 +2376,8 @@ subroutine setup_interactive(id) call prompt('Enter outward vfragout in m/s',vfragoutSI,1.) endif endif + call prompt('Enter porosity switch (0=off,1=on)',iporosity,0,1) + if (iporosity == 1) use_porosity = .true. endif endif diff --git a/src/setup/setup_shock.F90 b/src/setup/setup_shock.f90 similarity index 100% rename from src/setup/setup_shock.F90 rename to src/setup/setup_shock.f90 diff --git a/src/setup/setup_testparticles.F90 b/src/setup/setup_testparticles.f90 similarity index 100% rename from src/setup/setup_testparticles.F90 rename to src/setup/setup_testparticles.f90 diff --git a/src/tests/test_dust.F90 b/src/tests/test_dust.f90 similarity index 94% rename from src/tests/test_dust.F90 rename to src/tests/test_dust.f90 index f27bb670a..371059769 100644 --- a/src/tests/test_dust.F90 +++ b/src/tests/test_dust.f90 @@ -45,13 +45,13 @@ subroutine test_dust(ntests,npass) use physcon, only:solarm,au use units, only:set_units,unit_density,udist use eos, only:gamma - use dim, only:use_dust + use dim, only:use_dust,use_dustgrowth use mpiutils, only:barrier_mpi use options, only:use_dustfrac use table_utils, only:logspace use growth, only:init_growth integer, intent(inout) :: ntests,npass - integer :: nfailed(3),ierr,iregime + integer :: nfailed(3),ierr,iregime,j real :: rhoi,rhogasi,rhodusti,spsoundi,tsi,grainsizei,graindensi if (use_dust) then @@ -71,10 +71,10 @@ subroutine test_dust(ntests,npass) call init_drag(ierr) call checkval(ierr,0,0,nfailed(idrag),'drag initialisation') enddo -#ifdef DUSTGROWTH - call init_growth(ierr) - call checkval(ierr,0,0,nfailed(3),'growth initialisation') -#endif + if (use_dustgrowth) then + call init_growth(ierr) + call checkval(ierr,0,0,nfailed(3),'growth initialisation') + endif call update_test_scores(ntests,nfailed,npass) idrag = 1 @@ -94,36 +94,23 @@ subroutine test_dust(ntests,npass) call test_epsteinstokes(ntests,npass) call barrier_mpi() - if (id==master) write(*,"(/,a)") '--> testing drag with EXPLICIT scheme' use_dustfrac = .false. - ! - ! Test that drag conserves momentum and energy with explicit scheme - ! - drag_implicit = .false. - call test_drag(ntests,npass) - call barrier_mpi() - - ! - ! DUSTYBOX test with explicit scheme - ! + do j=1,2 + drag_implicit = .false. + if (j==2) drag_implicit = .true. + ! + ! Test that drag conserves momentum and energy with explicit/implicit scheme + ! + call test_drag(ntests,npass) + call barrier_mpi() + + ! + ! DUSTYBOX test with explicit/implicit scheme + ! + call test_dustybox(ntests,npass) + call barrier_mpi() + enddo drag_implicit = .false. - call test_dustybox(ntests,npass) - call barrier_mpi() - - if (id==master) write(*,"(/,a)") '--> testing DRAG with IMPLICIT scheme' - ! - ! Test that drag conserves momentum and energy with implicit scheme - ! - drag_implicit = .true. - call test_drag(ntests,npass) - call barrier_mpi() - - ! - ! DUSTYBOX test with explicit scheme - ! - drag_implicit = .true. - call test_dustybox(ntests,npass) - call barrier_mpi() ! ! DUSTYDIFFUSE test @@ -152,7 +139,7 @@ subroutine test_dustybox(ntests,npass) use energies, only:compute_energies,ekin use testutils, only:checkvalbuf,checkvalbuf_end use eos, only:ieos,polyk,gamma - use dust, only:K_code,idrag + use dust, only:K_code,idrag,drag_implicit use options, only:alpha,alphamax use unifdis, only:set_unifdis use dim, only:periodic,mhd,use_dust,use_dustgrowth @@ -173,6 +160,9 @@ subroutine test_dustybox(ntests,npass) real :: t, dt, dtext, dtnew real :: vg, vd, deltav, ekin_exact, fd real :: tol,tolvg,tolfg,tolfd + logical :: write_output = .false. + character(len=60) :: filename,string + integer, parameter :: lu = 36 if (index(kernelname,'quintic') /= 0) then tol = 1.e-5; tolvg = 2.5e-5; tolfg = 3.3e-4; tolfd = 3.3e-4 @@ -181,19 +171,19 @@ subroutine test_dustybox(ntests,npass) endif if (periodic .and. use_dust) then - if (id==master) write(*,"(/,a)") '--> testing DUSTYBOX' + string = '(explicit drag)' + if (drag_implicit) string = '(implicit drag)' + if (id==master) write(*,"(/,a)") '--> testing DUSTYBOX '//trim(string) else if (id==master) write(*,"(/,a)") '--> skipping DUSTYBOX (need -DPERIODIC and -DDUST)' return endif - if (use_dustgrowth .and. id==master) write(*,"(/,a)") '--> Adding dv interpolation test' - ! ! setup for dustybox problem ! call init_part() - nx = 32 + nx = 24 deltax = 1./nx dz = 2.*sqrt(6.)/nx call set_boundary(-0.5,0.5,-0.25,0.25,-dz,dz) @@ -272,6 +262,12 @@ subroutine test_dustybox(ntests,npass) vg = 0.5*(1. - deltav) vd = 0.5*(1. + deltav) fd = K_code(1)*(vg - vd) + if (write_output) then + write(filename,"(a,1pe8.2,a)") 'dustybox_t',t,'.out' + open(unit=lu,file=filename,status='replace') + print "(a)",' writing '//trim(filename) + endif + do j=1,npart if (iamdust(iphase(j))) then call checkvalbuf(vxyzu(1,j),vd,tol,'vd',nerr(1),ncheck(1),errmax(1)) @@ -282,8 +278,11 @@ subroutine test_dustybox(ntests,npass) else call checkvalbuf(vxyzu(1,j),vg,tolvg,'vg',nerr(3),ncheck(3),errmax(3)) call checkvalbuf(fxyzu(1,j),-fd,tolfg,'fg',nerr(4),ncheck(4),errmax(4)) + if (write_output) write(lu,*) vxyzu(1,j),fxyzu(1,j),vg,-fd endif enddo + if (write_output) close(lu) + !call checkval(npart/2-1,vxyzu(1,1:npart),vg,tolvg,nerr(2),'vg') ekin_exact = 0.5*totmass*(vd**2 + vg**2) !print*,' step ',i,'t = ',t,' ekin should be ',ekin_exact, ' got ',ekin,(ekin-ekin_exact)/ekin_exact @@ -507,7 +506,7 @@ subroutine test_drag(ntests,npass) use options, only:use_dustfrac use eos, only:polyk,ieos use kernel, only:hfact_default - use dust, only:K_code,idrag + use dust, only:K_code,idrag,drag_implicit use boundary, only:dxbound,dybound,dzbound,xmin,xmax,ymin,ymax,zmin,zmax,set_boundary use io, only:iverbose use unifdis, only:set_unifdis @@ -517,21 +516,25 @@ subroutine test_drag(ntests,npass) use vectorutils, only:cross_product3D use units, only:udist,unit_density use mpidomain, only:i_belong + use physcon, only:pi integer, intent(inout) :: ntests,npass integer(kind=8) :: npartoftypetot(maxtypes) integer :: nx,i,j,nfailed(7),itype,iseed,npart_previous,iu real :: da(3),dl(3),temp(3) real :: psep,time,rhozero,totmass,dekin,deint + character(len=10) :: string real, parameter :: tol_mom = 1.e-7 real, parameter :: tol_ang = 5.e-4 real, parameter :: tol_enj = 1.e-6 - if (id==master) write(*,"(/,a)") '--> testing DUST DRAG' + string = '(explicit)' + if (drag_implicit) string = '(implicit)' + if (id==master) write(*,"(/,a)") '--> testing DUST DRAG '//trim(string) ! ! set up particles in random distribution ! - nx = 50 + nx = 25 psep = 1./nx iseed= -14255 call init_part() @@ -583,7 +586,7 @@ subroutine test_drag(ntests,npass) if (use_dustgrowth) then dustprop(:,:) = 0. - dustprop(1,:) = grainsize(1) + dustprop(1,:) = 4./3.*pi*grainsize(1)**3*graindens(1) dustprop(2,:) = graindens(1) endif ! diff --git a/src/tests/test_growth.f90 b/src/tests/test_growth.f90 index 68dd5391a..01190a85b 100644 --- a/src/tests/test_growth.f90 +++ b/src/tests/test_growth.f90 @@ -104,14 +104,14 @@ 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,gsizemincgs + use growth, only:ifrag,init_growth,isnow,vfrag,gsizemincgs,get_size use options, only:alpha,alphamax,use_dustfrac use unifdis, only:set_unifdis use dim, only:periodic,mhd,use_dust,maxp,maxalpha use timestep, only:dtmax use io, only:iverbose use mpiutils, only:reduceall_mpi - use physcon, only:au,solarm,Ro,pi + use physcon, only:au,solarm,Ro,pi,fourpi use viscosity, only:shearparam use units, only:set_units,udist,unit_density!,unit_velocity use mpidomain, only:i_belong @@ -151,7 +151,7 @@ subroutine test_farmingbox(ntests,npass,frag,onefluid) dtgratio = 0.5 stringfrag = "fragmentation" else - sinit = 1.e-2/udist + sinit = 3.e-2/udist dtgratio = 1. stringfrag = "growth" endif @@ -219,7 +219,7 @@ subroutine test_farmingbox(ntests,npass,frag,onefluid) VrelVf(i) = 0. if (use_dustfrac) then dustfrac(1,i) = dtgratio - dustprop(1,i) = sinit + dustprop(1,i) = fourpi/3.*dens*sinit**3 dustprop(2,i) = dens else dustprop(:,i) = 0. @@ -245,7 +245,7 @@ subroutine test_farmingbox(ntests,npass,frag,onefluid) if (use_dust) then dustevol(:,i) = 0. dustfrac(:,i) = 0. - dustprop(1,i) = sinit + dustprop(1,i) = fourpi/3.*dens*sinit**3 dustprop(2,i) = dens dustgasprop(:,i) = 0. VrelVf(i) = 0. @@ -344,10 +344,10 @@ subroutine test_farmingbox(ntests,npass,frag,onefluid) s(j) = Stcomp(j)/(sqrt(pi*gamma/8)*dens/((rhog+rhod)*cscomp(j))*Omega_k(j)) if (onefluid) then call checkvalbuf(dustgasprop(3,j)/Stcomp(j),1.,tolst,'St',nerr(1),ncheck(1),errmax(1)) - call checkvalbuf(dustprop(1,j)/s(j),1.,tols,'size',nerr(2),ncheck(2),errmax(2)) + call checkvalbuf(get_size(dustprop(1,j),dustprop(2,j))/s(j),1.,tols,'size',nerr(2),ncheck(2),errmax(2)) else call checkvalbuf(dustgasprop(3,j)/Stcomp(j),1.,tolst,'St',nerr(1),ncheck(1),errmax(1)) - call checkvalbuf(dustprop(1,j)/s(j),1.,tols,'size',nerr(2),ncheck(2),errmax(2)) + call checkvalbuf(get_size(dustprop(1,j),dustprop(2,j))/s(j),1.,tols,'size',nerr(2),ncheck(2),errmax(2)) call checkvalbuf(dustgasprop(1,j)/cscomp(j),1.,tolcs,'csound',nerr(3),ncheck(3),errmax(3)) call checkvalbuf(dustgasprop(2,j)/rhozero,1.,tolrho,'rhogas',nerr(4),ncheck(4),errmax(4)) endif @@ -356,11 +356,11 @@ subroutine test_farmingbox(ntests,npass,frag,onefluid) enddo if (onefluid) then call checkvalbuf_end('Stokes number evaluation matches exact solution',ncheck(1),nerr(1),errmax(1),tolst) - call checkvalbuf_end('size evaluation matches exact solution',ncheck(2),nerr(2),errmax(2),tolcs) + call checkvalbuf_end('size evaluation matches exact solution',ncheck(2),nerr(2),errmax(2),tols) else call checkvalbuf_end('Stokes number interpolation matches exact solution',ncheck(1),nerr(1),errmax(1),tolst) - call checkvalbuf_end('size evaluation matches exact solution',ncheck(2),nerr(2),errmax(2),tolcs) - call checkvalbuf_end('sound speed interpolation matches exact number',ncheck(3),nerr(3),errmax(3),tols) + call checkvalbuf_end('size evaluation matches exact solution',ncheck(2),nerr(2),errmax(2),tols) + call checkvalbuf_end('sound speed interpolation matches exact number',ncheck(3),nerr(3),errmax(3),tolcs) call checkvalbuf_end('rhogas interpolation matches exact number',ncheck(4),nerr(4),errmax(4),tolrho) endif diff --git a/src/utils/interpolate3D.F90 b/src/utils/interpolate3D.f90 similarity index 98% rename from src/utils/interpolate3D.F90 rename to src/utils/interpolate3D.f90 index 1a6d0d75e..95fe2d7d6 100644 --- a/src/utils/interpolate3D.F90 +++ b/src/utils/interpolate3D.f90 @@ -17,7 +17,7 @@ module interpolations3D ! :Dependencies: einsteintk_utils, kernel ! use einsteintk_utils, only:exact_rendering - use kernel, only:radkern2,radkern,cnormk,wkern + use kernel, only:radkern2,radkern,cnormk,wkern implicit none integer, parameter :: doub_prec = kind(0.d0) real :: cnormk3D = cnormk @@ -990,12 +990,12 @@ pure elemental real function soft_func(x,eps) result(f) end function soft_func - !-------------------------------------------------------------------------- - ! - ! utility to wrap pixel index around periodic domain - ! indices that roll beyond the last position are re-introduced at the first - ! - !-------------------------------------------------------------------------- +!-------------------------------------------------------------------------- +! +! utility to wrap pixel index around periodic domain +! indices that roll beyond the last position are re-introduced at the first +! +!-------------------------------------------------------------------------- pure integer function iroll(i,n) integer, intent(in) :: i,n @@ -1008,5 +1008,5 @@ pure integer function iroll(i,n) endif end function iroll -end module interpolations3D +end module interpolations3D diff --git a/src/utils/moddump_LTE_to_rad.f90 b/src/utils/moddump_LTE_to_rad.f90 index 04b20e5c9..58e9a135a 100644 --- a/src/utils/moddump_LTE_to_rad.f90 +++ b/src/utils/moddump_LTE_to_rad.f90 @@ -14,7 +14,8 @@ module moddump ! ! :Runtime parameters: None ! -! :Dependencies: dim, eos, io, part +! :Dependencies: dim, eos, eos_idealplusrad, eos_mesa, io, +! mesa_microphysics, part, radiation_utils, units ! implicit none diff --git a/src/utils/moddump_dustadd.f90 b/src/utils/moddump_dustadd.f90 index 4df9d4c2a..588a74707 100644 --- a/src/utils/moddump_dustadd.f90 +++ b/src/utils/moddump_dustadd.f90 @@ -10,13 +10,17 @@ module moddump ! ! :References: None ! -! :Owner: Arnaud Vericel +! :Owner: Stephane Michoulier ! ! :Runtime parameters: None ! -! :Dependencies: dim, dust, growth, options, part, prompting, set_dust, -! table_utils, units +! :Dependencies: dim, dust, growth, options, part, porosity, prompting, +! set_dust, table_utils, units ! + + use part, only:delete_particles_outside_sphere,igas,idust + use prompting, only:prompt + implicit none contains @@ -27,8 +31,9 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) use part, only:igas,idust,set_particle_type,ndusttypes,ndustsmall,ndustlarge,& grainsize,graindens,dustfrac use set_dust, only:set_dustfrac,set_dustbinfrac - use options, only:use_dustfrac + use options, only:use_dustfrac,use_porosity use growth, only:set_dustprop,convert_to_twofluid + use porosity, only:iporosity use prompting, only:prompt use dust, only:grainsizecgs,graindenscgs use table_utils, only:logspace @@ -37,8 +42,13 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) integer, intent(inout) :: npartoftype(:) real, intent(inout) :: massoftype(:) real, intent(inout) :: xyzh(:,:),vxyzu(:,:) + real, dimension(3) :: incenter,outcenter integer :: i,j,itype,ipart,iloc,dust_method,np_ratio,np_gas,np_dust,maxdust real :: dust_to_gas,smincgs,smaxcgs,sindex,dustbinfrac(maxdusttypes),udens + integer :: iremoveparttype + real :: inradius,outradius,pwl_sizedistrib,R_ref,H_R_ref,q_index + logical :: icutinside,icutoutside,sizedistrib + if (.not. use_dust) then print*,' DOING NOTHING: COMPILE WITH DUST=yes' @@ -53,6 +63,19 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) smaxcgs = 1. sindex = 3.5 dustbinfrac = 0. + pwl_sizedistrib = -2 + R_ref = 100 + H_R_ref = 0.0895 + q_index = 0.25 + + icutinside = .false. + icutoutside = .false. + iremoveparttype = 0 + incenter(:) = 0. + outcenter(:) = 0. + inradius = 10. + outradius = 200. + !- grainsize and graindens already set if convert from one fluid to two fluid with growth if (.not. (use_dustfrac .and. use_dustgrowth)) then grainsize = 1. @@ -95,7 +118,20 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) graindens = graindens(1)/udens else if (use_dustgrowth) then - call prompt('Enter initial grain size in cm',grainsizecgs,0.) + call prompt('Use porosity ? (0=no,1=yes)',iporosity,0,1) + if (iporosity == 1) then + use_porosity = .true. + endif + call prompt('Set dust size via size distribution ?',sizedistrib) + if (sizedistrib) then + call prompt('Enter grain size in cm at Rref',grainsizecgs,0.) + call prompt('Enter power-law index ',pwl_sizedistrib) + call prompt('Enter R_ref ',R_ref,0.) + call prompt('Enter H/R at R_ref',H_R_ref,0.) + call prompt('Enter q index',q_index) + else + call prompt('Enter initial grain size in cm',grainsizecgs,0.) + endif else call prompt('Enter grain size in cm',grainsizecgs,0.) endif @@ -155,9 +191,62 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) enddo endif - if (use_dustgrowth) then - call set_dustprop(npart) + call set_dustprop(npart,xyzh,sizedistrib,pwl_sizedistrib,R_ref,H_R_ref,q_index) + endif + endif + !Delete particles if necessary + + ! + !--set the centers and the radius + ! + call prompt('Deleting particles inside a given radius ?',icutinside) + call prompt('Deleting particles outside a given radius ?',icutoutside) + if (icutinside) then + call prompt('Enter inward radius in au',inradius,0.) + call prompt('Enter x coordinate of the center of that sphere',incenter(1)) + call prompt('Enter y coordinate of the center of that sphere',incenter(2)) + call prompt('Enter z coordinate of the center of that sphere',incenter(3)) + endif + if (icutoutside) then + call prompt('Enter outward radius in au',outradius,0.) + call prompt('Enter x coordinate of the center of that sphere',outcenter(1)) + call prompt('Enter y coordinate of the center of that sphere',outcenter(2)) + call prompt('Enter z coordinate of the center of that sphere',outcenter(3)) + endif + + if (icutinside .or. icutoutside) then + call prompt('Deleting which particles (0=all, 1=gas only, 2=dust only)?', iremoveparttype) + ! add other types of particles here if needed + select case (iremoveparttype) + case (1) + iremoveparttype = igas + case (2) + iremoveparttype = idust + case default + iremoveparttype = 0 + end select + endif + + if (icutinside) then + print*,'Phantommoddump: Remove particles inside a particular radius' + print*,'Removing particles inside radius ',inradius + if (iremoveparttype > 0) then + print*,'Removing particles type ',iremoveparttype + call delete_particles_outside_sphere(incenter,inradius,npart,revert=.true.,mytype=iremoveparttype) + else + call delete_particles_outside_sphere(incenter,inradius,npart,revert=.true.) + endif + endif + + if (icutoutside) then + print*,'Phantommoddump: Remove particles outside a particular radius' + print*,'Removing particles outside radius ',outradius + if (iremoveparttype > 0) then + print*,'Removing particles type ',iremoveparttype + call delete_particles_outside_sphere(outcenter,outradius,npart,mytype=iremoveparttype) + else + call delete_particles_outside_sphere(outcenter,outradius,npart) endif endif diff --git a/src/utils/moddump_removeparticles_radius.f90 b/src/utils/moddump_removeparticles_radius.f90 index d123b0068..d9bbd3e94 100644 --- a/src/utils/moddump_removeparticles_radius.f90 +++ b/src/utils/moddump_removeparticles_radius.f90 @@ -10,14 +10,14 @@ module moddump ! ! :References: None ! -! :Owner: Arnaud Vericel +! :Owner: Daniel Mentiplay ! ! :Runtime parameters: None ! ! :Dependencies: part, prompting ! - use part, only:delete_particles_outside_sphere + use part, only:delete_particles_outside_sphere,igas,idust use prompting, only:prompt implicit none @@ -31,15 +31,10 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) real, dimension(:), intent(inout) :: massoftype real, dimension(:,:), intent(inout) :: xyzh,vxyzu real, dimension(3) :: incenter,outcenter + integer :: iremoveparttype real :: inradius,outradius logical :: icutinside,icutoutside - icutinside = .false. - icutoutside = .false. - incenter(:) = 0. - outcenter(:) = 0. - inradius = 10. - outradius = 200. ! !--set the centers and the radius @@ -59,16 +54,37 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) call prompt('Enter z coordinate of the center of that sphere',outcenter(3)) endif + call prompt('Deleting which particles (0=all, 1=gas only, 2=dust only)?', iremoveparttype) + ! add other types of particles here if needed + select case (iremoveparttype) + case (1) + iremoveparttype = igas + case (2) + iremoveparttype = idust + case default + iremoveparttype = 0 + end select + if (icutinside) then print*,'Phantommoddump: Remove particles inside a particular radius' print*,'Removing particles inside radius ',inradius - call delete_particles_outside_sphere(incenter,inradius,revert=.true.) + if (iremoveparttype > 0) then + print*,'Removing particles type ',iremoveparttype + call delete_particles_outside_sphere(incenter,inradius,npart,revert=.true.,mytype=iremoveparttype) + else + call delete_particles_outside_sphere(incenter,inradius,npart,revert=.true.) + endif endif if (icutoutside) then print*,'Phantommoddump: Remove particles outside a particular radius' print*,'Removing particles outside radius ',outradius - call delete_particles_outside_sphere(outcenter,outradius) + if (iremoveparttype > 0) then + print*,'Removing particles type ',iremoveparttype + call delete_particles_outside_sphere(outcenter,outradius,npart,mytype=iremoveparttype) + else + call delete_particles_outside_sphere(outcenter,outradius,npart) + endif endif end subroutine modify_dump diff --git a/src/utils/struct_part.f90 b/src/utils/struct_part.f90 index 99640148d..781a3c2fd 100644 --- a/src/utils/struct_part.f90 +++ b/src/utils/struct_part.f90 @@ -149,10 +149,10 @@ subroutine get_structure_fn(sf,nbins,norder,distmin,distmax,xbins,ncount,npart,x !$omp reduction(+:sf) do ipt=1,npts !$ if (.false.) then - if (mod(ipt,100)==0) then - call cpu_time(tcpu2) - print*,' ipt = ',ipt,tcpu2-tcpu1 - endif + if (mod(ipt,100)==0) then + call cpu_time(tcpu2) + print*,' ipt = ',ipt,tcpu2-tcpu1 + endif !$ endif i = list(ipt) xpt(1) = xyz(1,i)