Planetary wind

build/Makefile

@@ -552,8 +552,8 @@ SOURCES= physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 units.f90 \
          ${SRCINJECT_DEPS} wind_equations.f90 wind.F90  \
          ${SRCKROME} memory.f90 ${SRCREADWRITE_DUMPS} ${SRCINJECT} \
          H2regions.f90 subgroup.f90 \
-         quitdump.f90 ptmass.F90\
-         dens.F90 force.F90 deriv.F90 ${SRCAPR} readwrite_infile.F90 energies.F90 sort_particles.f90 \
+         quitdump.f90 dens.F90\
+         ptmass.F90 force.F90 deriv.F90 ${SRCAPR} readwrite_infile.F90 energies.F90 sort_particles.f90 \
          utils_shuffleparticles.F90 evwrite.f90 substepping.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

build/Makefile_setups

@@ -77,7 +77,7 @@ endif
 ifeq ($(SETUP), asteroidwind)
 #   asteroid emitting a wind (Trevascus et al. 2021)
     SETUPFILE=setup_asteroidwind.f90
-    SRCINJECT=utils_binary.f90 inject_randomwind.f90
+    SRCINJECT=utils_binary.f90 evolve_planet.f90 inject_randomwind.f90
     IND_TIMESTEPS=yes
     CONST_AV=yes
     ISOTHERMAL=yes
@@ -87,7 +87,7 @@ endif
 ifeq ($(SETUP), randomwind)
 #   asteroid emitting a wind (Trevascus et al. 2021)
     SETUPFILE=setup_disc.f90
-    SRCINJECT=utils_binary.f90 inject_randomwind.f90
+    SRCINJECT=utils_binary.f90 evolve_planet.f90 inject_randomwind.f90
     IND_TIMESTEPS=yes
     CONST_AV=yes
     ISOTHERMAL=yes
@@ -176,6 +176,17 @@ ifeq ($(SETUP), disc)
     IND_TIMESTEPS=yes
 endif
 
+ifeq ($(SETUP), boilingplanets)
+#   locally isothermal gas disc with planet evolution
+    DISC_VISCOSITY=yes
+    SETUPFILE= setup_disc.f90
+    ANALYSIS= analysis_disc.f90
+    ISOTHERMAL=yes
+    KNOWN_SETUP=yes
+    IND_TIMESTEPS=yes
+    SRCINJECT= utils_binary.f90 evolve_planet.f90 inject_randomwind.f90
+endif
+
 ifeq ($(SETUP), grtde)
 #   tidal disruption event in general relativity
     SETUPFILE= setup_grtde.f90

src/main/H2regions.f90

@@ -232,7 +232,7 @@ subroutine HII_feedback(nptmass,npart,xyzh,xyzmh_ptmass,vxyzu,isionised,dt)
           hcheck = Rmax
        endif
        do while(hcheck <= Rmax)
-          call getneigh_pos((/xi,yi,zi/),0.,hcheck,3,listneigh,nneigh,xyzh,xyzcache,maxcache,ifirstincell)
+          call getneigh_pos((/xi,yi,zi/),0.,hcheck,3,listneigh,nneigh,xyzcache,maxcache,ifirstincell)
           call set_r2func_origin(xi,yi,zi)
           call Knnfunc(nneigh,r2func_origin,xyzh,listneigh) !! Here still serial version of the quicksort. Parallel version in prep..
           if (nneigh > 0) exit

src/main/config.F90

@@ -42,7 +42,7 @@ module dim
 #else
  integer, parameter :: maxptmass = 1000
 #endif
- integer, parameter :: nsinkproperties = 22
+ integer, parameter :: nsinkproperties = 24
 
  logical :: store_sf_ptmass = .false.
 

src/main/dens.F90

@@ -27,10 +27,8 @@ module densityforce
  use timing,  only:getused,printused,print_time
 
  implicit none
- character(len=80), parameter, public :: &  ! module version
-    modid="$Id$"
 
- public :: densityiterate,get_neighbour_stats
+ public :: densityiterate,get_neighbour_stats,get_density_at_pos
 
  !--indexing for xpartveci array
  integer, parameter :: &
