Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added planetary destruction option and binding energy profile option to analysis_common_envelope.F90. #334

Merged
merged 1 commit into from
Nov 3, 2022
Merged

Added planetary destruction option and binding energy profile option to analysis_common_envelope.F90. #334

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
204 changes: 201 additions & 3 deletions src/utils/analysis_common_envelope.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)

!chose analysis type
if (dump_number==0) then
print "(29(a,/))", &
print "(31(a,/))", &
' 1) Sink separation', &
' 2) Bound and unbound quantities', &
' 3) Energies', &
Expand All @@ -96,6 +96,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
' 7) Simulation units and particle properties', &
' 8) Output .divv', &
' 9) EoS testing', &
'10) Planet destruction', &
'11) Profile of newly unbound particles', &
'12) Sink properties', &
'13) MESA EoS compute total entropy and other average td quantities', &
Expand All @@ -115,11 +116,12 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
'27) Print dumps number matching separation', &
'28) Companion mass coordinate vs. time', &
'29) Energy histogram',&
'30) Analyse disk'
'30) Analyse disk', &
'31) Binding energy profile'

analysis_to_perform = 1

call prompt('Choose analysis type ',analysis_to_perform,1,30)
call prompt('Choose analysis type ',analysis_to_perform,1,31)

endif

Expand Down Expand Up @@ -171,6 +173,8 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
call output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu)
case(9) !EoS testing
call eos_surfaces
case(10) !Planet destruction
call planet_destruction(time,npart,particlemass,xyzh,vxyzu)
case(11) !New unbound particle profiles in time
call unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu)
case(19) ! Rotation profile
Expand All @@ -197,6 +201,8 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
call energy_hist(time,npart,particlemass,xyzh,vxyzu)
case(30) ! Analyse disk around companion
call analyse_disk(num,npart,particlemass,xyzh,vxyzu)
case(31) !Binding energy profile
call create_bindingEnergy_profile(time,num,npart,particlemass,xyzh,vxyzu)
case(12) !sink properties
call sink_properties(time,npart,particlemass,xyzh,vxyzu)
case(13) !MESA EoS compute total entropy and other average thermodynamical quantities
Expand Down Expand Up @@ -2639,6 +2645,151 @@ subroutine J_E_plane(num,npart,particlemass,xyzh,vxyzu)

end subroutine J_E_plane

!-------------------------------------------------------------------
!+
! Planet destruction
!+
!-------------------------------------------------------------------
subroutine planet_destruction(time,npart,particlemass,xyzh,vxyzu)
use kernel, only:wkern
integer, intent(in) :: npart
real, intent(in) :: time,particlemass
real, intent(inout) :: xyzh(:,:),vxyzu(:,:)
character(len=17), allocatable :: columns(:)
character(len=18) :: filename
real, allocatable :: planetDestruction(:)
integer :: ncols,i,j
real, allocatable, save :: time_old
real, allocatable, save :: particleRho(:)
character(len=50) :: planetRadiusPromptString
real, allocatable, save :: planetRadii(:) !In units of Rsun

real, dimension(3) :: currentGasVel, currentVelContrast
real :: currentRho(1) !Is a one element array because sphInterpolation returns a 1 dimensional array.
real :: currentRhoScaled,currentVelContrastScaled,currentPlanetRhoScaled
real :: currentPlanetMassScaled,currentPlanetRadiusScaled
real, allocatable, save :: currentKhAblatedMass(:)

ncols=5
allocate(columns(ncols))
allocate(planetDestruction(ncols))
columns=(/" rhoGas", &
" kh_rhoCrit", &
" kh_lmax", &
" kh_mdot", &
" kh_ablatedM" /)

!Kelvin-Helmholtz instability planet destruction as described in "On the survival of brown dwarfs
!and planets by their giant host star" (https://arxiv.org/abs/1210.0879). Description of columns:
!rhoGas: background gas density at sink. In units of g/cm^3.
!kh_rhoCrit: paper equation 5. In units of g/cm^3.
!kh_lmax: paper equation 6. In units of Jupiter radii.
!kh_mdot: paper equation 7. In units of Jupiter mass/year.
!kh_ablatedM: kh_mdot integrated over time. In units of Jupiter masses.

