Skip to content
Permalink
Browse files

Merge pull request #28 from srichers/master

 fixed indexing bug in epannihil kernel calculation when using neutrino scheme 2 (separate heavies/anti-heavies)
  • Loading branch information...
evanoconnor committed Apr 6, 2019
2 parents f35bca3 + bb9fb5d commit c6a57caf0dd9f37c2bda1c73479f4192bee24c3e
Showing with 23 additions and 21 deletions.
  1. +19 −17 src/make_table_example.F90
  2. +4 −4 src/nulib.F90
@@ -273,7 +273,6 @@ program make_table_example
write(*,*) "Rho:", 100.0*dble(irho-1)/dble(final_table_size_rho),"%"
#endif
do itemp=1,final_table_size_temp
write(*,*) "Temp:", 100.0*dble(itemp-1)/dble(final_table_size_temp),"%"
do iye=1,final_table_size_ye

eos_variables = 0.0d0
@@ -311,10 +310,12 @@ program make_table_example
eos_variables(rhoindex),eos_variables(tempindex),eos_variables(yeindex),ns,ng
stop
endif
if (local_delta(ns,ng).ne.local_delta(ns,ng)) then
write(*,"(a,1P3E18.9,i6,i6)") "We have a NaN in scat delta", &
eos_variables(rhoindex),eos_variables(tempindex),eos_variables(yeindex),ns,ng
stop
if (.not. do_transport_opacities) then
if (local_delta(ns,ng).ne.local_delta(ns,ng)) then
write(*,"(a,1P3E18.9,i6,i6)") "We have a NaN in scat delta", &
eos_variables(rhoindex),eos_variables(tempindex),eos_variables(yeindex),ns,ng
stop
endif
endif

if (log10(local_emissivity(ns,ng)).ge.300.0d0) then
@@ -335,17 +336,19 @@ program make_table_example
eos_variables(tempindex),eos_variables(yeindex),ns,ng
stop
endif
if (local_delta(ns,ng).gt.1.0d0) then
write(*,"(a,1P4E18.9,i6,i6)") "delta > 1", &
local_delta(ns,ng),eos_variables(rhoindex), &
eos_variables(tempindex),eos_variables(yeindex),ns,ng
stop
endif
if (local_delta(ns,ng).lt.-1.0d0) then
write(*,"(a,1P4E18.9,i6,i6)") "delta < -1", &
local_delta(ns,ng),eos_variables(rhoindex), &
eos_variables(tempindex),eos_variables(yeindex),ns,ng
stop
if (.not. do_transport_opacities) then
if (local_delta(ns,ng).gt.1.0d0) then
write(*,"(a,1P4E18.9,i6,i6)") "delta > 1", &
local_delta(ns,ng),eos_variables(rhoindex), &
eos_variables(tempindex),eos_variables(yeindex),ns,ng
stop
endif
if (local_delta(ns,ng).lt.-1.0d0) then
write(*,"(a,1P4E18.9,i6,i6)") "delta < -1", &
local_delta(ns,ng),eos_variables(rhoindex), &
eos_variables(tempindex),eos_variables(yeindex),ns,ng
stop
endif
endif
enddo !do ng=1,mytable_number_groups
enddo !do ns=1,number_output_species
@@ -526,7 +529,6 @@ program make_table_example
#endif

do ieta=1,final_Itable_size_eta
write(*,*) "Eta:", 100.0*dble(ieta-1)/dble(final_Itable_size_eta),"%"
do iinE=final_Itable_size_inE,1,-1

#ifdef __MPI__
@@ -940,13 +940,13 @@ subroutine single_epannihil_kernel_point_return_all(iin,eta,temperature, &

!already ensure that amu and atau the same
if (add_anumu_kernel_epannihil) then
Phi0s(3,ng,1) = (epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,4,1,0) + &
Phi0s(4,ng,1) = (epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,4,1,0) + &
epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,6,1,0))/2.0d0 !production of Phi_0
Phi0s(3,ng,2) = (epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,4,2,0) + &
Phi0s(4,ng,2) = (epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,4,2,0) + &
epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,6,2,0))/2.0d0 !annihilation of Phi_0
Phi1s(3,ng,1) = (epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,4,1,1) + &
Phi1s(4,ng,1) = (epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,4,1,1) + &
epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,6,1,1))/2.0d0 !production of Phi_1
Phi1s(3,ng,2) = (epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,4,2,1) + &
Phi1s(4,ng,2) = (epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,4,2,1) + &
epannihil_Phi_Bruenn(nu_energy_x,nuother_energy_x,eta,6,2,1))/2.0d0 !annihilation of Phi_1

endif

0 comments on commit c6a57ca

Please sign in to comment.
You can’t perform that action at this time.