Improvements to windtunnel setup

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)

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)
 

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.