From 1bca84ca928dec833d69ae654259e5e116c71e2c Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Fri, 19 Jan 2024 00:21:18 +1100 Subject: [PATCH 1/8] (windtunnel) add option to have wind only, no star --- src/main/inject_windtunnel.f90 | 19 ++++++++++++++----- src/setup/setup_windtunnel.f90 | 32 +++++++++++++++++--------------- 2 files changed, 31 insertions(+), 20 deletions(-) diff --git a/src/main/inject_windtunnel.f90 b/src/main/inject_windtunnel.f90 index 7c304db07..f73f32bf2 100644 --- a/src/main/inject_windtunnel.f90 +++ b/src/main/inject_windtunnel.f90 @@ -27,10 +27,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 + set_default_options_inject,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. @@ -48,7 +51,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(:,:) @@ -72,6 +75,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 @@ -143,11 +152,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!') @@ -212,7 +221,7 @@ 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 diff --git a/src/setup/setup_windtunnel.f90 b/src/setup/setup_windtunnel.f90 index 2c8e20ba4..e5527dcbd 100644 --- a/src/setup/setup_windtunnel.f90 +++ b/src/setup/setup_windtunnel.f90 @@ -19,7 +19,7 @@ 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 @@ -132,22 +132,24 @@ 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) + ! 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 From 51bdb87019d583287689c9870c23f9b473164302 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Tue, 20 Feb 2024 10:59:05 +0100 Subject: [PATCH 2/8] (CE-analysis) calculate planet sep and vel wrt origin if no ptmasses are present --- src/utils/analysis_common_envelope.f90 | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 6cd1c6c27..d1c4a59bd 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -128,7 +128,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call prompt('Choose analysis type ',analysis_to_perform,1,41) 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) @@ -440,7 +439,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 @@ -470,6 +469,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 @@ -492,11 +499,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. From d49dd3142b24ab36aa87da55274a65162f8135a5 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Tue, 20 Feb 2024 21:01:15 +1100 Subject: [PATCH 3/8] (moddump_removeparticles) tweaks to prompt message --- src/utils/moddump_removeparticles_radius.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/utils/moddump_removeparticles_radius.f90 b/src/utils/moddump_removeparticles_radius.f90 index dd9ada106..057fbf6f8 100644 --- a/src/utils/moddump_removeparticles_radius.f90 +++ b/src/utils/moddump_removeparticles_radius.f90 @@ -17,7 +17,7 @@ module moddump ! :Dependencies: part, prompting ! - use part, only:delete_particles_outside_sphere + use part, only:delete_particles_outside_sphere,delete_particles_inside_radius use prompting, only:prompt implicit none @@ -47,13 +47,13 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) call prompt('Deleting particles inside a given radius ?',icutinside) call prompt('Deleting particles outside a given radius ?',icutoutside) if (icutinside) then - call prompt('Enter inward radius in au',inradius,0.) + call prompt('Enter inward radius in code units',inradius,0.) call prompt('Enter x coordinate of the center of that sphere',incenter(1)) call prompt('Enter y coordinate of the center of that sphere',incenter(2)) call prompt('Enter z coordinate of the center of that sphere',incenter(3)) endif if (icutoutside) then - call prompt('Enter outward radius in au',outradius,0.) + call prompt('Enter outward radius in code units',outradius,0.) call prompt('Enter x coordinate of the center of that sphere',outcenter(1)) call prompt('Enter y coordinate of the center of that sphere',outcenter(2)) call prompt('Enter z coordinate of the center of that sphere',outcenter(3)) @@ -62,13 +62,13 @@ subroutine modify_dump(npart,npartoftype,massoftype,xyzh,vxyzu) if (icutinside) then print*,'Phantommoddump: Remove particles inside a particular radius' print*,'Removing particles inside radius ',inradius - call delete_particles_outside_sphere(incenter,inradius,revert=.true.) + call delete_particles_inside_radius(incenter,inradius,npart,npartoftype) endif if (icutoutside) then print*,'Phantommoddump: Remove particles outside a particular radius' print*,'Removing particles outside radius ',outradius - call delete_particles_outside_sphere(outcenter,outradius) + call delete_particles_outside_sphere(outcenter,outradius,npart) endif end subroutine modify_dump From b3c219267aac46312c30f1ea943e132cecf942ba Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Tue, 20 Feb 2024 14:30:10 +0100 Subject: [PATCH 4/8] (windtunnel) add option to hold the gaseous sphere still inside the wind --- src/main/inject_windtunnel.f90 | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/src/main/inject_windtunnel.f90 b/src/main/inject_windtunnel.f90 index f73f32bf2..7c310811c 100644 --- a/src/main/inject_windtunnel.f90 +++ b/src/main/inject_windtunnel.f90 @@ -14,6 +14,7 @@ module inject ! :Runtime parameters: ! - lattice_type : *0: cubic distribution, 1: closepacked distribution* ! - handled_layers : *(integer) number of handled BHL wind layers* +! - hold_star : *1: subtract CM velocity of star particles at each timestep* ! - v_inf : *BHL wind speed* ! - Rstar : *BHL star radius (in accretion radii)* ! - BHL_radius : *radius of the wind cylinder (in star radii)* @@ -42,6 +43,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. @@ -228,6 +230,8 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& 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,vxyzu) + end subroutine inject_particles ! @@ -260,6 +264,26 @@ subroutine inject_or_update_particles(ifirst, n, position, velocity, h, u, bound end subroutine inject_or_update_particles +!----------------------------------------------------------------------- +!+ +! Subtracts centre-of-mass motion of star particles +! Assumes star particles have particle IDs 1 to nbulk +!+ +!----------------------------------------------------------------------- +subroutine subtract_star_vcom(nbulk,vxyzu) + integer, intent(in) :: nbulk + real, intent(inout) :: vxyzu(:,:) + real :: vstar(3) + integer :: i + + vstar = (/ sum(vxyzu(1,1:nbulk)), sum(vxyzu(2,1:nbulk)), sum(vxyzu(3,1:nbulk)) /) / real(nbulk) + do i=1,nbulk + vxyzu(1:3,i) = vxyzu(1:3,i) - vstar + enddo + +end subroutine subtract_star_vcom + + !----------------------------------------------------------------------- !+ ! Print summary of wind properties (assumes inputs are in code units) @@ -306,6 +330,7 @@ subroutine write_options_inject(iunit) call write_inopt(nstar,'nstar','No. of particles making up sphere',iunit) 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 +393,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) From 91867454b56413d74c0ac07a27cb87ec8fb2c05d Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Tue, 20 Feb 2024 14:30:28 +0100 Subject: [PATCH 5/8] (windtunnel) write number of star particles into infile --- src/setup/setup_windtunnel.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/setup/setup_windtunnel.f90 b/src/setup/setup_windtunnel.f90 index e5527dcbd..8bc49e23e 100644 --- a/src/setup/setup_windtunnel.f90 +++ b/src/setup/setup_windtunnel.f90 @@ -164,6 +164,7 @@ end subroutine setpart !----------------------------------------------------------------------- subroutine write_setupfile(filename) use infile_utils, only:write_inopt + use part, only:npart use dim, only:tagline use eos, only:gamma use setunits, only:write_options_units @@ -179,7 +180,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(npart,'nstar','number of particles resolving gas sphere',iunit) ! note: npart is output of set_sphere call write_inopt(Mstar,'Mstar','sphere mass in code units',iunit) call write_inopt(Rstar,'Rstar','sphere radius in code units',iunit) From 7d85592eb990dcca4b08a4083e8f433b1d45828d Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Wed, 21 Feb 2024 16:39:46 +0100 Subject: [PATCH 6/8] (windtunnel) subtract CM velocity of star excluding ablated particles --- src/main/inject_windtunnel.f90 | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/src/main/inject_windtunnel.f90 b/src/main/inject_windtunnel.f90 index 7c310811c..d69b2cc5f 100644 --- a/src/main/inject_windtunnel.f90 +++ b/src/main/inject_windtunnel.f90 @@ -230,7 +230,7 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& 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,vxyzu) + if ((hold_star>0) .and. (.not. windonly)) call subtract_star_vcom(nstarpart,xyzh,vxyzu) end subroutine inject_particles @@ -267,19 +267,32 @@ end subroutine inject_or_update_particles !----------------------------------------------------------------------- !+ ! Subtracts centre-of-mass motion of star particles -! Assumes star particles have particle IDs 1 to nbulk +! Assumes star particles have particle IDs 1 to nsphere !+ !----------------------------------------------------------------------- -subroutine subtract_star_vcom(nbulk,vxyzu) - integer, intent(in) :: nbulk +subroutine subtract_star_vcom(nsphere,xyzh,vxyzu) + integer, intent(in) :: nsphere + real, intent(in) :: xyzh(:,:) real, intent(inout) :: vxyzu(:,:) real :: vstar(3) - integer :: i - - vstar = (/ sum(vxyzu(1,1:nbulk)), sum(vxyzu(2,1:nbulk)), sum(vxyzu(3,1:nbulk)) /) / real(nbulk) - do i=1,nbulk - vxyzu(1:3,i) = vxyzu(1:3,i) - vstar + 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 From 69eb93e109d94ec2b8659fafc45684bea4193786 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Wed, 21 Feb 2024 18:07:50 +0100 Subject: [PATCH 7/8] (windtunnel) add call to relax_star when setting polytrope sphere --- src/setup/setup_windtunnel.f90 | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/setup/setup_windtunnel.f90 b/src/setup/setup_windtunnel.f90 index 8bc49e23e..dda592701 100644 --- a/src/setup/setup_windtunnel.f90 +++ b/src/setup/setup_windtunnel.f90 @@ -40,6 +40,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 @@ -57,9 +58,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 @@ -140,6 +141,10 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, 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 + + 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))) From ff2116db96c36a50fed22f2febd94806aaffb9a2 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Tue, 29 Oct 2024 12:28:28 +0100 Subject: [PATCH 8/8] fix unresolved merge conflict --- src/utils/moddump_removeparticles_radius.f90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/utils/moddump_removeparticles_radius.f90 b/src/utils/moddump_removeparticles_radius.f90 index 85d68818e..65bad1b90 100644 --- a/src/utils/moddump_removeparticles_radius.f90 +++ b/src/utils/moddump_removeparticles_radius.f90 @@ -17,11 +17,7 @@ module moddump ! :Dependencies: part, prompting ! -<<<<<<< HEAD - use part, only:delete_particles_outside_sphere,delete_particles_inside_radius -======= use part, only:delete_particles_outside_sphere,igas,idust ->>>>>>> e805ed68f91e4a807462a786fb1447b541deb978 use prompting, only:prompt implicit none