@@ -1692,4 +1690,70 @@ subroutine store_results(icall,cell,getdv,getdb,realviscosity,stressmax,xyzh,&
 
 end subroutine store_results
 
+subroutine get_density_at_pos(x,rho,itype)
+ use linklist, only:listneigh=>listneigh_global,getneigh_pos,ifirstincell
+ use kernel,   only:get_kernel,radkern2,cnormk
+ use boundary, only:dxbound,dybound,dzbound
+ use dim,      only:periodic,maxphase,maxp,use_apr
+ use part,     only:xyzh,iphase,iamtype,ibasetype,apr_level,massoftype,aprmassoftype
+ real, intent(in) :: x(3)
+ integer, intent(in) :: itype
+ real, intent(out) :: rho
+ integer, parameter :: maxcache = 12000
+ real, save :: xyzcache(maxcache,4)
+ integer :: n,j,iamtypej,nneigh
+ real :: dx,dy,dz,hj1,rij2,q2j,qj,pmassj,wabi,grkerni
+ logical :: same_type
+
+ call getneigh_pos(x,0.,0.,3,listneigh,nneigh,xyzcache,maxcache,ifirstincell,get_j=.true.)
+ same_type=.true.
+ rho = 0.
+ loop_over_neigh: do n=1,nneigh
+    j = listneigh(n)
+    if (n <=maxcache) then
+       ! positions from cache are already mod boundary
+       dx = x(1) - xyzcache(n,1)
+       dy = x(2) - xyzcache(n,2)
+       dz = x(3) - xyzcache(n,3)
+       hj1 = xyzcache(n,4)
+    else
+       dx = x(1) - xyzh(1,j)
+       dy = x(2) - xyzh(2,j)
+       dz = x(3) - xyzh(3,j)
+       hj1 = 1./xyzh(4,j)
+    endif
+    if (periodic) then
+       if (abs(dx) > 0.5*dxbound) dx = dx - dxbound*SIGN(1.0,dx)
+       if (abs(dy) > 0.5*dybound) dy = dy - dybound*SIGN(1.0,dy)
+       if (abs(dz) > 0.5*dzbound) dz = dz - dzbound*SIGN(1.0,dz)
+    endif
+    rij2 = dx*dx + dy*dy + dz*dz
+    q2j = rij2*hj1*hj1
+    if (q2j < radkern2) then
+       !
+       ! Density, gradh and div v are only computed using
+       ! neighbours of the same type
+       !
+       if (maxphase==maxp) then
+          iamtypej  = iamtype(iphase(j))
+          same_type = ((itype == iamtypej) .or. (ibasetype(iamtypej)==itype))
+       endif
+
+       ! adjust masses for apr
+       ! this defaults to massoftype if apr_level=1
+       if (use_apr) then
+          pmassj = aprmassoftype(iamtypej,apr_level(j))
+       else
+          pmassj = massoftype(iamtypej)
+       endif
+       if (same_type)  then
+          qj = sqrt(q2j)
+          call get_kernel(q2j,qj,wabi,grkerni)
+          rho = rho + wabi*pmassj*hj1*hj1*hj1*cnormk
+       endif
+    endif
+ enddo loop_over_neigh
+
+end subroutine get_density_at_pos
+
 end module densityforce

src/main/deriv.F90

@@ -43,7 +43,7 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
  use io,             only:iprint,fatal,error
  use linklist,       only:set_linklist
  use densityforce,   only:densityiterate
- use ptmass,         only:ipart_rhomax,ptmass_calc_enclosed_mass,ptmass_boundary_crossing
+ use ptmass,         only:ipart_rhomax,ptmass_calc_enclosed_mass,ptmass_boundary_crossing,get_pressure_on_sinks
  use externalforces, only:externalforce
  use part,           only:dustgasprop,dvdx,Bxyz,set_boundaries_to_active,&
                           nptmass,xyzmh_ptmass,sinks_have_heating,dust_temp,VrelVf,fxyz_drag
@@ -62,7 +62,7 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
  use cons2prim,      only:cons2primall,cons2prim_everything
  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_porosity
+ use options,        only:implicit_radiation,implicit_radiation_store_drad,use_porosity,need_pressure_on_sinks
  integer,      intent(in)    :: icall
  integer,      intent(inout) :: npart
  integer,      intent(in)    :: nactive
@@ -183,7 +183,10 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
     ! 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
-
+!
+! compute density and pressure at location of sink particles
+!
+ if (need_pressure_on_sinks) call get_pressure_on_sinks(nptmass,xyzmh_ptmass)
 !
 ! compute dust temperature
 !

src/main/evolve_planet.f90

@@ -0,0 +1,40 @@
+!--------------------------------------------------------------------------!
+! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. !
+! Copyright (c) 2007-2025 The Authors (see AUTHORS)                        !
+! See LICENCE file for usage and distribution conditions                   !
+! http://phantomsph.github.io/                                             !
+!--------------------------------------------------------------------------!
+module evolveplanet
+!
+!  This module evolve the embedded planet
+!   calculating the accretion and wind rate according
+!   to the pressure at the Bondi radius of the planet
+!
+! :References: Rogers et al. 2024 for the boil-off model
+!
+! :Owner: Cristiano Longarini
+!
+! :Runtime parameters: None
+!
+! :Dependencies: 
+!
+ implicit none
+ public :: evolve_planet
+
+ private
+
+contains
+
+subroutine evolve_planet(pbondi,rbondi,mdotacc,mdotwind)
+! this routine should call James' code to calculate mdotacc and mdotwind for boil-off
+ use units, only:umass,utime
+ use physcon, only:jupiterm,years
+ real, intent(in) :: pbondi,rbondi
+ real, intent(out) :: mdotacc,mdotwind
+
+ mdotacc = 0.
+ mdotwind = 1.e-3*jupiterm/umass/(years/utime)
+
+end subroutine evolve_planet
+
+end module evolveplanet

src/main/inject_randomwind.f90

@@ -66,9 +66,11 @@ module inject
 !+
 !-----------------------------------------------------------------------
 subroutine init_inject(ierr)
+ use options, only:need_pressure_on_sinks
  integer, intent(inout) :: ierr
 
  ierr = 0
+ if (wind_type==2) need_pressure_on_sinks= .true.
 
 end subroutine init_inject
 
@@ -80,7 +82,7 @@ end subroutine init_inject
 subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
                             npart,npart_old,npartoftype,dtinject)
  use io,            only:fatal
