Implemented GR BH-NS ejecta setup and r-process heating function

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

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!'

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.

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.

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