diff --git a/build/Makefile b/build/Makefile index c552669a5..dd066168e 100644 --- a/build/Makefile +++ b/build/Makefile @@ -272,6 +272,19 @@ 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 + GR=yes + METRIC=schwarzschild + GRAVITY=no + ISOTHERMAL=no + IND_TIMESTEPS=no + KNOWN_SETUP=yes + KERNEL=cubic +endif + ifeq ($(SETUP), srpolytrope) FPPFLAGS= -DPRIM2CONS_FIRST SETUPFILE= setup_srpolytrope.f90 @@ -539,6 +552,7 @@ ifeq ($(SETUP), srshock) CONST_AV=yes endif + ifeq ($(SETUP), testparticles) # test particles PERIODIC=no @@ -1223,6 +1237,28 @@ 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) + 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 + CC= cc +# CCFLAGS= -O3 -mcmodel=medium + CCFLAGS= -O3 + QSYS= slurm + OMPFLAGS= -qopenmp + NOMP=64 +# CRAYPE_LINK_TYPE=dynamic + CRAYPE_LINK_TYPE=static + MAXP=2000000 +# MAXP=8000000 +endif + ifeq ($(SYSTEM), ozstar) # ozstar facility include Makefile_defaults_ifort @@ -1233,6 +1269,8 @@ ifeq ($(SYSTEM), ozstar) WALLTIME='168:00:00' endif + + ifeq ($(SYSTEM), monarch) include Makefile_defaults_ifort OMPFLAGS=-qopenmp -qopt-report @@ -1720,7 +1758,7 @@ else endif SRCGR=inverse4x4.f90 $(SRCMETRIC) metric_tools.f90 utils_gr.f90 cons2primsolver.f90 -SRCCHEM= 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 @@ -2941,7 +2979,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/main/checksetup.F90 b/src/main/checksetup.F90 index 6f0570b90..9141ae9c7 100644 --- a/src/main/checksetup.F90 +++ b/src/main/checksetup.F90 @@ -440,7 +440,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,')' - if (.not.h2chemistry .and. maxvxyzu >= 4 .and. icooling >= 1 .and. iexternalforce/=iext_corotate) then +! (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 == 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..0410f43d8 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,45 @@ subroutine energ_cooling (xi, yi, zi, u, dudt, rho, dt, Trad, mu_in, K2, kappa) end subroutine energ_cooling +!----------------------------------------------------------------------- +! +! This function adds the rprocess specific heating rate q to du/dt at time t, all in code units +! +!----------------------------------------------------------------------- +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) :: 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, 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) + + dudt = dudt + q + +end subroutine get_rprocess_heating_rate + !----------------------------------------------------------------------- ! ! cooling using Townsend (2009), ApJS 181, 391-397 method with @@ -601,7 +657,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 +671,15 @@ 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(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) + 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 +734,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 >= 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. diff --git a/src/main/force.F90 b/src/main/force.F90 index b5f5428eb..8d9aabe56 100644 --- a/src/main/force.F90 +++ b/src/main/force.F90 @@ -2742,6 +2742,7 @@ 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 + #ifdef GR #ifdef ISENTROPIC fxyz4 = 0. diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 4db10c6df..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) @@ -405,6 +406,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,ieos) & +!$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) & @@ -706,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 @@ -716,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) @@ -780,7 +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(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) & @@ -803,6 +809,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) @@ -812,6 +826,23 @@ 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 + !----- (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) + 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 call equationofstate(ieos,pondensi,spsoundi,densi,xyz(1),xyz(2),xyz(3),uui) pri = pondensi*densi @@ -867,6 +898,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 diff --git a/src/setup/setup_grbhnsejecta.f90 b/src/setup/setup_grbhnsejecta.f90 new file mode 100644 index 000000000..316e38f93 --- /dev/null +++ b/src/setup/setup_grbhnsejecta.f90 @@ -0,0 +1,345 @@ +!--------------------------------------------------------------------------! +! 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, timestep, units +!+ +!-------------------------------------------------------------------------- +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(:), 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, densi + 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(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, dens_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 + dens_array(i) = dens_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 !----- 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) + !----- 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", 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. + 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.) / (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 + + 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(dens_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, 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(iunit, file=particle_file_name, status='old', action='read') + + do i=1,13 + read(iunit,*) + enddo + + i = 0 + do + 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' + exit + ! reached end of file + else if (iostatus < 0) then + exit + ! read successful + else + i = i + 1 + endif + enddo + npart = i + + close(iunit) + +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, 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(iunit, file=particle_file_name, status='old', action='read') + + do i=1,13 + read(iunit,*) + enddo + + i = 1 + do + 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' + exit + ! reached end of file + else if (iostatus < 0) then + exit + ! read successful + else + i = i + 1 + endif + enddo + + close(iunit) + +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..47dcefb0c --- /dev/null +++ b/src/utils/analysis_bhnsejecta.f90 @@ -0,0 +1,180 @@ +!--------------------------------------------------------------------------! +! 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,pxyzu,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(:,:) + 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, u + real :: m, densi + real :: T, s, y_e + real :: q + real :: time_output, time_output_seconds + character(len=120) :: particle_file_name, info_file_name + + + + particle_file_name = trim(dumpfile)//'.dat' + + open(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) + + 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) !----- 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 + 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.) + 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.)) + + !----- 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)) + + 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 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 + densi = densi * unit_density / ( solarm / (gg*solarm/(c**2.))**3. ) + + + + 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(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(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]' + + close(iunit) + + + +end subroutine do_analysis + +end module \ No newline at end of file