Merge pull request #582 from danieljprice/mega-merge

merge #578

AUTHORS

@@ -26,8 +26,8 @@ Terrence Tricco <[email protected]>
 Stephane Michoulier <[email protected]>
 Simone Ceppi <[email protected]>
 Spencer Magnall <[email protected]>
-Enrico Ragusa <[email protected]>
 Caitlyn Hardiman <[email protected]>
+Enrico Ragusa <[email protected]>
 Cristiano Longarini <[email protected]>
 Sergei Biriukov <[email protected]>
 Giovanni Dipierro <[email protected]>
@@ -38,38 +38,39 @@ Alison Young <[email protected]>
 Stephen Nielson <[email protected]>
 Martina Toscani <[email protected]>
 Benedetta Veronesi <[email protected]>
-Simon Glover <[email protected]>
 Sahl Rowther <[email protected]>
+Simon Glover <[email protected]>
 Thomas Reichardt <[email protected]>
 Jean-François Gonzalez <[email protected]>
-Christopher Russell <[email protected]>
 Madeline Overton <[email protected]>
-Alex Pettitt <[email protected]>
-Phantom benchmark bot <[email protected]>
+Christopher Russell <[email protected]>
 Jolien Malfait <[email protected]>
+Phantom benchmark bot <[email protected]>
 Alessia Franchini <[email protected]>
-Nicole Rodrigues <[email protected]>
+Alex Pettitt <[email protected]>
 Kieran Hirsh <[email protected]>
+Nicole Rodrigues <[email protected]>
 Farzana Meru <[email protected]>
-David Trevascus <[email protected]>
 Nicolás Cuello <[email protected]>
-Chris Nixon <[email protected]>
+David Trevascus <[email protected]>
 Miguel Gonzalez-Bolivar <[email protected]>
-Maxime Lombart <[email protected]>
-Zachary Pellow <[email protected]>
-Orsola De Marco <[email protected]>
+Chris Nixon <[email protected]>
+Shunquan Huang <[email protected]>
 Joe Fisher <[email protected]>
 Benoit Commercon <[email protected]>
+Zachary Pellow <[email protected]>
 Giulia Ballabio <[email protected]>
+Maxime Lombart <[email protected]>
+Orsola De Marco <[email protected]>
+Steven Rieder <[email protected]>
+Stéven Toupin <[email protected]>
+Taj Jankovič <[email protected]>
+Jeremy Smallwood <[email protected]>
+Ariel Chitan <[email protected]>
 Rebecca Martin <[email protected]>
 Jorge Cuadra <[email protected]>
-Hugh Griffiths <[email protected]>
-Jeremy Smallwood <[email protected]>
 David Bamba <[email protected]>
+Shunquan Huang <[email protected]>
 Cox, Samuel <[email protected]>
 Chunliang Mu <[email protected]>
-Shunquan Huang <[email protected]>
-Steven Rieder <[email protected]>
-Stéven Toupin <[email protected]>
-Taj Jankovič <[email protected]>
-Ariel Chitan <[email protected]>
+Hugh Griffiths <[email protected]>

src/main/geometry.f90

@@ -203,7 +203,7 @@ subroutine set_rotation_angles(a,b,sin_a,sin_b,cos_a,cos_b)
  endif
  if (present(cos_b)) then
     cosb = cos_b
-    if (.not.present(cos_b)) sinb = sqrt(1. - cosb**2)
+    if (.not.present(sin_b)) sinb = sqrt(1. - cosb**2)
  endif
 
 end subroutine set_rotation_angles

src/main/inject_randomwind.f90

@@ -12,23 +12,30 @@ module inject
 ! :References:
 !   Trevascus et al. (2021), MNRAS 505, L21-L25
 !
-! :Owner: David Liptai
+! :Owner: Shunquan Huang
 !
 ! :Runtime parameters:
-!   - delta_theta : *standard deviation for the gaussion distribution (random_type=1)*
+!   - delta_theta : *standard deviation for the gaussion distribution*
+!   - inject_pt   : *the particle that excites wind (when wind_type=1)*
 !   - mdot        : *mass injection rate with unit, e.g. 1e8*g/s, 1e-7M_s/yr*
