Skip to content

Commit

Permalink
Edits to analysis_disc_stresses.f90 for variable gamma (eos = 21) and…
Browse files Browse the repository at this point in the history 
… more accurate H
  • Loading branch information
@alisonkyoung1
alisonkyoung1 committed Apr 15, 2024
1 parent 1d0327b commit c67364c
Show file tree
Hide file tree
Showing 3 changed files with 75 additions and 27 deletions.
6 changes: 4 additions & 2 deletions src/main/cooling_stamatellos.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,7 +205,8 @@ subroutine write_options_cooling_stamatellos(iunit)

!N.B. Tfloor handled in cooling.F90
call write_inopt(eos_file,'EOS_file','File containing tabulated EOS values',iunit)
call write_inopt(od_method,'OD method','Method for estimating optical depth: (1) Stamatellos (2) Lombardi (3) combined (4) modified Lombardi',iunit)
call write_inopt(od_method,'OD method','Method for estimating optical depth: (1)'&
'Stamatellos (2) Lombardi (3) combined (4) modified Lombardi',iunit)
call write_inopt(Lstar,'Lstar','Luminosity of host star for calculating Tmin (Lsun)',iunit)
call write_inopt(FLD_opt,'do FLD','Do FLD? (1) yes (0) no',iunit)

Expand All @@ -229,7 +230,8 @@ subroutine read_options_cooling_stamatellos(name,valstring,imatch,igotallstam,ie
case('OD method')
read(valstring,*,iostat=ierr) od_method
if (od_method < 1 .or. od_method > 4) then
call fatal('cooling options','od_method must be 1, 2, 3 or 4',var='od_method',ival=od_method)
call fatal('cooling options','od_method must be 1, 2, 3 or 4', &
var='od_method',ival=od_method)
endif
ngot = ngot + 1
case('EOS_file')
Expand Down
4 changes: 2 additions & 2 deletions src/main/readwrite_dumps_fortran.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1286,9 +1286,9 @@ subroutine read_phantom_arrays(i1,i2,noffset,narraylengths,nums,npartread,nparto
if (eos_outputs_gasP(ieos) .or. eos_is_non_ideal(ieos)) then
call read_array(eos_vars(igasP,:),eos_vars_label(igasP),got_gasP,ik,i1,i2,noffset,idisk1,tag,match,ierr)
endif
if (eos_is_non_ideal(ieos)) then
! if (eos_is_non_ideal(ieos)) then
call read_array(eos_vars(itemp,:),eos_vars_label(itemp),got_temp,ik,i1,i2,noffset,idisk1,tag,match,ierr)
endif
! endif
call read_array(eos_vars(iX,:),eos_vars_label(iX),got_x,ik,i1,i2,noffset,idisk1,tag,match,ierr)
call read_array(eos_vars(iZ,:),eos_vars_label(iZ),got_z,ik,i1,i2,noffset,idisk1,tag,match,ierr)
call read_array(eos_vars(imu,:),eos_vars_label(imu),got_mu,ik,i1,i2,noffset,idisk1,tag,match,ierr)
Expand Down
92 changes: 69 additions & 23 deletions src/utils/analysis_disc_stresses.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,8 @@
!--------------------------------------------------------------------------!
module analysis
!
! Analysis routine for discs by DF, adapted from a routine by CJN
! Analysis routine for discs by DF, adapted from a routine by CJN.
! Edited for use with variable gammai and mui and more accurate alpha_ss calc by AKY
!
! :References: None
!
Expand All @@ -28,7 +29,7 @@ module analysis
real :: rin, rout,dr
integer, allocatable,dimension(:) :: ipartbin
real, allocatable,dimension(:) :: rad,ninbin,sigma,csbin,vrbin,vphibin, omega
real, allocatable,dimension(:) :: H, toomre_q,epicyc
real, allocatable,dimension(:) :: H, toomre_q,epicyc,part_scaleheight
real, allocatable,dimension(:) :: alpha_reyn,alpha_grav,alpha_mag,alpha_art
real, allocatable,dimension(:) :: rpart,phipart,vrpart,vphipart, gr,gphi,Br,Bphi
real, allocatable,dimension(:,:) :: gravxyz
Expand All @@ -42,7 +43,7 @@ module analysis

subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit)
use io, only:fatal
use part, only:gravity,mhd
use part, only:gravity,mhd,eos_vars

character(len=*), intent(in) :: dumpfile
real, intent(in) :: xyzh(:,:),vxyzu(:,:)
Expand Down Expand Up @@ -82,7 +83,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit)
call transform_to_cylindrical(npart,xyzh,vxyzu)

