diff --git a/src/setup/set_disc.F90 b/src/setup/set_disc.F90 index a8778cf66..59296fe1e 100644 --- a/src/setup/set_disc.F90 +++ b/src/setup/set_disc.F90 @@ -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 @@ -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 @@ -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(:,:) @@ -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') @@ -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) @@ -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) @@ -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 @@ -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 @@ -1161,15 +1184,14 @@ 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) @@ -1177,7 +1199,6 @@ subroutine get_disc_mass(disc_m,enc_m,rad,toomre_min,sigmaprofile,sigma_norm, & toomre_min = min(toomre_min,real(cs*kappa/(pi*G*sigma))) endif enddo - enc_m = enc_m + star_m end subroutine get_disc_mass diff --git a/src/setup/setup_disc.f90 b/src/setup/setup_disc.f90 index 2232ea75f..368208e45 100644 --- a/src/setup/setup_disc.f90 +++ b/src/setup/setup_disc.f90 @@ -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 @@ -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 @@ -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) @@ -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 @@ -1053,7 +1073,8 @@ 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), & @@ -1061,8 +1082,10 @@ subroutine calculate_disc_mass() 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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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