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

build/Makefile

@@ -272,6 +272,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
@@ -539,6 +556,7 @@ ifeq ($(SETUP), srshock)
     CONST_AV=yes
 endif
 
+
 ifeq ($(SETUP), testparticles)
 #   test particles
     PERIODIC=no
@@ -1223,6 +1241,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
@@ -1233,6 +1278,8 @@ ifeq ($(SYSTEM), ozstar)
     WALLTIME='168:00:00'
 endif
 
+
+
 ifeq ($(SYSTEM), monarch)
     include Makefile_defaults_ifort
     OMPFLAGS=-qopenmp -qopt-report
@@ -1720,7 +1767,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= 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_idealplusrad.f90 eos.F90
@@ -2941,7 +2988,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.