! Bin particles by radius
call radial_binning(npart,xyzh,vxyzu,pmass)
call radial_binning(npart,xyzh,vxyzu,pmass,eos_vars)

! Calculate stresses
call calc_stresses(npart,xyzh,vxyzu,pmass)
Expand Down Expand Up @@ -350,18 +351,21 @@ end subroutine transform_to_cylindrical
!+
!---------------------------------------------------------------

subroutine radial_binning(npart,xyzh,vxyzu,pmass)
use physcon, only:pi
use eos, only: gamma
subroutine radial_binning(npart,xyzh,vxyzu,pmass,eos_vars)
use physcon, only:pi,kb_on_mh
use eos, only:gamma,ieos
use part, only:itemp,imu
use units, only:udist

implicit none

integer,intent(in) :: npart
real,intent(in) :: pmass
real,intent(in) :: xyzh(:,:),vxyzu(:,:)
real,intent(in) :: xyzh(:,:),vxyzu(:,:),eos_vars(:,:)

integer :: ibin,ipart,nbinned
integer :: ibin,ipart,nbinned,iallocerr
real :: area
real,allocatable :: zsetgas(:,:)

print '(a,I4)', 'Carrying out radial binning, number of bins: ',nbins

Expand All @@ -373,6 +377,7 @@ subroutine radial_binning(npart,xyzh,vxyzu,pmass)
allocate(omega(nbins))
allocate(vrbin(nbins))
allocate(vphibin(nbins))
allocate(part_scaleheight(nbins))

ipartbin(:) = 0
ninbin(:) = 0.0
Expand All @@ -381,6 +386,15 @@ subroutine radial_binning(npart,xyzh,vxyzu,pmass)
omega(:) = 0.0
vrbin(:) = 0.0
vphibin(:) = 0.0
part_scaleheight(:) = 0.0

allocate(zsetgas(npart,nbins),stat=iallocerr)
! If you don't have enough memory to allocate zsetgas, then calculate H the slow way with less memory.
if (iallocerr/=0) then
write(*,'(/,a)') ' WARNING: Could not allocate memory for array zsetgas!'
write(*,'(a)') ' (It possibly requires too much memory)'
write(*,'(a,/)') ' Try calculate scaleheight the slow way.'
endif

! Set up radial bins

Expand Down Expand Up @@ -409,20 +423,28 @@ subroutine radial_binning(npart,xyzh,vxyzu,pmass)

ninbin(ibin) = ninbin(ibin) +1
ipartbin(ipart) = ibin

csbin(ibin) = csbin(ibin) + sqrt(gamma*(gamma-1)*vxyzu(4,ipart))
if (ieos==21) then
csbin(ibin) = csbin(ibin) + sqrt(kb_on_mh * eos_vars(itemp,ipart) / eos_vars(imu,ipart))
if (csbin(ibin) == 0) then
print *, eos_vars(itemp,ipart)
endif
else
csbin(ibin) = csbin(ibin) + sqrt(gamma*(gamma-1)*vxyzu(4,ipart))
endif

area = pi*((rad(ibin)+0.5*dr)**2-(rad(ibin)- 0.5*dr)**2)
sigma(ibin) = sigma(ibin) + pmass/area

vrbin(ibin) = vrbin(ibin) + vrpart(ipart)
vphibin(ibin) = vphibin(ibin) + vphipart(ipart)
omega(ibin) = omega(ibin) + vphipart(ipart)/rad(ibin)

zsetgas(int(ninbin(ibin)),ibin) = xyzh(3,ipart)
endif

enddo

call calculate_H(nbins,part_scaleheight,zsetgas,int(ninbin))
part_scaleheight(:) = part_scaleheight(:)
print*, nbinned, ' particles have been binned'

where(ninbin(:)/=0)
Expand All @@ -443,11 +465,11 @@ end subroutine radial_binning
!+
!--------------------------------------------------------------
subroutine calc_stresses(npart,xyzh,vxyzu,pmass)
use physcon, only: pi,gg
use physcon, only: pi,gg,kb_on_mh
use units, only: print_units, umass,udist,utime,unit_velocity,unit_density,unit_Bfield
use dim, only: gravity
use part, only: mhd,rhoh,alphaind
use eos, only: gamma
use part, only: mhd,rhoh,alphaind,eos_vars,imu,itemp
use eos, only: gamma,ieos