- use part,          only:nptmass,massoftype,igas,hfact,ihsoft
+ use part,          only:nptmass,massoftype,igas,hfact,ihsoft,ipbondi,irbondi
  use partinject,    only:add_or_update_particle
  use physcon,       only:twopi,gg,kboltz,mass_proton_cgs
  use random,        only:get_random_pos_on_sphere, get_gaussian_pos_on_sphere
@@ -89,6 +91,7 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
  use options,       only:iexternalforce
  use externalforces,only:mass1
  use binaryutils,   only:get_orbit_bits
+ use evolveplanet,  only:evolve_planet
  real,    intent(in)    :: time, dtlast
  real,    intent(inout) :: xyzh(:,:), vxyzu(:,:), xyzmh_ptmass(:,:), vxyz_ptmass(:,:)
  integer, intent(inout) :: npart, npart_old
@@ -98,8 +101,14 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
  real,    dimension(3)  :: xyz,vxyz,r1,r2,v2,vhat,v1
  integer :: i,ipart,npinject,seed,pt
  real    :: dmdt,rinject,h,u,speed,inject_this_step,m1,m2,r,dt
- real    :: dx(3), vecz(3), veczprime(3), rotaxis(3)
- real    :: theta_rad, phi_rad, cost, sint, cosp, sinp
+ real    :: dx(3), vecz(3), veczprime(3), rotaxis(3), cs
+ real    :: theta_rad,phi_rad,cost,sint,cosp,sinp,mdotacc,mdotwind
+
+ !
+ !-- no constraint on timestep
+ !
+ dtinject = huge(dtinject)
+ if (dtlast <= 0) return
 
  ! initialise some parameter to avoid warning...
  pt = 1
@@ -108,11 +117,19 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
 
  ! calculate the wind velocity and other quantities for different wind type
  select case (wind_type)
- case(1) ! set up random wind
+ case(1,2) ! set up random wind
     if (inject_pt > nptmass) call fatal('inject_randomwind', 'not enough point masses for inject target, check inject_pt')
+    if (wind_type==2) then
+       rinject = xyzmh_ptmass(irbondi,inject_pt)
+       cs = sqrt(xyzmh_ptmass(4,inject_pt)/rinject)
+       wind_speed = cs
+       call evolve_planet(xyzmh_ptmass(ipbondi,inject_pt),xyzmh_ptmass(irbondi,inject_pt),mdotacc,mdotwind)
+       if (rinject<=xyzmh_ptmass(ihsoft,inject_pt)) call fatal('inject_randomwind', 'Bondi radius not set for wind_type=2')
+    else
+       rinject = in_code_units(r_inject_str, ierr)
+    endif
     r2 = xyzmh_ptmass(1:3,inject_pt)
-    rinject   = in_code_units(r_inject_str, ierr)
-    v2        = vxyz_ptmass(1:3,pt)
+    v2        = vxyz_ptmass(1:3,inject_pt)
     wind_speed = wind_speed_factor*sqrt(xyzmh_ptmass(4, inject_pt)/rinject)
     u         = 0. ! setup is isothermal so utherm is not stored
     h         = hfact
@@ -147,9 +164,12 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
  !
  ! Add any dependency on radius to mass injection rate (and convert to code units)
  !
- mdot = in_code_units(mdot_str,ierr)
- dmdt = mdot*mdot_func(r,r_ref) ! r_ref is the radius for which mdot_fund = mdot
-
+ if (wind_type==2) then
+    dmdt = mdotwind
+ else
+    mdot = in_code_units(mdot_str,ierr)
+    dmdt = mdot*mdot_func(r,r_ref) ! r_ref is the radius for which mdot_fund = mdot
+ endif
  !
  !-- How many particles do we need to inject?
  !   (Seems to need at least eight gas particles to not crash) <-- This statement may or may not be true...
