Merge pull request #486 from lsiess/master

Improvements to cooling solver

.mailmap

@@ -73,10 +73,15 @@ Enrico Ragusa <enr.ragusa@gmail.com> <enrico.ragusa@gaspare.fisica.unimi.it>
 Enrico Ragusa <enr.ragusa@gmail.com> <er198@leicester.ac.uk>
 Kieran Hirsh <kieran.hirsh1@monash.edu> <kahir1@student.monash.edu>
 Giulia Ballabio <giulia.ballabio2@studenti.unimi.it> Giulia Ballabio <giuliaballabio@Giulias-MacBook-Pro.local>
+Mats Esseldeurs <mats.esseldeurs@kuleuven.be>
+Mats Esseldeurs <mats.esseldeurs@kuleuven.be> <matsesseldeurs@yahoo.com>
+Mats Esseldeurs <mats.esseldeurs@kuleuven.be> <matse@naos.ster.kuleuven.be>
 Lionel Siess <lionel.siess@astro.ulb.ac.be>
 Lionel Siess <lionel.siess@astro.ulb.ac.be> <Lionel.Siess@ulb.ac.be>
 Lionel Siess <lionel.siess@astro.ulb.ac.be> <Lionel.Siess@ulb.be>
 Lionel Siess <lionel.siess@astro.ulb.ac.be> <Lionel.Siess@astro.ulb.ac.be>
+Mats Esseldeurs <mats.esseldeurs@kuleuven.be>
+Mats Esseldeurs <mats.esseldeurs@kuleuven.be> <matsesseldeurs@yahoo.com>
 David Liptai <dliptai@swin.edu.au> <dliptai@hotmail.com>
 David Liptai <dliptai@swin.edu.au> <david.liptai@monash.edu>
 David Liptai <dliptai@swin.edu.au> <31463304+dliptai@users.noreply.github.com>

AUTHORS

@@ -16,13 +16,14 @@ Daniel Mentiplay <daniel.mentiplay@protonmail.com>
 Megha Sharma <msha0023@student.monash.edu>
 Arnaud Vericel <arnaud.vericel@univ-lyon1.fr>
 Mark Hutchison <markahutch@gmail.com>
-Mats Esseldeurs <matsesseldeurs@yahoo.com>
 Rebecca Nealon <rebecca.nealon@warwick.ac.uk>
 Elisabeth Borchert <elisabeth.borchert@monash.edu>
 Ward Homan <ward.homan@kuleuven.be>
 Christophe Pinte <christophe.pinte@univ-grenoble-alpes.fr>
 Terrence Tricco <tstricco@mun.ca>
 Simone Ceppi <simo.ceppi@gmail.com>
+Mats Esseldeurs <mats.esseldeurs@kuleuven.be>
+Mats Esseldeurs <matsesseldeurs@yahoo.com>
 Stephane Michoulier <stephane.michoulier@gmail.com>
 Spencer Magnall <spencermagnall@gmail.com>
 Caitlyn Hardiman <caitlyn.hardiman1@monash.edu>
@@ -31,36 +32,37 @@ Sergei Biriukov <svbiriukov@gmail.com>
 Cristiano Longarini <cristiano.longarini@unimi.it>
 Giovanni Dipierro <giovanni.dipierro@leicester.ac.uk>
 Roberto Iaconi <robertoiaconi1@gmail.com>
-Hauke Worpel <hworpel@aip.de>
 Amena Faruqi <Amena.Faruqi@warwick.ac.uk>
+Hauke Worpel <hworpel@aip.de>
 Alison Young <alison.young@ed.ac.uk>
 Stephen Neilson <36410751+s-neilson@users.noreply.github.com>
 Martina Toscani <mtoscani94@gmail.com>
 Benedetta Veronesi <benedetta.veronesi@unimi.it>
 Sahl Rowther <sahl.rowther@leicester.ac.uk>
-Thomas Reichardt <thomas.reichardt@students.mq.edu.au>
 Simon Glover <glover@uni-heidelberg.de>
+Thomas Reichardt <thomas.reichardt@students.mq.edu.au>
 Jean-François Gonzalez <Jean-Francois.Gonzalez@ens-lyon.fr>
 Christopher Russell <crussell@udel.edu>
