diff --git a/src/main/inject_windtunnel.f90 b/src/main/inject_windtunnel.f90 index 62cdf2f33..07f1fa166 100644 --- a/src/main/inject_windtunnel.f90 +++ b/src/main/inject_windtunnel.f90 @@ -17,6 +17,7 @@ module inject ! - BHL_radius : *radius of the wind cylinder (in star radii)* ! - Rstar : *sphere radius (code units)* ! - handled_layers : *(integer) number of handled BHL wind layers* +! - hold_star : *1: subtract CM velocity of star particles at each timestep* ! - lattice_type : *0: cubic distribution, 1: closepacked distribution* ! - nstar : *No. of particles making up sphere* ! - pres_inf : *ambient pressure (code units)* @@ -32,10 +33,13 @@ module inject character(len=*), parameter, public :: inject_type = 'windtunnel' public :: init_inject,inject_particles,write_options_inject,read_options_inject,& - set_default_options_inject,update_injected_par + set_default_options_inject,update_injected_par,windonly ! !--runtime settings for this module ! + + logical :: windonly = .false. + ! Main parameters: model MS6 from Ruffert & Arnett (1994) real, public :: v_inf = 1. real, public :: rho_inf = 1. @@ -44,6 +48,7 @@ module inject integer, public :: nstar = 0 ! Particle-related parameters + integer, public :: hold_star = 0 integer, public :: lattice_type = 1 integer, public :: handled_layers = 4 real, public :: wind_radius = 30. @@ -53,7 +58,7 @@ module inject private real :: wind_rad,wind_x,psep,distance_between_layers,& time_between_layers,h_inf,u_inf - integer :: max_layers,max_particles,nodd,neven + integer :: max_layers,max_particles,nodd,neven,nstarpart logical :: first_run = .true. real, allocatable :: layer_even(:,:),layer_odd(:,:) @@ -77,6 +82,12 @@ subroutine init_inject(ierr) ierr = 0 + if (windonly) then + nstarpart = 0 + else + nstarpart = nstar + endif + u_inf = pres_inf / (rho_inf*(gamma-1.)) cs_inf = sqrt(gamma*pres_inf/rho_inf) mach = v_inf/cs_inf @@ -148,11 +159,11 @@ subroutine init_inject(ierr) endif h_inf = hfact*(pmass/rho_inf)**(1./3.) max_layers = int(wind_length*Rstar/distance_between_layers) - max_particles = int(max_layers*(nodd+neven)/2) + nstar + max_particles = int(max_layers*(nodd+neven)/2) + nstarpart time_between_layers = distance_between_layers/v_inf call print_summary(v_inf,cs_inf,rho_inf,pres_inf,mach,pmass,distance_between_layers,& - time_between_layers,max_layers,nstar,max_particles) + time_between_layers,max_layers,nstarpart,max_particles) if (max_particles > maxp) call fatal('windtunnel', 'maxp too small for this simulation, please increase MAXP!') @@ -217,13 +228,15 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& endif print *, np, ' particles (npart=', npart, '/', max_particles, ')' endif - call inject_or_update_particles(i_part+nstar+1, np, xyz, vxyz, h, u, .false.) + call inject_or_update_particles(i_part+nstarpart+1, np, xyz, vxyz, h, u, .false.) deallocate(xyz, vxyz, h, u) enddo irrational_number_close_to_one = 3./pi dtinject = (irrational_number_close_to_one*time_between_layers)/utime + if ((hold_star>0) .and. (.not. windonly)) call subtract_star_vcom(nstarpart,xyzh,vxyzu) + end subroutine inject_particles ! @@ -260,6 +273,39 @@ subroutine update_injected_par ! -- does not do anything and will never be used end subroutine update_injected_par +!----------------------------------------------------------------------- +!+ +! Subtracts centre-of-mass motion of star particles +! Assumes star particles have particle IDs 1 to nsphere +!+ +!----------------------------------------------------------------------- +subroutine subtract_star_vcom(nsphere,xyzh,vxyzu) + integer, intent(in) :: nsphere + real, intent(in) :: xyzh(:,:) + real, intent(inout) :: vxyzu(:,:) + real :: vstar(3) + integer :: i,nbulk + +! vstar = (/ sum(vxyzu(1,1:nsphere)), sum(vxyzu(2,1:nsphere)), sum(vxyzu(3,1:nsphere)) /) / real(nsphere) + nbulk = 0 + vstar = 0. + do i=1,nsphere + if (xyzh(1,i) < 2.*Rstar) then + vstar = vstar + vxyzu(1:3,i) + nbulk = nbulk + 1 + endif + enddo + vstar = vstar/real(nbulk) + + do i=1,nsphere + if (xyzh(1,i) < 2.*Rstar) then + vxyzu(1:3,i) = vxyzu(1:3,i) - vstar + endif +enddo + +end subroutine subtract_star_vcom + + !----------------------------------------------------------------------- !+ ! Print summary of wind properties (assumes inputs are in code units) @@ -303,9 +349,10 @@ subroutine write_options_inject(iunit) call write_inopt(pres_inf,'pres_inf','ambient pressure (code units)',iunit) call write_inopt(rho_inf,'rho_inf','ambient density (code units)',iunit) call write_inopt(Rstar,'Rstar','sphere radius (code units)',iunit) - call write_inopt(nstar,'nstar','No. of particles making up sphere',iunit) + call write_inopt(nstar,'nstar','No. of particles making up sphere',iunit) ! need to write actual no. of particles, not nstar_in call write_inopt(lattice_type,'lattice_type','0: cubic distribution, 1: closepacked distribution',iunit) call write_inopt(handled_layers,'handled_layers','(integer) number of handled BHL wind layers',iunit) + call write_inopt(hold_star,'hold_star','1: subtract CM velocity of star particles at each timestep',iunit) call write_inopt(wind_radius,'BHL_radius','radius of the wind cylinder (in star radii)',iunit) call write_inopt(wind_injection_x,'wind_injection_x','x position of the wind injection boundary (in star radii)',iunit) call write_inopt(wind_length,'wind_length','crude wind length (in star radii)',iunit) @@ -368,9 +415,12 @@ subroutine read_options_inject(name,valstring,imatch,igotall,ierr) read(valstring,*,iostat=ierr) wind_length ngot = ngot + 1 if (wind_length <= 0.) call fatal(label,'wind_length must be positive') + case('hold_star') + read(valstring,*,iostat=ierr) hold_star + ngot = ngot + 1 end select - igotall = (ngot >= 10) + igotall = (ngot >= 11) end subroutine read_options_inject subroutine set_default_options_inject(flag) diff --git a/src/setup/setup_windtunnel.f90 b/src/setup/setup_windtunnel.f90 index f4df9bffb..0bbc21f56 100644 --- a/src/setup/setup_windtunnel.f90 +++ b/src/setup/setup_windtunnel.f90 @@ -33,12 +33,13 @@ module setup use io, only:master,fatal use inject, only:init_inject,nstar,Rstar,lattice_type,handled_layers,& wind_radius,wind_injection_x,wind_length,& - rho_inf,pres_inf,v_inf + rho_inf,pres_inf,v_inf,windonly implicit none public :: setpart real :: Mstar + integer :: nstar_in ! guess for no. of star prticles private @@ -54,6 +55,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, use eos, only:ieos,gmw use setstar_utils,only:set_star_density use rho_profile, only:rho_polytrope + use relaxstar, only:relax_star use extern_densprofile, only:nrhotab use physcon, only:solarm,solarr use units, only:udist,umass,utime,set_units,unit_velocity,unit_density,unit_pressure @@ -71,9 +73,9 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, real, intent(inout) :: time character(len=20), intent(in) :: fileprefix real :: rhocentre,rmin,pmass,densi,presi,ri - real, allocatable :: r(:),den(:),pres(:) - integer :: ierr,npts,np,i - logical :: use_exactN,setexists + real, allocatable :: r(:),den(:),pres(:),Xfrac(:),Yfrac(:),mu(:) + integer :: ierr,ierr_relax,npts,np,i + logical :: use_exactN,setexists,use_var_comp character(len=30) :: lattice character(len=120) :: setupfile @@ -100,7 +102,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! Star parameters Rstar = 0.1 Mstar = 1.e-3 - nstar = 1000 + nstar_in = 1000 lattice = 'closepacked' use_exactN = .true. @@ -131,7 +133,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, stop 'please check and edit .setup file and rerun phantomsetup' endif - pmass = Mstar / real(nstar) + pmass = Mstar / real(nstar_in) massoftype(igas) = pmass call check_setup(pmass,ierr) if (ierr /= 0) call fatal('windtunnel','errors in setup parameters') @@ -146,23 +148,29 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, vxyzu(:,:) = 0. ! Set polytropic star - allocate(r(nrhotab),den(nrhotab),pres(nrhotab)) - call rho_polytrope(gamma,polyk,Mstar,r,den,npts,rhocentre,set_polyk=.true.,Rstar=Rstar) - pres = polyk*den**gamma - rmin = r(1) - call set_star_density(lattice,id,master,rmin,Rstar,Mstar,hfact,& + if (.not. windonly) then + allocate(r(nrhotab),den(nrhotab),pres(nrhotab)) + call rho_polytrope(gamma,polyk,Mstar,r,den,npts,rhocentre,set_polyk=.true.,Rstar=Rstar) + pres = polyk*den**gamma + rmin = r(1) + call set_star_density(lattice,id,master,rmin,Rstar,Mstar,hfact,& npts,den,r,npart,npartoftype,massoftype,xyzh,& use_exactN,np,rhozero,npart_total,i_belong) ! Note: mass_is_set = .true., so np is not used - ! Set thermal energy - do i = 1,npart - ri = sqrt(dot_product(xyzh(1:3,i),xyzh(1:3,i))) - densi = yinterp(den(1:npts),r(1:npts),ri) - presi = yinterp(pres(1:npts),r(1:npts),ri) - vxyzu(4,i) = presi / ( (gamma-1.) * densi) - enddo - - deallocate(r,den,pres) - + nstar = npart + use_var_comp = .false. + call relax_star(npts,den,pres,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,ierr_relax) + + ! Set thermal energy + do i = 1,npart + ri = sqrt(dot_product(xyzh(1:3,i),xyzh(1:3,i))) + densi = yinterp(den(1:npts),r(1:npts),ri) + presi = yinterp(pres(1:npts),r(1:npts),ri) + vxyzu(4,i) = presi / ( (gamma-1.) * densi) + enddo + + deallocate(r,den,pres) + endif + print*, "udist = ", udist, "; umass = ", umass, "; utime = ", utime end subroutine setpart @@ -191,7 +199,7 @@ subroutine write_setupfile(filename) call write_options_units(iunit) write(iunit,"(/,a)") '# sphere settings' - call write_inopt(nstar,'nstar','number of particles resolving gas sphere',iunit) + call write_inopt(nstar_in,'nstar','number of particles resolving gas sphere',iunit) ! note: this is an estimate, actual no. of particles is npart outputted from set_sphere call write_inopt(Mstar,'Mstar','sphere mass in code units',iunit) call write_inopt(Rstar,'Rstar','sphere radius in code units',iunit) @@ -237,7 +245,7 @@ subroutine read_setupfile(filename,ierr) call read_options_and_set_units(db,nerr) - call read_inopt(nstar,'nstar',db,errcount=nerr) + call read_inopt(nstar_in,'nstar',db,errcount=nerr) call read_inopt(Mstar,'Mstar',db,errcount=nerr) call read_inopt(Rstar,'Rstar',db,errcount=nerr) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index a527f4b41..179b87a02 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -106,7 +106,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call prompt('Choose analysis type ',analysis_to_perform,1,38) endif - call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass) call adjust_corotating_velocities(npart,particlemass,xyzh,vxyzu,& xyzmh_ptmass,vxyz_ptmass,omega_corotate,dump_number) @@ -338,7 +337,7 @@ subroutine planet_rvm(time,particlemass,xyzh,vxyzu) real, intent(in) :: time,xyzh(:,:),vxyzu(:,:),particlemass character(len=17), allocatable :: columns(:) real, dimension(3) :: planet_com,planet_vel,sep,vel - real :: rhoi,rhoprev,sepi,si,smin,presi,Rthreshold + real :: rhoi,rhoprev,sepi,si,smin,presi,Rthreshold,xyz_origin(3),vxyz_origin(3) real, allocatable :: data_cols(:),mass(:),vthreshold(:) integer :: i,j,ncols,maxrho_ID,ientropy,Nmasks integer, save :: nplanet @@ -368,6 +367,14 @@ subroutine planet_rvm(time,particlemass,xyzh,vxyzu) if (dump_number == 0) call get_planetIDs(nplanet,planetIDs) isfulldump = (vxyzu(4,1) > 0.) + if (nptmass > 0) then + xyz_origin = xyzmh_ptmass(1:3,1) + vxyz_origin = vxyz_ptmass(1:3,1) + else + xyz_origin = (/0.,0.,0./) + vxyz_origin = (/0.,0.,0./) + endif + ! Find highest density and lowest entropy in planet rhoprev = 0. maxrho_ID = 1 @@ -390,11 +397,11 @@ subroutine planet_rvm(time,particlemass,xyzh,vxyzu) enddo planet_com = xyzh(1:3,maxrho_ID) - sep = planet_com - xyzmh_ptmass(1:3,1) + sep = planet_com - xyz_origin(1:3) if (isfulldump) then planet_vel = vxyzu(1:3,maxrho_ID) - vel = planet_vel - vxyz_ptmass(1:3,1) + vel = planet_vel - vxyz_origin(1:3) else vel = 0. smin = 0.