crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Remove most of remaining components of Material type on Fortran side

Browse source
This commit is contained in:
Paul Romano 2019-01-25 16:00:00 -06:00
parent 223a15b52f
commit 63aa3bf106
11 changed files with 225 additions and 510 deletions
Download patch file Download diff file Expand all files Collapse all files

9
include/openmc/material.h
View file

@ -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<bool> 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<int> mat_nuclide_index_;
// Thermal scattering tables
std::vector<ThermalTable> thermal_tables_;

4
src/initialize.cpp
View file

@ -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();

350
src/input_xml.F90
View file

@ -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 <nuclide> 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 <macroscopic>
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 <nuclide>
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

36
src/material.cpp
View file

@ -1,6 +1,6 @@
#include "openmc/material.h"
#include <algorithm> // for min, max, sort
#include <algorithm> // for min, max, sort, fill
#include <cmath>
#include <iterator>
#include <string>
@ -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);

90
src/material_header.F90
View file

@ -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
!===============================================================================

2
src/nuclide.cpp
View file

@ -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_; }

2
src/photon.cpp
View file

@ -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);

11
src/simulation.F90
View file

@ -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

9
src/simulation.cpp
View file

@ -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

202
src/tallies/tally.F90
View file

@ -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

20
src/volume_calc.F90
View file

@ -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