do i=1,nptmass
if (i==1) cycle !The first sink is assumed to be the core.

if ((dump_number==0) .and. (i==2)) then !This is only done once.
allocate(planetRadii(nptmass))
planetRadii=0.1
do j=2,nptmass
write(planetRadiusPromptString,"(A13,I0,A32)") "Enter planet ",j-1," radius in units of solar radii"
call prompt(planetRadiusPromptString,planetRadii(i),0.0,1.0)
enddo

allocate(time_old)
allocate(particleRho(npart))
allocate(currentKhAblatedMass(nptmass))

time_old=0.0
particleRho=getParticleRho(xyzh(4,:),particlemass)
currentKhAblatedMass=0.0
endif


currentRho=sphInterpolation(npart,particlemass,particleRho,xyzh,xyzmh_ptmass(1:3,i),reshape(particleRho,(/1,npart/)))
currentGasVel=sphInterpolation(npart,particlemass,particleRho,xyzh,xyzmh_ptmass(1:3,i),vxyzu(1:3,:))
currentVelContrast=vxyz_ptmass(1:3,i)-currentGasVel

currentPlanetRadiusScaled=planetRadii(i)/0.1 !In units of 0.1 Rsun.
currentPlanetMassScaled=xyzmh_ptmass(4,i)*104.74 !In units of 10 jupiter masses.
currentPlanetRhoScaled=(xyzmh_ptmass(4,i)/((4.0/3.0)*pi*(planetRadii(i)**3.0)))*0.44 !In units of 13.34 g/cm^3
currentRhoScaled=currentRho(1)*59000.0 !In units of 10^-4 g/cm^3.
currentVelContrastScaled=distance(currentVelContrast)*4.37 !In units of 100 km/s.

planetDestruction(1)=currentRho(1)*5.9
planetDestruction(2)=3.82*(currentPlanetRhoScaled**(4.0/3.0))*(currentPlanetMassScaled**(2.0/3.0))&
*(currentVelContrastScaled**(-2.0))
planetDestruction(3)=0.0000263*(currentVelContrastScaled**2.0)*currentRhoScaled*(currentPlanetRhoScaled**((-5.0)/3.0))&
*(currentPlanetMassScaled**((-1.0)/3.0))
planetDestruction(4)=11.0*currentVelContrastScaled*currentRhoScaled*(currentPlanetRadiusScaled**2.0)&
*(planetDestruction(3)/(currentPlanetRadiusScaled*0.973))

currentKhAblatedMass(i)=currentKhAblatedMass(i)+((time-time_old)*planetDestruction(4)*0.0000505)
planetDestruction(5)=currentKhAblatedMass(i)


write(filename, "(A17,I0)") "sink_destruction_",i
call write_time_file(filename, columns, time, planetDestruction, ncols, dump_number)
enddo

time_old=time
end subroutine planet_destruction

!-----------------------------------------------------------------------------------------
!+
!Binding energy profile
!+
!-----------------------------------------------------------------------------------------
subroutine create_bindingEnergy_profile(time,num,npart,particlemass,xyzh,vxyzu)
real, intent(in) :: time,particlemass
integer, intent(in) :: num,npart
real, intent(in) :: xyzh(4,npart),vxyzu(4,npart)

character(len=17), allocatable :: columns(:)
real, allocatable :: profile(:,:)
integer :: ncols,i,j,iorder(npart)
real :: currentInteriorMass,currentParticleGPE,currentCoreParticleSeparation
real :: previousBindingEnergy,previousBindingEnergyU

ncols=3
allocate(columns(ncols))
allocate(profile(ncols,npart))
columns=(/" radius",&
" bEnergy",& !Binding energy without internal energy.
" bEnergy (u)"/) !Binding energy with internal energy.


