From fe9f774d31fb271bb72bb87a691e65c3616d50b2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 22 Jan 2019 16:11:39 -0600 Subject: [PATCH] Add normalize_density, set_density, set_densities on C++ side --- include/openmc/material.h | 19 +- include/openmc/simulation.h | 2 +- src/input_xml.F90 | 18 +- src/material.cpp | 186 ++++++++++++++-- src/material_header.F90 | 417 ------------------------------------ src/simulation_header.F90 | 2 - 6 files changed, 191 insertions(+), 453 deletions(-) diff --git a/include/openmc/material.h b/include/openmc/material.h index 2ea9fdb86..eda102d22 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -48,12 +48,16 @@ public: // Methods void calculate_xs(const Particle& p) const; - //! Initialize bremsstrahlung data - void init_bremsstrahlung(); - - //! Initialize thermal scattering mapping + //! Assign thermal scattering tables to specific nuclides within the material + //! so the code knows when to apply bound thermal scattering data void init_thermal(); + //! Finalize the material, assigning tables, normalize density, etc. + void finalize(); + + //! Set total density of the material + int set_density(double density, std::string units); + // Data int32_t id_; //!< Unique ID std::string name_; //!< Name of material @@ -61,6 +65,7 @@ public: std::vector element_; //!< Indices in elements vector xt::xtensor atom_density_; //!< Nuclide atom density in [atom/b-cm] double density_; //!< Total atom density in [atom/b-cm] + double density_gpcc_; //!< Total atom density in [g/cm^3] double volume_ {-1.0}; //!< Volume in [cm^3] bool fissionable_ {false}; //!< Does this material contain fissionable nuclides bool depletable_ {false}; //!< Is the material depletable? @@ -78,6 +83,12 @@ public: std::unique_ptr ttb_; private: + //! Initialize bremsstrahlung data + void init_bremsstrahlung(); + + //! Normalize density + void normalize_density(); + void calculate_neutron_xs(const Particle& p) const; void calculate_photon_xs(const Particle& p) const; }; diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index cba0a2b71..f9e1f7c96 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -30,7 +30,7 @@ extern "C" double keff_std; //!< standard deviation of average k extern "C" double k_col_abs; //!< sum over batches of k_collision * k_absorption extern "C" double k_col_tra; //!< sum over batches of k_collision * k_tracklength extern "C" double k_abs_tra; //!< sum over batches of k_absorption * k_tracklength -extern "C" double log_spacing; //!< lethargy spacing for energy grid searches +extern double log_spacing; //!< lethargy spacing for energy grid searches extern "C" int n_lost_particles; //!< cumulative number of lost particles extern "C" bool need_depletion_rx; //!< need to calculate depletion rx? extern "C" int restart_batch; //!< batch at which a restart job resumed diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 5a2c4263b..bcca004cf 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2113,22 +2113,9 @@ contains ! Read multipole file into the appropriate entry on the nuclides array if (temperature_multipole) call read_multipole_data(i_nuclide) end if - - ! Check if material is fissionable - if (nuclides(materials(i) % nuclide(j)) % fissionable) then - call materials(i) % set_fissionable(.true.) - end if end do end do - call read_ce_cross_sections_c() - - ! Set up logarithmic grid for nuclides - do i = 1, size(nuclides) - call nuclides(i) % init_grid() - end do - log_spacing = log(energy_max(NEUTRON)/energy_min(NEUTRON)) / n_log_bins - do i = 1, size(materials) ! Skip materials with no S(a,b) tables if (.not. allocated(materials(i) % sab_names)) cycle @@ -2165,11 +2152,10 @@ contains call already_read % add(name) end if end do - - ! Associate S(a,b) tables with specific nuclides - call materials(i) % assign_sab_tables() end do + call read_ce_cross_sections_c() + ! Show which nuclide results in lowest energy for neutron transport do i = 1, size(nuclides) ! If a nuclide is present in a material that's not used in the model, its diff --git a/src/material.cpp b/src/material.cpp index 6557135c6..48026ca63 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -14,6 +14,7 @@ #include "openmc/container_util.h" #include "openmc/error.h" #include "openmc/math_functions.h" +#include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" #include "openmc/photon.h" #include "openmc/search.h" @@ -325,6 +326,77 @@ Material::Material(pugi::xml_node node) } } +void Material::finalize() +{ + // Set fissionable if any nuclide is fissionable + for (const auto& i_nuc : nuclide_) { + if (data::nuclides[i_nuc]->fissionable_) { + fissionable_ = true; + break; + } + } + + // Generate material bremsstrahlung data for electrons and positrons + if (settings::photon_transport && settings::electron_treatment == ELECTRON_TTB) { + this->init_bremsstrahlung(); + } + + // Assign thermal scattering tables + this->init_thermal(); + + // Normalize density + this->normalize_density(); +} + +void Material::normalize_density() +{ + bool percent_in_atom = (atom_density_(1) > 0.0); + bool density_in_atom = (density_ > 0.0); + + double sum_percent = 0.0; + for (int i = 0; i < nuclide_.size(); ++i) { + // determine atomic weight ratio + int i_nuc = nuclide_[i]; + double awr = settings::run_CE ? + data::nuclides[i_nuc]->awr_ : data::nuclides_MG[i_nuc].awr; + + // 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 (!percent_in_atom) atom_density_(i) = -atom_density_(i) / awr; + } + + // 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) + atom_density_ /= xt::sum(atom_density_)(); + + // 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 (!density_in_atom) { + sum_percent = 0.0; + for (int i = 0; i < nuclide_.size(); ++i) { + int i_nuc = nuclide_[i]; + double awr = settings::run_CE ? + data::nuclides[i_nuc]->awr_ : data::nuclides_MG[i_nuc].awr; + sum_percent += atom_density_(i)*awr; + } + sum_percent = 1.0 / sum_percent; + density_ *= -N_AVOGADRO / MASS_NEUTRON * sum_percent; + } + + // Calculate nuclide atom densities + atom_density_ *= density_; + + // Calculate density in g/cm^3. + double density_gpcc_ = 0.0; + for (int i = 0; i < nuclide_.size(); ++i) { + int i_nuc = nuclide_[i]; + double awr = settings::run_CE ? data::nuclides[i_nuc]->awr_ : 1.0; + density_gpcc_ += atom_density_(i) * awr * MASS_NEUTRON / N_AVOGADRO; + } +} + void Material::init_thermal() { std::vector tables; @@ -598,11 +670,11 @@ void Material::calculate_neutron_xs(const Particle& p) const // ====================================================================== // CHECK FOR S(A,B) TABLE - int i_sab = 0; + int i_sab = C_NONE; double sab_frac = 0.0; // Check if this nuclide matches one of the S(a,b) tables specified. - // This relies on i_sab_nuclides being in sorted order + // This relies on thermal_tables_ being sorted by .index_nuclide if (check_sab) { const auto& sab {thermal_tables_[j]}; if (i == sab.index_nuclide) { @@ -612,13 +684,13 @@ void Material::calculate_neutron_xs(const Particle& p) const // If particle energy is greater than the highest energy for the // S(a,b) table, then don't use the S(a,b) table - if (p.E > data::thermal_scatt[i_sab]->threshold()) i_sab = 0; + if (p.E > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE; // Increment position in i_sab_nuclides ++j; // Don't check for S(a,b) tables if there are no more left - if (j == thermal_tables_.size() - 1) check_sab = false; + if (j == thermal_tables_.size()) check_sab = false; } } @@ -689,6 +761,47 @@ void Material::calculate_photon_xs(const Particle& p) const } } +int Material::set_density(double density, std::string units) +{ + if (nuclide_.empty()) { + set_errmsg("No nuclides exist in material yet."); + return OPENMC_E_ALLOCATE; + } + + if (units == "atom/b-cm") { + // Set total density based on value provided + density_ = density; + + // Determine normalized atom percents + double sum_percent = xt::sum(atom_density_)(); + atom_density_ /= sum_percent; + + // Recalculate nuclide atom densities based on given density + atom_density_ *= density; + + // Calculate density in g/cm^3. + density_gpcc_ = 0.0; + for (int i = 0; i < nuclide_.size(); ++i) { + int i_nuc = nuclide_[i]; + double awr = data::nuclides[i_nuc]->awr_; + density_gpcc_ += atom_density_(i) * awr * MASS_NEUTRON / N_AVOGADRO; + } + } else if (units == "g/cm3" || units == "g/cc") { + // Determine factor by which to change densities + double previous_density_gpcc = density_gpcc_; + double f = density / previous_density_gpcc; + + // Update densities + density_gpcc_ = density; + density_ *= f; + atom_density_ *= f; + } else { + set_errmsg("Invalid units '" + units + "' specified."); + return OPENMC_E_INVALID_ARGUMENT; + } + return 0; +} + //============================================================================== // Non-method functions //============================================================================== @@ -719,12 +832,17 @@ read_materials(pugi::xml_node* node) extern "C" void read_ce_cross_sections_c() { for (auto& mat : model::materials) { - // Generate material bremsstrahlung data for electrons and positrons - if (settings::photon_transport && settings::electron_treatment == ELECTRON_TTB) { - mat->init_bremsstrahlung(); - } + mat->finalize(); } + // Set up logarithmic grid for nuclides + for (auto& nuc : data::nuclides) { + nuc->init_grid(); + } + int neutron = static_cast(ParticleType::neutron) - 1; + simulation::log_spacing = std::log(data::energy_max[neutron] / + data::energy_min[neutron]) / settings::n_log_bins; + if (settings::photon_transport && settings::electron_treatment == ELECTRON_TTB) { // Determine if minimum/maximum energy for bremsstrahlung is greater/less // than the current minimum/maximum @@ -765,6 +883,53 @@ openmc_material_get_volume(int32_t index, double* volume) } } +extern "C" int +openmc_material_set_density(int32_t index, double density, const char* units) +{ + if (index >= 1 && index <= model::materials.size()) { + return model::materials[index - 1]->set_density(density, units); + } else { + set_errmsg("Index in materials array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } +} + +extern "C" int +openmc_material_set_densities(int32_t index, int n, const char** name, const double* density) +{ + if (index >= 1 && index <= model::materials.size()) { + // TODO: off-by-one + auto& mat {model::materials[index - 1]}; + if (n != mat->nuclide_.size()) { + mat->nuclide_.resize(n); + mat->atom_density_ = xt::zeros({n}); + } + + double sum_density = 0.0; + for (int i = 0; i < n; ++i) { + std::string nuc {name[i]}; + if (data::nuclide_map.find(nuc) == data::nuclide_map.end()) { + int err = openmc_load_nuclide(nuc.c_str()); + if (err < 0) return err; + } + + mat->nuclide_[i] = data::nuclide_map[nuc]; + mat->atom_density_(i) = density[i]; + sum_density += density[i]; + } + + // Set total density to the sum of the vector + + int err = mat->set_density(sum_density, "atom/b-cm"); + + // Assign S(a,b) tables + mat->init_thermal(); + } else { + set_errmsg("Index in materials array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } +} + extern "C" int openmc_material_set_volume(int32_t index, double volume) { @@ -806,11 +971,6 @@ extern "C" { mat->fissionable_ = fissionable; } - void material_init_bremsstrahlung(Material* mat) - { - mat->init_bremsstrahlung(); - } - 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 75c3e3a7c..f1cc43602 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -9,7 +9,6 @@ module material_header use particle_header, only: Particle use photon_header use sab_header - use simulation_header, only: log_spacing use stl_vector, only: VectorReal, VectorInt use string, only: to_str, to_f_string, to_c_string @@ -23,8 +22,6 @@ module material_header public :: openmc_material_get_id public :: openmc_material_get_densities public :: openmc_material_get_volume - public :: openmc_material_set_density - public :: openmc_material_set_densities public :: openmc_material_set_id public :: material_pointer @@ -118,12 +115,8 @@ module material_header procedure :: set_id => material_set_id procedure :: fissionable => material_fissionable procedure :: set_fissionable => material_set_fissionable - procedure :: set_density => material_set_density procedure :: init_nuclide_index => material_init_nuclide_index - procedure :: assign_sab_tables => material_assign_sab_tables procedure :: calculate_xs => material_calculate_xs - procedure, private :: calculate_neutron_xs - procedure, private :: calculate_photon_xs end type Material integer(C_INT32_T), public, bind(C) :: n_materials ! # of materials @@ -135,10 +128,6 @@ module material_header contains -!=============================================================================== -! MATERIAL_SET_DENSITY sets the total density of a material in atom/b-cm. -!=============================================================================== - function material_id(this) result(id) class(Material), intent(in) :: this integer(C_INT32_T) :: id @@ -164,58 +153,6 @@ contains call material_set_fissionable_c(this % ptr, logical(fissionable, C_BOOL)) end subroutine material_set_fissionable - function material_set_density(this, density, units) result(err) - class(Material), intent(inout) :: this - real(8), intent(in) :: density - character(*), intent(in) :: units - integer :: err - - integer :: i - real(8) :: sum_percent - real(8) :: awr - real(8) :: previous_density_gpcc - real(8) :: f - - if (allocated(this % atom_density)) then - err = 0 - select case (units) - case ('atom/b-cm') - ! Set total density based on value provided - this % density = density - - ! Determine normalized atom percents - sum_percent = sum(this % atom_density) - this % atom_density(:) = this % atom_density / sum_percent - - ! Recalculate nuclide atom densities based on given density - this % atom_density(:) = density * this % atom_density - - ! Calculate density in g/cm^3. - this % density_gpcc = ZERO - do i = 1, this % n_nuclides - awr = nuclides(this % nuclide(i)) % awr - this % density_gpcc = this % density_gpcc & - + this % atom_density(i) * awr * MASS_NEUTRON / N_AVOGADRO - end do - case ('g/cm3', 'g/cc') - ! Determine factor by which to change densities - previous_density_gpcc = this % density_gpcc - f = density / previous_density_gpcc - - ! Update densities - this % density_gpcc = density - this % density = f * this % density - this % atom_density(:) = f * this % atom_density(:) - case default - err = E_INVALID_ARGUMENT - call set_errmsg("Invalid units '" // trim(units) // "' specified.") - end select - else - err = E_ALLOCATE - call set_errmsg("Material atom density array hasn't been allocated.") - end if - end function material_set_density - !=============================================================================== ! INIT_NUCLIDE_INDEX creates a mapping from indices in the global nuclides(:) ! array to the Material % nuclides array @@ -238,119 +175,6 @@ contains end do end subroutine material_init_nuclide_index -!=============================================================================== -! ASSIGN_SAB_TABLES assigns S(alpha,beta) tables to specific nuclides within -! materials so the code knows when to apply bound thermal scattering data -!=============================================================================== - - subroutine material_assign_sab_tables(this) - class(Material), intent(inout) :: this - - integer :: j ! index over nuclides in material - integer :: k ! index over S(a,b) tables in material - integer :: m ! position for sorting - integer :: i_sab - integer :: temp_nuclide ! temporary value for sorting - integer :: temp_table ! temporary value for sorting - real(8) :: temp_frac ! temporary value for sorting - logical :: found - type(VectorInt) :: i_sab_tables - type(VectorInt) :: i_sab_nuclides - type(VectorReal) :: sab_fracs - - if (.not. allocated(this % i_sab_tables)) return - - ASSIGN_SAB: do k = 1, size(this % i_sab_tables) - ! In order to know which nuclide the S(a,b) table applies to, we need - ! to search through the list of nuclides for one which has a matching - ! name - found = .false. - i_sab = this % i_sab_tables(k) - FIND_NUCLIDE: do j = 1, size(this % nuclide) - if (sab_has_nuclide(i_sab, to_c_string(nuclides(this % nuclide(j)) % name))) then - call i_sab_tables % push_back(i_sab) - call i_sab_nuclides % push_back(j) - call sab_fracs % push_back(this % sab_fracs(k)) - found = .true. - end if - end do FIND_NUCLIDE - - ! Check to make sure S(a,b) table matched a nuclide - if (.not. found) then - call fatal_error("S(a,b) table " // trim(this % & - sab_names(k)) // " did not match any nuclide on material " & - // trim(to_str(this % id()))) - end if - end do ASSIGN_SAB - - ! Make sure each nuclide only appears in one table. - do j = 1, i_sab_nuclides % size() - do k = j+1, i_sab_nuclides % size() - if (i_sab_nuclides % data(j) == i_sab_nuclides % data(k)) then - call fatal_error(trim( & - nuclides(this % nuclide(i_sab_nuclides % data(j))) % name) & - // " in material " // trim(to_str(this % id())) // " was found & - &in multiple S(a,b) tables. Each nuclide can only appear in & - &one S(a,b) table per material.") - end if - end do - end do - - ! Update i_sab_tables and i_sab_nuclides - deallocate(this % i_sab_tables) - deallocate(this % sab_fracs) - if (allocated(this % i_sab_nuclides)) deallocate(this % i_sab_nuclides) - m = i_sab_tables % size() - allocate(this % i_sab_tables(m)) - allocate(this % i_sab_nuclides(m)) - allocate(this % sab_fracs(m)) - this % i_sab_tables(:) = i_sab_tables % data(1:m) - this % i_sab_nuclides(:) = i_sab_nuclides % data(1:m) - this % sab_fracs(:) = sab_fracs % data(1:m) - - ! Clear entries in vectors for next material - call i_sab_tables % clear() - call i_sab_nuclides % clear() - call sab_fracs % clear() - - ! If there are multiple S(a,b) tables, we need to make sure that the - ! entries in i_sab_nuclides are sorted or else they won't be applied - ! correctly in the cross_section module. The algorithm here is a simple - ! insertion sort -- don't need anything fancy! - - if (size(this % i_sab_tables) > 1) then - SORT_SAB: do k = 2, size(this % i_sab_tables) - ! Save value to move - m = k - temp_nuclide = this % i_sab_nuclides(k) - temp_table = this % i_sab_tables(k) - temp_frac = this % i_sab_tables(k) - - MOVE_OVER: do - ! Check if insertion value is greater than (m-1)th value - if (temp_nuclide >= this % i_sab_nuclides(m-1)) exit - - ! Move values over until hitting one that's not larger - this % i_sab_nuclides(m) = this % i_sab_nuclides(m-1) - this % i_sab_tables(m) = this % i_sab_tables(m-1) - this % sab_fracs(m) = this % sab_fracs(m-1) - m = m - 1 - - ! Exit if we've reached the beginning of the list - if (m == 1) exit - end do MOVE_OVER - - ! Put the original value into its new position - this % i_sab_nuclides(m) = temp_nuclide - this % i_sab_tables(m) = temp_table - this % sab_fracs(m) = temp_frac - end do SORT_SAB - end if - - ! Deallocate temporary arrays for names of nuclides and S(a,b) tables - if (allocated(this % names)) deallocate(this % names) - end subroutine material_assign_sab_tables - !=============================================================================== ! MATERIAL_CALCULATE_XS determines the macroscopic cross sections for the ! material the particle is currently traveling through. @@ -371,172 +195,6 @@ contains call material_calculate_xs_c(this % ptr, p) end subroutine material_calculate_xs -!=============================================================================== -! CALCULATE_NEUTRON_XS determines the neutron cross section for the material the -! particle is traveling through -!=============================================================================== - - subroutine calculate_neutron_xs(this, p) - class(Material), intent(in) :: this - type(Particle), intent(in) :: p - - integer :: i ! loop index over nuclides - integer :: i_nuclide ! index into nuclides array - integer :: i_sab ! index into sab_tables array - integer :: j ! index in this % i_sab_nuclides - integer :: i_grid ! index into logarithmic mapping array or material - ! union grid - real(8) :: atom_density ! atom density of a nuclide - real(8) :: sab_frac ! fraction of atoms affected by S(a,b) - logical :: check_sab ! should we check for S(a,b) table? - - ! Find energy index on energy grid - i_grid = int(log(p % E/energy_min(NEUTRON))/log_spacing) - - ! Determine if this material has S(a,b) tables - check_sab = (this % n_sab > 0) - - ! Initialize position in i_sab_nuclides - j = 1 - - ! Add contribution from each nuclide in material - do i = 1, this % n_nuclides - ! ====================================================================== - ! CHECK FOR S(A,B) TABLE - - i_sab = 0 - sab_frac = ZERO - - ! Check if this nuclide matches one of the S(a,b) tables specified. - ! This relies on i_sab_nuclides being in sorted order - if (check_sab) then - if (i == this % i_sab_nuclides(j)) then - ! Get index in sab_tables - i_sab = this % i_sab_tables(j) - sab_frac = this % sab_fracs(j) - - ! If particle energy is greater than the highest energy for the - ! S(a,b) table, then don't use the S(a,b) table - if (p % E > sab_threshold(i_sab)) then - i_sab = 0 - end if - - ! Increment position in i_sab_nuclides - j = j + 1 - - ! Don't check for S(a,b) tables if there are no more left - if (j > size(this % i_sab_tables)) check_sab = .false. - end if - end if - - ! ====================================================================== - ! CALCULATE MICROSCOPIC CROSS SECTION - - ! Determine microscopic cross sections for this nuclide - i_nuclide = this % nuclide(i) - - ! Calculate microscopic cross section for this nuclide - if (p % E /= micro_xs(i_nuclide) % last_E & - .or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT & - .or. i_sab /= micro_xs(i_nuclide) % index_sab + 1 & - .or. sab_frac /= micro_xs(i_nuclide) % sab_frac) then - call nuclides(i_nuclide) % calculate_xs(i_sab, p % E, i_grid, & - p % sqrtkT, sab_frac) - end if - - ! ====================================================================== - ! ADD TO MACROSCOPIC CROSS SECTION - - ! Copy atom density of nuclide in material - atom_density = this % atom_density(i) - - ! Add contributions to material macroscopic total cross section - material_xs % total = material_xs % total + & - atom_density * micro_xs(i_nuclide) % total - - ! Add contributions to material macroscopic absorption cross section - material_xs % absorption = material_xs % absorption + & - atom_density * micro_xs(i_nuclide) % absorption - - ! Add contributions to material macroscopic fission cross section - material_xs % fission = material_xs % fission + & - atom_density * micro_xs(i_nuclide) % fission - - ! Add contributions to material macroscopic nu-fission cross section - material_xs % nu_fission = material_xs % nu_fission + & - atom_density * micro_xs(i_nuclide) % nu_fission - end do - - end subroutine calculate_neutron_xs - -!=============================================================================== -! CALCULATE_PHOTON_XS determines the macroscopic photon cross sections for the -! material the particle is currently traveling through. -!=============================================================================== - - subroutine calculate_photon_xs(this, p) - class(Material), intent(in) :: this - type(Particle), intent(in) :: p - - integer :: i ! loop index over nuclides - integer :: i_element ! index into elements array - real(8) :: atom_density ! atom density of a nuclide - - interface - subroutine photon_calculate_xs(i_element, E) bind(C) - import C_INT, C_DOUBLE - integer(C_INT), value :: i_element - real(C_DOUBLE), value :: E - end subroutine - end interface - - material_xs % coherent = ZERO - material_xs % incoherent = ZERO - material_xs % photoelectric = ZERO - material_xs % pair_production = ZERO - - ! Add contribution from each nuclide in material - do i = 1, this % n_nuclides - ! ======================================================================== - ! CALCULATE MICROSCOPIC CROSS SECTION - - ! Determine microscopic cross sections for this nuclide - i_element = this % element(i) - - ! Calculate microscopic cross section for this nuclide - if (p % E /= micro_photon_xs(i_element) % last_E) then - call photon_calculate_xs(i_element, p % E) - end if - - ! ======================================================================== - ! ADD TO MACROSCOPIC CROSS SECTION - - ! Copy atom density of nuclide in material - atom_density = this % atom_density(i) - - ! Add contributions to material macroscopic total cross section - material_xs % total = material_xs % total + & - atom_density * micro_photon_xs(i_element) % total - - ! Add contributions to material macroscopic coherent cross section - material_xs % coherent = material_xs % coherent + & - atom_density * micro_photon_xs(i_element) % coherent - - ! Add contributions to material macroscopic incoherent cross section - material_xs % incoherent = material_xs % incoherent + & - atom_density * micro_photon_xs(i_element) % incoherent - - ! Add contributions to material macroscopic photoelectric cross section - material_xs % photoelectric = material_xs % photoelectric + & - atom_density * micro_photon_xs(i_element) % photoelectric - - ! Add contributions to material macroscopic pair production cross section - material_xs % pair_production = material_xs % pair_production + & - atom_density * micro_photon_xs(i_element) % pair_production - end do - - end subroutine calculate_photon_xs - !=============================================================================== ! FREE_MEMORY_MATERIAL deallocates global arrays defined in this module !=============================================================================== @@ -757,81 +415,6 @@ contains end if end function openmc_material_set_id - - function openmc_material_set_density(index, density, units) result(err) bind(C) - ! Set the total density of a material - integer(C_INT32_T), value, intent(in) :: index - real(C_DOUBLE), value, intent(in) :: density - character(kind=C_CHAR), intent(in) :: units(*) - integer(C_INT) :: err - - character(:), allocatable :: units_ - - ! Convert C string to Fortran string - units_ = to_f_string(units) - - err = E_UNASSIGNED - if (index >= 1 .and. index <= size(materials)) then - associate (m => materials(index)) - err = m % set_density(density, units_) - end associate - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in materials array is out of bounds.") - end if - end function openmc_material_set_density - - - function openmc_material_set_densities(index, n, name, density) result(err) bind(C) - ! Sets the densities for a list of nuclides in a material. If the nuclides - ! don't already exist in the material, they will be added - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT), value, intent(in) :: n - type(C_PTR), intent(in) :: name(n) - real(C_DOUBLE), intent(in) :: density(n) - integer(C_INT) :: err - - integer :: i - integer :: stat - character(C_CHAR), pointer :: string(:) - character(len=:, kind=C_CHAR), allocatable :: name_ - - if (index >= 1 .and. index <= size(materials)) then - associate (m => materials(index)) - ! If nuclide/density arrays are not correct size, reallocate - if (n /= m % n_nuclides) then - deallocate(m % nuclide, m % atom_density, STAT=stat) - allocate(m % nuclide(n), m % atom_density(n)) - end if - - do i = 1, n - ! Convert C string to Fortran string - call c_f_pointer(name(i), string, [10]) - name_ = to_f_string(string) - - if (.not. nuclide_dict % has(name_)) then - err = openmc_load_nuclide(string) - if (err < 0) return - end if - - m % nuclide(i) = nuclide_dict % get(name_) - m % atom_density(i) = density(i) - end do - m % n_nuclides = n - - ! Set total density to the sum of the vector - err = m % set_density(sum(density), 'atom/b-cm') - - ! Assign S(a,b) tables - call m % assign_sab_tables() - end associate - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in materials array is out of bounds.") - end if - - end function openmc_material_set_densities - !=============================================================================== ! Fortran compatibility !=============================================================================== diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 1df4b0a4e..2643e3f28 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -15,8 +15,6 @@ module simulation_header ! Number of lost particles integer(C_INT), bind(C) :: n_lost_particles - real(C_DOUBLE), bind(C) :: log_spacing ! spacing on logarithmic grid - ! ============================================================================ ! SIMULATION VARIABLES