implicit none

Expand Down Expand Up @@ -481,7 +503,9 @@ subroutine calc_stresses(npart,xyzh,vxyzu,pmass)
call print_units

sigma(:) = sigma(:)*umass/(udist*udist)
csbin(:) = csbin(:)*unit_velocity
if (ieos /= 21) then
csbin(:) = csbin(:)*unit_velocity
endif
omega(:) = omega(:)/utime

Keplog = 1.5
Expand All @@ -500,16 +524,16 @@ subroutine calc_stresses(npart,xyzh,vxyzu,pmass)
do ipart=1,npart
ibin = ipartbin(ipart)

if (ibin<=0) cycle

cs2 = gamma*(gamma-1)*vxyzu(4,ipart)*unit_velocity*unit_velocity
if (ibin<=0) cycle


dvr = (vrpart(ipart) - vrbin(ibin))*unit_velocity
dvphi = (vphipart(ipart) -vphibin(ibin))*unit_velocity
rhopart = rhoh(xyzh(4,ipart),pmass)*unit_density

alpha_reyn(ibin) = alpha_reyn(ibin) + dvr*dvphi

! Handle constant alpha_sph
alpha_art(ibin) = alpha_art(ibin) + alphaind(1,ipart)*xyzh(4,ipart)*udist

if (gravity) alpha_grav(ibin) = alpha_grav(ibin) + gr(ipart)*gphi(ipart)/rhopart
Expand Down Expand Up @@ -587,7 +611,7 @@ subroutine write_radial_data(iunit,output,time)
print '(a,a)', 'Writing to file ',output
open(iunit,file=output)
write(iunit,'("# Disc Stress data at t = ",es20.12)') time
write(iunit,"('#',11(1x,'[',i2.2,1x,a11,']',2x))") &
write(iunit,"('#',12(1x,'[',i2.2,1x,a11,']',2x))") &
1,'radius (AU)', &
2,'sigma (cgs)', &
3,'cs (cgs)', &
Expand All @@ -598,12 +622,13 @@ subroutine write_radial_data(iunit,output,time)
8,'alpha_reyn',&
9,'alpha_grav',&
10,'alpha_mag',&
11,'alpha_art'
11,'alpha_art',&
12,'particle H (au)'

do ibin=1,nbins
write(iunit,'(11(es18.10,1X))') rad(ibin),sigma(ibin),csbin(ibin), &
write(iunit,'(12(es18.10,1X))') rad(ibin),sigma(ibin),csbin(ibin), &
omega(ibin),epicyc(ibin),H(ibin), abs(toomre_q(ibin)),alpha_reyn(ibin), &
alpha_grav(ibin),alpha_mag(ibin),alpha_art(ibin)
alpha_grav(ibin),alpha_mag(ibin),alpha_art(ibin),part_scaleheight(ibin)
enddo

close(iunit)
Expand All @@ -612,6 +637,26 @@ subroutine write_radial_data(iunit,output,time)

end subroutine write_radial_data

subroutine calculate_H(nbin,H,zsetgas,ninbin)
! copied from utils disc
integer, intent(in) :: nbin
real, intent(out) :: H(:)
real, intent(in) :: zsetgas(:,:)
integer, intent(in) :: ninbin(:)
integer :: ii
real :: meanzii

do ii = 1,nbin
if (ninbin(ii)==0) then
meanzii = 0.
else
meanzii = sum(zsetgas(1:ninbin(ii),ii))/real(ninbin(ii))
endif
H(ii) = sqrt(sum(((zsetgas(1:ninbin(ii),ii)-meanzii)**2)/(real(ninbin(ii)-1))))
enddo

end subroutine calculate_H

!--------------------------------------------------------
!+
! Deallocate arrays
Expand All @@ -628,6 +673,7 @@ subroutine deallocate_arrays
deallocate(gr,gphi,Br,Bphi,vrbin,vphibin)
deallocate(sigma,csbin,H,toomre_q,omega,epicyc)
deallocate(alpha_reyn,alpha_grav,alpha_mag,alpha_art)
deallocate(part_scaleheight)

end subroutine deallocate_arrays
!-------------------------------------------------------
Expand Down

0 comments on commit c67364c

Please sign in to comment.