-!   - mdot_type   : *injection rate (0=const, 1=cos(t), 2=r^(-2))*
+!   - mdot_type   : *injection rate (0=const, 2=r^(-2))*
+!   - phi         : *the tilt orientation of the star, (random_type=1)*
+!   - r_inject    : *inject radius with units, e.g. 1*AU, 1e8m, (when wind_type=1)*
+!   - r_ref       : *radius at whieh Mdot=mdot for 1/r^2 injection type*
 !   - random_type : *random position on the surface, 0 for random, 1 for gaussian*
+!   - theta       : *the tilt inclination of the star or planet (random_type=1)*
 !   - vlag        : *percentage lag in velocity of wind*
+!   - wind_type   : *wind setup (0=asteroidwind, 1=randomwind)*
 !
 ! :Dependencies: binaryutils, externalforces, infile_utils, io, options,
-!   part, partinject, physcon, random, units
+!   part, partinject, physcon, random, units, vectorutils
 !
  use io, only:error
  use physcon, only:pi
  implicit none
  character(len=*), parameter, public :: inject_type = 'randomwind'
  character(len=20), public :: mdot_str = "5.e8*g/s"
+ character(len=20), public :: r_inject_str = "0.5*au"
  real, public :: mdot = 1.e8     ! mass injection rate in grams/second
 
  public :: init_inject,inject_particles,write_options_inject,read_options_inject,&
@@ -37,13 +44,20 @@ module inject
  private
 
  real         :: npartperorbit = 1000.     ! particle injection rate in particles per orbit
+ integer      :: wind_type     = 0        ! wind setup (0=asteroidwind, 1=randomwind)
  real         :: vlag          = 0.0      ! percentage lag in velocity of wind
  integer      :: mdot_type     = 0        ! injection rate (0=const, 1=cos(t), 2=r^(-2))
  integer      :: random_type   = 0        ! random position on the surface, 0 for random, 1 for gaussian
  real         :: delta_theta   = 0.5      ! standard deviation for the gaussion distribution (random_type=1)
  real         :: have_injected = 0.
- real         :: t_old = 0.
- real         :: r_ref = 1.
+ real         :: t_old         = 0.
+ real         :: r_ref         = 1.       ! reference radius for mdot_type=2
+ real         :: theta         = 0.       ! the inclination of the star or planet
+ real         :: phi           = 0.       ! the orientation of the star
+ integer      :: inject_pt     = 2        ! the partical that produce wind (when wind_type=1)
+ real         :: wind_speed    = 1.0      ! wind speed in code unit (when wind_type=1)
+ real         :: wind_speed_factor = 1.2  ! factor to scale the wind speed based on the Keplerian speed at rinject
+ !real         :: rinject       = 1.0
 
 contains
 !-----------------------------------------------------------------------
@@ -71,6 +85,7 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
  use physcon,       only:twopi,gg,kboltz,mass_proton_cgs
  use random,        only:get_random_pos_on_sphere, get_gaussian_pos_on_sphere
  use units,         only:in_code_units
+ use vectorutils,   only:cross_product3D, rotatevec
  use options,       only:iexternalforce
  use externalforces,only:mass1
  use binaryutils,   only:get_orbit_bits
@@ -82,36 +97,52 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
  integer :: ierr
  real,    dimension(3)  :: xyz,vxyz,r1,r2,v2,vhat,v1
  integer :: i,ipart,npinject,seed,pt
