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

Jump to bottom

Compute total enclosed mass in self-gravitating disc initial conditions #427

Merged
merged 3 commits into from
Jun 6, 2023
Merged

Compute total enclosed mass in self-gravitating disc initial conditions #427

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
914c867
Trying to fix enclosed mass calculation
crislong Jun 6, 2023
bd5cba3
Merge https://github.com/danieljprice/phantom into dust_enclosed_mass
crislong Jun 6, 2023
413c837
bug fixed. now the enclosed mass of gas and dust is computed correctly
crislong Jun 6, 2023
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
43 changes: 32 additions & 11 deletions src/setup/set_disc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ module setdisc
public :: set_disc,set_incline_or_warp,get_disc_mass,scaled_sigma

private
integer, parameter :: maxbins = 4096
integer, parameter, public :: maxbins = 4096

contains

Expand All @@ -73,7 +73,7 @@ subroutine set_disc(id,master,mixture,nparttot,npart,npart_start,rmin,rmax, &
disc_mass,disc_massdust,sig_norm,star_mass,xyz_origin,vxyz_origin, &
particle_type,particle_mass,hfact,xyzh,vxyzu,polyk, &
position_angle,inclination,ismooth,alpha,rwarp,warp_smoothl, &
bh_spin,bh_spin_angle,rref,writefile,ierr,prefix,verbose)
bh_spin,bh_spin_angle,rref,enc_mass,r_grid,writefile,ierr,prefix,verbose)
use io, only:stdout
use part, only:maxp,idust,maxtypes
use centreofmass, only:get_total_angular_momentum
Expand All @@ -92,6 +92,7 @@ subroutine set_disc(id,master,mixture,nparttot,npart,npart_start,rmin,rmax, &
integer, optional, intent(in) :: particle_type
real, optional, intent(in) :: position_angle,inclination
real, optional, intent(in) :: rwarp,warp_smoothl,bh_spin,bh_spin_angle
real, optional, intent(in) :: enc_mass(maxbins),r_grid(maxbins)
logical, optional, intent(in) :: ismooth,mixture
real, intent(out) :: xyzh(:,:)
real, intent(out) :: vxyzu(:,:)
Expand Down Expand Up @@ -284,6 +285,15 @@ subroutine set_disc(id,master,mixture,nparttot,npart,npart_start,rmin,rmax, &
!
!--disc mass and sigma normalisation
!
if (present(r_grid)) then
rad = r_grid
else
do i=1,maxbins
rad(i) = R_in + (i-1) * (R_out-R_in)/real(maxbins-1)
!R_in + 0.5*(R_out-R_in)/maxbins
enddo
endif

if (present(sig_norm)) then
if (present(disc_mass)) then
call fatal('set_disc','cannot set disc_mass and sig_norm at same time')
Expand All @@ -306,12 +316,14 @@ subroutine set_disc(id,master,mixture,nparttot,npart,npart_start,rmin,rmax, &
sigma_norm = 0.
call fatal('set_disc','need to set disc mass directly or via sigma normalisation')
endif
enc_m = enc_m + star_m
!print*, 'encm in setdisc', enc_m(1:20)
!
!--dust mass
!
sigma_normdust = 1.d0
if (do_mixture) then
!--sigma_normdust set from dust disc mass
sigma_normdust = 1.d0
call get_disc_mass(disc_mdust,enc_m_tmp,rad_tmp,Q_tmp,sigmaprofiledust, &
sigma_normdust,star_m,p_indexdust,q_inddust, &
R_indust,R_outdust,R_ref,R_c_dust,H_Rdust)
Expand Down Expand Up @@ -345,14 +357,23 @@ subroutine set_disc(id,master,mixture,nparttot,npart,npart_start,rmin,rmax, &
R_indust,R_outdust,phi_min,phi_max,sigma_norm,sigma_normdust,&
sigmaprofile,sigmaprofiledust,R_c,R_c_dust,p_index,p_inddust,cs0,cs0dust,&
q_index,q_inddust,star_m,G,particle_mass,hfact,itype,xyzh,honH,do_verbose)
!
!--set particle velocities
!

if (present(inclination)) then
incl = inclination
else
incl = 0.
endif
!
!--override enclosed mass with the correct version
! for the case where multiple calls to set_disc are used
! e.g. for discs with gas and dust
!
if (present(enc_mass)) then
enc_m = enc_mass
endif
!
!--set particle velocities
!
call set_disc_velocities(npart_tot,npart_start_count,itype,G,star_m,aspin,aspin_angle, &
clight,cs0,exponential_taper,p_index,q_index,gamma,R_in, &
rad,enc_m,smooth_surface_density,xyzh,vxyzu,incl)
Expand Down Expand Up @@ -656,6 +677,7 @@ subroutine set_disc_velocities(npart_tot,npart_start_count,itype,G,star_m,aspin,

ierr = 0
ipart = npart_start_count - 1

do i=npart_start_count,npart_tot
if (i_belong_i4(i)) then
ipart = ipart + 1
Expand Down Expand Up @@ -1151,7 +1173,8 @@ subroutine get_disc_mass(disc_m,enc_m,rad,toomre_min,sigmaprofile,sigma_norm, &
real, intent(in) :: sigma_norm,star_m,pindex,qindex,R_in,R_out,R_ref,H_R
real, optional, intent(in) :: R_c
integer, intent(in) :: sigmaprofile
real, intent(out) :: disc_m,enc_m(:),rad(:),toomre_min
real, intent(in) :: rad(:)
real, intent(out) :: disc_m,enc_m(:),toomre_min

real :: dr,dM,R,sigma,cs0,cs,kappa,G
integer :: i
Expand All @@ -1161,23 +1184,21 @@ subroutine get_disc_mass(disc_m,enc_m,rad,toomre_min,sigmaprofile,sigma_norm, &
enc_m = 0.
toomre_min = huge(toomre_min)
disc_m = 0.
dR = (R_out-R_in)/real(maxbins-1)
dR = rad(2)-rad(1)
do i=1,maxbins
R = R_in + (i-1)*dR
R = rad(i)
sigma = sigma_norm * scaled_sigma(R,sigmaprofile,pindex,R_ref,R_in,R_out,R_c)
!--disc mass
dM = 2.*pi*R*sigma*dR
disc_m = disc_m + dM
enc_m(i) = disc_m
rad(i) = R + 0.5*dR
!--Toomre Q
cs = cs_func(cs0,R,qindex)
kappa = sqrt(G*star_m/R**3)
if (sigma > epsilon(sigma)) then
toomre_min = min(toomre_min,real(cs*kappa/(pi*G*sigma)))
endif
enddo
enc_m = enc_m + star_m

end subroutine get_disc_mass

Expand Down
54 changes: 37 additions & 17 deletions src/setup/setup_disc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,7 @@ module setup
ndustlarge,grainsize,graindens,nptmass,iamtype,dustgasprop,&
VrelVf,rad,radprop,ikappa,iradxi
use physcon, only:au,solarm,jupiterm,earthm,pi,years
use setdisc, only:scaled_sigma,get_disc_mass
use setdisc, only:scaled_sigma,get_disc_mass,maxbins
use set_dust_options, only:set_dust_default_options,dust_method,dust_to_gas,&
ndusttypesinp,ndustlargeinp,ndustsmallinp,isetdust,&
dustbinfrac,check_dust_method
Expand Down Expand Up @@ -198,6 +198,8 @@ module setup
real :: pindex_dust(maxdiscs,maxdusttypes),qindex_dust(maxdiscs,maxdusttypes)
real :: H_R_dust(maxdiscs,maxdusttypes)

real :: enc_mass(maxbins,maxdiscs)

!--planets
integer, parameter :: maxplanets = 9

Expand Down Expand Up @@ -1004,17 +1006,32 @@ end subroutine setup_dust_grain_distribution
subroutine calculate_disc_mass()

integer :: i,j
integer, parameter :: maxbins = 4096

real :: enc_m(maxbins),rad(maxbins)
real :: Q_mintmp,disc_mtmp,annulus_mtmp
real :: rgrid_min,rgrid_max,fac

totmass_gas = 0.
disc_mdust = 0.

do i=1,maxdiscs
if (iuse_disc(i)) then

!
!--set up a common radial grid for the enclosed mass including gas and dust
! even if the gas/dust discs have different radial extents
!
rgrid_min = R_in(i)
rgrid_max = R_out(i)
if (isetgas(i)==1) then
rgrid_min = min(rgrid_min,R_inann(i))
rgrid_max = max(rgrid_max,R_outann(i))
endif
if (use_dust) then
rgrid_min = min(rgrid_min,minval(R_indust_swap(i,1:ndusttypes)))
rgrid_max = min(rgrid_max,maxval(R_outdust_swap(i,1:ndusttypes)))
endif
do j=1,maxbins
rad(j) = rgrid_min + (j-1) * (rgrid_max-rgrid_min)/real(maxbins-1)
enddo
!--gas discs
select case(isetgas(i))
case (0)
Expand All @@ -1023,7 +1040,10 @@ subroutine calculate_disc_mass()
call get_disc_mass(disc_mtmp,enc_m,rad,Q_mintmp,sigmaprofilegas(i),sig_norm(i), &
star_m(i),pindex(i),qindex(i),R_in(i),R_out(i),R_ref(i),R_c(i), &
H_R(i))
sig_norm(i) = sig_norm(i) * disc_m(i) / disc_mtmp
fac = disc_m(i) / disc_mtmp
sig_norm(i) = sig_norm(i) * fac
enc_m = enc_m * fac

case (1)
!--set disc mass from annulus mass
sig_norm(i) = 1.d0
Expand Down Expand Up @@ -1053,16 +1073,19 @@ subroutine calculate_disc_mass()
call get_disc_mass(disc_mtmp,enc_m,rad,Q_mintmp,sigmaprofilegas(i),sig_norm(i), &
star_m(i),pindex(i),qindex(i),R_in(i),R_out(i),R_ref(i),R_c(i), &
H_R(i))
sig_norm(i) = sig_norm(i) * Q_mintmp / Q_min(i)
fac = Q_mintmp / Q_min(i)
sig_norm(i) = sig_norm(i) * fac
!--recompute actual disc mass and Toomre Q
call get_disc_mass(disc_m(i),enc_m,rad,Q_min(i),sigmaprofilegas(i),sig_norm(i), &
star_m(i),pindex(i),qindex(i),R_in(i),R_out(i),R_ref(i),R_c(i), &
H_R(i))
end select

totmass_gas = totmass_gas + disc_m(i)
enc_mass(:,i) = enc_m + star_m(i)

!--dust discs
print*,'dust'
if (use_dust) then
disc_mdust(i,:) = 0.
do j=1,ndusttypes
Expand All @@ -1071,7 +1094,9 @@ subroutine calculate_disc_mass()
call get_disc_mass(disc_mtmp,enc_m,rad,Q_mintmp,sigmaprofiledust(i,j), &
sig_normdust(i,j),star_m(i),pindex_dust(i,j),qindex_dust(i,j), &
R_indust_swap(i,j),R_outdust_swap(i,j),R_ref(i),R_c_dust(i,j),H_R_dust(i,j))
sig_normdust(i,j) = sig_normdust(i,j) * disc_mdust(i,j) / disc_mtmp
fac = disc_mdust(i,j) / disc_mtmp
sig_normdust(i,j) = sig_normdust(i,j) * fac
enc_mass(:,i) = enc_mass(:,i) + enc_m(:)*fac
enddo
endif
endif
Expand Down Expand Up @@ -1174,7 +1199,6 @@ subroutine setup_discs(id,fileprefix,hfact,gamma,npart,polyk,&
star_m(i) = m2

call get_hierarchical_level_com(disclabel, xorigini, vorigini, xyzmh_ptmass, vxyz_ptmass, fileprefix)
!print*,disclabel,' com_pos ', xorigini

endif

Expand Down Expand Up @@ -1237,6 +1261,7 @@ subroutine setup_discs(id,fileprefix,hfact,gamma,npart,polyk,&
rwarp = R_warp(i), &
warp_smoothl = H_warp(i), &
bh_spin = bhspin, &
enc_mass = enc_mass(:,i), &
prefix = prefix)

!--set dustfrac
Expand Down Expand Up @@ -1291,6 +1316,7 @@ subroutine setup_discs(id,fileprefix,hfact,gamma,npart,polyk,&
rwarp = R_warp(i), &
warp_smoothl = H_warp(i), &
bh_spin = bhspin, &
enc_mass = enc_mass(:,i), &
prefix = dustprefix(j))

npart = npart + npindustdisc
Expand Down Expand Up @@ -1332,6 +1358,7 @@ subroutine setup_discs(id,fileprefix,hfact,gamma,npart,polyk,&
rwarp = R_warp(i), &
warp_smoothl = H_warp(i), &
bh_spin = bhspin, &
enc_mass = enc_mass(:,i), &
prefix = prefix)

npart = npart + npingasdisc
Expand Down Expand Up @@ -1383,6 +1410,7 @@ subroutine setup_discs(id,fileprefix,hfact,gamma,npart,polyk,&
rwarp = R_warp(i), &
warp_smoothl = H_warp(i), &
bh_spin = bhspin, &
enc_mass = enc_mass(:,i), &
prefix = dustprefix(j))

npart = npart + npindustdisc
Expand Down Expand Up @@ -1967,10 +1995,7 @@ subroutine setup_interactive(id)
end select

case (5:)
!print "(/,a)",'================================'
!print "(a)", '+++ HIERARCHICAL SYSTEM +++'
!print "(a)", '================================'


call print_chess_logo()!id)

ibinary = 0
Expand Down Expand Up @@ -2183,19 +2208,14 @@ subroutine setup_interactive(id)
call prompt('Enter H/R of circum-'//trim(disclabel)//' at R_ref',H_R(higher_disc_index))

higher_mass = get_hier_level_mass(trim(disclabel))!, mass, sink_num, sink_labels)
!print*, disclabel, higher_mass
!return
do i=1,maxdiscs
if (iuse_disc(i) .and. i /= higher_disc_index) then
call get_hier_disc_label(i, disclabel)
current_mass = get_hier_level_mass(trim(disclabel))
!print*, R_ref(i), R_ref(i), &
! higher_mass, current_mass, qindex(higher_disc_index),&
! H_R(higher_disc_index)
H_R(i) = (R_ref(i)/R_ref(higher_disc_index) * &
higher_mass/current_mass)**(0.5-qindex(higher_disc_index)) * &
H_R(higher_disc_index)
!print*,'!!!!!!!!!!!!!!!!!!!!!!! ', H_R(i)
endif
enddo
endif
Expand Down