From 4748f43f9a78449bf618d863af7a098fee6b2568 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 3 Feb 2016 20:32:20 -0500 Subject: [PATCH] Switch temperature tag from materials to cells Note that this allows for distributed temperatures --- src/cross_section.F90 | 36 ++++++++++++------- src/geometry.F90 | 30 ++++++++++++++++ src/geometry_header.F90 | 3 ++ src/initialize.F90 | 78 +++++++++++++++++++++++++++++++++++++---- src/input_xml.F90 | 72 ++++++++++++++++++------------------- src/material_header.F90 | 1 - src/particle_header.F90 | 6 +++- 7 files changed, 168 insertions(+), 58 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index b62d9b389c..033e80b174 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -105,10 +105,13 @@ contains i_nuclide = mat % nuclide(i) ! Calculate microscopic cross section for this nuclide - if (p % E /= micro_xs(i_nuclide) % last_E .or. mat % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid, mat % sqrtkT) + if (p % E /= micro_xs(i_nuclide) % last_E & + .or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, & + i_grid, p % sqrtkT) else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then - call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid, mat % sqrtkT) + call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, & + i_grid, p % sqrtkT) end if ! ======================================================================== @@ -145,7 +148,8 @@ contains ! given index in the nuclides array at the energy of the given particle !=============================================================================== - subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, i_log_union, sqrtkT) + subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, & + i_log_union, sqrtkT) integer, intent(in) :: i_nuclide ! index into nuclides array integer, intent(in) :: i_sab ! index into sab_tables array real(8), intent(in) :: E ! energy @@ -590,14 +594,20 @@ contains ! sections in the resolved resonance regions !=============================================================================== - subroutine multipole_eval(multipole, Emev, sqrtkT, sigT, sigA, sigF) - type(MultipoleArray), intent(in) :: multipole ! The windowed multipole object to process. - real(8), intent(in) :: Emev ! The energy at which to evaluate the cross section in MeV - real(8), intent(in) :: sqrtkT ! The temperature in the form sqrt(kT (in eV)), at which to evaluate the cross section. - real(8), intent(out) :: sigT ! Total cross section - real(8), intent(out) :: sigA ! Absorption cross section - real(8), intent(out) :: sigF ! Fission cross section - complex(8) :: psi_ki ! The value of the psi-ki function for the asymptotic form + subroutine multipole_eval(multipole, Emev, sqrtkT_, sigT, sigA, sigF) + type(MultipoleArray), intent(in) :: multipole ! The windowed multipole + ! object to process. + real(8), intent(in) :: Emev ! The energy at which to + ! evaluate the cross section + ! in MeV + real(8), intent(in) :: sqrtkT_ ! The temperature in the form + ! sqrt(kT (in MeV)), at which + ! to evaluate the XS. + real(8), intent(out) :: sigT ! Total cross section + real(8), intent(out) :: sigA ! Absorption cross section + real(8), intent(out) :: sigF ! Fission cross section + complex(8) :: psi_ki ! The value of the psi-ki function for the asymptotic + ! form complex(8) :: c_temp ! complex temporary variable complex(8) :: w_val ! The faddeeva function evaluated at Z complex(8) :: Z ! sqrt(atomic weight ratio / kT) * (sqrt(E) - pole) @@ -607,6 +617,7 @@ contains real(8) :: dopp_ecoef ! sqrt(atomic weight ratio * pi / kT) / E real(8) :: temp ! real temporary value real(8) :: E ! energy, eV + real(8) :: sqrtkT ! sqrt(kT (in eV)) integer :: iP ! index of pole integer :: iC ! index of curvefit integer :: iW ! index of window @@ -617,6 +628,7 @@ contains ! Convert to eV E = Emev * 1.0e6_8 + sqrtkT = sqrtkT_ * 1.0e3_8 sqrtE = sqrt(E) invE = ONE/E diff --git a/src/geometry.F90 b/src/geometry.F90 index 8a38f982b5..5f37f1047b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -279,6 +279,36 @@ contains p % material = c % material(offset + 1) end if + ! Set the particle temperature + if (size(c % sqrtkT) == 1) then + ! Only one temperature for this cell; assign that one to the particle. + p % sqrtkT = c % sqrtkT(1) + else + ! Distributed instances of this cell have different temperatures. + ! Determine which instance this is and assign the matching temp. + distribcell_index = c % distribcell_index + offset = 0 + do k = 1, p % n_coord + if (cells(p % coord(k) % cell) % type == CELL_FILL) then + offset = offset + cells(p % coord(k) % cell) % & + offset(distribcell_index) + elseif (cells(p % coord(k) % cell) % type == CELL_LATTICE) then + if (lattices(p % coord(k + 1) % lattice) % obj & + % are_valid_indices([& + p % coord(k + 1) % lattice_x, & + p % coord(k + 1) % lattice_y, & + p % coord(k + 1) % lattice_z])) then + offset = offset + lattices(p % coord(k + 1) % lattice) % obj % & + offset(distribcell_index, & + p % coord(k + 1) % lattice_x, & + p % coord(k + 1) % lattice_y, & + p % coord(k + 1) % lattice_z) + end if + end if + end do + p % sqrtkT = c % sqrtkT(offset + 1) + end if + elseif (c % type == CELL_FILL) then CELL_TYPE ! ====================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 519b15bff4..0756e23dc4 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -139,6 +139,9 @@ module geometry_header ! only) integer :: distribcell_index ! Index corresponding to this cell in ! distribcell arrays + real(8), allocatable :: sqrtkT(:) ! Square root of k_Boltzmann * + ! temperature in MeV. Multiple for + ! distribcell ! Rotation matrix and translation vector real(8), allocatable :: translation(:) diff --git a/src/initialize.F90 b/src/initialize.F90 index de681fb949..610bfe538b 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -120,6 +120,9 @@ contains ! Create linked lists for multiple instances of the same nuclide call same_nuclide_list() + ! Set undefined cell temperatures to match the material data. + call lookup_material_temperatures() + ! Construct unionized or log energy grid for cross-sections select case (grid_method) case (GRID_NUCLIDE) @@ -967,10 +970,11 @@ contains end do end do - ! We also need distribcell if any distributed materials are present. + ! We also need distribcell if any distributed materials or distributed + ! temperatues are present. if (.not. distribcell_active) then do i = 1, n_cells - if (size(cells(i) % material) > 1) then + if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then distribcell_active = .true. exit end if @@ -995,8 +999,8 @@ contains end do end do - ! Make sure the number of materials matches the number of cell instances for - ! distributed materials. + ! Make sure the number of materials and temperatures matches the number of + ! cell instances. do i = 1, n_cells associate (c => cells(i)) if (size(c % material) > 1) then @@ -1008,6 +1012,15 @@ contains &equal one or the number of instances.") end if end if + if (size(c % sqrtkT) > 1) then + if (size(c % sqrtkT) /= c % instances) then + call fatal_error("Cell " // trim(to_str(c % id)) // " was & + &specified with " // trim(to_str(size(c % sqrtkT))) & + // " temperatures but has " // trim(to_str(c % instances)) & + // " distributed instances. The number of temperatures must & + &equal one or the number of instances.") + end if + end if end associate end do @@ -1049,9 +1062,9 @@ contains end do end do - ! List all cells with multiple (distributed) materials. + ! List all cells with multiple (distributed) materials or temperatures. do i = 1, n_cells - if (size(cells(i) % material) > 1) then + if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then call cell_list % add(i) end if end do @@ -1120,4 +1133,57 @@ contains end subroutine allocate_offsets +!=============================================================================== +! LOOKUP_MATERIAL_TEMPERATURES If any cells have undefined temperatures, try to +! find their temperatures from material data. +!=============================================================================== + + subroutine lookup_material_temperatures() + integer :: i, j, k + real(8) :: min_temp + logical :: warning_given + + warning_given = .false. + do i = 1, n_cells + ! Ignore non-normal cells and cells with defined temperature. + if (cells(i) % type /= CELL_NORMAL) cycle + if (cells(i) % sqrtkT(1) /= ERROR_REAL) cycle + + ! Set the number of temperatures equal to the number of materials. + deallocate(cells(i) % sqrtkT) + allocate(cells(i) % sqrtkT(size(cells(i) % material))) + + ! Check each of the cell materials for temperature data. + do j = 1, size(cells(i) % material) + ! Arbitrarily set void regions to 0K. + if (cells(i) % material(j) == MATERIAL_VOID) then + cells(i) % sqrtkT(j) = ZERO + cycle + end if + + associate (mat => materials(cells(i) % material(j))) + ! Find the temperature of the coldest nuclide. + min_temp = nuclides(mat % nuclide(1)) % kT + do k = 2, mat % n_nuclides + ! Warn the user if the nuclides don't have identical temperatues. + if (nuclides(mat % nuclide(k)) % kT /= min_temp & + .and. .not. warning_given) then + call warning("OpenMC cannot & + &identify the temperature of at least one cell. For the & + &purposes of multipole cross section evaluations, all cells & + &with unknown temperature will be set to the coldest & + &temperature found in the nuclear data for that cell's & + &material") + warning_given = .true. + end if + min_temp = min(min_temp, nuclides(mat % nuclide(k)) % kT) + end do + + ! Set the temperature for this cell instance. + cells(i) % sqrtkT(j) = sqrt(min_temp) + end associate + end do + end do + end subroutine lookup_material_temperatures + end module initialize diff --git a/src/input_xml.F90 b/src/input_xml.F90 index bbf4da4f9e..34cd0f1b95 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1297,6 +1297,40 @@ contains call get_node_array(node_cell, "translation", c % translation) end if + ! 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 (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. + if (c % material(1) == NONE) call fatal_error("Cell " & + // trim(to_str(c % id)) // " was specified with a temperature & + &but no material. Temperature specification is only valid for & + &cells filled with a material.") + + ! Copy in temperatures + allocate(c % sqrtkT(n)) + call get_node_array(node_cell, "temperature", c % sqrtkT) + + ! Make sure all temperatues are positive + do j = 1, size(c % sqrtkT) + if (c % sqrtkT(j) < ZERO) call fatal_error("Cell " & + // trim(to_str(c % id)) // " was specified with a negative & + &temperature. All cell temperatures must be non-negative.") + end do + + ! Convert to sqrt(kT) + c % sqrtkT(:) = sqrt(K_BOLTZMANN * c % sqrtkT(:)) + else + allocate(c % sqrtkT(1)) + c % sqrtkT(1) = ERROR_REAL + end if + else + allocate(c % sqrtkT(1)) + c % sqrtkT = ERROR_REAL + end if + ! Add cell to dictionary call cell_dict % add_key(c % id, i) @@ -1854,7 +1888,6 @@ contains real(8) :: temp_dble ! temporary double prec. real logical :: file_exists ! does materials.xml exist? logical :: sum_density ! density is taken to be sum of nuclide densities - logical :: temp_known ! Has the temperature yet been defined? character(12) :: name ! name of isotope, e.g. 92235.03c character(12) :: alias ! alias of nuclide, e.g. U-235.03c character(MAX_WORD_LEN) :: units ! units on density @@ -1938,17 +1971,6 @@ contains cycle end if - ! ======================================================================= - ! READ AND PARSE TAG - if (check_for_node(node_mat, "temperature")) then - call get_node_ptr(node_mat, "temperature", node_temp) - call get_node_value(node_temp, "value", temp_dble) - mat % sqrtkT = sqrt(temp_dble * K_BOLTZMANN * 1.0D6) - temp_known = .true. - else - temp_known = .false. - end if - ! ======================================================================= ! READ AND PARSE TAG @@ -2055,16 +2077,6 @@ contains call get_node_value(node_nuc, "xs", name) name = trim(temp_str) // "." // trim(name) - ! If needed, look up temperature - if (.not. temp_known) then - ! Find xs_listing and set the name/alias according to the listing - index_list = xs_listing_dict % get_key(to_lower(name)) - if(xs_listings(index_list) % kT /= 0.0_8) then - mat % sqrtkT = sqrt(xs_listings(index_list) % kT * 1.0D6) - temp_known = .true. - end if - end if - ! save name and density to list call list_names % append(name) @@ -2153,16 +2165,6 @@ contains temp_str = "data" end if - ! If still needed, look up temperature - if (.not. temp_known) then - ! Find xs_listing and set kT - index_list = xs_listing_dict % get_key(to_lower(list_names % tail % data)) - if(xs_listings(index_list) % kT /= 0.0_8) then - mat % sqrtkT = sqrt(xs_listings(index_list) % kT * 1.0D6) - temp_known = .true. - end if - end if - ! Set ace or iso-in-lab scattering for each nuclide in element do k = 1, n_nuc_ele if (adjustl(to_lower(temp_str)) == "iso-in-lab") then @@ -2177,12 +2179,6 @@ contains end do NATURAL_ELEMENTS - ! If still undefined, set the temperature to zero - if (.not. temp_known) then - mat % sqrtkT = 0.0_8 - temp_known = .true. - end if - ! ======================================================================== ! COPY NUCLIDES TO ARRAYS IN MATERIAL diff --git a/src/material_header.F90 b/src/material_header.F90 index 397fa0cd96..91c4cdfb82 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -13,7 +13,6 @@ module material_header integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm - real(8) :: sqrtkT ! sqrt(kT), kT in eV ! Energy grid information integer :: n_grid ! # of union material grid points diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 0b9b251ee5..f1d7375d7a 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -2,7 +2,7 @@ module particle_header use bank_header, only: Bank use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_SECONDARY, & - MAX_DELAYED_GROUPS + MAX_DELAYED_GROUPS, ERROR_REAL use error, only: fatal_error use geometry_header, only: BASE_UNIVERSE @@ -79,6 +79,9 @@ module particle_header integer :: material ! index for current material integer :: last_material ! index for last material + ! Temperature of the current cell + real(8) :: sqrtkT ! sqrt(k_Boltzmann * temperature) in MeV + ! Statistical data integer :: n_collision ! # of collisions @@ -124,6 +127,7 @@ contains this % absorb_wgt = ZERO this % n_bank = 0 this % wgt_bank = ZERO + this % sqrtkT = ERROR_REAL this % n_collision = 0 this % fission = .false. this % delayed_group = 0