- real    :: dmdt,rbody,h,u,speed,inject_this_step
- real    :: m1,m2,r
- real    :: dt
-
- if (nptmass < 1 .and. iexternalforce == 0) &
-    call fatal('inject_randomwind','not enough point masses for random wind injection')
- if (nptmass > 2) &
-    call fatal('inject_randomwind','too many point masses for random wind injection')
-
- if (nptmass == 2) then
-    pt = 2
-    r1 = xyzmh_ptmass(1:3,1)
-    m1 = xyzmh_ptmass(4,1)
-    v1 = vxyz_ptmass(1:3,1)
- else
-    pt = 1
-    r1 = 0.
-    m1 = mass1
-    v1 = 0.
- endif
-
- r2        = xyzmh_ptmass(1:3,pt)
- rbody     = xyzmh_ptmass(ihsoft,pt)
- m2        = xyzmh_ptmass(4,pt)
- v2        = vxyz_ptmass(1:3,pt)
-
- speed     = sqrt(dot_product(v2,v2))
- vhat      = v2/speed
-
- r         = sqrt(dot_product(r1-r2,r1-r2))
+ 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
+
+ ! initialise some parameter to avoid warning...
+ pt = 1
+ rinject = 1.0
+ r = 1.0
+
+ ! calculate the wind velocity and other quantities for different wind type
+ select case (wind_type)
+ case(1) ! set up random wind
+    if (inject_pt > nptmass) call fatal('inject_randomwind', 'not enough point masses for inject target, check inject_pt')
+    r2 = xyzmh_ptmass(1:3,inject_pt)
+    rinject   = in_code_units(r_inject_str, ierr)
+    v2        = vxyz_ptmass(1:3,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
+ case default ! set up asteroid wind
+    if (nptmass < 1 .and. iexternalforce == 0) &
+      call fatal('inject_asteroidwind','not enough point masses for asteroid wind injection')
+    if (nptmass > 2) &
+      call fatal('inject_asteroidwind','too many point masses for asteroid wind injection')
+    if (nptmass == 2) then
+       pt = 2
+       r1 = xyzmh_ptmass(1:3,1)
+       m1 = xyzmh_ptmass(4,1)
+       v1 = vxyz_ptmass(1:3,1)
+    else
+       pt = 1
+       r1 = 0.
+       m1 = mass1
+       v1 = 0.
+    endif
+    r2        = xyzmh_ptmass(1:3,pt)
+    rinject   = xyzmh_ptmass(ihsoft,pt)
+    m2        = xyzmh_ptmass(4,pt)
+    v2        = vxyz_ptmass(1:3,pt)
+    speed     = sqrt(dot_product(v2,v2))
+    vhat      = v2/speed
+    r         = sqrt(dot_product(r1-r2,r1-r2))
+    wind_speed  = (1.-vlag/100)*speed
+    u         = 0. ! setup is isothermal so utherm is not stored
+    h         = hfact*(rinject/2.)
+ end select
 
  !
  ! Add any dependency on radius to mass injection rate (and convert to code units)
@@ -129,22 +160,41 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
     ! Calculate how many extra particles from previous step to now
     dt = time - t_old
     inject_this_step = dt*dmdt/massoftype(igas)
-
     npinject = max(0, int(0.5 + have_injected + inject_this_step - npartoftype(igas) ))
-
     ! Save for next step (faster than integrating the whole thing each time)
     t_old = time
     have_injected = have_injected + inject_this_step
  endif
-
+ !
+ !-- set up the tilt of the star, and vectors for rotation
+ !
+ theta_rad = theta*pi/180.
+ phi_rad = phi*pi/180.
+ cost = cos(theta_rad)
+ sint = sin(theta_rad)
+ cosp = cos(phi_rad)
+ sinp = sin(phi_rad)
+ vecz = (/0.,0.,1./)
+ veczprime = (/sint*cosp,sint*sinp,cost/)
+ if (abs(theta_rad-0)<1e-6) then
+    rotaxis = (/0.,0.,1./)
+ else
+    call cross_product3D(vecz, veczprime, rotaxis)
+ endif
  !
  !-- Randomly inject particles around the body's outer 'radius'.
  !
  do i=1,npinject
-    xyz       = r2 + rbody*get_pos_on_sphere(seed, delta_theta)
-    vxyz      = (1.-vlag/100)*speed*vhat
-    u         = 0. ! setup is isothermal so utherm is not stored
-    h         = hfact*(rbody/2.)
+    select case (wind_type)
+    case (1)
+       dx = get_pos_on_sphere(seed, delta_theta)
+       call rotatevec(dx, rotaxis, theta_rad)
+       call cross_product3D(veczprime, dx, vhat)
+       vxyz      = v2 + wind_speed*vhat
+    case default
+       xyz       = r2 + rinject*get_pos_on_sphere(seed, delta_theta)
+       vxyz      = wind_speed*vhat
+    end select
     ipart     = npart + 1
     call add_or_update_particle(igas,xyz,vxyz,h,u,ipart,npart,npartoftype,xyzh,vxyzu)
  enddo
@@ -209,6 +259,7 @@ 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',&
                   'particle injection rate in particles/binary orbit',iunit)
