Merge pull request #427 from crislong/dust_enclosed_mass

Compute total enclosed mass in self-gravitating disc initial conditions

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,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
 

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,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
@@ -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