Update analysis_disc_stresses

build/Makefile_setups

@@ -923,10 +923,11 @@ ifeq ($(SETUP), isosgdisc)
     SETUPFILE= setup_disc.f90
     GRAVITY=yes
     IND_TIMESTEPS=yes
-    ANALYSIS=analysis_dustydisc.f90
+    #ANALYSIS=analysis_dustydisc.f90
     ISOTHERMAL=yes
     KNOWN_SETUP=yes
     SRCINJECT= inject_keplerian.f90
+    ANALYSIS = utils_getneighbours.F90 utils_omp.F90 analysis_disc_stresses.f90
 endif
 
 ifeq ($(SETUP), dustyisosgdisc)

src/utils/analysis_disc_stresses.f90

@@ -103,43 +103,47 @@ end subroutine do_analysis
 !+
 !-------------------------------------------
 subroutine read_analysis_options
-
- use prompting, only:prompt
-
- implicit none
-
+ use prompting,    only:prompt
+ use infile_utils, only:open_db_from_file,inopts,read_inopt,write_inopt,close_db
+ use io,           only:fatal
+ type(inopts), allocatable :: db(:)
+ integer :: ierr,nerr
  logical :: inputexist
  character(len=21) :: inputfile
+ integer, parameter :: iunit = 10
 
 ! Check for existence of input file
- inputfile = 'disc_stresses.options'
+ inputfile = trim(analysistype)//'.options'
  inquire(file=inputfile, exist=inputexist)
 
  if (inputexist) then
-
+    nerr = 0
     print '(a,a,a)', "Parameter file ",inputfile, " found: reading analysis options"
-
-    open(10,file=inputfile,form='formatted')
-    read(10,*) nbins
-    read(10,*) rin
-    read(10,*) rout
-    close(10)
+    call open_db_from_file(db,inputfile,iunit,ierr)
+    call read_inopt(nbins,'nbins',db,errcount=nerr)
+    call read_inopt(rin,'rin',db,errcount=nerr)
+    call read_inopt(rout,'rout',db,errcount=nerr)
+    call close_db(db)
+    if (nerr > 0) then
+      call fatal(trim(analysistype),'Error in reading '//trim(inputfile))
+    endif
 
  else
 
     print '(a,a,a)', "Parameter file ",inputfile, " NOT found"
-
+    nbins = 128; rin = 1.; rout = 100.
     call prompt('Enter the number of radial bins: ', nbins)
     call prompt('Enter the disc inner radius: ', rin)
     call prompt('Enter the disc outer radius: ', rout)
 
 ! Write choices to new inputfile
 
-    open(10,file=inputfile,status='new',form='formatted')
-    write(10,*) nbins, "      Number of radial bins"
-    write(10,*) rin,  "      Inner Disc Radius"
-    write(10,*) rout, "      Outer Disc Radius"
-    close(10)
+    open(unit=iunit,file=inputfile,status='new',form='formatted')
+    write(iunit,"(a)") '# parameter options for analysis of '//trim(analysistype)
+    call write_inopt(nbins,'nbins','Number of radial bins',iunit)
+    call write_inopt(rin,'rin','Inner Disc Radius',iunit)
+    call write_inopt(rout,'rout','Outer Disc Radius',iunit)
+    close(iunit)
  endif
 
 
@@ -162,8 +166,6 @@ subroutine calc_gravitational_forces(dumpfile,npart,xyzh,vxyzu)
  use part, only:poten,igas,iphase,maxphase,rhoh,massoftype,iamgas
  use kernel, only: get_kernel,get_kernel_grav1,cnormk
 
- implicit none
-
  character(len=*),intent(in) :: dumpfile
  real,intent(in) :: xyzh(:,:),vxyzu(:,:)
  integer,intent(in) :: npart
@@ -352,16 +354,15 @@ end subroutine transform_to_cylindrical
 
 subroutine radial_binning(npart,xyzh,vxyzu,pmass)
  use physcon, only:pi
- use eos, only: gamma
-
- implicit none
+ use eos,     only:get_spsound,ieos
+ use part,    only:rhoh,isdead_or_accreted
 
  integer,intent(in) :: npart
  real,intent(in) :: pmass
  real,intent(in) :: xyzh(:,:),vxyzu(:,:)
 
  integer :: ibin,ipart,nbinned
- real :: area
+ real :: area,csi
 
  print '(a,I4)', 'Carrying out radial binning, number of bins: ',nbins
 
@@ -395,7 +396,7 @@ subroutine radial_binning(npart,xyzh,vxyzu,pmass)
  do ipart=1,npart
 
     ! i refers to particle, ii refers to bin
-    if (xyzh(4,ipart)  >  tiny(xyzh)) then ! IF ACTIVE
+    if (.not.isdead_or_accreted(xyzh(4,ipart))) then ! IF ACTIVE
 
        ibin = int((rpart(ipart)-rad(1))/dr + 1)
 
@@ -410,7 +411,8 @@ 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))
+       csi = get_spsound(ieos,xyzh(1:3,ipart),rhoh(xyzh(4,ipart),pmass),vxyzu(:,ipart))
+       csbin(ibin) = csbin(ibin) + csi
 
        area = pi*((rad(ibin)+0.5*dr)**2-(rad(ibin)- 0.5*dr)**2)
        sigma(ibin) = sigma(ibin) + pmass/area