@@ -220,7 +271,15 @@ subroutine write_options_inject(iunit)
  call write_inopt(random_type, 'random_type', 'random position on the surface, 0 for random, 1 for gaussian', iunit)
  if (random_type==1) then
     call write_inopt(delta_theta, 'delta_theta', 'standard deviation for the gaussion distribution', iunit)
+    call write_inopt(theta, 'theta', 'the tilt inclination of the star or planet (random_type=1)', iunit)
+    call write_inopt(phi, 'phi', 'the tilt orientation of the star, (random_type=1)', iunit)
+ endif
+ if (wind_type==1) then
+    call write_inopt(inject_pt, 'inject_pt', 'the particle that excites wind (when wind_type=1)', iunit)
+    call write_inopt(r_inject_str, 'r_inject', 'inject radius with units, e.g. 1*AU, 1e8m, (when wind_type=1)', iunit)
  endif
+ call write_inopt(wind_speed_factor, &
+ & 'wind_speed_factor', 'factor to scale the wind speed based on the Keplerian speed at rinject', iunit)
 
 end subroutine write_options_inject
 
@@ -243,6 +302,9 @@ subroutine read_options_inject(name,valstring,imatch,igotall,ierr)
     read(valstring,'(A)',iostat=ierr) mdot_str
     ngot = ngot + 1
     ! if (mdot  <  0.) call fatal(label,'mdot < 0 in input options')
+ case('wind_type')
+    read(valstring,*,iostat=ierr) wind_type
+    ngot = ngot + 1
  case('npartperorbit')
     read(valstring,*,iostat=ierr) npartperorbit
     ngot = ngot + 1
@@ -262,6 +324,21 @@ subroutine read_options_inject(name,valstring,imatch,igotall,ierr)
  case('delta_theta')
     read(valstring,*,iostat=ierr) delta_theta
     ngot = ngot + 1
+ case('theta')
+    read(valstring,*,iostat=ierr) theta
+    ngot = ngot + 1
+ case('phi')
+    read(valstring,*,iostat=ierr) phi
+    ngot = ngot + 1
+ case('inject_pt')
+    read(valstring,*,iostat=ierr) inject_pt
+    ngot = ngot + 1
+ case('wind_speed_factor')
+    read(valstring,*,iostat=ierr) wind_speed_factor
+    ngot = ngot + 1
+ case('r_inject')
+    read(valstring,'(a)',iostat=ierr) r_inject_str
+    ngot = ngot + 1
  case default
     imatch = .false.
  end select

src/main/metric_tools.f90

@@ -19,7 +19,7 @@ module metric_tools
 !
 ! :Runtime parameters: None
 !
-! :Dependencies: fastmath, inverse4x4, metric
+! :Dependencies: inverse4x4, metric
 !
  use metric, only:imetric
  implicit none

src/main/ptmass.F90

@@ -28,7 +28,7 @@ module ptmass
 !   - h_soft_sinkgas  : *softening length for new sink particles*
 !   - h_soft_sinksink : *softening length between sink particles*
 !   - icreate_sinks   : *allow automatic sink particle creation*
-!   - isink_potential : *sink potential(0=1/r,1=surf)*
+!   - isink_potential : *sink potential (0=1/r,1=surf)*
 !   - r_crit          : *critical radius for point mass creation (no new sinks < r_crit from existing sink)*
 !   - r_merge_cond    : *sinks will merge if bound within this radius*
 !   - r_merge_uncond  : *sinks will unconditionally merge within this separation*

src/main/utils_gr.f90

@@ -14,7 +14,7 @@ module utils_gr
 !
 ! :Runtime parameters: None
 !
-! :Dependencies: fastmath, io, metric, metric_tools, part
+! :Dependencies: io, metric, metric_tools, part
 !
  implicit none