-Phantom benchmark bot <ubuntu@phantom-benchmarks.novalocal>
-Jolien Malfait <jolien.malfait@kuleuven.be>
 Alessia Franchini <alessia.franchini@unlv.edu>
 Alex Pettitt <alex@astro1.sci.hokudai.ac.jp>
-Nicole Rodrigues <nicole.rodrigues@monash.edu>
+Jolien Malfait <jolien.malfait@kuleuven.be>
+Phantom benchmark bot <ubuntu@phantom-benchmarks.novalocal>
 Kieran Hirsh <kieran.hirsh1@monash.edu>
-Nicolás Cuello <cuello.nicolas@gmail.com>
-Farzana Meru <f.meru@warwick.ac.uk>
+Nicole Rodrigues <nicole.rodrigues@monash.edu>
 David Trevascus <dtre10@student.monash.edu>
+Farzana Meru <f.meru@warwick.ac.uk>
+Nicolás Cuello <cuello.nicolas@gmail.com>
 Chris Nixon <cjn@leicester.ac.uk>
-Giulia Ballabio <giulia.ballabio2@studenti.unimi.it>
 Miguel Gonzalez-Bolivar <miguelgb@astro.unam.mx>
-Maxime Lombart <maxime.lombart@ens-lyon.fr>
 Benoit Commercon <benoit.commercon@gmail.com>
-Orsola De Marco <orsola.demarco@mq.edu.au>
+Giulia Ballabio <giulia.ballabio2@studenti.unimi.it>
 Joe Fisher <jwfis1@student.monash.edu>
-s-neilson <36410751+s-neilson@users.noreply.github.com>
+Maxime Lombart <maxime.lombart@ens-lyon.fr>
+Mike Lau <mike.lau@h-its.org>
+Orsola De Marco <orsola.demarco@mq.edu.au>
 Zachary Pellow <zpel1@student.monash.edu>
+s-neilson <36410751+s-neilson@users.noreply.github.com>
 Cox, Samuel <sc676@leicester.ac.uk>
+Jorge Cuadra <jcuadra@astro.puc.cl>
 Steven Rieder <steven@rieder.nl>
 Stéven Toupin <steven.toupin@gmail.com>
-Jorge Cuadra <jcuadra@astro.puc.cl>

build/Makefile

@@ -244,10 +244,6 @@ ifeq ($(NONIDEALMHD), yes)
     FPPFLAGS += -DNONIDEALMHD
 endif
 
-ifeq ($(H2CHEM), yes)
-    FPPFLAGS += -DH2CHEM
-endif
-
 ifeq ($(DISC_VISCOSITY), yes)
     FPPFLAGS += -DDISC_VISCOSITY
 endif
@@ -623,14 +619,14 @@ SRCDUMP= physcon.f90 ${CONFIG} ${SRCKERNEL} utils_omp.F90 io.F90 units.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 \
+         timestep.f90 ${SRCEOS} cullendehnen.f90 dust_formation.f90 \
          ${SRCGR} ${SRCPOT} \
          memory.F90 \
          utils_sphNG.f90 \
          setup_params.f90 ${SRCFASTMATH} checkoptions.F90 \
          viscosity.f90 damping.f90 options.f90 checkconserved.f90 prompting.f90 ${SRCDUST} \
          ${SRCREADWRITE_DUMPS} \
-         utils_sort.f90 sort_particles.F90
+         utils_sort.f90 sort_particles.f90
 OBJDUMP1= $(SRCDUMP:.f90=.o)
 OBJDUMP= $(OBJDUMP1:.F90=.o)
 
@@ -1328,7 +1324,7 @@ getdims:
 	@echo $(MAXP)
 
 get_setup_opts:
-	@echo "${GR:yes=GR} ${METRIC} ${MHD:yes=MHD} ${NONIDEALMHD:yes=non-ideal} ${DUST:yes=dust} ${GRAVITY:yes=self-gravity} ${RADIATION:yes=radiation}  ${H2CHEM:yes=H2_Chemistry} ${DISC_VISCOSITY:yes=disc_viscosity} ${ISOTHERMAL:yes=isothermal} ${PERIODIC:yes=periodic}" | xargs | sed -e 's/ /, /g' -e 's/_/ /g'
+	@echo "${GR:yes=GR} ${METRIC} ${MHD:yes=MHD} ${NONIDEALMHD:yes=non-ideal} ${DUST:yes=dust} ${GRAVITY:yes=self-gravity} ${RADIATION:yes=radiation} ${DISC_VISCOSITY:yes=disc_viscosity} ${ISOTHERMAL:yes=isothermal} ${PERIODIC:yes=periodic}" | xargs | sed -e 's/ /, /g' -e 's/_/ /g'
 
 get_setup_file:
 	@echo "$(SETUPFILE)"

