Permalink
Browse files

propagate changes for the inclusion of the scattering delta to the ma…

…ke_table_weakrates.F90 file
  • Loading branch information...
Evan O'Connor
Evan O'Connor committed Nov 8, 2017
1 parent 60da4cb commit bedf4c02fcabfa6739344ff7e51d375b2bf83f4d
Showing with 42 additions and 3 deletions.
  1. +42 −3 src/make_table_weakrates.F90
@@ -46,6 +46,7 @@ program make_table_example
real*8, allocatable,dimension(:,:,:,:,:) :: table_emission
real*8, allocatable,dimension(:,:,:,:,:) :: table_absopacity
real*8, allocatable,dimension(:,:,:,:,:) :: table_scatopacity
real*8, allocatable,dimension(:,:,:,:,:) :: table_delta
!final Itable parameters
@@ -73,6 +74,7 @@ program make_table_example
real*8, allocatable,dimension(:,:) :: local_emissivity
real*8, allocatable,dimension(:,:) :: local_absopacity
real*8, allocatable,dimension(:,:) :: local_scatopacity
real*8, allocatable,dimension(:,:) :: local_delta
real*8, allocatable,dimension(:,:) :: local_Phi0
real*8, allocatable,dimension(:,:) :: local_Phi1
real*8, allocatable,dimension(:,:,:) :: local_Phi0_epannihil
@@ -96,6 +98,7 @@ program make_table_example
real*8, allocatable,dimension(:,:,:,:,:) :: table_emission_node
real*8, allocatable,dimension(:,:,:,:,:) :: table_absopacity_node
real*8, allocatable,dimension(:,:,:,:,:) :: table_scatopacity_node
real*8, allocatable,dimension(:,:,:,:,:) :: table_delta_node
real*8, allocatable,dimension(:) :: Itable_temp_subset
real*8, allocatable,dimension(:,:,:,:,:) :: Itable_Phi0_node
real*8, allocatable,dimension(:,:,:,:,:) :: Itable_Phi1_node
@@ -210,6 +213,8 @@ program make_table_example
final_table_size_ye,number_output_species,mytable_number_groups))
allocate(table_scatopacity(final_table_size_rho,final_table_size_temp, &
final_table_size_ye,number_output_species,mytable_number_groups))
allocate(table_delta(final_table_size_rho,final_table_size_temp, &
final_table_size_ye,number_output_species,mytable_number_groups))
#ifdef __MPI__
!mpi node tables
@@ -219,6 +224,8 @@ program make_table_example
final_table_size_ye,number_output_species,mytable_number_groups))
allocate(table_scatopacity_node(final_table_size_rho,final_table_size_temp, &
final_table_size_ye,number_output_species,mytable_number_groups))
allocate(table_delta_node(final_table_size_rho,final_table_size_temp, &
final_table_size_ye,number_output_species,mytable_number_groups))
table_emission_node = 0.0d0
#endif
@@ -249,7 +256,7 @@ program make_table_example
!$OMP PARALLEL DO PRIVATE(itemp,iye,local_emissivity,local_absopacity,local_scatopacity, &
!$OMP ns,ng,eos_variables,keytemp,keyerr,matter_prs,matter_ent,matter_cs2,matter_dedt, &
!$OMP matter_dpderho,matter_dpdrhoe,hempel_lookup_table)
!$OMP matter_dpderho,matter_dpdrhoe,hempel_lookup_table,local_delta)
!loop over rho,temp,ye of table, do each point
#ifdef __MPI__
do irho=1,mpi_final_table_size_rho
@@ -260,6 +267,7 @@ program make_table_example
allocate(local_emissivity(number_output_species,mytable_number_groups))
allocate(local_absopacity(number_output_species,mytable_number_groups))
allocate(local_scatopacity(number_output_species,mytable_number_groups))
allocate(local_delta(number_output_species,mytable_number_groups))
allocate(eos_variables(total_eos_variables))
#ifdef __MPI__
write(*,*) "Rho:", 100.0*dble(displs(mpirank)+irho-1)/dble(final_table_size_rho),"%"
@@ -284,7 +292,7 @@ program make_table_example
!calculate the rho,temp,ye
call single_point_return_all(eos_variables, &
local_emissivity,local_absopacity,local_scatopacity, &
local_emissivity,local_absopacity,local_scatopacity,local_delta, &
mytable_neutrino_scheme)
!check that the number is not NaN or Inf (.gt.1.0d300)
@@ -305,6 +313,11 @@ 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
endif
if (log10(local_emissivity(ns,ng)).ge.300.0d0) then
write(*,"(a,1P4E18.9,i6,i6)") "We have a Inf in emissivity", &
@@ -324,6 +337,18 @@ 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
endif
enddo !do ng=1,mytable_number_groups
enddo !do ns=1,number_output_species
@@ -334,10 +359,12 @@ program make_table_example
table_emission_node(displs(mpirank)+irho,itemp,iye,ns,ng) = local_emissivity(ns,ng) !ergs/cm^3/s/MeV/srad
table_absopacity_node(displs(mpirank)+irho,itemp,iye,ns,ng) = local_absopacity(ns,ng) !cm^-1
table_scatopacity_node(displs(mpirank)+irho,itemp,iye,ns,ng) = local_scatopacity(ns,ng) !cm^-1
table_delta_node(displs(mpirank)+irho,itemp,iye,ns,ng) = local_delta(ns,ng) !dimensionless
#else
table_emission(irho,itemp,iye,ns,ng) = local_emissivity(ns,ng) !ergs/cm^3/s/MeV/srad
table_absopacity(irho,itemp,iye,ns,ng) = local_absopacity(ns,ng) !cm^-1
table_scatopacity(irho,itemp,iye,ns,ng) = local_scatopacity(ns,ng) !cm^-1
table_delta(irho,itemp,iye,ns,ng) = local_delta(ns,ng) !dimensionless
#endif
enddo !do ns=1,number_output_species
enddo !do ng=1,mytable_number_groups
@@ -348,6 +375,7 @@ program make_table_example
deallocate(local_emissivity)
deallocate(local_absopacity)
deallocate(local_scatopacity)
deallocate(local_delta)
deallocate(eos_variables)
enddo!do irho=1,final_table_size_rho
!$OMP END PARALLEL DO! end do
@@ -364,7 +392,8 @@ program make_table_example
table_size,mpi_double,mpi_sum,0,mpi_comm_world,ierror)
call mpi_reduce(table_scatopacity_node,table_scatopacity,&
table_size,mpi_double,mpi_sum,0,mpi_comm_world,ierror)
call mpi_reduce(table_delta_node,table_delta,&
table_size,mpi_double,mpi_sum,0,mpi_comm_world,ierror)
!begin inelastic, timing for mpi purposes
@@ -910,6 +939,16 @@ subroutine write_h5(filename,timestamp)
call h5sclose_f(dspace_id, error)
cerror = cerror + error
if(.not. do_transport_opacities) then
call h5screate_simple_f(rank, dims5, dspace_id, error)
call h5dcreate_f(file_id, "scattering_delta", H5T_NATIVE_DOUBLE, &
dspace_id, dset_id, error)
call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,table_delta, dims5, error)
call h5dclose_f(dset_id, error)
call h5sclose_f(dspace_id, error)
cerror = cerror + error
endif
if (doing_inelastic.or.doing_epannihil) then
rank = 1

0 comments on commit bedf4c0

Please sign in to comment.