crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Got the temperature dependent MGXS data hooked in to the rest of the code. On to debugging!

Browse source
This commit is contained in:
Adam Nelson 2016-09-04 06:34:33 -04:00
parent ab05e1b71b
commit aeb7f0c787
6 changed files with 186 additions and 135 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/cross_section.F90
View file

@ -161,6 +161,7 @@ contains
! If using multipole data but outside the RRR, pick the nearest
! temperature. Note that there is no tolerance here, so this
! temperature could be very far off!
kT = sqrtkT**2
i_temp = minloc(abs(nuclides(i_nuclide) % kTs - kT), dim=1)
end if
else

6
src/initialize.F90
View file

@ -113,12 +113,6 @@ contains
if (run_CE) then
! Construct log energy grid for cross-sections
call logarithmic_grid()
else
! Create material macroscopic data for MGXS
call time_read_xs%start()
call read_mgxs()
call create_macro_xs()
call time_read_xs%stop()
end if
! Allocate and setup tally stride, matching_bins, and tally maps

84
src/input_xml.F90
View file

@ -48,8 +48,8 @@ contains
call read_settings_xml()
call read_geometry_xml()
call read_materials()
call read_tallies_xml()
call read_materials()
if (cmfd_run) call configure_cmfd()
end subroutine read_input_xml
@ -1356,11 +1356,7 @@ contains
! Read cell temperatures. If the temperature is not specified, set it to
! ERROR_REAL for now. During initialization we'll replace ERROR_REAL with
! the temperature from the material data.
if (.not. run_CE) then
! Cell temperatures are not used for MG mode.
allocate(c % sqrtkT(1))
c % sqrtkT(1) = ZERO
else if (check_for_node(node_cell, "temperature")) then
if (check_for_node(node_cell, "temperature")) then
n = get_arraysize_double(node_cell, "temperature")
if (n > 0) then
! Make sure this is a "normal" cell.
@ -2038,7 +2034,7 @@ contains
if (run_CE) then
call read_ce_cross_sections_xml(libraries)
else
call read_mg_cross_sections_xml(libraries)
call read_mg_cross_sections_header(libraries)
end if
! Creating dictionary that maps the name of the material to the entry
@ -2069,10 +2065,16 @@ contains
call get_temperatures(nuc_temps, sab_temps)
! Read continuous-energy cross sections
if (run_CE .and. run_mode /= MODE_PLOTTING) then
call time_read_xs%start()
call read_ce_cross_sections(libraries, library_dict, nuc_temps, sab_temps)
call time_read_xs%stop()
if (run_mode /= MODE_PLOTTING) then
call time_read_xs % start()
if (run_CE) then
call read_ce_cross_sections(libraries, library_dict, nuc_temps, sab_temps)
else
! Create material macroscopic data for MGXS
call read_mgxs(nuc_temps)
call create_macro_xs()
end if
call time_read_xs % stop()
end if
! Normalize atom/weight percents
@ -4601,44 +4603,43 @@ contains
end subroutine read_ce_cross_sections_xml
subroutine read_mg_cross_sections_xml(libraries)
subroutine read_mg_cross_sections_header(libraries)
type(Library), allocatable, intent(out) :: libraries(:)
integer :: i ! loop index
integer :: n_libraries
logical :: file_exists ! does cross_sections.xml exist?
type(Node), pointer :: doc => null()
type(Node), pointer :: node_xsdata => null()
type(NodeList), pointer :: node_xsdata_list => null()
logical :: file_exists ! does mgxs.h5 exist?
integer(HID_T) :: file_id
real(8), allocatable :: rev_energy_bins(:)
character(len=255), allocatable :: names(:)
! Check if mgxs.xml exists
! Check if mgxs.h5 exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
! Could not find mgxs.xml file
call fatal_error("Cross sections XML file '" &
! Could not find mgxs.h5 file
call fatal_error("Cross sections HDF5 file '" &
// trim(path_cross_sections) // "' does not exist!")
end if
call write_message("Reading cross sections XML file...", 5)
call write_message("Reading cross sections HDF5 file...", 5)
! Parse mgxs.xml file
call open_xmldoc(doc, path_cross_sections)
! Open file for reading
file_id = file_open(path_cross_sections, 'r', parallel=.true.)
if (check_for_node(doc, "groups")) then
if (check_attribute(file_id, "groups")) then
! Get neutron group count
call get_node_value(doc, "groups", energy_groups)
call read_attribute(energy_groups, file_id, "groups")
else
call fatal_error("groups element must exist!")
call fatal_error("'groups' attribute must exist!")
end if
allocate(rev_energy_bins(energy_groups + 1))
allocate(energy_bins(energy_groups + 1))
if (check_for_node(doc, "group_structure")) then
if (check_attribute(file_id, "group structure")) then
! Get neutron group structure
call get_node_array(doc, "group_structure", energy_bins)
call read_attribute(energy_bins, file_id, "group structure")
else
call fatal_error("group_structures element must exist!")
call fatal_error("'group structure' attribute must exist!")
end if
! First reverse the order of energy_groups
@ -4650,9 +4651,9 @@ contains
end do
allocate(inverse_velocities(energy_groups))
if (check_for_node(doc, "inverse_velocities")) then
if (check_attribute(file_id, "inverse velocities")) then
! Get inverse velocities
call get_node_array(doc, "inverse_velocities", inverse_velocities)
call read_attribute(inverse_velocities, file_id, "inverse velocities")
else
! If not given, estimate them by using average energy in group which is
! assumed to be the midpoint
@ -4663,31 +4664,28 @@ contains
end do
end if
! Get node list of all <xsdata>
call get_node_list(doc, "xsdata", node_xsdata_list)
n_libraries = get_list_size(node_xsdata_list)
! Get the datasets present in the library
call get_groups(file_id, names)
n_libraries = size(names)
! Allocate xs_listings array
! Allocate libraries array
if (n_libraries == 0) then
call fatal_error("At least one <xsdata> element must be present in &
&mgxs.xml file!")
call fatal_error("At least one MGXS data set must be present in &
&mgxs library file!")
else
allocate(libraries(n_libraries))
end if
do i = 1, n_libraries
! Get pointer to xsdata table XML node
call get_list_item(node_xsdata_list, i, node_xsdata)
! Get name of material
allocate(libraries(i) % materials(1))
call get_node_value(node_xsdata, "name", libraries(i) % materials(1))
libraries(i) % materials(1) = names(i)
end do
! Close cross sections XML file
call close_xmldoc(doc)
! Close MGXS HDF file
call file_close(file_id)
end subroutine read_mg_cross_sections_xml
end subroutine read_mg_cross_sections_header
!===============================================================================
! EXPAND_NATURAL_ELEMENT converts natural elements specified using an <element>

117
src/mgxs_data.F90
View file

@ -3,13 +3,13 @@ module mgxs_data
use constants
use error, only: fatal_error
use global
use hdf5
use material_header, only: Material
use mgxs_header
use output, only: write_message
use set_header, only: SetChar
use stl_vector, only: VectorReal
use string, only: to_lower
use xml_interface
implicit none
contains
@ -19,21 +19,20 @@ contains
! nuclides and sab_tables arrays
!===============================================================================
subroutine read_mgxs()
subroutine read_mgxs(temps)
type(VectorReal), allocatable :: temps(:)
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: i_xsdata ! index in <xsdata> list
integer :: i_nuclide ! index in nuclides
character(20) :: name ! name of isotope, e.g. 92235.03c
integer :: i_nuclide ! index in nuclides array
character(20) :: name ! name of library to load
integer :: representation ! Data representation
type(Material), pointer :: mat
type(Material), pointer :: mat
type(SetChar) :: already_read
type(Node), pointer :: doc => null()
type(Node), pointer :: node_xsdata
type(NodeList), pointer :: node_xsdata_list => null()
integer(HID_T) :: file_id
integer(HID_T) :: xsdata_group
logical :: file_exists
character(MAX_LINE_LEN) :: temp_str
logical :: get_kfiss, get_fiss
integer :: l
type(DictCharInt) :: xsdata_dict
@ -42,30 +41,16 @@ contains
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
! Could not find cross_sections.xml file
call fatal_error("Cross sections XML file '" &
call fatal_error("Cross sections HDF5 file '" &
&// trim(path_cross_sections) // "' does not exist!")
end if
call write_message("Loading Cross Section Data...", 5)
! Parse cross_sections.xml file
call open_xmldoc(doc, path_cross_sections)
! Open file for reading
file_id = file_open(path_cross_sections, 'r', parallel=.true.)
! Get node list of all <xsdata>
call get_node_list(doc, "xsdata", node_xsdata_list)
! Build dictionary mapping nuclide names to an index in the <xsdata> node
! list
do i = 1, get_list_size(node_xsdata_list)
! Get pointer to xsdata table XML node
call get_list_item(node_xsdata_list, i, node_xsdata)
! Get name and create pair (name, i)
call get_node_value(node_xsdata, "name", name)
call xsdata_dict % add_key(to_lower(name), i)
end do
! allocate arrays for ACE table storage and cross section cache
! allocate arrays for MGXS storage and cross section cache
allocate(nuclides_MG(n_nuclides_total))
!$omp parallel
allocate(micro_xs(n_nuclides_total))
@ -102,20 +87,21 @@ contains
i_xsdata = xsdata_dict % get_key(to_lower(name))
i_nuclide = mat % nuclide(j)
! Get pointer to xsdata table XML node
call get_list_item(node_xsdata_list, i_xsdata, node_xsdata)
call write_message("Loading " // trim(name) // " Data...", 5)
! Check to make sure cross section set exists in the library
if (check_group(file_id, trim(name))) then
xsdata_group = open_group(file_id, trim(name))
else
call fatal_error("Data for '" // trim(name) // "' does not exist in "&
&// trim(path_cross_sections))
end if
! First find out the data representation
if (check_for_node(node_xsdata, "representation")) then
call get_node_value(node_xsdata, "representation", temp_str)
temp_str = trim(to_lower(temp_str))
if (temp_str == 'isotropic' .or. temp_str == 'iso') then
representation = MGXS_ISOTROPIC
else if (temp_str == 'angle') then
representation = MGXS_ANGLE
else
if (check_attribute(xsdata_group, "representation")) then
call read_attribute(representation, xsdata_group, "representation")
if (representation /= MGXS_ISOTROPIC .and. &
representation /= MGXS_ANGLE) then
call fatal_error("Invalid Data Representation!")
end if
else
@ -132,8 +118,9 @@ contains
end select
! Now read in the data specific to the type we just declared
call nuclides_MG(i_nuclide) % obj % init_file(node_xsdata, &
energy_groups, get_kfiss, get_fiss, max_order)
call nuclides_MG(i_nuclide) % obj % from_hdf5(xsdata_group, &
energy_groups, temps(i_nuclide), temperature_method, &
temperature_tolerance, get_kfiss, get_fiss, max_order)
! Add name to dictionary
call already_read % add(name)
@ -172,9 +159,13 @@ contains
subroutine create_macro_xs()
integer :: i_mat ! index in materials array
type(Material), pointer :: mat ! current material
type(VectorReal), allocatable :: kTs(:)
allocate(macro_xs(n_materials))
! Get temperatures to read for each material
call get_mat_kTs(kTs)
do i_mat = 1, n_materials
mat => materials(i_mat)
@ -187,9 +178,49 @@ contains
type is (MgxsAngle)
allocate(MgxsAngle :: macro_xs(i_mat) % obj)
end select
call macro_xs(i_mat) % obj % combine(mat, nuclides_MG, energy_groups, &
max_order)
call macro_xs(i_mat) % obj % combine(kTs(i_mat), mat, nuclides_MG, &
energy_groups, max_order)
end do
end subroutine create_macro_xs
!===============================================================================
! GET_MAT_kTs returns a list of temperatures (in MeV) that each
! material appears at in the model.
!===============================================================================
subroutine get_mat_kTs(kTs)
type(VectorReal), allocatable, intent(out) :: kTs(:)
integer :: i, j
integer :: i_material ! Index in materials array
real(8) :: kT ! temperature in MeV
allocate(kTs(size(materials)))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
! Get temperature of cell (rounding to nearest integer)
if (size(cells(i) % sqrtkT) > 1) then
kT = cells(i) % sqrtkT(j)**2
else
kT = cells(i) % sqrtkT(1)**2
end if
i_material = material_dict % get_key(cells(i) % material(j))
! Add temperature if it hasn't already been added
if (find(kTs(i_material), kT) == -1) then
call kTs(i_material) % push_back(kT)
end if
end do
end do
end subroutine get_mat_kTs
end module mgxs_data

109
src/mgxs_header.F90
View file

@ -57,6 +57,9 @@ module mgxs_header
logical :: fissionable ! mgxs object is fissionable?
integer :: scatter_type ! either legendre, histogram, or tabular.
! Caching information
integer :: index_temp ! temperature index for nuclide
contains
procedure(mgxs_from_hdf5_), deferred :: from_hdf5 ! Load the data
procedure(mgxs_combine_), deferred :: combine ! initializes object
@ -67,6 +70,7 @@ module mgxs_header
procedure(mgxs_sample_scatter_), deferred :: sample_scatter
! Calculate the material specific MGXS data from the nuclides
procedure(mgxs_calculate_xs_), deferred :: calculate_xs
procedure :: find_temperature => mgxs_find_temperature
end type Mgxs
!===============================================================================
@ -106,10 +110,9 @@ module mgxs_header
integer, intent(in) :: max_order ! Maximum requested order
end subroutine mgxs_combine_
pure function mgxs_get_xs_(this, t, xstype, gin, gout, uvw, mu) result(xs_val)
pure function mgxs_get_xs_(this, xstype, gin, gout, uvw, mu) result(xs_val)
import Mgxs
class(Mgxs), intent(in) :: this
integer, intent(in) :: t ! Temperature index of data
character(*), intent(in) :: xstype ! Cross Section Type
integer, intent(in) :: gin ! Incoming Energy group
integer, optional, intent(in) :: gout ! Outgoing Group
@ -118,20 +121,18 @@ module mgxs_header
real(8) :: xs_val ! Resultant xs
end function mgxs_get_xs_
function mgxs_sample_fission_(this, t, gin, uvw) result(gout)
function mgxs_sample_fission_(this, gin, uvw) result(gout)
import Mgxs
class(Mgxs), intent(in) :: this
integer, intent(in) :: t ! Temperature index
integer, intent(in) :: gin ! Incoming energy group
real(8), intent(in) :: uvw(3) ! Particle Direction
integer :: gout ! Sampled outgoing group
end function mgxs_sample_fission_
subroutine mgxs_sample_scatter_(this, t, uvw, gin, gout, mu, wgt)
subroutine mgxs_sample_scatter_(this, uvw, gin, gout, mu, wgt)
import Mgxs
class(Mgxs), intent(in) :: this
integer, intent(in) :: t ! Temperature index of data
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
@ -139,10 +140,9 @@ module mgxs_header
real(8), intent(inout) :: wgt ! Particle weight
end subroutine mgxs_sample_scatter_
subroutine mgxs_calculate_xs_(this, t, gin, uvw, xs)
subroutine mgxs_calculate_xs_(this, gin, uvw, xs)
import Mgxs, MaterialMacroXS
class(Mgxs), intent(in) :: this
integer, intent(in) :: t ! Temperature index
integer, intent(in) :: gin ! Incoming neutron group
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
type(MaterialMacroXS), intent(inout) :: xs ! Resultant Mgxs Data
@ -1006,7 +1006,7 @@ module mgxs_header
subroutine mgxsiso_combine(this, temps, mat, nuclides, groups, max_order)
class(MgxsIso), intent(inout) :: this ! The Mgxs to initialize
type(VectorReal), intent(in) :: temps ! Temperatures to obtain
type(VectorReal), intent(in) :: temps ! Temperatures to obtain [MeV]
type(Material), pointer, intent(in) :: mat ! base material
type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from
integer, intent(in) :: groups ! Number of E groups
@ -1406,15 +1406,17 @@ module mgxs_header
! MGXS*_GET_XS returns the requested data cross section data
!===============================================================================
pure function mgxsiso_get_xs(this, t, xstype, gin, gout, uvw, mu) result(xs)
pure function mgxsiso_get_xs(this, xstype, gin, gout, uvw, mu) result(xs)
class(MgxsIso), intent(in) :: this ! The Xs to get data from
integer, intent(in) :: t ! Temperature index of data
character(*) , intent(in) :: xstype ! Type of xs requested
integer, intent(in) :: gin ! Incoming Energy group
integer, optional, intent(in) :: gout ! Outgoing Energy group
real(8), optional, intent(in) :: uvw(3) ! Requested Angle
real(8), optional, intent(in) :: mu ! Change in angle
real(8) :: xs ! Requested x/s
integer :: t ! temperature index
t = this % index_temp
select case(xstype)
case('total')
@ -1494,9 +1496,8 @@ module mgxs_header
end select
end function mgxsiso_get_xs
pure function mgxsang_get_xs(this, t, xstype, gin, gout, uvw, mu) result(xs)
pure function mgxsang_get_xs(this, xstype, gin, gout, uvw, mu) result(xs)
class(MgxsAngle), intent(in) :: this ! The Mgxs to initialize
integer, intent(in) :: t ! Temperature index of data
character(*) , intent(in) :: xstype ! Type of xs requested
integer, intent(in) :: gin ! Incoming Energy group
integer, optional, intent(in) :: gout ! Outgoing Energy group
@ -1504,7 +1505,9 @@ module mgxs_header
real(8), optional, intent(in) :: mu ! Change in angle
real(8) :: xs ! Requested x/s
integer :: iazi, ipol
integer :: iazi, ipol, t
t = this % index_temp
if (present(uvw)) then
call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol)
@ -1594,9 +1597,8 @@ module mgxs_header
! MGXS*_SAMPLE_FISSION_ENERGY samples the outgoing energy from a fission event
!===============================================================================
function mgxsiso_sample_fission_energy(this, t, gin, uvw) result(gout)
function mgxsiso_sample_fission_energy(this, gin, uvw) result(gout)
class(MgxsIso), intent(in) :: this ! Data to work with
integer, intent(in) :: t ! Temperature index of data
integer, intent(in) :: gin ! Incoming energy group
real(8), intent(in) :: uvw(3) ! Particle Direction
integer :: gout ! Sampled outgoing group
@ -1605,18 +1607,17 @@ module mgxs_header
xi = prn()
gout = 1
prob = this % xs(t) % chi(gout,gin)
prob = this % xs(this % index_temp) % chi(gout,gin)
do while (prob < xi)
gout = gout + 1
prob = prob + this % xs(t) % chi(gout,gin)
prob = prob + this % xs(this % index_temp) % chi(gout,gin)
end do
end function mgxsiso_sample_fission_energy
function mgxsang_sample_fission_energy(this, t, gin, uvw) result(gout)
function mgxsang_sample_fission_energy(this, gin, uvw) result(gout)
class(MgxsAngle), intent(in) :: this ! Data to work with
integer, intent(in) :: t ! Temperature index of data
integer, intent(in) :: gin ! Incoming energy group
real(8), intent(in) :: uvw(3) ! Particle Direction
integer :: gout ! Sampled outgoing group
@ -1628,11 +1629,11 @@ module mgxs_header
xi = prn()
gout = 1
prob = this % xs(t) % chi(gout, gin, iazi, ipol)
prob = this % xs(this % index_temp) % chi(gout, gin, iazi, ipol)
do while (prob < xi)
gout = gout + 1
prob = prob + this % xs(t) % chi(gout, gin, iazi, ipol)
prob = prob + this % xs(this % index_temp) % chi(gout, gin, iazi, ipol)
end do
end function mgxsang_sample_fission_energy
@ -1641,22 +1642,20 @@ module mgxs_header
! MGXS*_SAMPLE_SCATTER Selects outgoing energy and angle after a scatter event
!===============================================================================
subroutine mgxsiso_sample_scatter(this, t, uvw, gin, gout, mu, wgt)
subroutine mgxsiso_sample_scatter(this, uvw, gin, gout, mu, wgt)
class(MgxsIso), intent(in) :: this
integer, intent(in) :: t ! Temperature index
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
real(8), intent(out) :: mu ! Sampled change in angle
real(8), intent(inout) :: wgt ! Particle weight
call this % xs(t) % scatter % sample(gin, gout, mu, wgt)
call this % xs(this % index_temp) % scatter % sample(gin, gout, mu, wgt)
end subroutine mgxsiso_sample_scatter
subroutine mgxsang_sample_scatter(this, t, uvw, gin, gout, mu, wgt)
subroutine mgxsang_sample_scatter(this, uvw, gin, gout, mu, wgt)
class(MgxsAngle), intent(in) :: this
integer, intent(in) :: t ! Temperature index
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
integer, intent(in) :: gin ! Incoming neutron group
integer, intent(out) :: gout ! Sampled outgoin group
@ -1666,7 +1665,8 @@ module mgxs_header
integer :: iazi, ipol ! Angular indices
call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol)
call this % xs(t) % scatter(iazi, ipol) % obj % sample(gin, gout, mu, wgt)
call this % xs(this % index_temp) % scatter(iazi, ipol) % obj % sample( &
gin, gout, mu, wgt)
end subroutine mgxsang_sample_scatter
@ -1675,24 +1675,22 @@ module mgxs_header
! for the material the particle is currently traveling through.
!===============================================================================
subroutine mgxsiso_calculate_xs(this, t, gin, uvw, xs)
subroutine mgxsiso_calculate_xs(this, gin, uvw, xs)
class(MgxsIso), intent(in) :: this
integer, intent(in) :: t ! Temperature index
integer, intent(in) :: gin ! Incoming neutron group
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
type(MaterialMacroXS), intent(inout) :: xs ! Resultant Mgxs Data
xs % total = this % xs(t) % total(gin)
xs % elastic = this % xs(t) % scatter % scattxs(gin)
xs % absorption = this % xs(t) % absorption(gin)
xs % fission = this % xs(t) % fission(gin)
xs % nu_fission = this % xs(t) % nu_fission(gin)
xs % total = this % xs(this % index_temp) % total(gin)
xs % elastic = this % xs(this % index_temp) % scatter % scattxs(gin)
xs % absorption = this % xs(this % index_temp) % absorption(gin)
xs % fission = this % xs(this % index_temp) % fission(gin)
xs % nu_fission = this % xs(this % index_temp) % nu_fission(gin)
end subroutine mgxsiso_calculate_xs
subroutine mgxsang_calculate_xs(this, t, gin, uvw, xs)
subroutine mgxsang_calculate_xs(this, gin, uvw, xs)
class(MgxsAngle), intent(in) :: this
integer, intent(in) :: t ! Temperature index
integer, intent(in) :: gin ! Incoming neutron group
real(8), intent(in) :: uvw(3) ! Incoming neutron direction
type(MaterialMacroXS), intent(inout) :: xs ! Resultant Mgxs Data
@ -1700,16 +1698,41 @@ module mgxs_header
integer :: iazi, ipol
call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol)
xs % total = this % xs(t) % total(gin, iazi, ipol)
xs % elastic = this % xs(t) % scatter(iazi, ipol) % obj % scattxs(gin)
xs % absorption = this % xs(t) % absorption(gin, iazi, ipol)
xs % fission = this % xs(t) % fission(gin, iazi, ipol)
xs % nu_fission = this % xs(t) % nu_fission(gin, iazi, ipol)
xs % total = this % xs(this % index_temp) % &
total(gin, iazi, ipol)
xs % elastic = this % xs(this % index_temp) % &
scatter(iazi, ipol) % obj % scattxs(gin)
xs % absorption = this % xs(this % index_temp) % &
absorption(gin, iazi, ipol)
xs % fission = this % xs(this % index_temp) % &
fission(gin, iazi, ipol)
xs % nu_fission = this % xs(this % index_temp) % &
nu_fission(gin, iazi, ipol)
end subroutine mgxsang_calculate_xs
!===============================================================================
! find_angle finds the closest angle on the data grid and returns that index
! MGXS_FIND_TEMPERATURE sets the temperature index for the given
! sqrt(temperature), (with temperature in units of MeV)
!===============================================================================
pure subroutine mgxs_find_temperature(this, sqrtkT, tolerance)
class(Mgxs), intent(in) :: this
real(8), intent(in) :: sqrtkT ! Temperature (in units of of MeV)
real(8), intent(in) :: tolerance ! Temperature tolerance
real(8) :: kT
integer :: i_temp
kT = sqrtkT**2
do i_temp = 1, size(this % kTs)
if (abs(this % kTs(i_temp) - kT) < K_BOLTZMANN * tolerance) exit
end do
end subroutine mgxs_find_temperature
!===============================================================================
! FIND_ANGLE finds the closest angle on the data grid and returns that index
!===============================================================================
pure subroutine find_angle(polar, azimuthal, uvw, i_azi, i_pol)

4
src/tracking.F90
View file

@ -92,6 +92,10 @@ contains
! Since the MGXS can be angle dependent, this needs to be done
! After every collision for the MGXS mode
if (p % material /= MATERIAL_VOID) then
! Update the temperature index
call macro_xs(p % material) % obj % find_temperature(p % sqrtkT, &
temperature_tolerance)
! Get the data
call macro_xs(p % material) % obj % calculate_xs(p % g, &
p % coord(p % n_coord) % uvw, material_xs)
else