build/Makefile_defaults_ifort

@@ -15,8 +15,8 @@ KNOWN_SYSTEM=yes
 
 # for ifort version 18+ -openmp flag is obsolete
 IFORT_VERSION_MAJOR=${shell ifort -v 2>&1 | head -1 | cut -d' ' -f 3 | cut -d'.' -f 1}
-ifeq ($(shell [ $(IFORT_VERSION_MAJOR) -gt 17 ] && echo true),true)
-    OMPFLAGS= -qopenmp
+ifeq ($(shell [ $(IFORT_VERSION_MAJOR) -lt 17 ] && echo true),true)
+    OMPFLAGS= -openmp
 else
-    OMPFLAGS = -openmp
+    OMPFLAGS = -qopenmp
 endif

build/Makefile_setups

@@ -775,7 +775,7 @@ ifeq ($(SETUP), dustystar)
     FPPFLAGS= -DDUST_NUCLEATION -DSTAR
     SETUPFILE= setup_star.f90
     MODFILE= utils_binary.f90 set_binary.f90 moddump_binary.f90
-    ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.f90 dust_formation.F90
+    ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.f90 dust_formation.f90
     KNOWN_SETUP=yes
     MAXP=10000000
     GRAVITY=yes

docs/examples.rst

@@ -17,3 +17,4 @@ This section contains some examples of physical problems that you can solve with
    density
    hierarchicalsystems
    selfgravity_gravitationalinstability
+   wind

docs/wind.rst

@@ -2,7 +2,7 @@
 Running a simulation with stellar wind and dust formation
 =========================================================
 
-The wind and dust formation algorithms are described in `Siess et al. (2022, in prep)`.
+The wind and dust formation algorithms are described in `Siess et al. (2022)`, and algortihms for the radiation field in `Esseldeurs et al. (2023)`
 
 If you find a bug, please send me an email at lionel.siess@ulb.be
 
@@ -50,12 +50,13 @@ Content of the .setup file
 
 The .setup file contains the stellar properties and sets the mass of the particle (see however  ``iwind_resolution``).
 Each star is considered as a sink particles and its properties, e.g. its luminosity, will be used to calculate the radiation pressure.
+Companions can be added using the icompanion_star parameter.
 
 Note also that 
 
-::
+.. math::
 
-      primary_lum = 4*pi*primary_Reff**2*sigma*primary_Teff**4 
+      \textrm{primary_lum} = 4\pi\times\textrm{primary_Reff}^2\times\sigma\times\textrm{primary_Teff}^4 
       
 so you only need to provide 2 out of these 3 variables. 
 
@@ -69,6 +70,9 @@ so you only need to provide 2 out of these 3 variables.
 Content of the .in file
 -----------------------
 
+Options controlling particle injection
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
 ::
 
   # options controlling particle injection
@@ -83,7 +87,7 @@ Content of the .in file
       iboundary_spheres =           5    ! number of boundary spheres (integer)
          outer_boundary =         50.    ! delete gas particles outside this radius (au)
 
