From 63aa3bf10648f1af1fe1afee1e458b5a593a23db Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 25 Jan 2019 16:00:00 -0600 Subject: [PATCH] Remove most of remaining components of Material type on Fortran side --- include/openmc/material.h | 9 + src/initialize.cpp | 4 - src/input_xml.F90 | 350 ++------------------------------------ src/material.cpp | 36 +++- src/material_header.F90 | 90 +++++----- src/nuclide.cpp | 2 + src/photon.cpp | 2 + src/simulation.F90 | 11 -- src/simulation.cpp | 9 + src/tallies/tally.F90 | 202 +++++++++++----------- src/volume_calc.F90 | 20 +-- 11 files changed, 225 insertions(+), 510 deletions(-) diff --git a/include/openmc/material.h b/include/openmc/material.h index 0443c3480..4792e22fc 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -53,6 +53,9 @@ public: //! so the code knows when to apply bound thermal scattering data void init_thermal(); + //! Set up mapping between global nuclides vector and indices in nuclide_ + void init_nuclide_index(); + //! Finalize the material, assigning tables, normalize density, etc. void finalize(); @@ -75,6 +78,12 @@ public: bool depletable_ {false}; //!< Is the material depletable? std::vector p0_; //!< Indicate which nuclides are to be treated with iso-in-lab scattering + // To improve performance of tallying, we store an array (direct address + // table) that indicates for each nuclide in data::nuclides the index of the + // corresponding nuclide in the nuclide_ vector. If it is not present in the + // material, the entry is set to -1. + std::vector mat_nuclide_index_; + // Thermal scattering tables std::vector thermal_tables_; diff --git a/src/initialize.cpp b/src/initialize.cpp index aaa318946..69986fdc0 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -30,7 +30,6 @@ #include "openmc/timer.h" // data/functions from Fortran side -extern "C" void normalize_ao(); extern "C" void print_usage(); extern "C" void print_version(); extern "C" void read_command_line(); @@ -296,9 +295,6 @@ void read_input_xml() if (mpi::master && settings::verbosity >= 5) print_plot(); } else { - // Normalize atom/weight percents - normalize_ao(); - // Write summary information if (mpi::master && settings::output_summary) write_summary(); diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c97642e1f..afbcc72fe 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -387,30 +387,24 @@ contains subroutine read_materials_xml() bind(C) integer :: i ! loop index for materials - integer :: j ! loop index for nuclides - integer :: n ! number of nuclides - integer :: index_nuclide ! index in nuclides - integer :: index_element ! index in elements - integer :: index_sab ! index in sab_tables logical :: file_exists ! does materials.xml exist? - character(20) :: name ! name of nuclide, e.g. U235 - character(MAX_WORD_LEN) :: units ! units on density character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml - real(8) :: val ! value entered for density - real(8) :: temp_dble ! temporary double prec. real - logical :: sum_density ! density is sum of nuclide densities - type(VectorChar) :: names ! temporary list of nuclide names - type(VectorReal) :: densities ! temporary list of nuclide densities type(Material), pointer :: mat => null() type(XMLDocument) :: doc type(XMLNode) :: root - type(XMLNode) :: node_mat - type(XMLNode) :: node_dens - type(XMLNode) :: node_nuc type(XMLNode), allocatable :: node_mat_list(:) - type(XMLNode), allocatable :: node_nuc_list(:) - type(XMLNode), allocatable :: node_ele_list(:) - type(XMLNode), allocatable :: node_macro_list(:) + + interface + function nuclides_size() bind(C) result(n) + import C_INT + integer(C_INT) :: n + end function + + function elements_size() bind(C) result(n) + import C_INT + integer(C_INT) :: n + end function + end interface ! Display output message call write_message("Reading materials XML file...", 5) @@ -436,247 +430,13 @@ contains n_materials = size(node_mat_list) allocate(materials(n_materials)) - ! Initialize count for number of nuclides/S(a,b) tables - index_nuclide = 0 - index_element = 0 - index_sab = 0 - do i = 1, n_materials - mat => materials(i) - - mat % ptr = material_pointer(i - 1) - - ! Get pointer to i-th material node - node_mat = node_mat_list(i) - - ! Check if material is depletable - if (check_for_node(node_mat, "depletable")) then - call get_node_value(node_mat, "depletable", mat % depletable) - end if - - ! Copy material name - if (check_for_node(node_mat, "name")) then - call get_node_value(node_mat, "name", mat % name) - end if - - ! Get pointer to density element - if (check_for_node(node_mat, "density")) then - node_dens = node_mat % child("density") - else - call fatal_error("Must specify density element in material " & - // trim(to_str(mat % id()))) - end if - - ! Copy units - call get_node_value(node_dens, "units", units) - - ! If the units is 'sum', then the total density of the material is taken - ! to be the sum of the atom fractions listed on the nuclides - if (units == 'sum') then - sum_density = .true. - - else if (units == 'macro') then - if (check_for_node(node_dens, "value")) then - call get_node_value(node_dens, "value", val) - else - val = ONE - end if - - ! Set density - mat % density = val - - sum_density = .false. - - else - call get_node_value(node_dens, "value", val) - - ! Check for erroneous density - sum_density = .false. - if (val <= ZERO) then - call fatal_error("Need to specify a positive density on material " & - // trim(to_str(mat % id())) // ".") - end if - - ! Adjust material density based on specified units - select case(to_lower(units)) - case ('g/cc', 'g/cm3') - mat % density = -val - case ('kg/m3') - mat % density = -0.001_8 * val - case ('atom/b-cm') - mat % density = val - case ('atom/cm3', 'atom/cc') - mat % density = 1.0e-24_8 * val - case default - call fatal_error("Unkwown units '" // trim(units) & - // "' specified on material " // trim(to_str(mat % id()))) - end select - end if - - ! Issue error if elements are provided - call get_node_list(node_mat, "element", node_ele_list) - - if (size(node_ele_list) > 0) then - call fatal_error("Unable to add an element to material " & - // trim(to_str(mat % id())) // " since the element option has & - &been removed from the xml input. Elements can only be added via & - &the Python API, which will expand elements into their natural & - &nuclides.") - end if - - ! ======================================================================= - ! READ AND PARSE TAGS - - ! Check to ensure material has at least one nuclide - if (.not. check_for_node(node_mat, "nuclide") .and. & - .not. check_for_node(node_mat, "macroscopic")) then - call fatal_error("No macroscopic data or nuclides specified on & - &material " // trim(to_str(mat % id()))) - end if - - ! Create list of macroscopic x/s based on those specified, just treat - ! them as nuclides. This is all really a facade so the user thinks they - ! are entering in macroscopic data but the code treats them the same - ! as nuclides internally. - ! Get pointer list of XML - call get_node_list(node_mat, "macroscopic", node_macro_list) - if (run_CE .and. (size(node_macro_list) > 0)) then - call fatal_error("Macroscopic can not be used in continuous-energy& - & mode!") - else if (size(node_macro_list) > 1) then - call fatal_error("Only one macroscopic object permitted per material, " & - // trim(to_str(mat % id()))) - else if (size(node_macro_list) == 1) then - - node_nuc = node_macro_list(1) - - ! Check for empty name on nuclide - if (.not. check_for_node(node_nuc, "name")) then - call fatal_error("No name specified on macroscopic data in material " & - // trim(to_str(mat % id()))) - end if - - ! store nuclide name - call get_node_value(node_nuc, "name", name) - name = trim(name) - - ! save name to list - call names % push_back(name) - - ! Set density for macroscopic data - if (units == 'macro') then - call densities % push_back(ONE) - else - call fatal_error("Units can only be macro for macroscopic data " & - // trim(name)) - end if - else - - ! Get pointer list of XML - call get_node_list(node_mat, "nuclide", node_nuc_list) - - ! Create list of nuclides based on those specified - INDIVIDUAL_NUCLIDES: do j = 1, size(node_nuc_list) - ! Combine nuclide identifier and cross section and copy into names - node_nuc = node_nuc_list(j) - - ! Check for empty name on nuclide - if (.not. check_for_node(node_nuc, "name")) then - call fatal_error("No name specified on nuclide in material " & - // trim(to_str(mat % id()))) - end if - - ! store nuclide name - call get_node_value(node_nuc, "name", name) - name = trim(name) - - ! save name to list - call names % push_back(name) - - ! Check if no atom/weight percents were specified or if both atom and - ! weight percents were specified - if (units == 'macro') then - call densities % push_back(ONE) - else - if (.not. check_for_node(node_nuc, "ao") .and. & - .not. check_for_node(node_nuc, "wo")) then - call fatal_error("No atom or weight percent specified for & - &nuclide" // trim(name)) - elseif (check_for_node(node_nuc, "ao") .and. & - check_for_node(node_nuc, "wo")) then - call fatal_error("Cannot specify both atom and weight percents & - &for a nuclide: " // trim(name)) - end if - - ! Copy atom/weight percents - if (check_for_node(node_nuc, "ao")) then - call get_node_value(node_nuc, "ao", temp_dble) - call densities % push_back(temp_dble) - else - call get_node_value(node_nuc, "wo", temp_dble) - call densities % push_back(-temp_dble) - end if - end if - end do INDIVIDUAL_NUCLIDES - end if - - ! ======================================================================== - ! COPY NUCLIDES TO ARRAYS IN MATERIAL - - ! allocate arrays in Material object - n = names % size() - mat % n_nuclides = n - allocate(mat % names(n)) - allocate(mat % nuclide(n)) - allocate(mat % atom_density(n)) - - ALL_NUCLIDES: do j = 1, mat % n_nuclides - ! Check that this nuclide is listed in the cross_sections.xml file - name = trim(names % data(j)) - if (.not. library_present(LIBRARY_NEUTRON, (name))) then - call fatal_error("Could not find nuclide " // trim(name) & - // " in cross_sections data file!") - end if - - ! If this nuclide hasn't been encountered yet, we need to add its name - ! and alias to the nuclide_dict - if (.not. nuclide_dict % has(name)) then - index_nuclide = index_nuclide + 1 - mat % nuclide(j) = index_nuclide - - call nuclide_dict % set(name, index_nuclide) - else - mat % nuclide(j) = nuclide_dict % get(name) - end if - - ! Copy name and atom/weight percent - mat % names(j) = name - mat % atom_density(j) = densities % data(j) - - end do ALL_NUCLIDES - - ! Check to make sure either all atom percents or all weight percents are - ! given - if (.not. (all(mat % atom_density >= ZERO) .or. & - all(mat % atom_density <= ZERO))) then - call fatal_error("Cannot mix atom and weight percents in material " & - // to_str(mat % id())) - end if - - ! Determine density if it is a sum value - if (sum_density) mat % density = sum(mat % atom_density) - - ! Clear lists - call names % clear() - call densities % clear() - - ! Add material to dictionary - call material_dict % set(mat % id(), i) + materials(i) % ptr = material_pointer(i - 1) end do - ! Set total number of nuclides and S(a,b) tables - n_nuclides = index_nuclide - n_elements = index_element + ! Set total number of nuclides and elements + n_nuclides = nuclides_size() + n_elements = elements_size() allocate(nuclides(n_nuclides)) ! Close materials XML file @@ -1721,82 +1481,4 @@ contains end subroutine read_mg_cross_sections_header -!=============================================================================== -! NORMALIZE_AO Normalize the nuclide atom percents -!=============================================================================== - - subroutine normalize_ao() bind(C) - integer :: i ! index in materials array - integer :: j ! index over nuclides in material - real(8) :: sum_percent ! summation - real(8) :: awr ! atomic weight ratio - real(8) :: x ! atom percent - logical :: percent_in_atom ! nuclides specified in atom percent? - logical :: density_in_atom ! density specified in atom/b-cm? - - do i = 1, size(materials) - associate (mat => materials(i)) - percent_in_atom = (mat % atom_density(1) > ZERO) - density_in_atom = (mat % density > ZERO) - - sum_percent = ZERO - do j = 1, size(mat % nuclide) - ! determine atomic weight ratio - if (run_CE) then - awr = nuclide_awr(mat % nuclide(j)) - else - awr = get_awr_c(mat % nuclide(j)) - end if - - ! if given weight percent, convert all values so that they are divided - ! by awr. thus, when a sum is done over the values, it's actually - ! sum(w/awr) - if (.not. percent_in_atom) then - mat % atom_density(j) = -mat % atom_density(j) / awr - end if - end do - - ! determine normalized atom percents. if given atom percents, this is - ! straightforward. if given weight percents, the value is w/awr and is - ! divided by sum(w/awr) - sum_percent = sum(mat % atom_density) - mat % atom_density = mat % atom_density / sum_percent - - ! Change density in g/cm^3 to atom/b-cm. Since all values are now in - ! atom percent, the sum needs to be re-evaluated as 1/sum(x*awr) - if (.not. density_in_atom) then - sum_percent = ZERO - do j = 1, mat % n_nuclides - if (run_CE) then - awr = nuclide_awr(mat % nuclide(j)) - else - awr = get_awr_c(mat % nuclide(j)) - end if - x = mat % atom_density(j) - sum_percent = sum_percent + x*awr - end do - sum_percent = ONE / sum_percent - mat%density = -mat % density * N_AVOGADRO & - / MASS_NEUTRON * sum_percent - end if - - ! Calculate nuclide atom densities - mat % atom_density = mat % density * mat % atom_density - - ! Calculate density in g/cm^3. - mat % density_gpcc = ZERO - do j = 1, mat % n_nuclides - if (run_CE) then - awr = nuclide_awr(mat % nuclide(j)) - else - awr = ONE - end if - mat % density_gpcc = mat % density_gpcc & - + mat % atom_density(j) * awr * MASS_NEUTRON / N_AVOGADRO - end do - end associate - end do - - end subroutine normalize_ao - end module input_xml diff --git a/src/material.cpp b/src/material.cpp index c877868f1..ade0f0dc0 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1,6 +1,6 @@ #include "openmc/material.h" -#include // for min, max, sort +#include // for min, max, sort, fill #include #include #include @@ -630,6 +630,15 @@ void Material::init_bremsstrahlung() } } +void Material::init_nuclide_index() +{ + mat_nuclide_index_.resize(data::nuclides.size()); + std::fill(mat_nuclide_index_.begin(), mat_nuclide_index_.end(), -1); + for (int i = 0; i < nuclide_.size(); ++i) { + mat_nuclide_index_[nuclide_[i]] = i; + } +} + void Material::calculate_xs(const Particle& p) const { // Set all material macroscopic cross sections to zero @@ -1092,6 +1101,31 @@ extern "C" { int32_t material_id(Material* mat) {return mat->id_;} + int material_nuclide(int32_t i_mat, int idx) + { + return model::materials[i_mat - 1]->nuclide_[idx - 1] + 1; + } + + int material_nuclide_size(int32_t i_mat) + { + return model::materials[i_mat - 1]->nuclide_.size(); + } + + double material_atom_density(int32_t i_mat, int idx) + { + return model::materials[i_mat - 1]->atom_density_(idx - 1); + } + + double material_density_gpcc(int32_t i_mat) + { + return model::materials[i_mat - 1]->density_gpcc_; + } + + int material_nuclide_index(int32_t i_mat, int i_nuc) + { + return model::materials[i_mat - 1]->mat_nuclide_index_[i_nuc - 1] + 1; + } + void material_calculate_xs_c(Material* mat, const Particle* p) { mat->calculate_xs(*p); diff --git a/src/material_header.F90 b/src/material_header.F90 index 103cd1b5b..844e0d731 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -16,6 +16,11 @@ module material_header private public :: free_memory_material + public :: material_nuclide + public :: material_nuclide_size + public :: material_nuclide_index + public :: material_atom_density + public :: material_density_gpcc public :: openmc_extend_materials public :: openmc_material_get_volume public :: material_pointer @@ -33,6 +38,39 @@ module material_header integer(C_INT32_T) :: id end function material_id_c + function material_nuclide(i_mat, idx) bind(C) result(nuc) + import C_INT32_T, C_INT + integer(C_INT32_T), value :: i_mat + integer(C_INT), value :: idx + integer(C_INT) :: nuc + end function + + function material_nuclide_size(i_mat) bind(C) result(n) + import C_INT32_T, C_INT + integer(C_INT32_T), value :: i_mat + integer(C_INT) :: n + end function + + function material_nuclide_index(i_mat, i_nuc) bind(C) result(idx) + import C_INT32_T, C_INT + integer(C_INT32_T), value :: i_mat + integer(C_INT), value :: i_nuc + integer(C_INT) :: idx + end function + + function material_atom_density(i_mat, idx) bind(C) result(density) + import C_INT32_T, C_INT, C_DOUBLE + integer(C_INT32_T), value :: i_mat + integer(C_INT), value :: idx + real(C_DOUBLE) :: density + end function + + function material_density_gpcc(i_mat) bind(C) result(density) + import C_INT32_T, C_DOUBLE + integer(C_INT32_T), value :: i_mat + real(C_DOUBLE) :: density + end function + subroutine extend_materials_c(n) bind(C) import C_INT32_T integer(C_INT32_T), intent(in), value :: n @@ -52,34 +90,8 @@ module material_header type, public :: Material type(C_PTR) :: ptr - character(len=104) :: name = "" ! User-defined name - integer :: n_nuclides = 0 ! number of nuclides - integer, allocatable :: nuclide(:) ! index in nuclides array - real(8) :: density ! total atom density in atom/b-cm - real(C_DOUBLE), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm - real(8) :: density_gpcc ! total density in g/cm^3 - - ! To improve performance of tallying, we store an array (direct address - ! table) that indicates for each nuclide in the global nuclides(:) array the - ! index of the corresponding nuclide in the Material % nuclide(:) array. If - ! it is not present in the material, the entry is set to zero. - integer, allocatable :: mat_nuclide_index(:) - - ! S(a,b) data - integer :: n_sab = 0 ! number of S(a,b) tables - integer, allocatable :: i_sab_nuclides(:) ! index of corresponding nuclide - integer, allocatable :: i_sab_tables(:) ! index in sab_tables - real(8), allocatable :: sab_fracs(:) ! how often to use S(a,b) - - ! Temporary names read during initialization - character(20), allocatable :: names(:) ! isotope names - character(20), allocatable :: sab_names(:) ! name of S(a,b) table - - ! Does this material contain fissionable nuclides? Is it depletable? - logical :: depletable = .false. contains procedure :: id => material_id - procedure :: init_nuclide_index => material_init_nuclide_index procedure :: calculate_xs => material_calculate_xs end type Material @@ -87,9 +99,6 @@ module material_header type(Material), public, allocatable, target :: materials(:) - ! Dictionary that maps user IDs to indices in 'materials' - type(DictIntInt), public :: material_dict - contains function material_id(this) result(id) @@ -98,28 +107,6 @@ contains id = material_id_c(this % ptr) end function material_id -!=============================================================================== -! INIT_NUCLIDE_INDEX creates a mapping from indices in the global nuclides(:) -! array to the Material % nuclides array -!=============================================================================== - - subroutine material_init_nuclide_index(this) - class(Material), intent(inout) :: this - - integer :: i - - ! Allocate nuclide index array and set to zeros - if (allocated(this % mat_nuclide_index)) & - deallocate(this % mat_nuclide_index) - allocate(this % mat_nuclide_index(n_nuclides)) - this % mat_nuclide_index(:) = 0 - - ! Assign entries in the index array - do i = 1, this % n_nuclides - this % mat_nuclide_index(this % nuclide(i)) = i - end do - end subroutine material_init_nuclide_index - !=============================================================================== ! MATERIAL_CALCULATE_XS determines the macroscopic cross sections for the ! material the particle is currently traveling through. @@ -152,7 +139,6 @@ contains call free_memory_material_c() n_materials = 0 if (allocated(materials)) deallocate(materials) - call material_dict % clear() end subroutine free_memory_material !=============================================================================== diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 0e4d326d5..de3d94007 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -943,6 +943,8 @@ set_particle_energy_bounds(int particle, double E_min, double E_max) data::energy_max[particle - 1] = E_max; } +extern "C" int nuclides_size() { return data::nuclide_map.size(); } + extern "C" void nuclide_init_grid_c(Nuclide* nuc) { nuc->init_grid(); } extern "C" double nuclide_awr(int i_nuc) { return data::nuclides[i_nuc - 1]->awr_; } diff --git a/src/photon.cpp b/src/photon.cpp index d799e1f84..52d3118dc 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -792,6 +792,8 @@ extern "C" void photon_from_hdf5(hid_t group) } } +extern "C" int elements_size() { return data::element_map.size(); } + extern "C" void photon_calculate_xs(int i_element, double E) { data::elements[i_element - 1].calculate_xs(E); diff --git a/src/simulation.F90 b/src/simulation.F90 index f00b58ed1..3de296f9d 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -2,7 +2,6 @@ module simulation use, intrinsic :: ISO_C_BINDING - use material_header, only: n_materials, materials use nuclide_header, only: micro_xs, n_nuclides use photon_header, only: micro_photon_xs, n_elements use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer @@ -20,11 +19,6 @@ contains integer :: i - ! Set up material nuclide index mapping - do i = 1, n_materials - call materials(i) % init_nuclide_index() - end do - !$omp parallel ! Allocate array for microscopic cross section cache allocate(micro_xs(n_nuclides)) @@ -46,12 +40,7 @@ contains subroutine simulation_finalize_f() bind(C) - integer :: i ! loop index - ! Free up simulation-specific memory - do i = 1, n_materials - deallocate(materials(i) % mat_nuclide_index) - end do !$omp parallel deallocate(micro_xs, micro_photon_xs, filter_matches) !$omp end parallel diff --git a/src/simulation.cpp b/src/simulation.cpp index 3591953cb..dd1a2c2a7 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -5,6 +5,7 @@ #include "openmc/container_util.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" +#include "openmc/material.h" #include "openmc/message_passing.h" #include "openmc/nuclide.h" #include "openmc/output.h" @@ -91,6 +92,11 @@ int openmc_simulation_init() // Allocate tally results arrays if they're not allocated yet allocate_tally_results(); + // Set up material nuclide index mapping + for (auto& mat : model::materials) { + mat->init_nuclide_index(); + } + // Call Fortran initialization simulation_init_f(); set_micro_xs(); @@ -144,6 +150,9 @@ int openmc_simulation_finalize() } // Deallocate Fortran variables, set tallies to inactive + for (auto& mat : model::materials) { + mat->mat_nuclide_index_.clear(); + } simulation_finalize_f(); // Increment total number of generations diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index bca3420e9..066d715e3 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -8,6 +8,7 @@ module tally use dict_header, only: EMPTY use error, only: fatal_error use geometry_header + use material_header use math, only: t_percentile use message_passing use mgxs_interface @@ -385,13 +386,13 @@ contains ! Loop over all nuclides in the current material if (p % material /= MATERIAL_VOID) then - do l = 1, materials(p % material) % n_nuclides + do l = 1, material_nuclide_size(p % material) ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) + atom_density_ = material_atom_density(p % material, l) ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) + i_nuc = material_nuclide(p % material, l) ! Accumulate the contribution from each nuclide score = score + micro_xs(i_nuc) % fission * nuclides(i_nuc) % & @@ -544,13 +545,13 @@ contains ! Loop over all nuclides in the current material if (p % material /= MATERIAL_VOID) then - do l = 1, materials(p % material) % n_nuclides + do l = 1, material_nuclide_size(p % material) ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) + atom_density_ = material_atom_density(p % material, l) ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) + i_nuc = material_nuclide(p % material, l) ! Loop over all delayed group bins and tally to them ! individually @@ -578,13 +579,13 @@ contains ! Loop over all nuclides in the current material if (p % material /= MATERIAL_VOID) then - do l = 1, materials(p % material) % n_nuclides + do l = 1, material_nuclide_size(p % material) ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) + atom_density_ = material_atom_density(p % material, l) ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) + i_nuc = material_nuclide(p % material, l) ! Accumulate the contribution from each nuclide score = score + micro_xs(i_nuc) % fission * nuclides(i_nuc) %& @@ -796,13 +797,13 @@ contains ! Loop over all nuclides in the current material if (p % material /= MATERIAL_VOID) then - do l = 1, materials(p % material) % n_nuclides + do l = 1, material_nuclide_size(p % material) ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) + atom_density_ = material_atom_density(p % material, l) ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) + i_nuc = material_nuclide(p % material, l) if (nuclides(i_nuc) % fissionable) then @@ -841,13 +842,13 @@ contains ! Loop over all nuclides in the current material if (p % material /= MATERIAL_VOID) then - do l = 1, materials(p % material) % n_nuclides + do l = 1, material_nuclide_size(p % material) ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) + atom_density_ = material_atom_density(p % material, l) ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) + i_nuc = material_nuclide(p % material, l) if (nuclides(i_nuc) % fissionable) then @@ -927,10 +928,10 @@ contains if (p % material == MATERIAL_VOID) then score = ZERO else - do l = 1, materials(p % material) % n_nuclides + do l = 1, material_nuclide_size(p % material) ! Determine atom density and index of nuclide - atom_density_ = materials(p % material) % atom_density(l) - i_nuc = materials(p % material) % nuclide(l) + atom_density_ = material_atom_density(p % material, l) + i_nuc = material_nuclide(p % material, l) ! If nuclide is fissionable, accumulate kappa fission associate(nuc => nuclides(i_nuc)) @@ -964,12 +965,12 @@ contains else score = ZERO if (p % material /= MATERIAL_VOID) then - do l = 1, materials(p % material) % n_nuclides + do l = 1, material_nuclide_size(p % material) ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) + atom_density_ = material_atom_density(p % material, l) ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) + i_nuc = material_nuclide(p % material, l) if (micro_xs(i_nuc) % elastic == CACHE_INVALID) then call nuclides(i_nuc) % calculate_elastic_xs() end if @@ -1038,9 +1039,9 @@ contains end associate else if (p % material /= MATERIAL_VOID) then - do l = 1, materials(p % material) % n_nuclides - atom_density_ = materials(p % material) % atom_density(l) - i_nuc = materials(p % material) % nuclide(l) + do l = 1, material_nuclide_size(p % material) + atom_density_ = material_atom_density(p % material, l) + i_nuc = material_nuclide(p % material, l) associate (nuc => nuclides(i_nuc)) if (score_bin == SCORE_FISS_Q_PROMPT) then @@ -1086,9 +1087,9 @@ contains score = ZERO if (p % material /= MATERIAL_VOID) then associate (mat => materials(p % material)) - do l = 1, materials(p % material) % n_nuclides - i_nuc = mat % nuclide(l) - atom_density_ = mat % atom_density(l) + do l = 1, material_nuclide_size(p % material) + i_nuc = material_nuclide(p % material, l) + atom_density_ = material_atom_density(p % material, l) score = score + micro_xs(i_nuc) % reaction(m) * atom_density_ * flux end do end associate @@ -1139,12 +1140,12 @@ contains else if (p % material /= MATERIAL_VOID) then - do l = 1, materials(p % material) % n_nuclides + do l = 1, material_nuclide_size(p % material) ! Get atom density - atom_density_ = materials(p % material) % atom_density(l) + atom_density_ = material_atom_density(p % material, l) ! Get index in nuclides array - i_nuc = materials(p % material) % nuclide(l) + i_nuc = material_nuclide(p % material, l) m = nuclides(i_nuc) % reaction_index(score_bin) if (m /= 0) then @@ -2013,12 +2014,12 @@ contains ! ========================================================================== ! SCORE ALL INDIVIDUAL NUCLIDE REACTION RATES - NUCLIDE_LOOP: do i = 1, mat % n_nuclides + NUCLIDE_LOOP: do i = 1, material_nuclide_size(p % material) ! Determine index in nuclides array and atom density for i-th nuclide in ! current material - i_nuclide = mat % nuclide(i) - atom_density = mat % atom_density(i) + i_nuclide = material_nuclide(p % material, i) + atom_density = material_atom_density(p % material, i) ! Determine score for each bin call score_general(p, t, (i_nuclide-1)*t % n_score_bins, filter_index, & @@ -2264,11 +2265,11 @@ contains mat => materials(p % material) ! Determine index of nuclide in Material % atom_density array - j = mat % mat_nuclide_index(i_nuclide) + j = material_nuclide_index(p % material, i_nuclide) if (j == 0) cycle NUCLIDE_LOOP ! Copy corresponding atom density - atom_density = mat % atom_density(j) + atom_density = material_atom_density(p % material, j) end if else atom_density = ZERO @@ -2634,11 +2635,11 @@ contains mat => materials(p % material) ! Determine index of nuclide in Material % atom_density array - j = mat % mat_nuclide_index(i_nuclide) + j = material_nuclide_index(p % material, i_nuclide) if (j == 0) cycle NUCLIDE_BIN_LOOP ! Copy corresponding atom density - atom_density = mat % atom_density(j) + atom_density = material_atom_density(p % material, j) else atom_density = ZERO end if @@ -2791,11 +2792,11 @@ contains mat => materials(p % material) ! Determine index of nuclide in Material % atom_density array - j = mat % mat_nuclide_index(i_nuclide) + j = material_nuclide_index(p % material, i_nuclide) if (j == 0) cycle NUCLIDE_BIN_LOOP ! Copy corresponding atom density - atom_density = mat % atom_density(j) + atom_density = material_atom_density(p % material, j) else atom_density = ZERO end if @@ -2993,6 +2994,7 @@ contains real(8), intent(inout) :: score integer :: l + integer :: i_nuc logical :: scoring_diff_nuclide real(8) :: flux_deriv real(8) :: dsig_s, dsig_a, dsig_f, cum_dsig @@ -3031,7 +3033,7 @@ contains SCORE_NU_FISSION) if (materials(p % material) % id() == deriv % diff_material) then score = score * (flux_deriv + ONE & - / materials(p % material) % density_gpcc) + / material_density_gpcc(p % material)) else score = score * flux_deriv end if @@ -3052,7 +3054,7 @@ contains SCORE_NU_FISSION) if (materials(p % material) % id() == deriv % diff_material) then score = score * (flux_deriv + ONE & - / materials(p % material) % density_gpcc) + / material_density_gpcc(p % material)) else score = score * flux_deriv end if @@ -3096,12 +3098,12 @@ contains .and. p % event_nuclide == deriv % diff_nuclide) then associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. - do l = 1, mat % n_nuclides - if (mat % nuclide(l) == deriv % diff_nuclide) exit + do l = 1, material_nuclide_size(p % material) + if (material_nuclide(p % material, l) == deriv % diff_nuclide) exit end do score = score * (flux_deriv & - + ONE / mat % atom_density(l)) + + ONE / material_atom_density(p % material, l)) end associate else score = score * flux_deriv @@ -3232,8 +3234,8 @@ contains micro_xs(p % event_nuclide) % total > ZERO) then associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. - do l = 1, mat % n_nuclides - if (mat % nuclide(l) == p % event_nuclide) exit + do l = 1, material_nuclide_size(p % material) + if (material_nuclide(p % material, l) == p % event_nuclide) exit end do dsig_s = ZERO @@ -3245,7 +3247,7 @@ contains end if end associate score = score * (flux_deriv & - + (dsig_s + dsig_a) * mat % atom_density(l) & + + (dsig_s + dsig_a) * material_atom_density(p % material, l) & / material_xs % total) end associate else @@ -3258,8 +3260,8 @@ contains - micro_xs(p % event_nuclide) % absorption) > ZERO) then associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. - do l = 1, mat % n_nuclides - if (mat % nuclide(l) == p % event_nuclide) exit + do l = 1, material_nuclide_size(p % material) + if (material_nuclide(p % material, l) == p % event_nuclide) exit end do dsig_s = ZERO @@ -3269,7 +3271,7 @@ contains p % sqrtkT, dsig_s, dsig_a, dsig_f) end if end associate - score = score * (flux_deriv + dsig_s * mat % atom_density(l) / & + score = score * (flux_deriv + dsig_s * material_atom_density(p % material, l) / & (material_xs % total - material_xs % absorption)) end associate else @@ -3281,8 +3283,8 @@ contains micro_xs(p % event_nuclide) % absorption > ZERO) then associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. - do l = 1, mat % n_nuclides - if (mat % nuclide(l) == p % event_nuclide) exit + do l = 1, material_nuclide_size(p % material) + if (material_nuclide(p % material, l) == p % event_nuclide) exit end do dsig_a = ZERO @@ -3292,7 +3294,7 @@ contains p % sqrtkT, dsig_s, dsig_a, dsig_f) end if end associate - score = score * (flux_deriv + dsig_a * mat % atom_density(l) & + score = score * (flux_deriv + dsig_a * material_atom_density(p % material, l) & / material_xs % absorption) end associate else @@ -3304,8 +3306,8 @@ contains micro_xs(p % event_nuclide) % fission > ZERO) then associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. - do l = 1, mat % n_nuclides - if (mat % nuclide(l) == p % event_nuclide) exit + do l = 1, material_nuclide_size(p % material) + if (material_nuclide(p % material, l) == p % event_nuclide) exit end do dsig_f = ZERO @@ -3316,7 +3318,7 @@ contains end if end associate score = score * (flux_deriv & - + dsig_f * mat % atom_density(l) / material_xs % fission) + + dsig_f * material_atom_density(p % material, l) / material_xs % fission) end associate else score = score * flux_deriv @@ -3327,8 +3329,8 @@ contains micro_xs(p % event_nuclide) % nu_fission > ZERO) then associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. - do l = 1, mat % n_nuclides - if (mat % nuclide(l) == p % event_nuclide) exit + do l = 1, material_nuclide_size(p % material) + if (material_nuclide(p % material, l) == p % event_nuclide) exit end do dsig_f = ZERO @@ -3339,7 +3341,7 @@ contains end if end associate score = score * (flux_deriv & - + dsig_f * mat % atom_density(l) / material_xs % nu_fission& + + dsig_f * material_atom_density(p % material, l) / material_xs % nu_fission& * micro_xs(p % event_nuclide) % nu_fission & / micro_xs(p % event_nuclide) % fission) end associate @@ -3365,14 +3367,15 @@ contains material_xs % total > ZERO) then cum_dsig = ZERO associate(mat => materials(p % material)) - do l = 1, mat % n_nuclides - associate (nuc => nuclides(mat % nuclide(l))) + do l = 1, material_nuclide_size(p % material) + i_nuc = material_nuclide(p % material, l) + associate (nuc => nuclides(i_nuc)) if (multipole_in_range(nuc % ptr, p % last_E) .and. & - micro_xs(mat % nuclide(l)) % total > ZERO) then + micro_xs(i_nuc) % total > ZERO) then call multipole_deriv_eval(nuc % ptr, p % last_E, & p % sqrtkT, dsig_s, dsig_a, dsig_f) cum_dsig = cum_dsig + (dsig_s + dsig_a) & - * mat % atom_density(l) + * material_atom_density(p % material, l) end if end associate end do @@ -3401,14 +3404,15 @@ contains (material_xs % total - material_xs % absorption) > ZERO) then cum_dsig = ZERO associate(mat => materials(p % material)) - do l = 1, mat % n_nuclides - associate (nuc => nuclides(mat % nuclide(l))) + do l = 1, material_nuclide_size(p % material) + i_nuc = material_nuclide(p % material, l) + associate (nuc => nuclides(i_nuc)) if (multipole_in_range(nuc % ptr, p % last_E) .and. & - (micro_xs(mat % nuclide(l)) % total & - - micro_xs(mat % nuclide(l)) % absorption) > ZERO) then + (micro_xs(i_nuc) % total & + - micro_xs(i_nuc) % absorption) > ZERO) then call multipole_deriv_eval(nuc % ptr, p % last_E, & p % sqrtkT, dsig_s, dsig_a, dsig_f) - cum_dsig = cum_dsig + dsig_s * mat % atom_density(l) + cum_dsig = cum_dsig + dsig_s * material_atom_density(p % material, l) end if end associate end do @@ -3438,13 +3442,14 @@ contains material_xs % absorption > ZERO) then cum_dsig = ZERO associate(mat => materials(p % material)) - do l = 1, mat % n_nuclides - associate (nuc => nuclides(mat % nuclide(l))) + do l = 1, material_nuclide_size(p % material) + i_nuc = material_nuclide(p % material, l) + associate (nuc => nuclides(i_nuc)) if (multipole_in_range(nuc % ptr, p % last_E) .and. & - micro_xs(mat % nuclide(l)) % absorption > ZERO) then + micro_xs(i_nuc) % absorption > ZERO) then call multipole_deriv_eval(nuc % ptr, p % last_E, & p % sqrtkT, dsig_s, dsig_a, dsig_f) - cum_dsig = cum_dsig + dsig_a * mat % atom_density(l) + cum_dsig = cum_dsig + dsig_a * material_atom_density(p % material, l) end if end associate end do @@ -3472,13 +3477,14 @@ contains material_xs % fission > ZERO) then cum_dsig = ZERO associate(mat => materials(p % material)) - do l = 1, mat % n_nuclides - associate (nuc => nuclides(mat % nuclide(l))) + do l = 1, material_nuclide_size(p % material) + i_nuc = material_nuclide(p % material, l) + associate (nuc => nuclides(i_nuc)) if (multipole_in_range(nuc % ptr, p % last_E) .and. & - micro_xs(mat % nuclide(l)) % fission > ZERO) then + micro_xs(i_nuc) % fission > ZERO) then call multipole_deriv_eval(nuc % ptr, p % last_E, & p % sqrtkT, dsig_s, dsig_a, dsig_f) - cum_dsig = cum_dsig + dsig_f * mat % atom_density(l) + cum_dsig = cum_dsig + dsig_f * material_atom_density(p % material, l) end if end associate end do @@ -3506,15 +3512,16 @@ contains material_xs % nu_fission > ZERO) then cum_dsig = ZERO associate(mat => materials(p % material)) - do l = 1, mat % n_nuclides - associate (nuc => nuclides(mat % nuclide(l))) + do l = 1, material_nuclide_size(p % material) + i_nuc = material_nuclide(p % material, l) + associate (nuc => nuclides(i_nuc)) if (multipole_in_range(nuc % ptr, p % last_E) .and. & - micro_xs(mat % nuclide(l)) % nu_fission > ZERO) then + micro_xs(i_nuc) % nu_fission > ZERO) then call multipole_deriv_eval(nuc % ptr, p % last_E, & p % sqrtkT, dsig_s, dsig_a, dsig_f) - cum_dsig = cum_dsig + dsig_f * mat % atom_density(l) & - * micro_xs(mat % nuclide(l)) % nu_fission & - / micro_xs(mat % nuclide(l)) % fission + cum_dsig = cum_dsig + dsig_f * material_atom_density(p % material, l) & + * micro_xs(i_nuc) % nu_fission & + / micro_xs(i_nuc) % fission end if end associate end do @@ -3575,7 +3582,7 @@ contains ! (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist ! (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist deriv % flux_deriv = deriv % flux_deriv & - - distance * material_xs % total / mat % density_gpcc + - distance * material_xs % total / material_density_gpcc(p % material) end if end associate @@ -3593,8 +3600,8 @@ contains case (DIFF_TEMPERATURE) associate (mat => materials(p % material)) if (mat % id() == deriv % diff_material) then - do l=1, mat % n_nuclides - associate (nuc => nuclides(mat % nuclide(l))) + do l=1, material_nuclide_size(p % material) + associate (nuc => nuclides(material_nuclide(p % material, l))) if (multipole_in_range(nuc % ptr, p % E)) then ! phi is proportional to e^(-Sigma_tot * dist) ! (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist @@ -3602,7 +3609,7 @@ contains call multipole_deriv_eval(nuc % ptr, p % E, & p % sqrtkT, dsig_s, dsig_a, dsig_f) deriv % flux_deriv = deriv % flux_deriv & - - distance * (dsig_s + dsig_a) * mat % atom_density(l) + - distance * (dsig_s + dsig_a) * material_atom_density(p % material, l) end if end associate end do @@ -3632,7 +3639,7 @@ contains subroutine score_collision_derivative(p) type(Particle), intent(in) :: p - integer :: i, j, l + integer :: i, j, l, i_nuc real(8) :: dsig_s, dsig_a, dsig_f ! A void material cannot be perturbed so it will not affect flux derivatives @@ -3649,7 +3656,7 @@ contains ! (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s ! (1 / phi) * (d_phi / d_rho) = 1 / rho deriv % flux_deriv = deriv % flux_deriv & - + ONE / mat % density_gpcc + + ONE / material_density_gpcc(p % material) end if end associate @@ -3658,11 +3665,11 @@ contains if (mat % id() == deriv % diff_material & .and. p % event_nuclide == deriv % diff_nuclide) then ! Find the index in this material for the diff_nuclide. - do j = 1, mat % n_nuclides - if (mat % nuclide(j) == deriv % diff_nuclide) exit + do j = 1, material_nuclide_size(p % material) + if (material_nuclide(p % material, j) == deriv % diff_nuclide) exit end do ! Make sure we found the nuclide. - if (mat % nuclide(j) /= deriv % diff_nuclide) then + if (material_nuclide(p % material, j) /= deriv % diff_nuclide) then call fatal_error("Couldn't find the right nuclide.") end if ! phi is proportional to Sigma_s @@ -3670,16 +3677,17 @@ contains ! (1 / phi) * (d_phi / d_N) = sigma_s / Sigma_s ! (1 / phi) * (d_phi / d_N) = 1 / N deriv % flux_deriv = deriv % flux_deriv & - + ONE / mat % atom_density(j) + + ONE / material_atom_density(p % material, j) end if end associate case (DIFF_TEMPERATURE) associate (mat => materials(p % material)) if (mat % id() == deriv % diff_material) then - do l=1, mat % n_nuclides - associate (nuc => nuclides(mat % nuclide(l))) - if (mat % nuclide(l) == p % event_nuclide .and. & + do l=1, material_nuclide_size(p % material) + i_nuc = material_nuclide(p % material, l) + associate (nuc => nuclides(i_nuc)) + if (i_nuc == p % event_nuclide .and. & multipole_in_range(nuc % ptr, p % last_E)) then ! phi is proportional to Sigma_s ! (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s @@ -3687,8 +3695,8 @@ contains call multipole_deriv_eval(nuc % ptr, p % last_E, & p % sqrtkT, dsig_s, dsig_a, dsig_f) deriv % flux_deriv = deriv % flux_deriv + dsig_s& - / (micro_xs(mat % nuclide(l)) % total & - - micro_xs(mat % nuclide(l)) % absorption) + / (micro_xs(i_nuc) % total & + - micro_xs(i_nuc) % absorption) ! Note that this is an approximation! The real scattering ! cross section is Sigma_s(E'->E, uvw'->uvw) = ! Sigma_s(E') * P(E'->E, uvw'->uvw). We are assuming that diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index 5379a3b2e..ff079c658 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -13,7 +13,7 @@ module volume_calc use hdf5_interface, only: file_open, file_close, write_attribute, & create_group, close_group, write_dataset, HID_T use output, only: header, time_stamp - use material_header, only: materials + use material_header use message_passing use nuclide_header, only: nuclides use particle_header @@ -328,16 +328,14 @@ contains i_material = master_indices(i_domain) % data(j) if (i_material == MATERIAL_VOID) cycle - associate (mat => materials(i_material)) - do k = 1, size(mat % nuclide) - ! Accumulate nuclide density - i_nuclide = mat % nuclide(k) - atoms(1, i_nuclide) = atoms(1, i_nuclide) + & - mat % atom_density(k) * f - atoms(2, i_nuclide) = atoms(2, i_nuclide) + & - mat % atom_density(k)**2 * var_f - end do - end associate + do k = 1, material_nuclide_size(i_material) + ! Accumulate nuclide density + i_nuclide = material_nuclide(i_material, k) + atoms(1, i_nuclide) = atoms(1, i_nuclide) + & + material_atom_density(i_material, k) * f + atoms(2, i_nuclide) = atoms(2, i_nuclide) + & + material_atom_density(i_material, k)**2 * var_f + end do end do ! Determine volume