Permalink
Browse files

add in separate extraction and evolution radii

  • Loading branch information...
Evan O'Connor
Evan O'Connor committed Jul 13, 2017
1 parent db0f20b commit 4d11160b4f93f71b03e4694d6bb8cf26e0662531
Showing with 31 additions and 21 deletions.
  1. +2 −0 src/GR1D_module.F90
  2. +6 −1 src/M1/M1_init.F90
  3. +2 −0 src/initialize_vars.F90
  4. +2 −1 src/input_parser.F90
  5. +19 −19 src/output.F90
View
@@ -154,7 +154,9 @@ module GR1D_module
integer :: v_order
real*8 :: M1_maxradii
real*8 :: M1_extractradii
integer :: M1_imaxradii
integer :: M1_iextractradii
integer :: number_species
integer :: number_species_to_evolve
integer :: number_groups
View
@@ -8,7 +8,7 @@ subroutine M1_init
temp_mev_to_kelvin,number_species,volume,pi,M1_moment_to_distro,clite, &
hbarc_mevcm,M1_testcase_number,v_order,include_nes_kernels, &
M1_moment_to_distro_inverse,nulib_kernel_gf,number_species_to_evolve, &
include_epannihil_kernels
include_epannihil_kernels,M1_extractradii,M1_iextractradii
use nulibtable
implicit none
@@ -60,6 +60,11 @@ subroutine M1_init
if (x1(i)/length_gf.lt.M1_maxradii) M1_imaxradii = i
enddo
!find zone to extract
do i=1,n1
if (x1(i)/length_gf.lt.M1_extractradii) M1_iextractradii = i
enddo
if (v_order.eq.0.or.v_order.eq.-1) then
write(*,*) "Velocity order is:",v_order," (-1 for all orders)"
else
View
@@ -41,7 +41,9 @@ subroutine initialize_vars
M1_do_backwardfix = 0
v_order = -1
M1_maxradii = 0.0d0
M1_extractradii = 0.0d0
M1_imaxradii = 0
M1_iextractradii = 0
number_species = 0
number_species = -1
number_species_to_evolve = -1
View
@@ -144,7 +144,8 @@ subroutine input_parser
if(do_M1) then
call get_integer_parameter('v_order',v_order)
if (fake_neutrinos) stop "M1 neutrino's are not fake, how dare you (fake_neutrino needs to be 0 for M1)"
call get_double_parameter('extraction_radii',M1_maxradii)
call get_double_parameter('evolution_radii',M1_maxradii)
call get_double_parameter('extraction_radii',M1_extractradii)
call get_integer_parameter('number_species',number_species)
call get_integer_parameter('number_groups',number_groups)
call get_integer_parameter('number_eas',number_eas)
View
@@ -421,19 +421,19 @@ subroutine output_all(modeflag)
average_energy_fluid = 0.0d0
rms_energy_fluid = 0.0d0
do k=1,number_species
luminosity(k) = sum(q_M1(M1_imaxradii,k,:,2))&
* (4.0d0*pi)**2*x1(M1_imaxradii)**2 * (clite**5/ggrav)
luminosity_fluid(k) = sum(q_M1_fluid(M1_imaxradii,k,:,2))&
* (4.0d0*pi)**2*x1(M1_imaxradii)**2 * (clite**5/ggrav)
num_luminosity_fluid(k) = sum(q_M1_fluid(M1_imaxradii,k,:,2)/nulibtable_energies(:))&
* (4.0d0*pi)**2*x1(M1_imaxradii)**2 * (clite**3/ggrav*mass_gf)
average_energy_fluid(k) = sum(q_M1_fluid(M1_imaxradii,k,:,2))&
/ sum(q_M1_fluid(M1_imaxradii,k,:,2)/nulibtable_energies(:)) / (mev_to_erg*energy_gf)
rms_energy_fluid(k) = sqrt( sum(q_M1_fluid(M1_imaxradii,k,:,2)*nulibtable_energies(:))&
/ sum(q_M1_fluid(M1_imaxradii,k,:,2)/nulibtable_energies(:)) ) / (mev_to_erg*energy_gf)
average_energy(k) = average_energy_fluid(k)/sqrt(1.0d0-v(M1_imaxradii)**2)*(1.0d0+v(M1_imaxradii))
rms_energy(k) = rms_energy_fluid(k)/sqrt(1.0d0-v(M1_imaxradii)**2)*(1.0d0+v(M1_imaxradii))
do k=1,number_species
luminosity(k) = sum(q_M1(M1_iextractradii,k,:,2))&
* (4.0d0*pi)**2*x1(M1_iextractradii)**2 * (clite**5/ggrav)
luminosity_fluid(k) = sum(q_M1_fluid(M1_iextractradii,k,:,2))&
* (4.0d0*pi)**2*x1(M1_iextractradii)**2 * (clite**5/ggrav)
num_luminosity_fluid(k) = sum(q_M1_fluid(M1_iextractradii,k,:,2)/nulibtable_energies(:))&
* (4.0d0*pi)**2*x1(M1_iextractradii)**2 * (clite**3/ggrav*mass_gf)
average_energy_fluid(k) = sum(q_M1_fluid(M1_iextractradii,k,:,2))&
/ sum(q_M1_fluid(M1_iextractradii,k,:,2)/nulibtable_energies(:)) / (mev_to_erg*energy_gf)
rms_energy_fluid(k) = sqrt( sum(q_M1_fluid(M1_iextractradii,k,:,2)*nulibtable_energies(:))&
/ sum(q_M1_fluid(M1_iextractradii,k,:,2)/nulibtable_energies(:)) ) / (mev_to_erg*energy_gf)
average_energy(k) = average_energy_fluid(k)/sqrt(1.0d0-v(M1_iextractradii)**2)*(1.0d0+v(M1_iextractradii))
rms_energy(k) = rms_energy_fluid(k)/sqrt(1.0d0-v(M1_iextractradii)**2)*(1.0d0+v(M1_iextractradii))
num_luminosity(k) = luminosity(k)/(average_energy(k)*mev_to_erg)
enddo
@@ -587,27 +587,27 @@ subroutine output_all(modeflag)
call output_many_scalars(scalars,nscalars0,nscalars,filename)
filename = trim(adjustl(outdir))//"/M1_nue_fluxspectra_out.xg"
spectrum = q_M1(M1_imaxradii,1,:,2)*M1_moment_to_distro(:)
spectrum = q_M1(M1_iextractradii,1,:,2)*M1_moment_to_distro(:)
call output_spectra(spectrum,filename)
filename = trim(adjustl(outdir))//"/M1_anue_fluxspectra_out.xg"
spectrum = q_M1(M1_imaxradii,2,:,2)*M1_moment_to_distro(:)
spectrum = q_M1(M1_iextractradii,2,:,2)*M1_moment_to_distro(:)
call output_spectra(spectrum,filename)
filename = trim(adjustl(outdir))//"/M1_nux_fluxspectra_out.xg"
spectrum = q_M1(M1_imaxradii,3,:,2)*M1_moment_to_distro(:)
spectrum = q_M1(M1_iextractradii,3,:,2)*M1_moment_to_distro(:)
call output_spectra(spectrum,filename)
filename = trim(adjustl(outdir))//"/M1_nue_enspectra_out.xg"
spectrum = q_M1_fluid(M1_imaxradii,1,:,1)*M1_moment_to_distro(:)
spectrum = q_M1_fluid(M1_iextractradii,1,:,1)*M1_moment_to_distro(:)
call output_spectra(spectrum,filename)
filename = trim(adjustl(outdir))//"/M1_anue_enspectra_out.xg"
spectrum = q_M1_fluid(M1_imaxradii,2,:,1)*M1_moment_to_distro(:)
spectrum = q_M1_fluid(M1_iextractradii,2,:,1)*M1_moment_to_distro(:)
call output_spectra(spectrum,filename)
filename = trim(adjustl(outdir))//"/M1_nux_enspectra_out.xg"
spectrum = q_M1_fluid(M1_imaxradii,3,:,1)*M1_moment_to_distro(:)
spectrum = q_M1_fluid(M1_iextractradii,3,:,1)*M1_moment_to_distro(:)
call output_spectra(spectrum,filename)
filename = trim(adjustl(outdir))//"/M1_nue_fluxspectra_cen.xg"

0 comments on commit 4d11160

Please sign in to comment.