@@ -193,16 +213,12 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
        vxyz      = v2 + wind_speed*vhat
     case default
        xyz       = r2 + rinject*get_pos_on_sphere(seed, delta_theta)
-       vxyz      = wind_speed*vhat
+       vxyz      = v2 + wind_speed*vhat
     end select
     ipart     = npart + 1
     call add_or_update_particle(igas,xyz,vxyz,h,u,ipart,npart,npartoftype,xyzh,vxyzu)
  enddo
 
- !
- !-- no constraint on timestep
- !
- dtinject = huge(dtinject)
 
 end subroutine inject_particles
 
@@ -259,14 +275,16 @@ subroutine write_options_inject(iunit)
  use infile_utils, only:write_inopt
  integer, intent(in) :: iunit
 
- call write_inopt(wind_type, 'wind_type', 'wind setup (0=asteroidwind, 1=randomwind)', iunit)
- call write_inopt(mdot_str,'mdot','mass injection rate with unit, e.g. 1e8*g/s, 1e-7M_s/yr',iunit)
- call write_inopt(npartperorbit,'npartperorbit',&
+ call write_inopt(wind_type, 'wind_type', 'wind setup (0=asteroidwind, 1=randomwind, 2=boil-off)', iunit)
+ if (wind_type /= 2) then
+    call write_inopt(mdot_str,'mdot','mass injection rate with unit, e.g. 1e8*g/s, 1e-7M_s/yr',iunit)
+    call write_inopt(npartperorbit,'npartperorbit',&
                   'particle injection rate in particles/binary orbit',iunit)
- call write_inopt(vlag,'vlag','percentage lag in velocity of wind',iunit)
- call write_inopt(mdot_type,'mdot_type','injection rate (0=const, 2=r^(-2))',iunit)
- if (mdot_type==2) then
-    call write_inopt(r_ref,'r_ref','radius at whieh Mdot=mdot for 1/r^2 injection type',iunit)
+    call write_inopt(vlag,'vlag','percentage lag in velocity of wind',iunit)
+    call write_inopt(mdot_type,'mdot_type','injection rate (0=const, 2=r^(-2))',iunit)
+    if (mdot_type==2) then
+       call write_inopt(r_ref,'r_ref','radius at whieh Mdot=mdot for 1/r^2 injection type',iunit)
+    endif
  endif
  call write_inopt(random_type, 'random_type', 'random position on the surface, 0 for random, 1 for gaussian', iunit)
  if (random_type==1) then

src/main/linklist_kdtree.F90

@@ -272,19 +272,22 @@ subroutine get_neighbour_list(inode,mylistneigh,nneigh,xyzh,xyzcache,ixyzcachesi
 
 end subroutine get_neighbour_list
 
-subroutine getneigh_pos(xpos,xsizei,rcuti,ndim,mylistneigh,nneigh,xyzh,xyzcache,ixyzcachesize,ifirstincell)
+subroutine getneigh_pos(xpos,xsizei,rcuti,ndim,mylistneigh,nneigh,xyzcache,ixyzcachesize,ifirstincell,get_j)
  use kdtree, only:getneigh
  integer, intent(in)  :: ndim,ixyzcachesize
  real,    intent(in)  :: xpos(ndim)
  real,    intent(in)  :: xsizei,rcuti
  integer, intent(out) :: mylistneigh(:)
  integer, intent(out) :: nneigh
- real,    intent(in)  :: xyzh(:,:)
  real,    intent(out) :: xyzcache(:,:)
  integer, intent(in)  :: ifirstincell(:) !ncellsmax+1)
+ logical, intent(in), optional :: get_j
+ logical :: getj
 
+ getj = .false.
+ if (present(get_j)) getj=get_j
  call getneigh(node,xpos,xsizei,rcuti,ndim,mylistneigh,nneigh,xyzcache,ixyzcachesize, &
-               ifirstincell,.false.,.false.)
+               ifirstincell,get_j,.false.)
 
 end subroutine getneigh_pos
 

src/main/options.f90

@@ -55,6 +55,9 @@ module options
  real(kind=sp), public :: mcfost_keep_part
  character(len=80), public :: Voronoi_limits_file
 
+ ! pressure on sinks
+ logical, public :: need_pressure_on_sinks
+
  ! radiation
  logical, public :: exchange_radiation_energy, limit_radiation_flux, implicit_radiation
  logical, public :: implicit_radiation_store_drad
@@ -169,6 +172,8 @@ subroutine set_default_options
  ! variable composition
  use_var_comp = .false.
 
+ need_pressure_on_sinks = .false.
+
 end subroutine set_default_options
 
 end module options