From 4d61c3a0d681bca6560d5249d53f5296365637f1 Mon Sep 17 00:00:00 2001 From: Shunquan Huang Date: Wed, 31 Jul 2024 18:42:14 -0700 Subject: [PATCH 1/6] new features for randomwind --- src/main/geometry.f90 | 2 +- src/main/inject_randomwind.f90 | 138 ++++++++++++++++++++++++------- src/setup/setup_asteroidwind.f90 | 18 ++-- 3 files changed, 117 insertions(+), 41 deletions(-) diff --git a/src/main/geometry.f90 b/src/main/geometry.f90 index e0fcf88d2..d90782a91 100644 --- a/src/main/geometry.f90 +++ b/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 diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 index e0f78cf44..69f6c235c 100644 --- a/src/main/inject_randomwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -16,11 +16,20 @@ module inject ! ! :Runtime parameters: ! - mdot_str : *mdot with unit* +! - wind_type : wind setup (0=asteroidwind, 1=randomwind) ! - mdot : *mass injection rate in grams/second* ! - mdot_type : *injection rate (0=const, 1=cos(t), 2=r^(-2))* -! - vlag : *percentage lag in velocity of wind* +! - vlag : *percentage lag in velocity of wind* ! - random_type : random position on the surface, 0 for random, 1 for gaussian ! - delta_theta : standard deviation for the gaussion distribution (random_type=1) +! - theta : the inclination of the star or planet (random_type=1, +! if theta = 90, more particles are injected in z direction) +! - phi : the orientation of the star, (random_type=1, +! if theta=90 and phi=90 more particles are injected in x-z plane) +! - inject_pt : the partical that produce wind (when wind_type=1) +! - wind_speed_factor : factor to scale the wind speed based on the Keplerian speed at rinject (when wind_type=1) +! - wind_speed : wind speed +! - r_inject_str : rinject with unit (when wind_type=1) ! ! :Dependencies: binaryutils, externalforces, infile_utils, io, options, ! part, partinject, physcon, random, units @@ -30,6 +39,7 @@ module inject implicit none character(len=*), parameter, public :: inject_type = 'asteroidwind' 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,& @@ -38,10 +48,17 @@ 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 :: 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 !----------------------------------------------------------------------- @@ -68,7 +85,8 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& 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 - use units, only:umass, utime, in_code_units + 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 @@ -80,45 +98,62 @@ 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,rasteroid,h,u,speed,inject_this_step - real :: m1,m2,r - real :: dt + 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, save :: have_injected,t_old real, save :: semia - - if (nptmass < 2 .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 + + ! initial 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(0) ! set up asteroid wind + if (nptmass < 2 .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 + else pt = 1 r1 = 0. m1 = mass1 v1 = 0. - endif - - r2 = xyzmh_ptmass(1:3,pt) - rasteroid = xyzmh_ptmass(ihsoft,pt) - m2 = xyzmh_ptmass(4,pt) - v2 = vxyz_ptmass(1:3,pt) + 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.) + + 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 - speed = sqrt(dot_product(v2,v2)) - vhat = v2/speed - - r = sqrt(dot_product(r1-r2,r1-r2)) + end select ! ! 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,semia) ! Use semi-major axis as r_ref - ! !-- 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... @@ -129,23 +164,38 @@ 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/) + call cross_product3D(vecz, veczprime, rotaxis) ! !-- Randomly inject particles around the asteroids outer 'radius'. ! Only inject them on the side that is facing the central sink ! do i=1,npinject - xyz = r2 + rasteroid*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*(rasteroid/2.) + select case (wind_type) + case (0) + xyz = r2 + rinject*get_pos_on_sphere(seed, delta_theta) + vxyz = wind_speed*vhat + 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 + 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) @@ -216,6 +267,13 @@ subroutine write_options_inject(iunit) call write_inopt(mdot_type,'mdot_type','injection rate (0=const, 1=cos(t), 2=r^(-2))',iunit) call write_inopt(random_type, 'random_type', 'random position on the surface, 0 for random, 1 for gaussian', iunit) call write_inopt(delta_theta, 'delta_theta', 'standard deviation for the gaussion distribution (random_type=1)', 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) + call write_inopt(inject_pt, 'inject_pt', 'the partical that produce wind (when wind_type=1)', iunit) + + call write_inopt(wind_speed_factor, & + & 'wind_speed_factor', 'factor to scale the wind speed based on the Keplerian speed at rinject', iunit) + call write_inopt(r_inject_str, 'r_inject', 'inject radius with units, e.g. 1*AU, 1e8m, (when wind_type=1)', iunit) end subroutine write_options_inject @@ -238,6 +296,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 @@ -254,6 +315,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 diff --git a/src/setup/setup_asteroidwind.f90 b/src/setup/setup_asteroidwind.f90 index 0a2af47e0..47b021259 100644 --- a/src/setup/setup_asteroidwind.f90 +++ b/src/setup/setup_asteroidwind.f90 @@ -26,7 +26,7 @@ module setup ! - mdot_str : *mdot with unit* ! - norbits : *number of orbits* ! - npart_at_end : *number of particles injected after norbits* -! - rasteroid : *radius of asteroid (km)* +! - rinject : *radius of asteroid (km)* ! - semia : *semi-major axis (solar radii)* ! ! :Dependencies: eos, extern_lensethirring, externalforces, infile_utils, @@ -38,7 +38,7 @@ module setup implicit none public :: setpart - real :: m1,m2,ecc,semia,hacc1,rasteroid,norbits,gastemp + real :: m1,m2,ecc,semia,hacc1,rinject,norbits,gastemp integer :: npart_at_end,dumpsperorbit,ipot private @@ -49,7 +49,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, use part, only:nptmass,xyzmh_ptmass,vxyz_ptmass,ihacc,ihsoft,idust,set_particle_type,igas use setbinary, only:set_binary,get_a_from_period use spherical, only:set_sphere - use units, only:set_units,umass,udist,utime,unit_velocity,in_code_units + use units, only:set_units,umass,udist,unit_velocity,in_code_units use physcon, only:solarm,au,pi,solarr,ceresm,km,kboltz,mass_proton_cgs use externalforces, only:iext_binary, iext_einsteinprec, update_externalforce, & mass1,accradius1 @@ -83,7 +83,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ecc = 0.54 ! (eccentricity) semia = 0.73 ! (solar radii) hacc1 = 0.1679 ! (solar radii) - rasteroid = 2338.3 ! (km) + rinject = 2338.3 ! (km) gastemp = 5000. ! (K) norbits = 1000. mdot = 5.e8 ! Mass injection rate (will change later by the mdot_str) @@ -125,7 +125,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, semia = semia*solarr/udist hacc1 = hacc1*solarr/udist - rasteroid = rasteroid*km/udist + rinject = rinject*km/udist ! !--general parameters @@ -159,7 +159,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, !--Set a binary orbit given the desired orbital parameters ! call set_binary(m1,m2,semia,ecc,hacc1,hacc2,xyzmh_ptmass,vxyz_ptmass,nptmass,ierr) - xyzmh_ptmass(ihsoft,2) = rasteroid ! Asteroid radius softening + xyzmh_ptmass(ihsoft,2) = rinject ! Asteroid radius softening else @@ -179,7 +179,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, xyzmh_ptmass(4,1) = m2 xyzmh_ptmass(ihacc,1) = hacc2 ! asteroid should not accrete - xyzmh_ptmass(ihsoft,1) = rasteroid ! asteroid radius softening + xyzmh_ptmass(ihsoft,1) = rinject ! asteroid radius softening endif ! we use the estimated injection rate and the final time to set the particle mass @@ -218,7 +218,7 @@ subroutine write_setupfile(filename) call write_inopt(ecc, 'ecc', 'eccentricity', iunit) call write_inopt(semia, 'semia', 'semi-major axis (solar radii)', iunit) call write_inopt(hacc1, 'hacc1', 'white dwarf (sink) accretion radius (solar radii)',iunit) - call write_inopt(rasteroid, 'rasteroid', 'radius of asteroid (km)', iunit) + call write_inopt(rinject, 'rinject', 'radius of asteroid (km)', iunit) call write_inopt(gastemp, 'gastemp', 'gas temperature in K', iunit) call write_inopt(norbits, 'norbits', 'number of orbits', iunit) call write_inopt(dumpsperorbit,'dumpsperorbit','number of dumps per orbit', iunit) @@ -248,7 +248,7 @@ subroutine read_setupfile(filename,ierr) call read_inopt(ecc, 'ecc', db,min=0.,errcount=nerr) call read_inopt(semia, 'semia', db,min=0.,errcount=nerr) call read_inopt(hacc1, 'hacc1', db,min=0.,errcount=nerr) - call read_inopt(rasteroid, 'rasteroid', db,min=0.,errcount=nerr) + call read_inopt(rinject, 'rinject', db,min=0.,errcount=nerr) call read_inopt(gastemp, 'gastemp', db,min=0.,errcount=nerr) call read_inopt(norbits, 'norbits', db,min=0.,errcount=nerr) call read_inopt(dumpsperorbit,'dumpsperorbit',db,min=0 ,errcount=nerr) From ea350d88e466ddf10a9326f167041e3e000ea2ad Mon Sep 17 00:00:00 2001 From: Shunquan Huang Date: Wed, 31 Jul 2024 21:59:09 -0700 Subject: [PATCH 2/6] fix a bug when theta=0 and some typo --- src/main/inject_randomwind.f90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 index 69f6c235c..00cb7abfb 100644 --- a/src/main/inject_randomwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -180,7 +180,11 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& sinp = sin(phi_rad) vecz = (/0.,0.,1./) veczprime = (/sint*cosp,sint*sinp,cost/) - call cross_product3D(vecz, veczprime, rotaxis) + 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 asteroids outer 'radius'. ! Only inject them on the side that is facing the central sink @@ -269,8 +273,7 @@ subroutine write_options_inject(iunit) call write_inopt(delta_theta, 'delta_theta', 'standard deviation for the gaussion distribution (random_type=1)', 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) - call write_inopt(inject_pt, 'inject_pt', 'the partical that produce wind (when wind_type=1)', iunit) - + call write_inopt(inject_pt, 'inject_pt', 'the particle that excites wind (when wind_type=1)', iunit) call write_inopt(wind_speed_factor, & & 'wind_speed_factor', 'factor to scale the wind speed based on the Keplerian speed at rinject', iunit) call write_inopt(r_inject_str, 'r_inject', 'inject radius with units, e.g. 1*AU, 1e8m, (when wind_type=1)', iunit) From aa28b8b7fd71216f4ed2eed00716caf47f2d1145 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Tue, 6 Aug 2024 10:58:22 +1000 Subject: [PATCH 3/6] [header-bot] updated file headers --- src/main/inject_randomwind.f90 | 30 +++++++++++++----------------- src/main/metric_tools.f90 | 2 +- src/main/ptmass.F90 | 2 +- src/main/utils_gr.f90 | 2 +- src/setup/setup_asteroidwind.f90 | 2 +- 5 files changed, 17 insertions(+), 21 deletions(-) diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 index 4570d778c..0830b37cd 100644 --- a/src/main/inject_randomwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -12,27 +12,23 @@ module inject ! :References: ! Trevascus et al. (2021), MNRAS 505, L21-L25 ! -! :Owner: David Liptai +! :Owner: Shunquan Huang ! ! :Runtime parameters: -! - mdot_str : *mdot with unit* -! - wind_type : wind setup (0=asteroidwind, 1=randomwind) -! - mdot : *mass injection rate in grams/second* -! - mdot_type : *injection rate (0=const, 1=cos(t), 2=r^(-2))* -! - vlag : *percentage lag in velocity of wind* -! - random_type : random position on the surface, 0 for random, 1 for gaussian -! - delta_theta : standard deviation for the gaussion distribution (random_type=1) -! - theta : the inclination of the star or planet (random_type=1, -! if theta = 90, more particles are injected in z direction) -! - phi : the orientation of the star, (random_type=1, -! if theta=90 and phi=90 more particles are injected in x-z plane) -! - inject_pt : the partical that produce wind (when wind_type=1) -! - wind_speed_factor : factor to scale the wind speed based on the Keplerian speed at rinject (when wind_type=1) -! - wind_speed : wind speed -! - r_inject_str : rinject with unit (when wind_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, 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 diff --git a/src/main/metric_tools.f90 b/src/main/metric_tools.f90 index b0cf56c3f..05b485fcb 100644 --- a/src/main/metric_tools.f90 +++ b/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 diff --git a/src/main/ptmass.F90 b/src/main/ptmass.F90 index 2f79f2738..505e32623 100644 --- a/src/main/ptmass.F90 +++ b/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* diff --git a/src/main/utils_gr.f90 b/src/main/utils_gr.f90 index 1308aeef8..6ecc4be43 100644 --- a/src/main/utils_gr.f90 +++ b/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 diff --git a/src/setup/setup_asteroidwind.f90 b/src/setup/setup_asteroidwind.f90 index 717f596ae..6a34c93b4 100644 --- a/src/setup/setup_asteroidwind.f90 +++ b/src/setup/setup_asteroidwind.f90 @@ -25,7 +25,7 @@ module setup ! - mdot : *mass injection rate with unit, e.g. 1e8*g/s, 1e-7M_s/yr (from setup)* ! - norbits : *number of orbits* ! - npart_at_end : *number of particles injected after norbits* -! - rinject : *radius of asteroid (km)* +! - rinject : *radius of asteroid (km)* ! - semia : *semi-major axis (solar radii)* ! ! :Dependencies: eos, extern_lensethirring, externalforces, infile_utils, From d28ed996414cd47f64e9755ffca60cae263890f6 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Tue, 6 Aug 2024 10:58:50 +1000 Subject: [PATCH 4/6] [space-bot] whitespace at end of lines removed --- src/main/inject_randomwind.f90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 index 0830b37cd..54c6f77b9 100644 --- a/src/main/inject_randomwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -51,13 +51,13 @@ module inject 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. ! reference radius for mdot_type=2 + 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 + !real :: rinject = 1.0 contains !----------------------------------------------------------------------- @@ -85,7 +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 vectorutils, only:cross_product3D, rotatevec use options, only:iexternalforce use externalforces,only:mass1 use binaryutils, only:get_orbit_bits @@ -100,12 +100,12 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& 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 @@ -174,13 +174,13 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& sint = sin(theta_rad) cosp = cos(phi_rad) sinp = sin(phi_rad) - vecz = (/0.,0.,1./) + vecz = (/0.,0.,1./) veczprime = (/sint*cosp,sint*sinp,cost/) if (abs(theta_rad-0)<1e-6) then - rotaxis = (/0.,0.,1./) + rotaxis = (/0.,0.,1./) else call cross_product3D(vecz, veczprime, rotaxis) - endif + endif ! !-- Randomly inject particles around the body's outer 'radius'. ! @@ -278,7 +278,7 @@ subroutine write_options_inject(iunit) 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, & + 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 From ef312b07edb206c0192114d23d9c782f2e7e705d Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Tue, 6 Aug 2024 10:58:50 +1000 Subject: [PATCH 5/6] [author-bot] updated AUTHORS file --- AUTHORS | 37 +++++++++++++++++++------------------ 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/AUTHORS b/AUTHORS index b3a77dfde..9f3013a96 100644 --- a/AUTHORS +++ b/AUTHORS @@ -26,8 +26,8 @@ Terrence Tricco Stephane Michoulier Simone Ceppi Spencer Magnall -Enrico Ragusa Caitlyn Hardiman +Enrico Ragusa Cristiano Longarini Sergei Biriukov Giovanni Dipierro @@ -38,38 +38,39 @@ Alison Young Stephen Nielson Martina Toscani Benedetta Veronesi -Simon Glover Sahl Rowther +Simon Glover Thomas Reichardt Jean-François Gonzalez -Christopher Russell Madeline Overton -Alex Pettitt -Phantom benchmark bot +Christopher Russell Jolien Malfait +Phantom benchmark bot Alessia Franchini -Nicole Rodrigues +Alex Pettitt Kieran Hirsh +Nicole Rodrigues Farzana Meru -David Trevascus Nicolás Cuello -Chris Nixon +David Trevascus Miguel Gonzalez-Bolivar -Maxime Lombart -Zachary Pellow -Orsola De Marco +Chris Nixon +Shunquan Huang Joe Fisher Benoit Commercon +Zachary Pellow Giulia Ballabio +Maxime Lombart +Orsola De Marco +Steven Rieder +Stéven Toupin +Taj Jankovič +Jeremy Smallwood +Ariel Chitan Rebecca Martin Jorge Cuadra -Hugh Griffiths -Jeremy Smallwood David Bamba +Shunquan Huang Cox, Samuel Chunliang Mu -Shunquan Huang -Steven Rieder -Stéven Toupin -Taj Jankovič -Ariel Chitan +Hugh Griffiths From 1bb03aeed5a614c791c16486d967d0541dba14aa Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Tue, 6 Aug 2024 10:59:32 +1000 Subject: [PATCH 6/6] [indent-bot] standardised indentation --- src/main/inject_randomwind.f90 | 76 +++++++++++++++++----------------- 1 file changed, 38 insertions(+), 38 deletions(-) diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 index 54c6f77b9..7c830098c 100644 --- a/src/main/inject_randomwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -109,39 +109,39 @@ 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 - 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 + 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) & + if (nptmass < 1 .and. iexternalforce == 0) & call fatal('inject_asteroidwind','not enough point masses for asteroid wind injection') - if (nptmass > 2) & + 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.) + 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 ! @@ -177,7 +177,7 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& vecz = (/0.,0.,1./) veczprime = (/sint*cosp,sint*sinp,cost/) if (abs(theta_rad-0)<1e-6) then - rotaxis = (/0.,0.,1./) + rotaxis = (/0.,0.,1./) else call cross_product3D(vecz, veczprime, rotaxis) endif @@ -187,13 +187,13 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& do i=1,npinject 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 + 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 + 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) @@ -279,7 +279,7 @@ subroutine write_options_inject(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) + & 'wind_speed_factor', 'factor to scale the wind speed based on the Keplerian speed at rinject', iunit) end subroutine write_options_inject