-Here’s a brief description of each of them (remember that technical details can be found in `Siess et al. (in prep)
+Here’s a brief description of each of them (remember that technical details can be found in `Siess et al. (2023)`
 
 ::
 
@@ -150,6 +154,10 @@ set the number of shells that serve as inner boundary condition for the wind
 To limit the number of particles, delete from the memory the particles that go beyond ``outer_boundary`` (in astronomical unit).
 This option is slightly different from ``rkill`` where in this case the particles are declared dead and remained allocated.
 
+
+Options controlling dust
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
 ::
 
    # options controlling dust
@@ -175,12 +183,17 @@ default gas opacity. Only activated if ``idust_opacity > 0``
 
 set the C/O ratio of the ejected wind material. For the moment only C-rich chemistry (C/O > 1) is implemented. Option only available with ``idust_opacity = 2``
 
+
+Options controlling radiation pressure from sink particles
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
 ::
 
    # options controling radiation pressure from sink particles
         isink_radiation =           3    ! sink radiation pressure method (0=off,1=alpha,2=dust,3=alpha+dust)
               alpha_rad =       1.000    ! fraction of the gravitational acceleration imparted to the gas
-             iget_tdust =           1    ! dust temperature (0:Tdust=Tgas 1:T(r) 2:Lucy (devel)
+             iget_tdust =           1    ! dust temperature (0:Tdust=Tgas 1:T(r) 2:Flux dilution 3:Lucy 4:MCfost)
+        iray_resolution =          -1    ! set the number of rays to 12*4**iray_resolution (deactivated if <0)
               tdust_exp =         0.5    ! exponent of the dust temperature profile
 
 ::
@@ -189,10 +202,12 @@ set the C/O ratio of the ejected wind material. For the moment only C-rich chemi
 
 set how radiation pressure is accounted for. The star's effective gravity is given by
 
-              g = Gm/r**2 *(1-alpha_rad-Gamma)
+.. math::
+
+              g_\mathrm{eff} = \frac{Gm}{r^2} \times (1-\alpha_\mathrm{rad}-\Gamma)
 
 alpha is an ad-hoc parameter that allows the launching of the wind in case of a cool wind for example when dust is not accounted for.
-Gamma = is the Eddington factor that depends on the dust opacity. gamma is therefore <> 0 only when nucleation is activated (``idust_opacity = 2``)
+Gamma is the Eddington factor that depends on the dust opacity. gamma is therefore <> 0 only when dust is activated (``idust_opacity > 0``)
 
 ::
 
@@ -202,19 +217,30 @@ parameter entering in the above equation for the effective gravity
 
 ::
 
-             iget_tdust =           1    ! dust temperature (0:Tdust=Tgas 1:T(r) 2:Lucy (devel))
+             iget_tdust =           1    ! dust temperature (0:Tdust=Tgas 1:T(r) 2:Flux dilution 3:Lucy 4:MCfost)
 
-defines how the dust temperature is calculated. By default one assumes Tdust = Tgas but option (1, under development!) should be available soon.
+defines how the dust temperature is calculated. By default one assumes Tdust = Tgas but other options are availabe as well. 
+Options 1-3 use analytical prescriptions, and option 4 uses full 3D RT using the MCfost code (under development!)
+
+::
+
+        iray_resolution =          -1    ! set the number of rays to 12*4**iray_resolution (deactivated if <0)
+
+If ``iget_tdust = 1-3``, the dust temperature profile is then given by an analytical prescription.
+In these prescriptions (see `Esseldeurs et al. (2023)`), there is directional dependance, where the resolution of this directional dependance is set by iray_resolution.
 
 ::
 
               tdust_exp =         0.5    ! exponent of the dust temperature profile
 
 If ``iget_tdust = 1``, the dust temperature profile is then given by
 
-              Tdust(r) = T_star*(R_star/r)**tdust_exp
+.. math::
+
+              T_\mathrm{dust}(r) = T_\mathrm{star}*(R_\mathrm{star}/r)^\textrm{tdust_exp}
 
 where T_star and R_star are the stellar (effective) temperature and radius as defined in the .setup file
 
 
+
 **Have fun :)**

src/main/checkoptions.F90

@@ -30,8 +30,8 @@ module checkoptions
 !
 !-------------------------------------------------------------------
 subroutine check_compile_time_settings(ierr)
- use part,     only:mhd,gravity,ngradh,h2chemistry,maxvxyzu,use_dust,gr
- use dim,      only:use_dustgrowth,maxtypes,mpi,inject_parts
+ use part,     only:mhd,gravity,ngradh,maxvxyzu,use_dust,gr
+ use dim,      only:use_dustgrowth,maxtypes,mpi,inject_parts,h2chemistry
  use io,       only:error,id,master,fatal,warning
  use mpiutils, only:barrier_mpi
 #ifdef GR

src/main/checksetup.F90

@@ -37,10 +37,10 @@ module checksetup
 !+
 !------------------------------------------------------------------
 subroutine check_setup(nerror,nwarn,restart)
- use dim,  only:maxp,maxvxyzu,periodic,use_dust,ndim,mhd,use_dustgrowth, &
+ use dim,  only:maxp,maxvxyzu,periodic,use_dust,ndim,mhd,use_dustgrowth,h2chemistry, &
                 do_radiation,n_nden_phantom,mhd_nonideal,do_nucleation,use_krome
  use part, only:xyzh,massoftype,hfact,vxyzu,npart,npartoftype,nptmass,gravity, &
-                iphase,maxphase,isetphase,labeltype,igas,h2chemistry,maxtypes,&
+                iphase,maxphase,isetphase,labeltype,igas,maxtypes,&
                 idust,xyzmh_ptmass,vxyz_ptmass,iboundary,isdeadh,ll,ideadhead,&
                 kill_particle,shuffle_part,iamtype,iamdust,Bxyz,rad,radprop, &
                 remove_particle_from_npartoftype,ien_type,ien_etotal,gr
@@ -104,7 +104,7 @@ subroutine check_setup(nerror,nwarn,restart)
        nerror = nerror + 1
     endif
  else
-    if (polyk < tiny(0.) .and. ieos /= 2) then
+    if (polyk < tiny(0.) .and. ieos /= 2 .and. ieos /= 5) then
        print*,'WARNING! polyk = ',polyk,' in setup, speed of sound will be zero in equation of state'
        nwarn = nwarn + 1
     endif
@@ -238,7 +238,7 @@ subroutine check_setup(nerror,nwarn,restart)
        nerror = nerror + 1
     endif
  else
-    if (abs(gamma-1.) > tiny(gamma) .and. (ieos /= 2 .and. ieos /=9)) then
+    if (abs(gamma-1.) > tiny(gamma) .and. (ieos /= 2 .and. ieos /= 5 .and. ieos /=9)) then
        print*,'*** ERROR: using isothermal EOS, but gamma = ',gamma
        gamma = 1.
        print*,'*** Resetting gamma to 1, gamma = ',gamma
@@ -620,7 +620,7 @@ subroutine check_setup_ptmass(nerror,nwarn,hmin)
  !
  !  check that radiation properties are sensible
  !
- if (isink_radiation > 1 .and. xyzmh_ptmass(ilum,1) < 1e-10) then
+ if (isink_radiation > 1 .and. xyzmh_ptmass(ilum,1) < 1e-15) then
     nerror = nerror + 1
     print*,'ERROR: isink_radiation > 1 and sink particle has no luminosity'
     return

src/main/config.F90

@@ -241,12 +241,8 @@ module dim
 ! H2 Chemistry
 !--------------------
  integer :: maxp_h2 = 0
-#ifdef H2CHEM
- logical, parameter :: h2chemistry = .true.
-#else
- logical, parameter :: h2chemistry = .false.
-#endif
  integer, parameter :: nabundances = 5
+ logical :: h2chemistry = .false.
 
 !--------------------
 ! Self-gravity
@@ -287,10 +283,11 @@ module dim
 !--------------------
 ! Dust formation
 !--------------------
- logical :: do_nucleation = .false.
- integer :: itau_alloc    = 0
- integer :: itauL_alloc   = 0
- integer :: inucleation   = 0
+ logical :: do_nucleation  = .false.
+ logical :: update_muGamma = .false.
+ integer :: itau_alloc     = 0
+ integer :: itauL_alloc    = 0
+ integer :: inucleation    = 0
  !number of elements considered in the nucleation chemical network
  integer, parameter :: nElements = 10
 #ifdef DUST_NUCLEATION
@@ -365,9 +362,9 @@ subroutine update_max_sizes(n,ntot)
 
  maxp = n
 
-#ifdef KROME
- maxp_krome = maxp
-#endif
+ if (use_krome) maxp_krome = maxp
+
+ if (h2chemistry) maxp_h2 = maxp
 
 #ifdef SINK_RADIATION
  store_dust_temperature = .true.
@@ -416,10 +413,6 @@ subroutine update_max_sizes(n,ntot)
 #endif
 #endif
 
-#ifdef H2CHEM
- maxp_h2 = maxp
-#endif
-
 #ifdef GRAVITY
  maxgrav = maxp
 #endif