call set_r2func_origin(xyzmh_ptmass(1,1),xyzmh_ptmass(2,1),xyzmh_ptmass(3,1))
call indexxfunc(npart,r2func_origin,xyzh,iorder)
currentInteriorMass=xyzmh_ptmass(4,1)+(npart*particlemass) !Initally set to the entire mass of the star.

do i=npart,1,-1 !Loops over all particles from outer to inner.
j=iorder(i)
currentInteriorMass=currentInteriorMass-particlemass
currentCoreParticleSeparation=separation(xyzmh_ptmass(1:3,1),xyzh(1:3,j))
currentParticleGPE=(currentInteriorMass*particlemass)/currentCoreParticleSeparation

!The binding energy at a particular radius is the sum of the gravitational potential energies
!(and internal energies in the case of the third column) of all particles beyond that radius.
if (i==npart) then
previousBindingEnergy=0.0
previousBindingEnergyU=0.0
else
previousBindingEnergy=profile(2,i+1)
previousBindingEnergyU=profile(3,i+1)
endif

profile(1,i)=currentCoreParticleSeparation
profile(2,i)=previousBindingEnergy+currentParticleGPE
profile(3,i)=previousBindingEnergyU+currentParticleGPE-(vxyzu(4,j)*particlemass)
enddo


call write_file('bEnergyProfile','bEnergyProfiles',columns,profile,npart,ncols,num)
end subroutine create_bindingEnergy_profile


subroutine get_core_gas_com(time,npart,xyzh,vxyzu)
use sortutils, only:set_r2func_origin,r2func_origin,indexxfunc
Expand Down Expand Up @@ -3459,6 +3610,53 @@ real function separation(a,b)
separation = distance(a - b)
end function separation

!Creates an array of SPH particle densities for each value of h.
elemental real function getParticleRho(h,particlemass)
real, intent(in) :: h,particlemass
getParticleRho=rhoh(h,particlemass)
end function getParticleRho

!Performs SPH interpolation on the SPH particle property toInterpolate at the location interpolateXyz.
!The smoothing length used is the smoothing length of the closest SPH particle to interpolateXyz.
function sphInterpolation(npart,particlemass,particleRho,particleXyzh,interpolateXyz,toInterpolate) result(interpolatedData)
use kernel, only:wkern
integer, intent(in) :: npart
real, intent(in) :: particlemass
real, intent(in) :: particleRho(npart)
real, intent(in) :: particleXyzh(4,npart)
real, intent(in) :: interpolateXyz(3)
real, intent(in) :: toInterpolate(:,:)
real :: interpolatedData(size(toInterpolate,1))

integer :: i,j,iorder(npart)
real :: currentR,currentQ,currentQ2
real :: nearestSphH
real :: currentParticleRho,currentSphSummandFactor

interpolatedData=0.0
call set_r2func_origin(interpolateXyz(1),interpolateXyz(2),interpolateXyz(3))
call indexxfunc(npart,r2func_origin,particleXyzh,iorder) !Gets the order of SPH particles from the interpolation point.
nearestSphH=particleXyzh(4,iorder(1)) !The smoothing length of the nearest SPH particle to the ineterpolation point.

do i=1,npart
j=iorder(i)

currentR=separation(interpolateXyz,particleXyzh(1:3,j))
currentQ=currentR/nearestSphH !currentR is scaled in units of nearestSphH
currentQ2=currentQ**2.0

!All SPH particles beyond 2 smoothing lengths are ignored.
if (currentQ>2) then
exit
endif

!SPH interpolation is done below.
currentParticleRho=particleRho(j)
currentSphSummandFactor=(particlemass/currentParticleRho)*((1.0/((nearestSphH**3.0)*pi))*wkern(currentQ2,currentQ))
interpolatedData=interpolatedData+(currentSphSummandFactor*toInterpolate(:,j))
enddo
end function sphInterpolation

!Sorting routines
recursive subroutine quicksort(a, first, last, ncols, sortcol)
integer, intent(in) :: first, last, ncols, sortcol
Expand Down