From 77654a107baf93b24998a7911868af38240a86a6 Mon Sep 17 00:00:00 2001 From: sdarbha Date: Sat, 20 Jun 2020 21:28:48 -0700 Subject: [PATCH 01/15] Initial commit. Updated Makefile to compile on NERSC Cori. Created setup and analysis files for BH-NS tidal ejecta --- build/Makefile | 53 ++++- src/setup/setup_grbhnsejecta.f90 | 338 ++++++++++++++++++++++++++++++ src/utils/analysis_bhnsejecta.f90 | 170 +++++++++++++++ 3 files changed, 559 insertions(+), 2 deletions(-) create mode 100644 src/setup/setup_grbhnsejecta.f90 create mode 100644 src/utils/analysis_bhnsejecta.f90 diff --git a/build/Makefile b/build/Makefile index eca738a7b..9153876cd 100644 --- a/build/Makefile +++ b/build/Makefile @@ -271,6 +271,23 @@ ifeq ($(SETUP), grtde) ANALYSIS=analysis_tde.f90 endif +# I (SD) added this setup for the ejecta from a BH-NS merger +ifeq ($(SETUP), grbhnsejecta) + SETUPFILE=setup_grbhnsejecta.f90 + ANALYSIS=analysis_bhnsejecta.f90 +# SRCINJECT=inject_bhnsejecta_particles.f90 + GR=yes + METRIC=schwarzschild + GRAVITY=no +# GRAVITY=yes + ISOTHERMAL=no + IND_TIMESTEPS=no + KNOWN_SETUP=yes + KERNEL=cubic +# KERNEL=quintic +# KERNEL=WendlandC4 +endif + ifeq ($(SETUP), srpolytrope) FPPFLAGS= -DPRIM2CONS_FIRST SETUPFILE= setup_srpolytrope.f90 @@ -538,6 +555,7 @@ ifeq ($(SETUP), srshock) CONST_AV=yes endif + ifeq ($(SETUP), testparticles) # test particles PERIODIC=no @@ -1217,6 +1235,33 @@ ifeq ($(SYSTEM), g2) MPIEXEC='mpiexec -npernode 1' endif +# I (SD) added this system to be able to compile Phantom on NERSC Cori +# On Cori, you load programming environments, which use wrappers for the compilers, you do not use the vendor specific compiler names directly +# I (SD) just copied the parameters from Makefile_defaults_ifort and ozstar below and made several replacements +# Most importantly, I (SD) made the replacements ifort -> ftn and icc -> cc, and added CRAYPE_LINK_TYPE to tell Cori how to link +ifeq ($(SYSTEM), cori) +#Cori (NERSC machine) + FC= ftn +# FFLAGS= -O3 -inline-factor=500 -mcmodel=medium -shared-intel -warn uninitialized -warn unused -warn truncated_source + FFLAGS= -O3 -inline-factor=500 -shared-intel -warn uninitialized -warn unused -warn truncated_source + DBLFLAG= -r8 + DEBUGFLAG= -check all -WB -traceback -g -debug -all + ENDIANFLAGBIG= -convert big_endian + ENDIANFLAGLITTLE= -convert little_endian + CC= cc +# CCFLAGS= -O3 -mcmodel=medium + CCFLAGS= -O3 + LIBCXX = -cxxlib + QSYS= slurm + OMPFLAGS= -qopenmp + NOMP=64 + KNOWN_SYSTEM=yes +# CRAYPE_LINK_TYPE=dynamic + CRAYPE_LINK_TYPE=static + MAXP=2000000 +# MAXP=8000000 +endif + ifeq ($(SYSTEM), ozstar) # ozstar facility include Makefile_defaults_ifort @@ -1227,6 +1272,8 @@ ifeq ($(SYSTEM), ozstar) WALLTIME='168:00:00' endif + + ifeq ($(SYSTEM), monarch) include Makefile_defaults_ifort OMPFLAGS=-qopenmp -qopt-report @@ -1931,7 +1978,7 @@ endif OBJMOD1 = utils_omp.F90 utils_summary.f90 utils_indtimesteps.F90 \ utils_sort.f90 checksetup.f90 set_Bfield.f90 \ partinject.F90 random.f90 set_disc.F90 set_dust.F90 set_binary.f90 \ - icosahedron.f90 utils_inject.f90 ${SRCINJECT} \ + icosahedron.f90 utils_inject.f90 ${SRCINJECT} \ ${SRCTURB} ${SRCNIMHD} ${SRCCHEM} \ density_profiles.f90 ptmass.F90 damping.f90 readwrite_infile.f90 ${MODFILE:.f90=.o} OBJMOD2 = ${OBJMOD1:.F90=.o} @@ -1985,6 +2032,8 @@ OBJLIB1=${SRCLIB:.f90=.o} OBJLIB=${OBJLIB1:.F90=.o} libphantom: checksystem checkparams phantom ${OBJLIB} ar rc $(BINDIR)/libphantom.a ${OBJLIB} ${OBJECTS} $(LDFLAGS) +# $(LIBTOOL) -static ${OBJLIB} ${OBJECTS} $(LDFLAGS) -o $(BINDIR)/libphantom.a +# $(LIBTOOL) --tag=FC --mode=link ${OBJLIB} ${OBJECTS} $(LDFLAGS) -o $(BINDIR)/libphantom.a $(FC) -shared $(FFLAGS) $(FPPFLAGS) $(DBLFLAG) ${OBJLIB} ${OBJECTS} $(LDFLAGS) -o $(BINDIR)/libphantom.so cleanlibphantom: @@ -2852,7 +2901,7 @@ checkparams: ifneq ($(FPPFLAGS),$(LASTFPPFLAGS)) @echo 'pre-processor flags changed from "'${LASTFPPFLAGS}'" to "'${FPPFLAGS}'"' @${MAKE} clean; - #for x in ../src/*/*.F90; do y=`basename $$x`; rm -f $${y/.F90/.o}; done +# for x in ../src/*/*.F90; do y=`basename $$x`; rm -f $${y/.F90/.o}; done endif @echo "Preprocessor flags are "${FPPFLAGS} @echo "${FPPFLAGS}" > .make_lastfppflags diff --git a/src/setup/setup_grbhnsejecta.f90 b/src/setup/setup_grbhnsejecta.f90 new file mode 100644 index 000000000..19a8a0927 --- /dev/null +++ b/src/setup/setup_grbhnsejecta.f90 @@ -0,0 +1,338 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2020 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +!+ +! MODULE: setup +! +! DESCRIPTION: Read particle file from SpEC BH-NS merger simulation and set up SPH particles +! +! REFERENCES: None +! +! OWNER: Siva Darbha +! +! RUNTIME PARAMETERS: +! mhole -- mass of black hole (solar mass) +! particle_file_name -- name of SpEC particle file +! +! DEPENDENCIES: eos, externalforces, infile_utils, io, metric, part, +! physcon, rho_profile, timestep, units, vectorutils +!+ +!-------------------------------------------------------------------------- +module setup + implicit none + public :: setpart + + real :: mhole + character(len=120) :: particle_file_name + + private + +contains + +!---------------------------------------------------------------- +!+ +! setup for sink particle binary simulation (no gas) +!+ +!---------------------------------------------------------------- +subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,time,fileprefix) + use part, only:nptmass,ihacc,ihsoft,igas,set_particle_type,rhoh,gravity + use units, only:set_units,umass,udist,unit_velocity,unit_energ,unit_density,unit_ergg + use physcon, only:solarm,pi,solarr,gg,c,eV,radconst,kboltz + use io, only:master,fatal,warning + use timestep, only:tmax,dtmax + use metric, only:mass1 + use eos, only:ieos + use externalforces,only:accradius1,accradius1_hard + use allocutils, only:allocate_array + integer, intent(in) :: id + integer, intent(inout) :: npart + integer, intent(out) :: npartoftype(:) + !----- xyzh(:,i) is the array (x,y,z,h) for the ith SPH particle, where (x,y,z) are the coordinates and h is the smoothing length. + !----- In Price et al PASA (2018) and in Liptai & Price (2019), h is given by the symbols h and h_a. + real, intent(out) :: xyzh(:,:) + real, intent(out) :: massoftype(:) + real, intent(out) :: polyk,gamma,hfact + real, intent(inout) :: time + character(len=20), intent(in) :: fileprefix + !----- vxyzu(:,i) is the array (vx,vy,vz,u) for the ith SPH particle, where (vx,vy,vz) are the coordinate velocities and u is the specific internal energy in the rest frame. + !----- In Price et al PASA (2018) and in Liptai & Price (2019), u is given the symbol u. + real, intent(out) :: vxyzu(:,:) + character(len=120) :: filename + integer, parameter :: ntab=5000 + integer :: ierr,i + logical :: iexist + real :: x,y,z + real, allocatable :: m_array(:), r_array(:), theta_array(:), phi_array(:), e_array(:), l_array(:), u_r_array(:), u_theta_array(:), y_e_array(:), s_array(:), T_array(:), rho_array(:) + real :: m_ejecta + real :: r, theta, phi, e, l, u_r, u_theta, dt_dtau, dr_dtau, dtheta_dtau, dphi_dtau, dr_dt, dtheta_dt, dphi_dt + real :: vx, vy, vz + real :: m, rho, delta + real :: T, y_e + +! +!-- general parameters +! + time = 0. + ieos = 2 !----- flag for EOS type. The value ieos=2 chooses an adiabatic EOS + gamma = 4./3. !----- gamma in adiabatic gamma-law EOS + polyk = 0. +! if (.not.gravity) call fatal('setup','recompile with GRAVITY=yes') + +! +!-- space available for injected gas particles +! + npart = 0 + npartoftype(:) = 0 + xyzh(:,:) = 0. + vxyzu(:,:) = 0. + nptmass = 0 + +! +!-- Default runtime parameters +! +! + mhole = 0. ! (solar masses) + particle_file_name = 'particle_file_name.dat' + +! +!-- Read runtime parameters from setup file +! + if (id==master) print "(/,65('-'),1(/,a),/,65('-'),/)",' BH-NS merger ejecta in GR' + filename = trim(fileprefix)//'.setup' + inquire(file=filename,exist=iexist) + if (iexist) call read_setupfile(filename,ierr) + if (.not. iexist .or. ierr /= 0) then + if (id==master) then + call write_setupfile(filename) + print*,' Edit '//trim(filename)//' and rerun phantomsetup' + endif + stop + endif + +! +!-- Convert to code untis +! + mhole = mhole*solarm + call set_units(mass=mhole,c=1.,G=1.) !--Set central mass to M=1 in code units + + accradius1_hard = 2.02*mass1 + accradius1 = accradius1_hard + + call get_num_particles(npart) + + allocate(m_array(npart)) + allocate(r_array(npart)) + allocate(theta_array(npart)) + allocate(phi_array(npart)) + allocate(e_array(npart)) + allocate(l_array(npart)) + allocate(u_r_array(npart)) + allocate(u_theta_array(npart)) + allocate(y_e_array(npart)) + allocate(s_array(npart)) + allocate(T_array(npart)) + allocate(rho_array(npart)) + + call read_particle_data(npart, m_array, r_array, theta_array, phi_array, e_array, l_array, u_r_array, u_theta_array, y_e_array, s_array, T_array, rho_array) + + m_ejecta = 0 + do i=1,npart + !----- Convert to code units, i.e. first convert to cgs then normalize to code units + !----- Note: the particle file is in units of G = c = M_Sun = 1 (for most of the variables), whereas in set_units above we set the code units to be G = c = M_bh = 1 + !----- Note: the particle file gives temperature in units of MeV, which we have to convert here + m_array(i) = m_array(i) * solarm / umass + r_array(i) = r_array(i) * (gg*solarm/(c**2.)) / udist + e_array(i) = e_array(i) * (c**2.) / (unit_velocity**2.) + l_array(i) = l_array(i) * (gg*solarm/c) / (udist*unit_velocity) + u_r_array(i) = u_r_array(i) * c / unit_velocity + u_theta_array(i) = u_theta_array(i) * (gg*solarm/c) / (udist*unit_velocity) + ! y_e_array(i) = y_e_array(i) !----- Y_e is dimensionless, so we don't need to do any unit conversion + ! s_array(i) = s_array(i)* !----- I'm not sure what the units are for the entropy in the particle file, I have to look into this + T_array(i) = T_array(i) * (1.e6)*eV / unit_energ + rho_array(i) = rho_array(i) * ((c**6.)/((gg**3.)*(solarm**2.))) / unit_density + + m_ejecta = m_ejecta + m_array(i) + enddo + + npartoftype(igas) = npart + !----------- The SPH gas particles must all have the same mass, since massoftype(igas) can only take one value ------------- + massoftype(igas) = m_array(1) !----- set the mass of the SPH gas particles + do i=1,npart + call set_particle_type(i,igas) + r = r_array(i) + theta = theta_array(i) + phi = phi_array(i) + x = r*sin(theta)*cos(phi) + y = r*sin(theta)*sin(phi) + z = r*cos(theta) + e = e_array(i) + u_r = u_r_array(i) + u_theta = u_theta_array(i) + l = l_array(i) + dt_dtau = e / (1. - (2.*mass1/r)) + dr_dtau = (1. - (2.*mass1/r)) * u_r + dtheta_dtau = u_theta / (r**2.) + dphi_dtau = l / ((r**2.)*(sin(theta)**2.)) + dr_dt = dr_dtau / dt_dtau + dtheta_dt = dtheta_dtau / dt_dtau + dphi_dt = dphi_dtau / dt_dtau + vx = sin(theta)*cos(phi)*dr_dt + r*cos(theta)*cos(phi)*dtheta_dt - r*sin(theta)*sin(phi)*dphi_dt + vy = sin(theta)*sin(phi)*dr_dt + r*cos(theta)*sin(phi)*dtheta_dt + r*sin(theta)*cos(phi)*dphi_dt + vz = cos(theta)*dr_dt - r*sin(theta)*dtheta_dt + !----- xyzh(:,i) is the array (x,y,z,h) for the ith SPH particle, where (x,y,z) are the coordinates and h is the smoothing length. + !----- In Price et al PASA (2018) and in Liptai & Price (2019), h is given by the symbols h and h_a. + m = m_array(i) + rho = rho_array(i) + delta = (m/rho)**(1./3.) !----- local interparticle spacing + xyzh(4,i) = hfact*delta !----- smoothing length + xyzh(1:3,i) = (/x,y,z/) + !----- vxyzu(:,i) is the array (vx,vy,vz,u) for the ith SPH particle, where (vx,vy,vz) are the coordinate velocities and u is the specific internal energy in the rest frame. + !----- In Price et al PASA (2018) and in Liptai & Price (2019), u is given the symbol u. + T = T_array(i) + y_e = y_e_array(i) + vxyzu(4,i) = ( radconst * ((T*unit_energ/kboltz)**4.) / (rho*unit_density) ) / unit_ergg !----- the specific internal energy in the rest frame, u = a * T^4 / rho + vxyzu(1:3,i) = (/vx,vy,vz/) + enddo + + deallocate(m_array) + deallocate(r_array) + deallocate(theta_array) + deallocate(phi_array) + deallocate(e_array) + deallocate(l_array) + deallocate(u_r_array) + deallocate(u_theta_array) + deallocate(y_e_array) + deallocate(s_array) + deallocate(T_array) + deallocate(rho_array) + + if (id==master) then + print "(/,a)", ' EJECTA SETUP:' + print "(a,1f10.3)" ,' EOS gamma = ', gamma + endif + + if (id==master) print "(/,a,i10,/)",' Number of particles setup = ',npart + + if (npart == 0) call fatal('setup','no particles setup') + if (ierr /= 0) call fatal('setup','ERROR during setup') + +end subroutine setpart + +! +!---Read/write setup file-------------------------------------------------- +! +subroutine write_setupfile(filename) + use infile_utils, only:write_inopt + character(len=*), intent(in) :: filename + integer, parameter :: iunit = 20 + + print "(a)",' writing setup options file '//trim(filename) + open(unit=iunit,file=filename,status='replace',form='formatted') + write(iunit,"(a)") '# input file for binary setup routines' + call write_inopt(mhole, 'mhole', 'mass of black hole (solar mass)', iunit) + call write_inopt(particle_file_name,'particle_file_name','name of SpEC particle file', iunit) + close(iunit) + +end subroutine write_setupfile + +subroutine read_setupfile(filename,ierr) + use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db + use io, only:error + character(len=*), intent(in) :: filename + integer, intent(out) :: ierr + integer, parameter :: iunit = 21 + integer :: nerr + type(inopts), allocatable :: db(:) + + print "(a)",'reading setup options from '//trim(filename) + nerr = 0 + ierr = 0 + call open_db_from_file(db,filename,iunit,ierr) + call read_inopt(mhole, 'mhole', db,min=0.,errcount=nerr) + call read_inopt(particle_file_name,'particle_file_name',db,errcount=nerr) + call close_db(db) + if (nerr > 0) then + print "(1x,i2,a)",nerr,' error(s) during read of setup file: re-writing...' + ierr = nerr + endif + +end subroutine read_setupfile + +! +!---Get number of particles in particle data-------------------------------------------------- +! +subroutine get_num_particles(npart) + integer :: iostatus, i + integer, intent(out) :: npart + real :: m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, rho + + print *, 'getting number of particles from '//trim(particle_file_name) + + open(1, file=particle_file_name, status='old', action='read') + + do i=1,13 + read(1,*) + enddo + + i = 0 + do + read(1,*,IOSTAT=iostatus) m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, rho + ! read failed + if (iostatus > 0) then + print *, 'read failed, input data not in correct format' + exit + ! reached end of file + else if (iostatus < 0) then + exit + ! read successful + else + i = i + 1 + endif + enddo + npart = i + + close(1) + +end subroutine get_num_particles + +! +!---Read particle data-------------------------------------------------- +! +subroutine read_particle_data(npart, m_array, r_array, theta_array, phi_array, e_array, l_array, u_r_array, u_theta_array, y_e_array, s_array, T_array, rho_array) + integer :: iostatus, i + integer, intent(in) :: npart + real, intent(out) :: m_array(:), r_array(:), theta_array(:), phi_array(:), e_array(:), l_array(:), u_r_array(:), u_theta_array(:), y_e_array(:), s_array(:), T_array(:), rho_array(:) + + print *, 'reading particle data from '//trim(particle_file_name) + + open(1, file=particle_file_name, status='old', action='read') + + do i=1,13 + read(1,*) + enddo + + i = 1 + do + read(1,*,IOSTAT=iostatus) m_array(i), r_array(i), theta_array(i), phi_array(i), e_array(i), l_array(i), u_r_array(i), u_theta_array(i), y_e_array(i), s_array(i), T_array(i), rho_array(i) + ! read failed + if (iostatus > 0) then + print *, 'read failed, input data not in correct format' + exit + ! reached end of file + else if (iostatus < 0) then + exit + ! read successful + else + i = i + 1 + endif + enddo + + close(1) + +end subroutine read_particle_data + +end module setup \ No newline at end of file diff --git a/src/utils/analysis_bhnsejecta.f90 b/src/utils/analysis_bhnsejecta.f90 new file mode 100644 index 000000000..8d289adbc --- /dev/null +++ b/src/utils/analysis_bhnsejecta.f90 @@ -0,0 +1,170 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2020 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +!+ +! MODULE: analysis +! +! DESCRIPTION: +! Analysis routine for BH-NS merger simulations +! +! REFERENCES: None +! +! OWNER: Siva Darbha +! +! RUNTIME PARAMETERS: None +! +! DEPENDENCIES: None +!+ +!-------------------------------------------------------------------------- +module analysis + implicit none + character(len=20), parameter, public :: analysistype = 'bhnsejecta' + public :: do_analysis + + private + +contains + +subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) + use metric, only:mass1 + use units, only:umass,udist,utime,unit_velocity,unit_energ,unit_density,unit_ergg + use physcon, only:solarm,pi,solarr,gg,c,eV,radconst,kboltz + use part, only:hfact + character(len=*), intent(in) :: dumpfile + integer, intent(in) :: numfile,npart,iunit + real, intent(in) :: xyzh(:,:),vxyzu(:,:) + real, intent(in) :: pmass,time + integer :: iostatus, i + real :: x,y,z,h + real :: r, theta, phi, e, l, u_r, u_theta, dt_dtau, dr_dtau, dtheta_dtau, dphi_dtau, dr_dt, dtheta_dt, dphi_dt + real :: vx, vy, vz, uint + real :: m, rho, delta + real :: T, s, y_e + real :: q + real :: time_output + character(len=120) :: particle_file_name, info_file_name + + + + particle_file_name = trim(dumpfile)//'.dat' + + open(1, file=particle_file_name, action='write', status='replace') + + write(1,'(A)',IOSTAT=iostatus) '# Particle data' + write(1,'(A)',IOSTAT=iostatus) '# [1] = Mass' + write(1,'(A)',IOSTAT=iostatus) '# [2] = Radius' + write(1,'(A)',IOSTAT=iostatus) '# [3] = Theta' + write(1,'(A)',IOSTAT=iostatus) '# [4] = Phi' + write(1,'(A)',IOSTAT=iostatus) '# [5] = e (-u_t)' + write(1,'(A)',IOSTAT=iostatus) '# [6] = l (u_phi)' + write(1,'(A)',IOSTAT=iostatus) '# [7] = u_r' + write(1,'(A)',IOSTAT=iostatus) '# [8] = u_theta' + write(1,'(A)',IOSTAT=iostatus) '# [9] = Y_e' + write(1,'(A)',IOSTAT=iostatus) '# [10]= s' + write(1,'(A)',IOSTAT=iostatus) '# [11]= T' + write(1,'(A)',IOSTAT=iostatus) '# [12]= rho' + + do i=1,npart + + m = pmass + + !----- xyzh(:,i) is the array (x,y,z,h) for the ith SPH particle, where (x,y,z) are the coordinates and h is the smoothing length. + !----- In Price et al PASA (2018) and in Liptai & Price (2019), h is given by the symbols h and h_a. + x = xyzh(1,i) + y = xyzh(2,i) + z = xyzh(3,i) + h = xyzh(4,i) !----- the smoothing length + + !----- vxyzu(:,i) is the array (vx,vy,vz,u) for the ith SPH particle, where (vx,vy,vz) are the coordinate velocities and u is the specific internal energy in the rest frame. + !----- In Price et al PASA (2018) and in Liptai & Price (2019), u is given the symbol u. + vx = vxyzu(1,i) + vy = vxyzu(2,i) + vz = vxyzu(3,i) + uint = vxyzu(4,i) !----- the specific internal energy in the rest frame + + r = (x**2. + y**2. + z**2.)**(1./2.) + q = (x**2. + y**2.)**(1./2.) + theta = atan2(q,z) + phi = atan2(y,x) + if (phi < 0) then + phi = phi + 2.*pi + endif + + dr_dt = (x/r)*vx + (y/r)*vy + (z/r)*vz + dtheta_dt = (z/(r**2.))*(x/q)*vx + (z/(r**2.))*(y/q)*vy - (q/(r**2.))*vz + dphi_dt = (-y/(q**2.))*vx + (x/(q**2.))*vy + + dt_dtau = ( -1. * ( -1.*(1. - (2.*mass1/r)) + ((1. - (2.*mass1/r))**(-1.)) * (dr_dt**2.) + (r**2.) * (dtheta_dt**2.) + ((r**2.)*(sin(theta)**2.)) * (dphi_dt**2.) ) )**(-1./2.) + + dr_dtau = dr_dt * dt_dtau + dtheta_dtau = dtheta_dt * dt_dtau + dphi_dtau = dphi_dt * dt_dtau + + e = dt_dtau * (1. - (2.*mass1/r)) + u_r = dr_dtau / (1. - (2.*mass1/r)) + u_theta = dtheta_dtau * (r**2.) + l = dphi_dtau * ((r**2.)*(sin(theta)**2.)) + + rho = m * (hfact/h)**3. + + T = ( (uint*unit_ergg) * (rho*unit_density) / radconst )**(1./4.) !----- temperature in K + + y_e = 0. !----- The SPH gas particles do not store the electron fraction Y_e since the code doesn't evolve the composition, so I just set it to zero here + + s = 0. !----- I'm not sure how to calculate the entropy of the SPH particles, or what the units are for the entropy in the particle file, so I just set it to zero here + + + + !----- Convert from code units to units of particle file, i.e. first convert to cgs then normalize to units of particle file + !----- Note: the particle file is in units of G = c = M_Sun = 1 (for most of the variables), whereas in set_units in setup_grbhnsejecta.f90 we set the code units to be G = c = M_bh = 1 + !----- Note: the particle file gives temperature in units of MeV, which we have to convert here + m = m * umass / solarm + r = r * udist / (gg*solarm/(c**2.)) + e = e * (unit_velocity**2.) / (c**2.) + l = l * (udist*unit_velocity) / (gg*solarm/c) + u_r = u_r * unit_velocity / c + u_theta = u_theta * (udist*unit_velocity) / (gg*solarm/c) + T = kboltz * T / ( (1.e6)*eV ) !----- temperature in MeV + rho = rho * unit_density / ( solarm / (gg*solarm/(c**2.))**3. ) + + + + write(1,'(ES11.4,1X)',advance='no',IOSTAT=iostatus) m + write(1,'(F10.3,1X)',advance='no',IOSTAT=iostatus) r + write(1,'(F9.5,1X)',advance='no',IOSTAT=iostatus) theta + write(1,'(F10.6,1X)',advance='no',IOSTAT=iostatus) phi + write(1,'(F10.5,1X)',advance='no',IOSTAT=iostatus) e + write(1,'(F10.4,1X)',advance='no',IOSTAT=iostatus) l + write(1,'(F12.6,1X)',advance='no',IOSTAT=iostatus) u_r + write(1,'(F10.4,1X)',advance='no',IOSTAT=iostatus) u_theta + write(1,'(F10.7,1X)',advance='no',IOSTAT=iostatus) y_e + write(1,'(F10.7,1X)',advance='no',IOSTAT=iostatus) s + write(1,'(F10.6,1X)',advance='no',IOSTAT=iostatus) T + write(1,'(ES12.5)',IOSTAT=iostatus) rho + + enddo + + close(1) + + + + !----- Convert time from code units to units of particle file, i.e. first convert to cgs then normalize to units of particle file + !----- Note: the particle file is in units of G = c = M_Sun = 1 (for most of the variables), whereas in set_units in setup_grbhnsejecta.f90 we set the code units to be G = c = M_bh = 1 + time_output = time * utime / (gg*solarm/(c**3.)) + + info_file_name = trim(dumpfile)//'.info' + + open(2, file=info_file_name, action='write', status='replace') + + write(2,'(A,1X,ES10.3,1X,A)',IOSTAT=iostatus) 'time =', time_output, '[G = c = M_Sun = 1]' + + close(2) + + + +end subroutine do_analysis + +end module \ No newline at end of file From feb527838dc5c3f315eabcc9a758fe91c742cfb0 Mon Sep 17 00:00:00 2001 From: sdarbha Date: Tue, 4 Aug 2020 06:59:30 -0700 Subject: [PATCH 02/15] Added rprocess heating function for BH-NS tidal ejecta --- build/Makefile | 2 +- src/main/checksetup.F90 | 3 +- src/main/cooling.f90 | 11 ++- src/main/force.F90 | 10 ++ src/main/rprocess_heating.f90 | 179 ++++++++++++++++++++++++++++++++++ src/main/step_leapfrog.F90 | 8 +- 6 files changed, 208 insertions(+), 5 deletions(-) create mode 100644 src/main/rprocess_heating.f90 diff --git a/build/Makefile b/build/Makefile index 9153876cd..52bcce106 100644 --- a/build/Makefile +++ b/build/Makefile @@ -1777,7 +1777,7 @@ endif SRCTEST = $(SRCTEST2:test_corotate.f90=) test_gr.f90 endif -SRCCHEM= coolfunc.f90 fs_data.f90 mol_data.f90 utils_spline.f90 h2cooling.f90 h2chem.f90 cooling.f90 +SRCCHEM= coolfunc.f90 rprocess_heating.f90 fs_data.f90 mol_data.f90 utils_spline.f90 h2cooling.f90 h2chem.f90 cooling.f90 SRCMESA= eos_mesa_microphysics.F90 eos_mesa.f90 SRCEOS = ${SRCMESA} eos_shen.f90 eos_helmholtz.f90 eos.F90 diff --git a/src/main/checksetup.F90 b/src/main/checksetup.F90 index 084170971..25bb4cabb 100644 --- a/src/main/checksetup.F90 +++ b/src/main/checksetup.F90 @@ -410,7 +410,8 @@ subroutine check_setup(nerror,nwarn,restart) if (id==master) & write(*,"(a,2(es10.3,', '),es10.3,a)") ' Centre of mass is at (x,y,z) = (',xcom,')' - if (.not.h2chemistry .and. maxvxyzu >= 4 .and. icooling >= 1 .and. iexternalforce/=iext_corotate) then +! if (.not.h2chemistry .and. maxvxyzu >= 4 .and. icooling >= 1 .and. iexternalforce/=iext_corotate) then + if (.not.h2chemistry .and. maxvxyzu >= 4 .and. icooling == 1 .and. iexternalforce/=iext_corotate) then if (dot_product(xcom,xcom) > 1.e-2) then print*,'Error in setup: Gammie (2001) cooling (icooling=1) assumes Omega = 1./r^1.5' print*,' but the centre of mass is not at the origin!' diff --git a/src/main/cooling.f90 b/src/main/cooling.f90 index e601567d2..00e491523 100644 --- a/src/main/cooling.f90 +++ b/src/main/cooling.f90 @@ -22,7 +22,7 @@ ! beta_cool -- beta factor in Gammie (2001) cooling ! icooling -- cooling function (0=off, 1=Gammie cooling 2=SD93 3=cooling function) ! -! DEPENDENCIES: coolfunc, dim, eos, h2cooling, infile_utils, io, options, +! DEPENDENCIES: coolfunc, rprocess_heating, dim, eos, h2cooling, infile_utils, io, options, ! part, physcon, timestep, units !+ !-------------------------------------------------------------------------- @@ -149,6 +149,7 @@ subroutine write_options_cooling(iunit) use infile_utils, only:write_inopt use h2cooling, only:write_options_h2cooling use coolfunc, only:write_options_coolfunc + use rprocess_heating, only:write_options_rprocess integer, intent(in) :: iunit write(iunit,"(/,a)") '# options controlling cooling' @@ -166,6 +167,8 @@ subroutine write_options_cooling(iunit) call write_inopt(beta_cool,'beta_cool','beta factor in Gammie (2001) cooling',iunit) elseif (icooling == 3) then call write_options_coolfunc(iunit) + elseif (icooling == 4) then + call write_options_rprocess(iunit) endif end subroutine write_options_cooling @@ -178,12 +181,13 @@ end subroutine write_options_cooling subroutine read_options_cooling(name,valstring,imatch,igotall,ierr) use h2cooling, only:read_options_h2cooling use coolfunc, only:read_options_coolfunc + use rprocess_heating, only:read_options_rprocess use io, only:fatal character(len=*), intent(in) :: name,valstring logical, intent(out) :: imatch,igotall integer, intent(out) :: ierr integer, save :: ngot = 0 - logical :: igotallh2,igotallcf + logical :: igotallh2,igotallcf,igotallrp imatch = .true. igotall = .false. ! cooling options are compulsory @@ -206,6 +210,7 @@ subroutine read_options_cooling(name,valstring,imatch,igotall,ierr) imatch = .false. if (h2chemistry) call read_options_h2cooling(name,valstring,imatch,igotallh2,ierr) if (icooling==3) call read_options_coolfunc(name,valstring,imatch,igotallcf,ierr) + if (icooling==4) call read_options_rprocess(name,valstring,imatch,igotallrp,ierr) end select if (igotallh2 .and. ngot >= 1) igotall = .true. @@ -214,6 +219,8 @@ subroutine read_options_cooling(name,valstring,imatch,igotall,ierr) if (.not.igotallcf) igotall = .false. + if (.not.igotallrp) igotall = .false. + if (.not.h2chemistry .and. (icooling == 1 .and. ngot < 3)) igotall= .false. end subroutine read_options_cooling diff --git a/src/main/force.F90 b/src/main/force.F90 index 4a72d968e..de1f6099e 100644 --- a/src/main/force.F90 +++ b/src/main/force.F90 @@ -2493,6 +2493,8 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv use nicil, only:nimhd_get_dudt,nimhd_get_dt use cooling, only:energ_cooling use chem, only:energ_h2cooling + use rprocess_heating, only:get_rprocess_heating_rate + use timestep, only:time use timestep, only:C_cour,C_cool,C_force,bignumber,dtmax use timestep_sts, only:use_sts use units, only:unit_ergg,unit_density @@ -2586,6 +2588,7 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv real :: posi(3),veli(3),gcov(0:3,0:3),metrici(0:3,0:3,2) integer :: ii,ia,ib,ic,ierror #endif + real :: q eni = 0. realviscosity = (irealvisc > 0) @@ -2774,6 +2777,13 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv if (gr) then fxyz4 = fxyz4 + (gamma - 1.)*densi**(1.-gamma)*u0i*fsum(idendtdissi) endif + + ! add r-process heating, if option selected + !if (gr) then + ! if (icooling==4) call get_rprocess_heating_rate(q,time) + ! fxyz4 = fxyz4 + (gamma - 1.) * densi**(1. - gamma) * q + !endif + #ifdef GR #ifdef ISENTROPIC fxyz4 = 0. diff --git a/src/main/rprocess_heating.f90 b/src/main/rprocess_heating.f90 new file mode 100644 index 000000000..daee2c1c9 --- /dev/null +++ b/src/main/rprocess_heating.f90 @@ -0,0 +1,179 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2020 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +!+ +! MODULE: rprocess_heating +! +! DESCRIPTION: +! Implements r-process heating. +! +! OWNER: Siva Darbha +! +! $Id$ +! +! RUNTIME PARAMETERS: +! +! q_0_cgs -- heating rate coefficient, in [ergs/s/g] +! t_b1_seconds -- time of break 1 in heating rate, in [s]' +! exp_1 -- exponent 1 in radioactive power component +! t_b2_seconds -- time of break 2 in heating rate, in [s]' +! exp_2 -- exponent 2 in radioactive power component +! t_start_seconds -- time after merger at which the simulation begins [s] +! +! DEPENDENCIES: infile_utils, io, physcon, units +!+ +!-------------------------------------------------------------------------- +module rprocess_heating + implicit none + integer, parameter :: maxt = 1000 + real :: temper(maxt),lambda(maxt),slope(maxt),yfunc(maxt) + integer :: nt + ! + ! set default values for input parameters + ! + real :: q_0_cgs = 0 ! heating rate coefficient [ergs/s/g] + real :: t_b1_seconds = 0 ! time of break 1 in heating rate [s] + real :: exp_1 = 0 ! exponent 1 in heating rate + real :: t_b2_seconds = 0 ! time of break 2 in heating rate [s] + real :: exp_2 = 0 ! exponent 2 in heating rate + real :: t_start_seconds = 0 ! time after merger at which the simulation begins [s] + + public :: write_options_rprocess,read_options_rprocess + public :: get_rprocess_heating_rate, energ_rprocess + + private + +contains + +!----------------------------------------------------------------------- +!+ +! get the r-process specific heating rate q at time t +!+ +!----------------------------------------------------------------------- +subroutine get_rprocess_heating_rate(q,t) + use physcon, only:days + use units, only:unit_ergg,utime + real, intent(in) :: t + real, intent(out) :: q + real :: t_seconds, t_days + real :: t_b1_days, t_b2_days + real :: t_start_days, t_new_days + real :: q_cgs + + t_seconds = t * utime + t_days = t_seconds / days + + t_b1_days = t_b1_seconds / days + t_b2_days = t_b2_seconds / days + + t_start_days = t_start_seconds / days + t_new_days = t_days + t_start_days + + !----- Heating rate in [ergs/s/g] + if (t_new_days < t_b1_days) then + q_cgs = q_0_cgs + else if ( (t_b1_days <= t_new_days) .and. (t_new_days < t_b2_days) ) then + q_cgs = q_0_cgs * (t_days/t_b1_days)**exp_1 + else + q_cgs = q_0_cgs * (t_b2_days/t_b1_days)**exp_1 * (t_days/t_b2_days)**exp_2 + endif + + q = q_cgs / (unit_ergg/utime) + +end subroutine get_rprocess_heating_rate + +!----------------------------------------------------------------------- +!+ +! get the change in the entropy variable K in the timestep t -> t+dt +!+ +!----------------------------------------------------------------------- +subroutine energ_rprocess(entrop_var,rho,t,dt) + use eos, only:gamma + use physcon, only:days + use units, only:unit_ergg,utime + real, intent(in) :: rho,t,dt + real, intent(inout) :: entrop_var + real :: dt_seconds + real :: q, q_cgs + real :: delta_u, delta_u_cgs + real :: delta_entrop_var + + dt_seconds = dt * utime + + call get_rprocess_heating_rate(q, t) + + q_cgs = q * (unit_ergg/utime) + + !----- Change in the specific internal energy [ergs/g] during the time step + delta_u_cgs = q_cgs * dt_seconds + delta_u = delta_u_cgs / unit_ergg + + !----- Change in the entropy variable K during the time step + !----- In Liptai & Price (2019), K is given the symbol K. ----- + delta_entrop_var = (gamma - 1.) * rho**(1. - gamma) * delta_u + entrop_var = entrop_var + delta_entrop_var + +end subroutine energ_rprocess + +!----------------------------------------------------------------------- +!+ +! writes input options to the input file +!+ +!----------------------------------------------------------------------- +subroutine write_options_rprocess(iunit) + use infile_utils, only:write_inopt + integer, intent(in) :: iunit + + call write_inopt(q_0_cgs,'q_0_cgs','heating rate coefficient [ergs/s/g]',iunit) + call write_inopt(t_b1_seconds,'t_b1_seconds','time of break 1 in heating rate [s]',iunit) + call write_inopt(exp_1,'exp_1','exponent 1 in heating rate',iunit) + call write_inopt(t_b2_seconds,'t_b2_seconds','time of break 2 in heating rate [s]',iunit) + call write_inopt(exp_2,'exp_2','exponent 2 in heating rate',iunit) + call write_inopt(t_start_seconds,'t_start_seconds','time after merger at which the simulation begins [s]',iunit) + +end subroutine write_options_rprocess + +!----------------------------------------------------------------------- +!+ +! reads input options from the input file +!+ +!----------------------------------------------------------------------- +subroutine read_options_rprocess(name,valstring,imatch,igotall,ierr) + use io, only:fatal + character(len=*), intent(in) :: name,valstring + logical, intent(out) :: imatch,igotall + integer, intent(out) :: ierr + integer, save :: ngot = 0 + + imatch = .true. + igotall = .false. ! cooling options are compulsory + select case(trim(name)) + case('q_0_cgs') + read(valstring,*,iostat=ierr) q_0_cgs + ngot = ngot + 1 + case('t_b1_seconds') + read(valstring,*,iostat=ierr) t_b1_seconds + ngot = ngot + 1 + case('exp_1') + read(valstring,*,iostat=ierr) exp_1 + ngot = ngot + 1 + case('t_b2_seconds') + read(valstring,*,iostat=ierr) t_b2_seconds + ngot = ngot + 1 + case('exp_2') + read(valstring,*,iostat=ierr) exp_2 + ngot = ngot + 1 + case('t_start_seconds') + read(valstring,*,iostat=ierr) t_start_seconds + ngot = ngot + 1 + case default + imatch = .false. + end select + if (ngot >= 6) igotall = .true. + +end subroutine read_options_rprocess + +end module rprocess_heating \ No newline at end of file diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 30611f97a..2416c8f14 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -26,7 +26,7 @@ ! ! RUNTIME PARAMETERS: None ! -! DEPENDENCIES: bowen_dust, chem, cons2prim, cons2primsolver, coolfunc, +! DEPENDENCIES: bowen_dust, chem, cons2prim, cons2primsolver, coolfunc, rprocess_heating, ! damping, deriv, derivutils, dim, eos, extern_gr, externalforces, ! growth, io, io_summary, metric_tools, mpiutils, options, part, ptmass, ! timestep, timestep_ind, timestep_sts, timing @@ -116,6 +116,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) use part, only:nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,ibin_wake,this_is_a_test use io_summary, only:summary_printout,summary_variable,iosumtvi,iowake use coolfunc, only:energ_coolfunc + use rprocess_heating, only:energ_rprocess #ifdef IND_TIMESTEPS use timestep, only:dtmax,dtmax_ifactor,dtdiff use timestep_ind, only:get_dt,nbinmax,decrease_dtmax,ibinnow @@ -414,6 +415,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) !$omp shared(dustprop,ddustprop,dustproppred) & !$omp shared(xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,nptmass,massoftype) & !$omp shared(dtsph,icooling) & +!$omp shared(timei) & #ifdef IND_TIMESTEPS !$omp shared(ibin,ibin_old,ibin_sts,twas,timei,use_sts,dtsph_next,ibin_wake,sts_it_n) & !$omp shared(ibin_dts,nbinmax,ibinnow) & @@ -519,6 +521,10 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) pxyzu(2,i) = pyi pxyzu(3,i) = pzi pxyzu(4,i) = eni + + !----- pxyzu(:,i) is the array (px,py,pz,K) for the ith SPH particle, where (px,py,pz) are the conserved coordinate momenta and K is the entropy variable in the rest frame. + !----- In Liptai & Price (2019), K is given the symbol K. ----- + if (icooling==4) call energ_rprocess(pxyzu(4,i),rhoh(xyzh(4,i),massoftype(itype)),timei,dtsph) else vxi = vxyzu(1,i) + hdtsph*fxyzu(1,i) vyi = vxyzu(2,i) + hdtsph*fxyzu(2,i) From 226d057a6ccee468f26080cf8818dd9e8bb89d99 Mon Sep 17 00:00:00 2001 From: sdarbha Date: Tue, 4 Aug 2020 10:04:48 -0700 Subject: [PATCH 03/15] Minor fix to cooling.f90 --- src/main/cooling.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/main/cooling.f90 b/src/main/cooling.f90 index 00e491523..20b76a5ae 100644 --- a/src/main/cooling.f90 +++ b/src/main/cooling.f90 @@ -193,6 +193,7 @@ subroutine read_options_cooling(name,valstring,imatch,igotall,ierr) igotall = .false. ! cooling options are compulsory igotallh2 = .true. igotallcf = .true. + igotallrp = .true. if (maxvxyzu < 4) igotall = .true. ! options unnecessary if isothermal select case(trim(name)) From 215dac0441c6e2689480b2083125afcb03278cae Mon Sep 17 00:00:00 2001 From: sdarbha Date: Tue, 18 Aug 2020 07:55:38 -0700 Subject: [PATCH 04/15] Updated print format in analysis_bhnsejecta.f90 --- src/utils/analysis_bhnsejecta.f90 | 33 ++++++++++++++++++++----------- 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/src/utils/analysis_bhnsejecta.f90 b/src/utils/analysis_bhnsejecta.f90 index 8d289adbc..1c7ee9d52 100644 --- a/src/utils/analysis_bhnsejecta.f90 +++ b/src/utils/analysis_bhnsejecta.f90 @@ -32,7 +32,9 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) use metric, only:mass1 use units, only:umass,udist,utime,unit_velocity,unit_energ,unit_density,unit_ergg use physcon, only:solarm,pi,solarr,gg,c,eV,radconst,kboltz - use part, only:hfact + use part, only:hfact,pxyzu,dens,metrics,metricderivs + use metric_tools, only:init_metric + use utils_gr, only:h2dens character(len=*), intent(in) :: dumpfile integer, intent(in) :: numfile,npart,iunit real, intent(in) :: xyzh(:,:),vxyzu(:,:) @@ -41,7 +43,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) real :: x,y,z,h real :: r, theta, phi, e, l, u_r, u_theta, dt_dtau, dr_dtau, dtheta_dtau, dphi_dtau, dr_dt, dtheta_dt, dphi_dt real :: vx, vy, vz, uint - real :: m, rho, delta + real :: m, rho real :: T, s, y_e real :: q real :: time_output @@ -67,6 +69,8 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) write(1,'(A)',IOSTAT=iostatus) '# [11]= T' write(1,'(A)',IOSTAT=iostatus) '# [12]= rho' + call init_metric(npart,xyzh,metrics,metricderivs) + do i=1,npart m = pmass @@ -80,9 +84,9 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) !----- vxyzu(:,i) is the array (vx,vy,vz,u) for the ith SPH particle, where (vx,vy,vz) are the coordinate velocities and u is the specific internal energy in the rest frame. !----- In Price et al PASA (2018) and in Liptai & Price (2019), u is given the symbol u. - vx = vxyzu(1,i) - vy = vxyzu(2,i) - vz = vxyzu(3,i) + vx = vxyzu(1,i) !----- vx here is defined as vx = dx/dt + vy = vxyzu(2,i) !----- vy here is defined as vy = dy/dt + vz = vxyzu(3,i) !----- vz here is defined as vz = dz/dt uint = vxyzu(4,i) !----- the specific internal energy in the rest frame r = (x**2. + y**2. + z**2.)**(1./2.) @@ -108,7 +112,12 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) u_theta = dtheta_dtau * (r**2.) l = dphi_dtau * ((r**2.)*(sin(theta)**2.)) - rho = m * (hfact/h)**3. + !----- rho here is the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho + !----- Note: + !----- the routines in GR Phantom use a different naming convention: + !----- they use the variable name 'dens' for the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho + !----- they use the variable name 'rho' for the "relativistic rest mass density." In Liptai & Price (2019), the "relativistic rest mass density" is given the symbol rho^* + call h2dens(rho, xyzh(:,i), metrics(:,:,:,i), vxyzu(1:3,i)) T = ( (uint*unit_ergg) * (rho*unit_density) / radconst )**(1./4.) !----- temperature in K @@ -133,16 +142,16 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) write(1,'(ES11.4,1X)',advance='no',IOSTAT=iostatus) m - write(1,'(F10.3,1X)',advance='no',IOSTAT=iostatus) r + write(1,'(ES12.5,1X)',advance='no',IOSTAT=iostatus) r write(1,'(F9.5,1X)',advance='no',IOSTAT=iostatus) theta - write(1,'(F10.6,1X)',advance='no',IOSTAT=iostatus) phi + write(1,'(F9.5,1X)',advance='no',IOSTAT=iostatus) phi write(1,'(F10.5,1X)',advance='no',IOSTAT=iostatus) e - write(1,'(F10.4,1X)',advance='no',IOSTAT=iostatus) l - write(1,'(F12.6,1X)',advance='no',IOSTAT=iostatus) u_r - write(1,'(F10.4,1X)',advance='no',IOSTAT=iostatus) u_theta + write(1,'(ES12.5,1X)',advance='no',IOSTAT=iostatus) l + write(1,'(ES12.5,1X)',advance='no',IOSTAT=iostatus) u_r + write(1,'(ES12.5,1X)',advance='no',IOSTAT=iostatus) u_theta write(1,'(F10.7,1X)',advance='no',IOSTAT=iostatus) y_e write(1,'(F10.7,1X)',advance='no',IOSTAT=iostatus) s - write(1,'(F10.6,1X)',advance='no',IOSTAT=iostatus) T + write(1,'(ES12.4,1X)',advance='no',IOSTAT=iostatus) T write(1,'(ES12.5)',IOSTAT=iostatus) rho enddo From 876e9cf8e35181fd44dcec003c4f97cef4e97ee2 Mon Sep 17 00:00:00 2001 From: sdarbha Date: Thu, 27 Aug 2020 16:02:45 -0700 Subject: [PATCH 05/15] Fixed mistakes related to using conserved density 'rho' vs rest mass density 'dens' --- src/main/rprocess_heating.f90 | 6 ++--- src/main/step_leapfrog.F90 | 7 +++-- src/setup/setup_grbhnsejecta.f90 | 45 +++++++++++++++++-------------- src/utils/analysis_bhnsejecta.f90 | 26 +++++++++--------- 4 files changed, 46 insertions(+), 38 deletions(-) diff --git a/src/main/rprocess_heating.f90 b/src/main/rprocess_heating.f90 index daee2c1c9..b2ff1465e 100644 --- a/src/main/rprocess_heating.f90 +++ b/src/main/rprocess_heating.f90 @@ -90,11 +90,11 @@ end subroutine get_rprocess_heating_rate ! get the change in the entropy variable K in the timestep t -> t+dt !+ !----------------------------------------------------------------------- -subroutine energ_rprocess(entrop_var,rho,t,dt) +subroutine energ_rprocess(entrop_var,dens,t,dt) use eos, only:gamma use physcon, only:days use units, only:unit_ergg,utime - real, intent(in) :: rho,t,dt + real, intent(in) :: dens,t,dt real, intent(inout) :: entrop_var real :: dt_seconds real :: q, q_cgs @@ -113,7 +113,7 @@ subroutine energ_rprocess(entrop_var,rho,t,dt) !----- Change in the entropy variable K during the time step !----- In Liptai & Price (2019), K is given the symbol K. ----- - delta_entrop_var = (gamma - 1.) * rho**(1. - gamma) * delta_u + delta_entrop_var = (gamma - 1.) * delta_u / dens**(gamma - 1.) entrop_var = entrop_var + delta_entrop_var end subroutine energ_rprocess diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 2416c8f14..4a1c4bdcf 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -415,7 +415,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) !$omp shared(dustprop,ddustprop,dustproppred) & !$omp shared(xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,nptmass,massoftype) & !$omp shared(dtsph,icooling) & -!$omp shared(timei) & +!$omp shared(dens,timei) & !----- This line is needed to call the function energ_rprocess below #ifdef IND_TIMESTEPS !$omp shared(ibin,ibin_old,ibin_sts,twas,timei,use_sts,dtsph_next,ibin_wake,sts_it_n) & !$omp shared(ibin_dts,nbinmax,ibinnow) & @@ -524,7 +524,10 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) !----- pxyzu(:,i) is the array (px,py,pz,K) for the ith SPH particle, where (px,py,pz) are the conserved coordinate momenta and K is the entropy variable in the rest frame. !----- In Liptai & Price (2019), K is given the symbol K. ----- - if (icooling==4) call energ_rprocess(pxyzu(4,i),rhoh(xyzh(4,i),massoftype(itype)),timei,dtsph) + ! if (icooling==4) call energ_rprocess(pxyzu(4,i), rhoh(xyzh(4,i),massoftype(itype)), timei, dtsph) + if (icooling==4) then + call energ_rprocess(pxyzu(4,i), dens(i), timei, dtsph) + endif else vxi = vxyzu(1,i) + hdtsph*fxyzu(1,i) vyi = vxyzu(2,i) + hdtsph*fxyzu(2,i) diff --git a/src/setup/setup_grbhnsejecta.f90 b/src/setup/setup_grbhnsejecta.f90 index 19a8a0927..80adc26c2 100644 --- a/src/setup/setup_grbhnsejecta.f90 +++ b/src/setup/setup_grbhnsejecta.f90 @@ -18,7 +18,7 @@ ! particle_file_name -- name of SpEC particle file ! ! DEPENDENCIES: eos, externalforces, infile_utils, io, metric, part, -! physcon, rho_profile, timestep, units, vectorutils +! physcon, timestep, units !+ !-------------------------------------------------------------------------- module setup @@ -65,11 +65,11 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, integer :: ierr,i logical :: iexist real :: x,y,z - real, allocatable :: m_array(:), r_array(:), theta_array(:), phi_array(:), e_array(:), l_array(:), u_r_array(:), u_theta_array(:), y_e_array(:), s_array(:), T_array(:), rho_array(:) + real, allocatable :: m_array(:), r_array(:), theta_array(:), phi_array(:), e_array(:), l_array(:), u_r_array(:), u_theta_array(:), y_e_array(:), s_array(:), T_array(:), dens_array(:) real :: m_ejecta real :: r, theta, phi, e, l, u_r, u_theta, dt_dtau, dr_dtau, dtheta_dtau, dphi_dtau, dr_dt, dtheta_dt, dphi_dt real :: vx, vy, vz - real :: m, rho, delta + real :: m, dens real :: T, y_e ! @@ -134,9 +134,9 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, allocate(y_e_array(npart)) allocate(s_array(npart)) allocate(T_array(npart)) - allocate(rho_array(npart)) + allocate(dens_array(npart)) - call read_particle_data(npart, m_array, r_array, theta_array, phi_array, e_array, l_array, u_r_array, u_theta_array, y_e_array, s_array, T_array, rho_array) + call read_particle_data(npart, m_array, r_array, theta_array, phi_array, e_array, l_array, u_r_array, u_theta_array, y_e_array, s_array, T_array, dens_array) m_ejecta = 0 do i=1,npart @@ -152,7 +152,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! y_e_array(i) = y_e_array(i) !----- Y_e is dimensionless, so we don't need to do any unit conversion ! s_array(i) = s_array(i)* !----- I'm not sure what the units are for the entropy in the particle file, I have to look into this T_array(i) = T_array(i) * (1.e6)*eV / unit_energ - rho_array(i) = rho_array(i) * ((c**6.)/((gg**3.)*(solarm**2.))) / unit_density + dens_array(i) = dens_array(i) * ((c**6.)/((gg**3.)*(solarm**2.))) / unit_density m_ejecta = m_ejecta + m_array(i) enddo @@ -179,21 +179,26 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, dr_dt = dr_dtau / dt_dtau dtheta_dt = dtheta_dtau / dt_dtau dphi_dt = dphi_dtau / dt_dtau - vx = sin(theta)*cos(phi)*dr_dt + r*cos(theta)*cos(phi)*dtheta_dt - r*sin(theta)*sin(phi)*dphi_dt - vy = sin(theta)*sin(phi)*dr_dt + r*cos(theta)*sin(phi)*dtheta_dt + r*sin(theta)*cos(phi)*dphi_dt - vz = cos(theta)*dr_dt - r*sin(theta)*dtheta_dt + vx = sin(theta)*cos(phi)*dr_dt + r*cos(theta)*cos(phi)*dtheta_dt - r*sin(theta)*sin(phi)*dphi_dt !----- vx here is defined as vx = dx/dt + vy = sin(theta)*sin(phi)*dr_dt + r*cos(theta)*sin(phi)*dtheta_dt + r*sin(theta)*cos(phi)*dphi_dt !----- vy here is defined as vy = dy/dt + vz = cos(theta)*dr_dt - r*sin(theta)*dtheta_dt !----- vz here is defined as vz = dz/dt + m = m_array(i) + !----- dens is the rest mass density + !----- The routines in GR Phantom use the following naming convention: + !----- The variable name 'dens' is the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho + !----- The variable name 'rho' is the "relativistic rest-mass density"/"conserved density". In Liptai & Price (2019), the "relativistic rest-mass density" is given the symbol rho^* + !----- The initial value for dens provided here is irrelevant, the code will calculate the relativistic rest-mass density rho from the particle distribution. + dens = dens_array(i) !----- xyzh(:,i) is the array (x,y,z,h) for the ith SPH particle, where (x,y,z) are the coordinates and h is the smoothing length. !----- In Price et al PASA (2018) and in Liptai & Price (2019), h is given by the symbols h and h_a. - m = m_array(i) - rho = rho_array(i) - delta = (m/rho)**(1./3.) !----- local interparticle spacing - xyzh(4,i) = hfact*delta !----- smoothing length + !----- The initial value for the smoothing length h provided here is irrelevant, the code will calculate it from the particle distribution. We simply set it equal to the proportionality factor hfact. + xyzh(4,i) = hfact xyzh(1:3,i) = (/x,y,z/) !----- vxyzu(:,i) is the array (vx,vy,vz,u) for the ith SPH particle, where (vx,vy,vz) are the coordinate velocities and u is the specific internal energy in the rest frame. !----- In Price et al PASA (2018) and in Liptai & Price (2019), u is given the symbol u. T = T_array(i) y_e = y_e_array(i) - vxyzu(4,i) = ( radconst * ((T*unit_energ/kboltz)**4.) / (rho*unit_density) ) / unit_ergg !----- the specific internal energy in the rest frame, u = a * T^4 / rho + vxyzu(4,i) = ( radconst * ((T*unit_energ/kboltz)**4.) / (dens*unit_density) ) / unit_ergg !----- the specific internal energy in the rest frame, u = a * T^4 / dens vxyzu(1:3,i) = (/vx,vy,vz/) enddo @@ -208,7 +213,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, deallocate(y_e_array) deallocate(s_array) deallocate(T_array) - deallocate(rho_array) + deallocate(dens_array) if (id==master) then print "(/,a)", ' EJECTA SETUP:' @@ -268,7 +273,7 @@ end subroutine read_setupfile subroutine get_num_particles(npart) integer :: iostatus, i integer, intent(out) :: npart - real :: m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, rho + real :: m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, dens print *, 'getting number of particles from '//trim(particle_file_name) @@ -280,7 +285,7 @@ subroutine get_num_particles(npart) i = 0 do - read(1,*,IOSTAT=iostatus) m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, rho + read(1,*,IOSTAT=iostatus) m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, dens ! read failed if (iostatus > 0) then print *, 'read failed, input data not in correct format' @@ -302,10 +307,10 @@ end subroutine get_num_particles ! !---Read particle data-------------------------------------------------- ! -subroutine read_particle_data(npart, m_array, r_array, theta_array, phi_array, e_array, l_array, u_r_array, u_theta_array, y_e_array, s_array, T_array, rho_array) +subroutine read_particle_data(npart, m_array, r_array, theta_array, phi_array, e_array, l_array, u_r_array, u_theta_array, y_e_array, s_array, T_array, dens_array) integer :: iostatus, i integer, intent(in) :: npart - real, intent(out) :: m_array(:), r_array(:), theta_array(:), phi_array(:), e_array(:), l_array(:), u_r_array(:), u_theta_array(:), y_e_array(:), s_array(:), T_array(:), rho_array(:) + real, intent(out) :: m_array(:), r_array(:), theta_array(:), phi_array(:), e_array(:), l_array(:), u_r_array(:), u_theta_array(:), y_e_array(:), s_array(:), T_array(:), dens_array(:) print *, 'reading particle data from '//trim(particle_file_name) @@ -317,7 +322,7 @@ subroutine read_particle_data(npart, m_array, r_array, theta_array, phi_array, e i = 1 do - read(1,*,IOSTAT=iostatus) m_array(i), r_array(i), theta_array(i), phi_array(i), e_array(i), l_array(i), u_r_array(i), u_theta_array(i), y_e_array(i), s_array(i), T_array(i), rho_array(i) + read(1,*,IOSTAT=iostatus) m_array(i), r_array(i), theta_array(i), phi_array(i), e_array(i), l_array(i), u_r_array(i), u_theta_array(i), y_e_array(i), s_array(i), T_array(i), dens_array(i) ! read failed if (iostatus > 0) then print *, 'read failed, input data not in correct format' diff --git a/src/utils/analysis_bhnsejecta.f90 b/src/utils/analysis_bhnsejecta.f90 index 1c7ee9d52..9998d4df4 100644 --- a/src/utils/analysis_bhnsejecta.f90 +++ b/src/utils/analysis_bhnsejecta.f90 @@ -42,8 +42,8 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) integer :: iostatus, i real :: x,y,z,h real :: r, theta, phi, e, l, u_r, u_theta, dt_dtau, dr_dtau, dtheta_dtau, dphi_dtau, dr_dt, dtheta_dt, dphi_dt - real :: vx, vy, vz, uint - real :: m, rho + real :: vx, vy, vz, u + real :: m, densi real :: T, s, y_e real :: q real :: time_output @@ -87,7 +87,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) vx = vxyzu(1,i) !----- vx here is defined as vx = dx/dt vy = vxyzu(2,i) !----- vy here is defined as vy = dy/dt vz = vxyzu(3,i) !----- vz here is defined as vz = dz/dt - uint = vxyzu(4,i) !----- the specific internal energy in the rest frame + u = vxyzu(4,i) !----- the specific internal energy in the rest frame r = (x**2. + y**2. + z**2.)**(1./2.) q = (x**2. + y**2.)**(1./2.) @@ -112,18 +112,18 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) u_theta = dtheta_dtau * (r**2.) l = dphi_dtau * ((r**2.)*(sin(theta)**2.)) - !----- rho here is the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho - !----- Note: - !----- the routines in GR Phantom use a different naming convention: - !----- they use the variable name 'dens' for the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho - !----- they use the variable name 'rho' for the "relativistic rest mass density." In Liptai & Price (2019), the "relativistic rest mass density" is given the symbol rho^* - call h2dens(rho, xyzh(:,i), metrics(:,:,:,i), vxyzu(1:3,i)) + !----- dens(:) is the array of rest mass densities of the SPH particles, where dens(i) is the rest mass density of the ith SPH particle + !----- The routines in GR Phantom use the following naming convention: + !----- The variable name 'dens' is the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho + !----- The variable name 'rho' is the "relativistic rest-mass density"/"conserved density". In Liptai & Price (2019), the "relativistic rest-mass density" is given the symbol rho^* + !call h2dens(densi, xyzh(:,i), metrics(:,:,:,i), vxyzu(1:3,i)) + densi = dens(i) - T = ( (uint*unit_ergg) * (rho*unit_density) / radconst )**(1./4.) !----- temperature in K + T = ( (u*unit_ergg) * (densi*unit_density) / radconst )**(1./4.) !----- temperature in K y_e = 0. !----- The SPH gas particles do not store the electron fraction Y_e since the code doesn't evolve the composition, so I just set it to zero here - s = 0. !----- I'm not sure how to calculate the entropy of the SPH particles, or what the units are for the entropy in the particle file, so I just set it to zero here + s = 0. !----- I'm not sure what the units are for the entropy in the particle file, so I just set it to zero here @@ -137,7 +137,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) u_r = u_r * unit_velocity / c u_theta = u_theta * (udist*unit_velocity) / (gg*solarm/c) T = kboltz * T / ( (1.e6)*eV ) !----- temperature in MeV - rho = rho * unit_density / ( solarm / (gg*solarm/(c**2.))**3. ) + densi = densi * unit_density / ( solarm / (gg*solarm/(c**2.))**3. ) @@ -152,7 +152,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) write(1,'(F10.7,1X)',advance='no',IOSTAT=iostatus) y_e write(1,'(F10.7,1X)',advance='no',IOSTAT=iostatus) s write(1,'(ES12.4,1X)',advance='no',IOSTAT=iostatus) T - write(1,'(ES12.5)',IOSTAT=iostatus) rho + write(1,'(ES12.5)',IOSTAT=iostatus) densi enddo From 14d8f55840aeca7509e780abc9805262883d94b1 Mon Sep 17 00:00:00 2001 From: sdarbha Date: Thu, 27 Aug 2020 17:30:18 -0700 Subject: [PATCH 06/15] Updated a variable name and some comments --- src/setup/setup_grbhnsejecta.f90 | 16 ++++++++-------- src/utils/analysis_bhnsejecta.f90 | 1 + 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/setup/setup_grbhnsejecta.f90 b/src/setup/setup_grbhnsejecta.f90 index 80adc26c2..43862e018 100644 --- a/src/setup/setup_grbhnsejecta.f90 +++ b/src/setup/setup_grbhnsejecta.f90 @@ -69,7 +69,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, real :: m_ejecta real :: r, theta, phi, e, l, u_r, u_theta, dt_dtau, dr_dtau, dtheta_dtau, dphi_dtau, dr_dt, dtheta_dt, dphi_dt real :: vx, vy, vz - real :: m, dens + real :: m, densi real :: T, y_e ! @@ -183,12 +183,12 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, vy = sin(theta)*sin(phi)*dr_dt + r*cos(theta)*sin(phi)*dtheta_dt + r*sin(theta)*cos(phi)*dphi_dt !----- vy here is defined as vy = dy/dt vz = cos(theta)*dr_dt - r*sin(theta)*dtheta_dt !----- vz here is defined as vz = dz/dt m = m_array(i) - !----- dens is the rest mass density + !----- densi is the rest mass density of the ith SPH particle !----- The routines in GR Phantom use the following naming convention: !----- The variable name 'dens' is the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho - !----- The variable name 'rho' is the "relativistic rest-mass density"/"conserved density". In Liptai & Price (2019), the "relativistic rest-mass density" is given the symbol rho^* - !----- The initial value for dens provided here is irrelevant, the code will calculate the relativistic rest-mass density rho from the particle distribution. - dens = dens_array(i) + !----- The variable name 'rho' is the "relativistic rest-mass density", aka the "conserved density". In Liptai & Price (2019), the "relativistic rest-mass density" is given the symbol rho^* + !----- The initial value for densi provided here is irrelevant, the code will calculate the relativistic rest-mass density rho from the particle distribution. + densi = dens_array(i) !----- xyzh(:,i) is the array (x,y,z,h) for the ith SPH particle, where (x,y,z) are the coordinates and h is the smoothing length. !----- In Price et al PASA (2018) and in Liptai & Price (2019), h is given by the symbols h and h_a. !----- The initial value for the smoothing length h provided here is irrelevant, the code will calculate it from the particle distribution. We simply set it equal to the proportionality factor hfact. @@ -198,7 +198,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, !----- In Price et al PASA (2018) and in Liptai & Price (2019), u is given the symbol u. T = T_array(i) y_e = y_e_array(i) - vxyzu(4,i) = ( radconst * ((T*unit_energ/kboltz)**4.) / (dens*unit_density) ) / unit_ergg !----- the specific internal energy in the rest frame, u = a * T^4 / dens + vxyzu(4,i) = ( radconst * ((T*unit_energ/kboltz)**4.) / (densi*unit_density) ) / unit_ergg !----- the specific internal energy in the rest frame, u = a * T^4 / densi vxyzu(1:3,i) = (/vx,vy,vz/) enddo @@ -273,7 +273,7 @@ end subroutine read_setupfile subroutine get_num_particles(npart) integer :: iostatus, i integer, intent(out) :: npart - real :: m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, dens + real :: m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, densi print *, 'getting number of particles from '//trim(particle_file_name) @@ -285,7 +285,7 @@ subroutine get_num_particles(npart) i = 0 do - read(1,*,IOSTAT=iostatus) m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, dens + read(1,*,IOSTAT=iostatus) m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, densi ! read failed if (iostatus > 0) then print *, 'read failed, input data not in correct format' diff --git a/src/utils/analysis_bhnsejecta.f90 b/src/utils/analysis_bhnsejecta.f90 index 9998d4df4..905b4e5c4 100644 --- a/src/utils/analysis_bhnsejecta.f90 +++ b/src/utils/analysis_bhnsejecta.f90 @@ -113,6 +113,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) l = dphi_dtau * ((r**2.)*(sin(theta)**2.)) !----- dens(:) is the array of rest mass densities of the SPH particles, where dens(i) is the rest mass density of the ith SPH particle + !----- The dens(:) array resides in the part module in the file part.F90, and is updated when the code calls the conservative to primitive variable transformation function !----- The routines in GR Phantom use the following naming convention: !----- The variable name 'dens' is the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho !----- The variable name 'rho' is the "relativistic rest-mass density"/"conserved density". In Liptai & Price (2019), the "relativistic rest-mass density" is given the symbol rho^* From e093aa19e0290a1fbaaa4d9dca72aa79dc72d7eb Mon Sep 17 00:00:00 2001 From: sdarbha Date: Mon, 31 Aug 2020 16:39:08 -0700 Subject: [PATCH 07/15] Fixed mistake in analysis_bhnsejecta.f90, minor change to step_leapfrog.F90 --- src/main/step_leapfrog.F90 | 9 ++++++--- src/utils/analysis_bhnsejecta.f90 | 12 ++++++------ 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 4a1c4bdcf..4dd35f73c 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -524,10 +524,13 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) !----- pxyzu(:,i) is the array (px,py,pz,K) for the ith SPH particle, where (px,py,pz) are the conserved coordinate momenta and K is the entropy variable in the rest frame. !----- In Liptai & Price (2019), K is given the symbol K. ----- + !----- dens(:) is the array of rest mass densities of the SPH particles, where dens(i) is the rest mass density of the ith SPH particle + !----- The dens(:) array resides in the part module in the file part.F90, and is updated when the code calls the conservative to primitive variable transformation function + !----- The routines in GR Phantom use the following naming convention: + !----- The variable name 'dens' is the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho + !----- The variable name 'rho' is the "relativistic rest-mass density"/"conserved density". In Liptai & Price (2019), the "relativistic rest-mass density" is given the symbol rho^* ! if (icooling==4) call energ_rprocess(pxyzu(4,i), rhoh(xyzh(4,i),massoftype(itype)), timei, dtsph) - if (icooling==4) then - call energ_rprocess(pxyzu(4,i), dens(i), timei, dtsph) - endif + if (icooling==4) call energ_rprocess(pxyzu(4,i), dens(i), timei, dtsph) else vxi = vxyzu(1,i) + hdtsph*fxyzu(1,i) vyi = vxyzu(2,i) + hdtsph*fxyzu(2,i) diff --git a/src/utils/analysis_bhnsejecta.f90 b/src/utils/analysis_bhnsejecta.f90 index 905b4e5c4..231ef2116 100644 --- a/src/utils/analysis_bhnsejecta.f90 +++ b/src/utils/analysis_bhnsejecta.f90 @@ -32,7 +32,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) use metric, only:mass1 use units, only:umass,udist,utime,unit_velocity,unit_energ,unit_density,unit_ergg use physcon, only:solarm,pi,solarr,gg,c,eV,radconst,kboltz - use part, only:hfact,pxyzu,dens,metrics,metricderivs + use part, only:hfact,pxyzu,metrics,metricderivs use metric_tools, only:init_metric use utils_gr, only:h2dens character(len=*), intent(in) :: dumpfile @@ -46,7 +46,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) real :: m, densi real :: T, s, y_e real :: q - real :: time_output + real :: time_output, time_output_seconds character(len=120) :: particle_file_name, info_file_name @@ -112,13 +112,11 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) u_theta = dtheta_dtau * (r**2.) l = dphi_dtau * ((r**2.)*(sin(theta)**2.)) - !----- dens(:) is the array of rest mass densities of the SPH particles, where dens(i) is the rest mass density of the ith SPH particle - !----- The dens(:) array resides in the part module in the file part.F90, and is updated when the code calls the conservative to primitive variable transformation function + !----- densi is the rest mass density of the ith SPH particle !----- The routines in GR Phantom use the following naming convention: !----- The variable name 'dens' is the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho !----- The variable name 'rho' is the "relativistic rest-mass density"/"conserved density". In Liptai & Price (2019), the "relativistic rest-mass density" is given the symbol rho^* - !call h2dens(densi, xyzh(:,i), metrics(:,:,:,i), vxyzu(1:3,i)) - densi = dens(i) + call h2dens(densi, xyzh(:,i), metrics(:,:,:,i), vxyzu(1:3,i)) T = ( (u*unit_ergg) * (densi*unit_density) / radconst )**(1./4.) !----- temperature in K @@ -164,12 +162,14 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) !----- Convert time from code units to units of particle file, i.e. first convert to cgs then normalize to units of particle file !----- Note: the particle file is in units of G = c = M_Sun = 1 (for most of the variables), whereas in set_units in setup_grbhnsejecta.f90 we set the code units to be G = c = M_bh = 1 time_output = time * utime / (gg*solarm/(c**3.)) + time_output_seconds = time * utime info_file_name = trim(dumpfile)//'.info' open(2, file=info_file_name, action='write', status='replace') write(2,'(A,1X,ES10.3,1X,A)',IOSTAT=iostatus) 'time =', time_output, '[G = c = M_Sun = 1]' + write(2,'(A,1X,ES10.3,1X,A)',IOSTAT=iostatus) 'time_seconds =', time_output_seconds, '[s]' close(2) From 87ffa23330ff74ec809e72eff0e6e62409f69c6b Mon Sep 17 00:00:00 2001 From: sdarbha Date: Wed, 2 Sep 2020 09:45:42 -0700 Subject: [PATCH 08/15] Updated rprocess heating to align with new Phantom changes --- src/main/checksetup.F90 | 3 +- src/main/cooling.F90 | 91 ++++++++++++++++++++++++++++++++++++-- src/main/force.F90 | 9 ---- src/main/step_leapfrog.F90 | 43 +++++++++++++----- 4 files changed, 120 insertions(+), 26 deletions(-) diff --git a/src/main/checksetup.F90 b/src/main/checksetup.F90 index 6fb97f6ea..7afdf6879 100644 --- a/src/main/checksetup.F90 +++ b/src/main/checksetup.F90 @@ -419,8 +419,9 @@ subroutine check_setup(nerror,nwarn,restart) if (id==master) & write(*,"(a,2(es10.3,', '),es10.3,a)") ' Centre of mass is at (x,y,z) = (',xcom,')' +! (Siva Darbha): I modified this condition, I think the change is correct ! if (.not.h2chemistry .and. maxvxyzu >= 4 .and. icooling >= 1 .and. iexternalforce/=iext_corotate) then - if (.not.h2chemistry .and. maxvxyzu >= 4 .and. icooling == 1 .and. iexternalforce/=iext_corotate) then + if (.not.h2chemistry .and. maxvxyzu >= 4 .and. icooling == 3 .and. iexternalforce/=iext_corotate) then if (dot_product(xcom,xcom) > 1.e-2) then print*,'Error in setup: Gammie (2001) cooling (icooling=1) assumes Omega = 1./r^1.5' print*,' but the centre of mass is not at the origin!' diff --git a/src/main/cooling.F90 b/src/main/cooling.F90 index 38da61709..35d858035 100644 --- a/src/main/cooling.F90 +++ b/src/main/cooling.F90 @@ -24,8 +24,15 @@ module cooling ! - icool_radiation_H0 : *H0 cooling on/off* ! - icool_relax_bowen : *Bowen (diffusive) relaxation on/off* ! - icool_relax_stefan : *radiative relaxation on/off* -! - icooling : *cooling function (0=off, 1=physics, 2=cooling table, 3=gammie)* +! - icooling : *cooling function (0=off, 1=physics, 2=cooling table, 3=gammie, 4=rprocess heating rate (only works with GR Phantom for now))* ! - temp_floor : *Minimum allowed temperature in K* +! ----- Input parameters for the rprocess heating rate function: +! - q_0_cgs : heating rate coefficient, in [ergs/s/g] +! - t_b1_seconds : time of break 1 in heating rate, in [s]' +! - exp_1 : exponent 1 in radioactive power component +! - t_b2_seconds : time of break 2 in heating rate, in [s]' +! - exp_2 : exponent 2 in radioactive power component +! - t_start_seconds : time after merger at which the simulation begins [s] ! ! :Dependencies: datafiles, eos, h2cooling, infile_utils, io, options, ! part, physcon, timestep, units @@ -55,6 +62,13 @@ module cooling real :: Tgrid(nTg) integer :: icool_radiation_H0 = 0, icool_relax_Bowen = 0, icool_dust_collision = 0, icool_relax_Stefan = 0 character(len=120) :: cooltable = 'cooltable.dat' + !----- Input parameters for the rprocess heating rate function: + real :: q_0_cgs = 0 ! heating rate coefficient [ergs/s/g] + real :: t_b1_seconds = 0 ! time of break 1 in heating rate [s] + real :: exp_1 = 0 ! exponent 1 in heating rate + real :: t_b2_seconds = 0 ! time of break 2 in heating rate [s] + real :: exp_2 = 0 ! exponent 2 in heating rate + real :: t_start_seconds = 0 ! time after merger at which the simulation begins [s] contains @@ -396,12 +410,15 @@ end subroutine implicit_cooling ! this routine returns the effective cooling rate du/dt ! !----------------------------------------------------------------------- -subroutine energ_cooling (xi, yi, zi, u, dudt, rho, dt, Trad, mu_in, K2, kappa) - real, intent(in) :: xi, yi, zi, u, rho, dt !in code unit +subroutine energ_cooling (xi, yi, zi, u, dudt, rho, dt, Trad, mu_in, K2, kappa, t) + real, intent(in) :: xi, yi, zi, u, rho, dt ! in code unit real, intent(in), optional :: Trad, mu_in, K2, kappa ! in cgs!!! + real, intent(in), optional :: t ! in code units real, intent(inout) :: dudt select case (icooling) + case (4) + call get_rprocess_heating_rate(dudt, t) case (3) call cooling_Gammie(xi, yi, zi, u, dudt) case (2) @@ -414,6 +431,43 @@ subroutine energ_cooling (xi, yi, zi, u, dudt, rho, dt, Trad, mu_in, K2, kappa) end subroutine energ_cooling +!----------------------------------------------------------------------- +! +! This function returns the rprocess specific heating rate q at time t, both in code units +! +!----------------------------------------------------------------------- +subroutine get_rprocess_heating_rate(q,t) + use physcon, only:days + use units, only:unit_ergg,utime + real, intent(in) :: t !----- time, in code units + real, intent(out) :: q !----- specific heating rate, in code units + real :: t_seconds, t_days + real :: t_b1_days, t_b2_days + real :: t_start_days, t_new_days + real :: q_cgs + + t_seconds = t * utime + t_days = t_seconds / days + + t_b1_days = t_b1_seconds / days + t_b2_days = t_b2_seconds / days + + t_start_days = t_start_seconds / days + t_new_days = t_days + t_start_days + + !----- Heating rate in [ergs/s/g] + if (t_new_days < t_b1_days) then + q_cgs = q_0_cgs + else if ( (t_b1_days <= t_new_days) .and. (t_new_days < t_b2_days) ) then + q_cgs = q_0_cgs * (t_days/t_b1_days)**exp_1 + else + q_cgs = q_0_cgs * (t_b2_days/t_b1_days)**exp_1 * (t_days/t_b2_days)**exp_2 + endif + + q = q_cgs / (unit_ergg/utime) + +end subroutine get_rprocess_heating_rate + !----------------------------------------------------------------------- ! ! cooling using Townsend (2009), ApJS 181, 391-397 method with @@ -601,7 +655,7 @@ subroutine write_options_cooling(iunit) call write_options_h2cooling(iunit) endif else - call write_inopt(icooling,'icooling','cooling function (0=off, 1=physics, 2=cooling table, 3=gammie)',iunit) + call write_inopt(icooling,'icooling','cooling function (0=off, 1=physics, 2=cooling table, 3=gammie, 4=rprocess heating rate (only works with GR Phantom for now))',iunit) select case(icooling) case(1) call write_inopt(icool_radiation_H0,'icool_radiation_H0','H0 cooling on/off',iunit) @@ -615,6 +669,14 @@ subroutine write_options_cooling(iunit) call write_inopt(temp_floor,'temp_floor','Minimum allowed temperature in K',iunit) case(3) call write_inopt(beta_cool,'beta_cool','beta factor in Gammie (2001) cooling',iunit) + !----- Input parameters for the rprocess heating rate function: + case(4) + call write_inopt(q_0_cgs,'q_0_cgs','heating rate coefficient [ergs/s/g]',iunit) + call write_inopt(t_b1_seconds,'t_b1_seconds','time of break 1 in heating rate [s]',iunit) + call write_inopt(exp_1,'exp_1','exponent 1 in heating rate',iunit) + call write_inopt(t_b2_seconds,'t_b2_seconds','time of break 2 in heating rate [s]',iunit) + call write_inopt(exp_2,'exp_2','exponent 2 in heating rate',iunit) + call write_inopt(t_start_seconds,'t_start_seconds','time after merger at which the simulation begins [s]',iunit) end select endif @@ -669,12 +731,33 @@ subroutine read_options_cooling(name,valstring,imatch,igotall,ierr) read(valstring,*,iostat=ierr) beta_cool ngot = ngot + 1 if (beta_cool < 1.) call fatal('read_options','beta_cool must be >= 1') + !----- Input parameters for the rprocess heating rate function: + case('q_0_cgs') + read(valstring,*,iostat=ierr) q_0_cgs + ngot = ngot + 1 + case('t_b1_seconds') + read(valstring,*,iostat=ierr) t_b1_seconds + ngot = ngot + 1 + case('exp_1') + read(valstring,*,iostat=ierr) exp_1 + ngot = ngot + 1 + case('t_b2_seconds') + read(valstring,*,iostat=ierr) t_b2_seconds + ngot = ngot + 1 + case('exp_2') + read(valstring,*,iostat=ierr) exp_2 + ngot = ngot + 1 + case('t_start_seconds') + read(valstring,*,iostat=ierr) t_start_seconds + ngot = ngot + 1 + !----- End of input parameters for rprocess heating rate function case default imatch = .false. if (h2chemistry) then call read_options_h2cooling(name,valstring,imatch,igotallh2,ierr) endif end select + if (icooling == 4 .and. ngot >= 6) igotall = .true. if (icooling == 3 .and. ngot >= 1) igotall = .true. if (icooling == 2 .and. ngot >= 3) igotall = .true. if (icooling == 1 .and. ngot >= 5) igotall = .true. diff --git a/src/main/force.F90 b/src/main/force.F90 index f65136c34..361b6dbdd 100644 --- a/src/main/force.F90 +++ b/src/main/force.F90 @@ -2438,8 +2438,6 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv use linklist, only:get_distance_from_centre_of_mass use kdtree, only:expand_fgrav_in_taylor_series use nicil, only:nimhd_get_dudt,nimhd_get_dt - use rprocess_heating, only:get_rprocess_heating_rate !----- This line is needed to call the function get_rprocess_heating_rate below - use timestep, only:time !----- This line is needed to call the function get_rprocess_heating_rate below use timestep, only:C_cour,C_cool,C_force,bignumber,dtmax use timestep_sts, only:use_sts use units, only:unit_ergg,unit_density @@ -2529,7 +2527,6 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv real :: densi, vxi,vyi,vzi,u0i real :: posi(3),veli(3),gcov(0:3,0:3),metrici(0:3,0:3,2) integer :: ii,ia,ib,ic,ierror - real :: q !----- This line is needed to call the function get_rprocess_heating_rate below eni = 0. realviscosity = (irealvisc > 0) @@ -2737,12 +2734,6 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv fxyz4 = fxyz4 + (gamma - 1.)*densi**(1.-gamma)*u0i*fsum(idendtdissi) endif - ! add r-process heating, if option selected - !if (gr) then - ! if (icooling==4) call get_rprocess_heating_rate(q,time) - ! fxyz4 = fxyz4 + q * (gamma - 1.) / densi**(gamma - 1.) - !endif - #ifdef GR #ifdef ISENTROPIC fxyz4 = 0. diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index cf8a3212d..2f54add58 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -25,7 +25,7 @@ module step_lf_global ! :Dependencies: chem, cons2prim, cons2primsolver, cooling, damping, deriv, ! derivutils, dim, dust_formation, eos, extern_gr, externalforces, ! growth, h2cooling, io, io_summary, krome_interface, metric_tools, -! mpiutils, options, part, ptmass, ptmass_radiation, rprocess_heating, timestep, +! mpiutils, options, part, ptmass, ptmass_radiation, timestep, ! timestep_ind, timestep_sts, timing ! use dim, only:maxp,maxvxyzu,do_radiation @@ -111,7 +111,6 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) use mpiutils, only:reduceall_mpi use part, only:nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,ibin_wake use io_summary, only:summary_printout,summary_variable,iosumtvi,iowake - use rprocess_heating, only:energ_rprocess #ifdef KROME use part, only:gamma_chem #endif @@ -512,16 +511,6 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) pxyzu(2,i) = pyi pxyzu(3,i) = pzi pxyzu(4,i) = eni - - !----- pxyzu(:,i) is the array (px,py,pz,K) for the ith SPH particle, where (px,py,pz) are the conserved coordinate momenta and K is the entropy variable in the rest frame. - !----- In Liptai & Price (2019), K is given the symbol K. ----- - !----- dens(:) is the array of rest mass densities of the SPH particles, where dens(i) is the rest mass density of the ith SPH particle - !----- The dens(:) array resides in the part module in the file part.F90, and is updated when the code calls the conservative to primitive variable transformation function - !----- The routines in GR Phantom use the following naming convention: - !----- The variable name 'dens' is the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho - !----- The variable name 'rho' is the "relativistic rest-mass density"/"conserved density". In Liptai & Price (2019), the "relativistic rest-mass density" is given the symbol rho^* - ! if (icooling==4) call energ_rprocess(pxyzu(4,i), rhoh(xyzh(4,i),massoftype(itype)), timei, dtsph) - if (icooling==4) call energ_rprocess(pxyzu(4,i), dens(i), timei, dtsph) else vxi = vxyzu(1,i) + hdtsph*fxyzu(1,i) vyi = vxyzu(2,i) + hdtsph*fxyzu(2,i) @@ -792,6 +781,7 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me !$omp shared(npart,xyzh,vxyzu,fext,iphase,ntypes,massoftype) & !$omp shared(maxphase,maxp) & !$omp shared(dt,hdt,xtol,ptol) & + !$omp shared(timei,dKdtcool) & !----- This line is needed to call energ_cooling below, which calls the rprocess heating function when icooling==4 !$omp shared(ieos,gamma,pxyzu,dens,metrics,metricderivs) & !$omp private(i,its,pondensi,spsoundi,rhoi,hi,eni,uui,densi) & !$omp private(converged,pmom_err,x_err,pri,ierr) & @@ -815,6 +805,14 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me ! ! make local copies of array quantities ! + !----- (Siva Darbha): Some notes for myself about the variable names: + !----- pxyzu(:,i) is the array (px,py,pz,K) for the ith SPH particle, where (px,py,pz) are the conserved coordinate momenta and K is the entropy variable in the rest frame. + !----- In Liptai & Price (2019), K is given the symbol K. ----- + !----- dens(:) is the array of rest mass densities of the SPH particles, where dens(i) is the rest mass density of the ith SPH particle + !----- The dens(:) array resides in the part module in the file part.F90, and is updated when the code calls the conservative to primitive variable transformation function + !----- The routines in GR Phantom use the following naming convention: + !----- The variable name 'dens' is the rest mass density. In Liptai & Price (2019), the rest mass density is given the symbol rho + !----- The variable name 'rho' is the "relativistic rest-mass density"/"conserved density". In Liptai & Price (2019), the "relativistic rest-mass density" is given the symbol rho^* pxyz(1:3) = pxyzu(1:3,i) eni = pxyzu(4,i) vxyz(1:3) = vxyzu(1:3,i) @@ -824,6 +822,26 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me pxyz = pxyz + hdt*fexti + !----- (Siva Darbha): I call the heating/cooling function in the predictor step, in analogy with nonrelativistic SPH + if (maxvxyzu >= 4) then + dudtcool = 0. + ! + ! COOLING + ! + !----- Calls the rprocess heating rate + if (icooling == 4) then + + !----- The rprocess heating rate only requires the arguments dudtcool and timei, so I pass zeros to the rest + call energ_cooling(0, 0, 0, 0, dudtcool, 0, 0, 0, 0, 0, 0, timei) + !call energ_cooling(xyzh(1,i), xyzh(2,i), xyzh(3,i), vxyzu(4,i), dudtcool, rhoh(xyzh(4,i), pmassi), dt, 0, 0, 0, 0, timei) + + !----- Updates the entropy variable + dKdtcool = dudtcool * (gamma - 1) / densi**(gamma - 1) + eni = eni + dt * dKdtcool !----- (Siva Darbha): I used dt here, but I might have to use hdt, I'm not sure + + endif + endif + !-- Compute pressure for the first guess in cons2prim call equationofstate(ieos,pondensi,spsoundi,densi,xyz(1),xyz(2),xyz(3),uui) pri = pondensi*densi @@ -879,6 +897,7 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me ! re-pack arrays back where they belong xyzh(1:3,i) = xyz(1:3) pxyzu(1:3,i) = pxyz(1:3) + pxyzu(4,i) = eni !----- This line is needed if you call energ_cooling above, which calls the rprocess heating function when icooling==4, which updates the entropy variable vxyzu(1:3,i) = vxyz(1:3) vxyzu(4,i) = uui fext(:,i) = fexti From 296ae1eab49266067a04d7f99f9a61d95c3e616b Mon Sep 17 00:00:00 2001 From: sdarbha Date: Wed, 2 Sep 2020 11:30:57 -0700 Subject: [PATCH 09/15] Fixed mistakes in updated rprocess heating treatment --- src/main/cooling.F90 | 10 ++++++---- src/main/step_leapfrog.F90 | 23 ++++++++++++----------- 2 files changed, 18 insertions(+), 15 deletions(-) diff --git a/src/main/cooling.F90 b/src/main/cooling.F90 index 35d858035..24361b277 100644 --- a/src/main/cooling.F90 +++ b/src/main/cooling.F90 @@ -433,18 +433,18 @@ end subroutine energ_cooling !----------------------------------------------------------------------- ! -! This function returns the rprocess specific heating rate q at time t, both in code units +! This function adds the rprocess specific heating rate q to du/dt at time t, all in code units ! !----------------------------------------------------------------------- -subroutine get_rprocess_heating_rate(q,t) +subroutine get_rprocess_heating_rate(dudt,t) use physcon, only:days use units, only:unit_ergg,utime real, intent(in) :: t !----- time, in code units - real, intent(out) :: q !----- specific heating rate, in code units + real, intent(out) :: dudt !----- specific heating/cooling rate, in code units real :: t_seconds, t_days real :: t_b1_days, t_b2_days real :: t_start_days, t_new_days - real :: q_cgs + real :: q, q_cgs t_seconds = t * utime t_days = t_seconds / days @@ -466,6 +466,8 @@ subroutine get_rprocess_heating_rate(q,t) q = q_cgs / (unit_ergg/utime) + dudt = dudt + q + end subroutine get_rprocess_heating_rate !----------------------------------------------------------------------- diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 2f54add58..9f6fb2714 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -233,7 +233,8 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) !---------------------------------------------------------------------- call get_timings(t1,tcpu1) #ifdef GR - if ((iexternalforce > 0 .and. imetric /= imet_minkowski) .or. idamp > 0) then + ! This now includes a condition on icooling. If icooling == 4, then step_extern_gr calls the function energ_cooling, which calls the rprocess heating function + if ((iexternalforce > 0 .and. imetric /= imet_minkowski) .or. idamp > 0 .or. icooling == 4) then call cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars) call get_grforce_all(npart,xyzh,metrics,metricderivs,vxyzu,dens,fext,dtextforce) call step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,metrics,metricderivs,fext,t) @@ -707,7 +708,7 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me use dim, only:maxptmass,maxp,maxvxyzu use io, only:iverbose,id,master,iprint,warning use externalforces, only:externalforce,accrete_particles,update_externalforce - use options, only:iexternalforce,idamp + use options, only:iexternalforce,idamp,icooling !----- icooling is needed on this line to call energ_cooling below use part, only:maxphase,isdead_or_accreted,iamboundary,igas,iphase,iamtype,massoftype,rhoh use io_summary, only:summary_variable,iosumextsr,iosumextst,iosumexter,iosumextet,iosumextr,iosumextt, & summary_accrete,summary_accrete_fail @@ -717,12 +718,14 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me use extern_gr, only:get_grforce use metric_tools, only:pack_metric,pack_metricderivs use damping, only:calc_damp,apply_damp + use cooling, only:energ_cooling !----- this line is needed to call energ_cooling below integer, intent(in) :: npart,ntypes real, intent(in) :: dtsph,time real, intent(inout) :: dtextforce real, intent(inout) :: xyzh(:,:),vxyzu(:,:),fext(:,:),pxyzu(:,:),dens(:),metrics(:,:,:,:),metricderivs(:,:,:,:) integer :: i,itype,nsubsteps,naccreted,its,ierr real :: timei,t_end_step,hdt,pmassi + real :: dudtcool,dKdtcool !----- these variables are needed to use energ_cooling below real :: dt,dtf,dtextforcenew,dtextforce_min real :: pri,spsoundi,pondensi real, save :: pprev(3),xyz_prev(3),fstar(3),vxyz_star(3),xyz(3),pxyz(3),vxyz(3),fexti(3) @@ -781,8 +784,9 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me !$omp shared(npart,xyzh,vxyzu,fext,iphase,ntypes,massoftype) & !$omp shared(maxphase,maxp) & !$omp shared(dt,hdt,xtol,ptol) & - !$omp shared(timei,dKdtcool) & !----- This line is needed to call energ_cooling below, which calls the rprocess heating function when icooling==4 + !$omp shared(icooling,timei) & !----- This line is needed to call energ_cooling below, which calls the rprocess heating function when icooling==4 !$omp shared(ieos,gamma,pxyzu,dens,metrics,metricderivs) & + !$omp private(dudtcool,dKdtcool) & !----- This line is needed to call energ_cooling below, which calls the rprocess heating function when icooling==4 !$omp private(i,its,pondensi,spsoundi,rhoi,hi,eni,uui,densi) & !$omp private(converged,pmom_err,x_err,pri,ierr) & !$omp firstprivate(pmassi,itype) & @@ -830,16 +834,13 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me ! !----- Calls the rprocess heating rate if (icooling == 4) then - - !----- The rprocess heating rate only requires the arguments dudtcool and timei, so I pass zeros to the rest - call energ_cooling(0, 0, 0, 0, dudtcool, 0, 0, 0, 0, 0, 0, timei) + !----- (Siva Darbha): The rprocess heating rate only requires the arguments dudtcool and timei, so I pass zeros to the rest + call energ_cooling(0., 0., 0., 0., dudtcool, 0., 0., 0., 0., 0., 0., timei) !call energ_cooling(xyzh(1,i), xyzh(2,i), xyzh(3,i), vxyzu(4,i), dudtcool, rhoh(xyzh(4,i), pmassi), dt, 0, 0, 0, 0, timei) - - !----- Updates the entropy variable - dKdtcool = dudtcool * (gamma - 1) / densi**(gamma - 1) - eni = eni + dt * dKdtcool !----- (Siva Darbha): I used dt here, but I might have to use hdt, I'm not sure - endif + !----- Updates the entropy variable + dKdtcool = dudtcool * (gamma - 1) / densi**(gamma - 1) + eni = eni + dt * dKdtcool !----- (Siva Darbha): I used dt here, but I might have to use hdt, I'm not sure endif !-- Compute pressure for the first guess in cons2prim From c8d53ab798e3ad721efc13aaea0aa056b7b54a7c Mon Sep 17 00:00:00 2001 From: sdarbha Date: Thu, 24 Sep 2020 11:09:51 -0700 Subject: [PATCH 10/15] Added C_cool as input parameter for rprocess heating option --- src/main/cooling.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/main/cooling.F90 b/src/main/cooling.F90 index 24361b277..0410f43d8 100644 --- a/src/main/cooling.F90 +++ b/src/main/cooling.F90 @@ -673,6 +673,7 @@ subroutine write_options_cooling(iunit) call write_inopt(beta_cool,'beta_cool','beta factor in Gammie (2001) cooling',iunit) !----- Input parameters for the rprocess heating rate function: case(4) + call write_inopt(C_cool,'C_cool','factor controlling cooling timestep',iunit) call write_inopt(q_0_cgs,'q_0_cgs','heating rate coefficient [ergs/s/g]',iunit) call write_inopt(t_b1_seconds,'t_b1_seconds','time of break 1 in heating rate [s]',iunit) call write_inopt(exp_1,'exp_1','exponent 1 in heating rate',iunit) @@ -759,7 +760,7 @@ subroutine read_options_cooling(name,valstring,imatch,igotall,ierr) call read_options_h2cooling(name,valstring,imatch,igotallh2,ierr) endif end select - if (icooling == 4 .and. ngot >= 6) igotall = .true. + if (icooling == 4 .and. ngot >= 7) igotall = .true. if (icooling == 3 .and. ngot >= 1) igotall = .true. if (icooling == 2 .and. ngot >= 3) igotall = .true. if (icooling == 1 .and. ngot >= 5) igotall = .true. From 211b57f1ec15975df35cef17b03003bc1ff38204 Mon Sep 17 00:00:00 2001 From: sdarbha Date: Mon, 7 Dec 2020 14:04:38 -0800 Subject: [PATCH 11/15] Removed vestigial file rprocess_heating.f90 --- build/Makefile | 2 +- src/main/rprocess_heating.f90 | 179 ---------------------------------- 2 files changed, 1 insertion(+), 180 deletions(-) delete mode 100644 src/main/rprocess_heating.f90 diff --git a/build/Makefile b/build/Makefile index f9b741943..8e93e04c2 100644 --- a/build/Makefile +++ b/build/Makefile @@ -1767,7 +1767,7 @@ else endif SRCGR=inverse4x4.f90 $(SRCMETRIC) metric_tools.f90 utils_gr.f90 cons2primsolver.f90 -SRCCHEM= rprocess_heating.f90 fs_data.f90 mol_data.f90 utils_spline.f90 h2cooling.f90 h2chem.f90 cooling.f90 +SRCCHEM= fs_data.f90 mol_data.f90 utils_spline.f90 h2cooling.f90 h2chem.f90 cooling.f90 SRCMESA= eos_mesa_microphysics.f90 eos_mesa.f90 SRCEOS = ${SRCMESA} eos_shen.f90 eos_helmholtz.f90 eos_idealplusrad.f90 eos.F90 diff --git a/src/main/rprocess_heating.f90 b/src/main/rprocess_heating.f90 deleted file mode 100644 index b2ff1465e..000000000 --- a/src/main/rprocess_heating.f90 +++ /dev/null @@ -1,179 +0,0 @@ -!--------------------------------------------------------------------------! -! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! -! Copyright (c) 2007-2020 The Authors (see AUTHORS) ! -! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! -!--------------------------------------------------------------------------! -!+ -! MODULE: rprocess_heating -! -! DESCRIPTION: -! Implements r-process heating. -! -! OWNER: Siva Darbha -! -! $Id$ -! -! RUNTIME PARAMETERS: -! -! q_0_cgs -- heating rate coefficient, in [ergs/s/g] -! t_b1_seconds -- time of break 1 in heating rate, in [s]' -! exp_1 -- exponent 1 in radioactive power component -! t_b2_seconds -- time of break 2 in heating rate, in [s]' -! exp_2 -- exponent 2 in radioactive power component -! t_start_seconds -- time after merger at which the simulation begins [s] -! -! DEPENDENCIES: infile_utils, io, physcon, units -!+ -!-------------------------------------------------------------------------- -module rprocess_heating - implicit none - integer, parameter :: maxt = 1000 - real :: temper(maxt),lambda(maxt),slope(maxt),yfunc(maxt) - integer :: nt - ! - ! set default values for input parameters - ! - real :: q_0_cgs = 0 ! heating rate coefficient [ergs/s/g] - real :: t_b1_seconds = 0 ! time of break 1 in heating rate [s] - real :: exp_1 = 0 ! exponent 1 in heating rate - real :: t_b2_seconds = 0 ! time of break 2 in heating rate [s] - real :: exp_2 = 0 ! exponent 2 in heating rate - real :: t_start_seconds = 0 ! time after merger at which the simulation begins [s] - - public :: write_options_rprocess,read_options_rprocess - public :: get_rprocess_heating_rate, energ_rprocess - - private - -contains - -!----------------------------------------------------------------------- -!+ -! get the r-process specific heating rate q at time t -!+ -!----------------------------------------------------------------------- -subroutine get_rprocess_heating_rate(q,t) - use physcon, only:days - use units, only:unit_ergg,utime - real, intent(in) :: t - real, intent(out) :: q - real :: t_seconds, t_days - real :: t_b1_days, t_b2_days - real :: t_start_days, t_new_days - real :: q_cgs - - t_seconds = t * utime - t_days = t_seconds / days - - t_b1_days = t_b1_seconds / days - t_b2_days = t_b2_seconds / days - - t_start_days = t_start_seconds / days - t_new_days = t_days + t_start_days - - !----- Heating rate in [ergs/s/g] - if (t_new_days < t_b1_days) then - q_cgs = q_0_cgs - else if ( (t_b1_days <= t_new_days) .and. (t_new_days < t_b2_days) ) then - q_cgs = q_0_cgs * (t_days/t_b1_days)**exp_1 - else - q_cgs = q_0_cgs * (t_b2_days/t_b1_days)**exp_1 * (t_days/t_b2_days)**exp_2 - endif - - q = q_cgs / (unit_ergg/utime) - -end subroutine get_rprocess_heating_rate - -!----------------------------------------------------------------------- -!+ -! get the change in the entropy variable K in the timestep t -> t+dt -!+ -!----------------------------------------------------------------------- -subroutine energ_rprocess(entrop_var,dens,t,dt) - use eos, only:gamma - use physcon, only:days - use units, only:unit_ergg,utime - real, intent(in) :: dens,t,dt - real, intent(inout) :: entrop_var - real :: dt_seconds - real :: q, q_cgs - real :: delta_u, delta_u_cgs - real :: delta_entrop_var - - dt_seconds = dt * utime - - call get_rprocess_heating_rate(q, t) - - q_cgs = q * (unit_ergg/utime) - - !----- Change in the specific internal energy [ergs/g] during the time step - delta_u_cgs = q_cgs * dt_seconds - delta_u = delta_u_cgs / unit_ergg - - !----- Change in the entropy variable K during the time step - !----- In Liptai & Price (2019), K is given the symbol K. ----- - delta_entrop_var = (gamma - 1.) * delta_u / dens**(gamma - 1.) - entrop_var = entrop_var + delta_entrop_var - -end subroutine energ_rprocess - -!----------------------------------------------------------------------- -!+ -! writes input options to the input file -!+ -!----------------------------------------------------------------------- -subroutine write_options_rprocess(iunit) - use infile_utils, only:write_inopt - integer, intent(in) :: iunit - - call write_inopt(q_0_cgs,'q_0_cgs','heating rate coefficient [ergs/s/g]',iunit) - call write_inopt(t_b1_seconds,'t_b1_seconds','time of break 1 in heating rate [s]',iunit) - call write_inopt(exp_1,'exp_1','exponent 1 in heating rate',iunit) - call write_inopt(t_b2_seconds,'t_b2_seconds','time of break 2 in heating rate [s]',iunit) - call write_inopt(exp_2,'exp_2','exponent 2 in heating rate',iunit) - call write_inopt(t_start_seconds,'t_start_seconds','time after merger at which the simulation begins [s]',iunit) - -end subroutine write_options_rprocess - -!----------------------------------------------------------------------- -!+ -! reads input options from the input file -!+ -!----------------------------------------------------------------------- -subroutine read_options_rprocess(name,valstring,imatch,igotall,ierr) - use io, only:fatal - character(len=*), intent(in) :: name,valstring - logical, intent(out) :: imatch,igotall - integer, intent(out) :: ierr - integer, save :: ngot = 0 - - imatch = .true. - igotall = .false. ! cooling options are compulsory - select case(trim(name)) - case('q_0_cgs') - read(valstring,*,iostat=ierr) q_0_cgs - ngot = ngot + 1 - case('t_b1_seconds') - read(valstring,*,iostat=ierr) t_b1_seconds - ngot = ngot + 1 - case('exp_1') - read(valstring,*,iostat=ierr) exp_1 - ngot = ngot + 1 - case('t_b2_seconds') - read(valstring,*,iostat=ierr) t_b2_seconds - ngot = ngot + 1 - case('exp_2') - read(valstring,*,iostat=ierr) exp_2 - ngot = ngot + 1 - case('t_start_seconds') - read(valstring,*,iostat=ierr) t_start_seconds - ngot = ngot + 1 - case default - imatch = .false. - end select - if (ngot >= 6) igotall = .true. - -end subroutine read_options_rprocess - -end module rprocess_heating \ No newline at end of file From 8b56f63f143dc9ed142ee6caaf80ee51308ae344 Mon Sep 17 00:00:00 2001 From: sdarbha Date: Mon, 7 Dec 2020 14:21:20 -0800 Subject: [PATCH 12/15] Updated formatting in analysis_bhnsejecta.f90 --- src/utils/analysis_bhnsejecta.f90 | 64 +++++++++++++++---------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/src/utils/analysis_bhnsejecta.f90 b/src/utils/analysis_bhnsejecta.f90 index 231ef2116..d0c8f0ec0 100644 --- a/src/utils/analysis_bhnsejecta.f90 +++ b/src/utils/analysis_bhnsejecta.f90 @@ -53,21 +53,21 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) particle_file_name = trim(dumpfile)//'.dat' - open(1, file=particle_file_name, action='write', status='replace') - - write(1,'(A)',IOSTAT=iostatus) '# Particle data' - write(1,'(A)',IOSTAT=iostatus) '# [1] = Mass' - write(1,'(A)',IOSTAT=iostatus) '# [2] = Radius' - write(1,'(A)',IOSTAT=iostatus) '# [3] = Theta' - write(1,'(A)',IOSTAT=iostatus) '# [4] = Phi' - write(1,'(A)',IOSTAT=iostatus) '# [5] = e (-u_t)' - write(1,'(A)',IOSTAT=iostatus) '# [6] = l (u_phi)' - write(1,'(A)',IOSTAT=iostatus) '# [7] = u_r' - write(1,'(A)',IOSTAT=iostatus) '# [8] = u_theta' - write(1,'(A)',IOSTAT=iostatus) '# [9] = Y_e' - write(1,'(A)',IOSTAT=iostatus) '# [10]= s' - write(1,'(A)',IOSTAT=iostatus) '# [11]= T' - write(1,'(A)',IOSTAT=iostatus) '# [12]= rho' + open(newunit=iunit, file=particle_file_name, action='write', status='replace') + + write(iunit,'(a)',iostat=iostatus) '# Particle data' + write(iunit,'(a)',iostat=iostatus) '# [1] = Mass' + write(iunit,'(a)',iostat=iostatus) '# [2] = Radius' + write(iunit,'(a)',iostat=iostatus) '# [3] = Theta' + write(iunit,'(a)',iostat=iostatus) '# [4] = Phi' + write(iunit,'(a)',iostat=iostatus) '# [5] = e (-u_t)' + write(iunit,'(a)',iostat=iostatus) '# [6] = l (u_phi)' + write(iunit,'(a)',iostat=iostatus) '# [7] = u_r' + write(iunit,'(a)',iostat=iostatus) '# [8] = u_theta' + write(iunit,'(a)',iostat=iostatus) '# [9] = Y_e' + write(iunit,'(a)',iostat=iostatus) '# [10]= s' + write(iunit,'(a)',iostat=iostatus) '# [11]= T' + write(iunit,'(a)',iostat=iostatus) '# [12]= rho' call init_metric(npart,xyzh,metrics,metricderivs) @@ -140,22 +140,22 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) - write(1,'(ES11.4,1X)',advance='no',IOSTAT=iostatus) m - write(1,'(ES12.5,1X)',advance='no',IOSTAT=iostatus) r - write(1,'(F9.5,1X)',advance='no',IOSTAT=iostatus) theta - write(1,'(F9.5,1X)',advance='no',IOSTAT=iostatus) phi - write(1,'(F10.5,1X)',advance='no',IOSTAT=iostatus) e - write(1,'(ES12.5,1X)',advance='no',IOSTAT=iostatus) l - write(1,'(ES12.5,1X)',advance='no',IOSTAT=iostatus) u_r - write(1,'(ES12.5,1X)',advance='no',IOSTAT=iostatus) u_theta - write(1,'(F10.7,1X)',advance='no',IOSTAT=iostatus) y_e - write(1,'(F10.7,1X)',advance='no',IOSTAT=iostatus) s - write(1,'(ES12.4,1X)',advance='no',IOSTAT=iostatus) T - write(1,'(ES12.5)',IOSTAT=iostatus) densi + write(iunit,'(es11.4,1x)',advance='no',iostat=iostatus) m + write(iunit,'(es12.5,1x)',advance='no',iostat=iostatus) r + write(iunit,'(f9.5,1x)',advance='no',iostat=iostatus) theta + write(iunit,'(f9.5,1x)',advance='no',iostat=iostatus) phi + write(iunit,'(f10.5,1x)',advance='no',iostat=iostatus) e + write(iunit,'(es12.5,1x)',advance='no',iostat=iostatus) l + write(iunit,'(es12.5,1x)',advance='no',iostat=iostatus) u_r + write(iunit,'(es12.5,1x)',advance='no',iostat=iostatus) u_theta + write(iunit,'(f10.7,1x)',advance='no',iostat=iostatus) y_e + write(iunit,'(f10.7,1x)',advance='no',iostat=iostatus) s + write(iunit,'(es12.4,1x)',advance='no',iostat=iostatus) T + write(iunit,'(es12.5)',iostat=iostatus) densi enddo - close(1) + close(iunit) @@ -166,12 +166,12 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) info_file_name = trim(dumpfile)//'.info' - open(2, file=info_file_name, action='write', status='replace') + open(newunit=iunit, file=info_file_name, action='write', status='replace') - write(2,'(A,1X,ES10.3,1X,A)',IOSTAT=iostatus) 'time =', time_output, '[G = c = M_Sun = 1]' - write(2,'(A,1X,ES10.3,1X,A)',IOSTAT=iostatus) 'time_seconds =', time_output_seconds, '[s]' + write(iunit,'(a,1x,es10.3,1x,a)',iostat=iostatus) 'time =', time_output, '[G = c = M_Sun = 1]' + write(iunit,'(a,1x,es10.3,1x,a)',iostat=iostatus) 'time_seconds =', time_output_seconds, '[s]' - close(2) + close(iunit) From 3694759efcf6bd29ddf2b23c631de306d491418b Mon Sep 17 00:00:00 2001 From: sdarbha Date: Mon, 7 Dec 2020 14:30:03 -0800 Subject: [PATCH 13/15] Fixed mistake in formatting in analysis_bhnsejecta.f90 --- src/utils/analysis_bhnsejecta.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utils/analysis_bhnsejecta.f90 b/src/utils/analysis_bhnsejecta.f90 index d0c8f0ec0..47dcefb0c 100644 --- a/src/utils/analysis_bhnsejecta.f90 +++ b/src/utils/analysis_bhnsejecta.f90 @@ -53,7 +53,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) particle_file_name = trim(dumpfile)//'.dat' - open(newunit=iunit, file=particle_file_name, action='write', status='replace') + open(iunit, file=particle_file_name, action='write', status='replace') write(iunit,'(a)',iostat=iostatus) '# Particle data' write(iunit,'(a)',iostat=iostatus) '# [1] = Mass' @@ -166,7 +166,7 @@ subroutine do_analysis(dumpfile,numfile,xyzh,vxyzu,pmass,npart,time,iunit) info_file_name = trim(dumpfile)//'.info' - open(newunit=iunit, file=info_file_name, action='write', status='replace') + open(iunit, file=info_file_name, action='write', status='replace') write(iunit,'(a,1x,es10.3,1x,a)',iostat=iostatus) 'time =', time_output, '[G = c = M_Sun = 1]' write(iunit,'(a,1x,es10.3,1x,a)',iostat=iostatus) 'time_seconds =', time_output_seconds, '[s]' From 0c45ae8e5bdad6f676f9e695dddb10a3976b61ec Mon Sep 17 00:00:00 2001 From: sdarbha Date: Mon, 7 Dec 2020 14:36:39 -0800 Subject: [PATCH 14/15] Updated formatting of setup_grbhnsejecta.f90 --- src/setup/setup_grbhnsejecta.f90 | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/setup/setup_grbhnsejecta.f90 b/src/setup/setup_grbhnsejecta.f90 index 43862e018..316e38f93 100644 --- a/src/setup/setup_grbhnsejecta.f90 +++ b/src/setup/setup_grbhnsejecta.f90 @@ -272,20 +272,21 @@ end subroutine read_setupfile ! subroutine get_num_particles(npart) integer :: iostatus, i + integer, parameter :: iunit = 22 integer, intent(out) :: npart real :: m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, densi print *, 'getting number of particles from '//trim(particle_file_name) - open(1, file=particle_file_name, status='old', action='read') + open(iunit, file=particle_file_name, status='old', action='read') do i=1,13 - read(1,*) + read(iunit,*) enddo i = 0 do - read(1,*,IOSTAT=iostatus) m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, densi + read(iunit,*,iostat=iostatus) m, r, theta, phi, e, l, u_r, u_theta, y_e, s, T, densi ! read failed if (iostatus > 0) then print *, 'read failed, input data not in correct format' @@ -300,7 +301,7 @@ subroutine get_num_particles(npart) enddo npart = i - close(1) + close(iunit) end subroutine get_num_particles @@ -309,20 +310,21 @@ end subroutine get_num_particles ! subroutine read_particle_data(npart, m_array, r_array, theta_array, phi_array, e_array, l_array, u_r_array, u_theta_array, y_e_array, s_array, T_array, dens_array) integer :: iostatus, i + integer, parameter :: iunit = 23 integer, intent(in) :: npart real, intent(out) :: m_array(:), r_array(:), theta_array(:), phi_array(:), e_array(:), l_array(:), u_r_array(:), u_theta_array(:), y_e_array(:), s_array(:), T_array(:), dens_array(:) print *, 'reading particle data from '//trim(particle_file_name) - open(1, file=particle_file_name, status='old', action='read') + open(iunit, file=particle_file_name, status='old', action='read') do i=1,13 - read(1,*) + read(iunit,*) enddo i = 1 do - read(1,*,IOSTAT=iostatus) m_array(i), r_array(i), theta_array(i), phi_array(i), e_array(i), l_array(i), u_r_array(i), u_theta_array(i), y_e_array(i), s_array(i), T_array(i), dens_array(i) + read(iunit,*,iostat=iostatus) m_array(i), r_array(i), theta_array(i), phi_array(i), e_array(i), l_array(i), u_r_array(i), u_theta_array(i), y_e_array(i), s_array(i), T_array(i), dens_array(i) ! read failed if (iostatus > 0) then print *, 'read failed, input data not in correct format' @@ -336,7 +338,7 @@ subroutine read_particle_data(npart, m_array, r_array, theta_array, phi_array, e endif enddo - close(1) + close(iunit) end subroutine read_particle_data From 0a75a3b1c8170fe67cbb99b7d4d0af4949354ed8 Mon Sep 17 00:00:00 2001 From: sdarbha Date: Mon, 7 Dec 2020 16:12:20 -0800 Subject: [PATCH 15/15] Updated Makefile for SYSTEM=cori and SETUP=grbhnsejecta --- build/Makefile | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/build/Makefile b/build/Makefile index 8e93e04c2..dd066168e 100644 --- a/build/Makefile +++ b/build/Makefile @@ -276,17 +276,13 @@ endif ifeq ($(SETUP), grbhnsejecta) SETUPFILE=setup_grbhnsejecta.f90 ANALYSIS=analysis_bhnsejecta.f90 -# SRCINJECT=inject_bhnsejecta_particles.f90 GR=yes METRIC=schwarzschild GRAVITY=no -# GRAVITY=yes ISOTHERMAL=no IND_TIMESTEPS=no KNOWN_SETUP=yes KERNEL=cubic -# KERNEL=quintic -# KERNEL=WendlandC4 endif ifeq ($(SETUP), srpolytrope) @@ -1247,21 +1243,16 @@ endif # Most importantly, I (SD) made the replacements ifort -> ftn and icc -> cc, and added CRAYPE_LINK_TYPE to tell Cori how to link ifeq ($(SYSTEM), cori) #Cori (NERSC machine) + include Makefile_defaults_ifort FC= ftn # FFLAGS= -O3 -inline-factor=500 -mcmodel=medium -shared-intel -warn uninitialized -warn unused -warn truncated_source FFLAGS= -O3 -inline-factor=500 -shared-intel -warn uninitialized -warn unused -warn truncated_source - DBLFLAG= -r8 - DEBUGFLAG= -check all -WB -traceback -g -debug -all - ENDIANFLAGBIG= -convert big_endian - ENDIANFLAGLITTLE= -convert little_endian CC= cc # CCFLAGS= -O3 -mcmodel=medium CCFLAGS= -O3 - LIBCXX = -cxxlib QSYS= slurm OMPFLAGS= -qopenmp NOMP=64 - KNOWN_SYSTEM=yes # CRAYPE_LINK_TYPE=dynamic CRAYPE_LINK_TYPE=static MAXP=2000000