Merge pull request #526 from danieljprice/StephaneMichoulier-master

porosity evolution module

AUTHORS

@@ -22,17 +22,17 @@ Elisabeth Borchert <[email protected]>
 Ward Homan <[email protected]>
 Christophe Pinte <[email protected]>
 Terrence Tricco <[email protected]>
-Simone Ceppi <[email protected]>
 Stephane Michoulier <[email protected]>
+Simone Ceppi <[email protected]>
 Spencer Magnall <[email protected]>
 Caitlyn Hardiman <[email protected]>
 Enrico Ragusa <[email protected]>
 Sergei Biriukov <[email protected]>
 Cristiano Longarini <[email protected]>
 Giovanni Dipierro <[email protected]>
 Roberto Iaconi <[email protected]>
-Amena Faruqi <[email protected]>
 Hauke Worpel <[email protected]>
+Amena Faruqi <[email protected]>
 Alison Young <[email protected]>
 Stephen Neilson <[email protected]>
 Martina Toscani <[email protected]>
@@ -49,21 +49,23 @@ Phantom benchmark bot <[email protected]>
 Kieran Hirsh <[email protected]>
 Nicole Rodrigues <[email protected]>
 David Trevascus <[email protected]>
-Farzana Meru <[email protected]>
 Nicolás Cuello <[email protected]>
+Farzana Meru <[email protected]>
+Mike Lau <[email protected]>
 Chris Nixon <[email protected]>
 Miguel Gonzalez-Bolivar <[email protected]>
-Benoit Commercon <[email protected]>
-Giulia Ballabio <[email protected]>
-Joe Fisher <[email protected]>
-Maxime Lombart <[email protected]>
-Mike Lau <[email protected]>
 Orsola De Marco <[email protected]>
+Maxime Lombart <[email protected]>
+Joe Fisher <[email protected]>
 Zachary Pellow <[email protected]>
+Benoit Commercon <[email protected]>
+Giulia Ballabio <[email protected]>
 s-neilson <[email protected]>
-Cox, Samuel <[email protected]>
+MICHOULIER Stephane <smichoulier@sidus11>
+Steven Rieder <[email protected]>
 Jeremy Smallwood <[email protected]>
+Cox, Samuel <[email protected]>
 Jorge Cuadra <[email protected]>
-Steven Rieder <[email protected]>
 Stéven Toupin <[email protected]>
 Taj Jankovič <[email protected]>
+Chunliang Mu <[email protected]>

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 \

build/Makefile_setups

@@ -420,15 +420,15 @@ endif
 ifeq ($(SETUP), shock)
 #   shock tube tests
     PERIODIC=yes
-    SETUPFILE= setup_shock.F90
+    SETUPFILE= setup_shock.f90
     KERNEL=quintic
     KNOWN_SETUP=yes
 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,15 +460,15 @@ 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
 
 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

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

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

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