From e3e388bb6d16ef6c88f1ee50290a72c720565889 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Fri, 12 Oct 2018 10:40:35 -0400 Subject: [PATCH] Completed replacement of physics_mg.F90 with C++ --- CMakeLists.txt | 1 - include/openmc/material.h | 1 + include/openmc/physics_mg.h | 24 +++++-- src/input_xml.F90 | 2 +- src/material.cpp | 7 ++ src/material_header.F90 | 31 ++++++++- src/mgxs_data.F90 | 3 +- src/physics_mg.F90 | 126 ------------------------------------ src/physics_mg.cpp | 67 ++++++++++++++++++- src/tracking.F90 | 18 +++++- 10 files changed, 140 insertions(+), 140 deletions(-) delete mode 100644 src/physics_mg.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 135a724311..b72bb57eb9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -345,7 +345,6 @@ add_library(libopenmc SHARED src/photon_physics.F90 src/physics_common.F90 src/physics.F90 - src/physics_mg.F90 src/plot.F90 src/plot_header.F90 src/progress_header.F90 diff --git a/include/openmc/material.h b/include/openmc/material.h index 1af6f85976..909215ee37 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -26,6 +26,7 @@ class Material public: int32_t id; //!< Unique ID double volume_ {-1.0}; //!< Volume in [cm^3] + bool fissionable {false}; //!< Does this material contain fissionable nuclides //! \brief Default temperature for cells containing this material. //! diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index 7fed561f34..fdfebf71e4 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -14,19 +14,33 @@ namespace openmc { // SCATTER //============================================================================== -//! \brief Samples the scattering event +//! \brief samples particle behavior after a collision event. extern "C" void +collision_mg(Particle* p, Bank* fission_bank, const int64_t fission_bank_size, + const double* energy_bin_avg, const MaterialMacroXS& material_xs); + +//! \brief samples a reaction type. +//! +//! Note that there is special logic when suvival biasing is turned on since +//! fission and disappearance are treated implicitly. +void +sample_reaction(Particle* p, Bank* fission_bank, + const int64_t fission_bank_size, const double* energy_bin_avg, + const MaterialMacroXS& material_xs); + +//! \brief Samples the scattering event +void scatter(Particle* p, const double* energy_bin_avg); //! \brief Determines the average total, prompt and delayed neutrons produced //! from fission and creates the appropriate bank sites. -extern "C" void +void create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, - int64_t& bank_array_size, MaterialMacroXS& material_xs); + const int64_t& bank_array_size, const MaterialMacroXS& material_xs); //! \brief Handles an absorption event -extern "C" void -absorption(Particle* p, MaterialMacroXS& material_xs); +void +absorption(Particle* p, const MaterialMacroXS& material_xs); } // namespace openmc #endif // OPENMC_PHYSICS_MG_H diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 420ae5bbb8..12b97d9407 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2986,7 +2986,7 @@ contains ! Check if material is fissionable if (nuclides(materials(i) % nuclide(j)) % fissionable) then - materials(i) % fissionable = .true. + call materials(i) % set_fissionable(logical(.true., C_BOOL)) end if end do diff --git a/src/material.cpp b/src/material.cpp index 1a0c2d2078..2bd9228753 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -122,6 +122,13 @@ extern "C" { material_map[id] = index - 1; } + bool material_fissionable(Material* mat) {return mat->fissionable;} + + void material_set_fissionable(Material* mat, bool fissionable, int32_t index) + { + mat->fissionable = fissionable; + } + void extend_materials_c(int32_t n) { materials.reserve(materials.size() + n); diff --git a/src/material_header.F90 b/src/material_header.F90 index 71e6ec5f1a..a067269ce2 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -51,6 +51,20 @@ module material_header integer(C_INT32_T), intent(in), value :: index end subroutine material_set_id_c + function material_fissionable_c(mat_ptr) & + bind(C, name='material_fissionable') result(fissionable) + import C_PTR, C_BOOL + type(C_PTR), intent(in), value :: mat_ptr + logical(C_BOOL) :: fissionable + end function material_fissionable_c + + subroutine material_set_fissionable_c(mat_ptr, fissionable) & + bind(C, name='material_set_fissionable') + import C_PTR, C_BOOL + type(C_PTR), intent(in), value :: mat_ptr + logical(C_BOOL), intent(in), value :: fissionable + end subroutine material_set_fissionable_c + subroutine extend_materials_c(n) bind(C) import C_INT32_T integer(C_INT32_T), intent(in), value :: n @@ -95,7 +109,6 @@ module material_header character(20), allocatable :: sab_names(:) ! name of S(a,b) table ! Does this material contain fissionable nuclides? Is it depletable? - logical :: fissionable = .false. logical :: depletable = .false. ! enforce isotropic scattering in lab for specific nuclides @@ -105,6 +118,8 @@ module material_header contains procedure :: id => material_id 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 @@ -139,6 +154,18 @@ contains call material_set_id_c(this % ptr, id, index) end subroutine material_set_id + function material_fissionable(this) result(fissionable) + class(Material), intent(in) :: this + logical(C_BOOL) :: fissionable + fissionable = material_fissionable_c(this % ptr) + end function material_fissionable + + subroutine material_set_fissionable(this, fissionable) + class(Material), intent(in) :: this + logical(C_BOOL), intent(in) :: fissionable + call material_set_fissionable_c(this % ptr, fissionable) + end subroutine material_set_fissionable + function material_set_density(this, density) result(err) class(Material), intent(inout) :: this real(8), intent(in) :: density @@ -685,7 +712,7 @@ contains integer(C_INT) :: err if (index >= 1 .and. index <= size(materials)) then - fissionable = materials(index) % fissionable + fissionable = materials(index) % fissionable() err = 0 else err = E_OUT_OF_BOUNDS diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index 03a0e9402f..e6eb5a3da5 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -90,7 +90,8 @@ contains end if end do NUCLIDE_LOOP - mat % fissionable = query_fissionable_c(mat % n_nuclides, mat % nuclide) + call mat % set_fissionable(query_fissionable_c(mat % n_nuclides, & + mat % nuclide)) end do MATERIAL_LOOP diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 deleted file mode 100644 index d0a69c8631..0000000000 --- a/src/physics_mg.F90 +++ /dev/null @@ -1,126 +0,0 @@ -module physics_mg - ! This module contains the multi-group specific physics routines so as to not - ! hinder performance of the CE versions with multiple if-thens. - - use bank_header - use constants - use error, only: fatal_error, warning, write_message - use material_header, only: Material, materials - use math, only: rotate_angle - use mgxs_interface - use message_passing - use nuclide_header, only: MaterialMacroXS, material_xs - use particle_header - use physics_common - use random_lcg, only: prn - use settings - use simulation_header - use string, only: to_str - use tally_header - - implicit none - - interface - subroutine scatter(p, energy_bin_avg) bind(C) - import Particle, C_DOUBLE - type(Particle), intent(inout) :: p - real(C_DOUBLE), intent(in) :: energy_bin_avg(*) - end subroutine scatter - - subroutine create_fission_sites(p, bank_array, size_bank, bank_array_size, & - material_xs) bind(C) - import Particle, Bank, C_INT64_T, MaterialMacroXS - type(Particle), intent(inout) :: p - type(Bank), intent(inout) :: bank_array(*) - integer(C_INT64_T), intent(inout) :: size_bank - integer(C_INT64_T), intent(in) :: bank_array_size - type(MaterialMacroXS), intent(in) :: material_xs - end subroutine create_fission_sites - - subroutine absorption(p, material_xs) bind(C) - import Particle, MaterialMacroXS - type(Particle), intent(inout) :: p - type(MaterialMacroXS), intent(in) :: material_xs - end subroutine absorption - end interface - -contains - -!=============================================================================== -! COLLISION_MG samples a nuclide and reaction and then calls the appropriate -! routine for that reaction -!=============================================================================== - - subroutine collision_mg(p) - - type(Particle), intent(inout) :: p - - ! Add to collision counter for particle - p % n_collision = p % n_collision + 1 - - ! Sample nuclide/reaction for the material the particle is in - call sample_reaction(p) - - ! Display information about collision - if (verbosity >= 10 .or. trace) then - call write_message(" " // "Energy Group = " // trim(to_str(p % g))) - end if - - end subroutine collision_mg - -!=============================================================================== -! SAMPLE_REACTION samples a nuclide based on the macroscopic cross sections for -! each nuclide within a material and then samples a reaction for that nuclide -! and calls the appropriate routine to process the physics. Note that there is -! special logic when suvival biasing is turned on since fission and -! disappearance are treated implicitly. -!=============================================================================== - - subroutine sample_reaction(p) - - type(Particle), intent(inout) :: p - - type(Material), pointer :: mat - - mat => materials(p % material) - - ! Create fission bank sites. Note that while a fission reaction is sampled, - ! it never actually "happens", i.e. the weight of the particle does not - ! change when sampling fission sites. The following block handles all - ! absorption (including fission) - - if (mat % fissionable) then - if (run_mode == MODE_EIGENVALUE) then - call create_fission_sites(p, fission_bank, n_bank, & - size(fission_bank, KIND=C_INT64_T), & - material_xs) - elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then - call create_fission_sites(p, p % secondary_bank, p % n_secondary, & - size(p % secondary_bank, KIND=C_INT64_T), & - material_xs) - end if - end if - - ! If survival biasing is being used, the following subroutine adjusts the - ! weight of the particle. Otherwise, it checks to see if absorption occurs - - if (material_xs % absorption > ZERO) then - call absorption(p, material_xs) - else - p % absorb_wgt = ZERO - end if - if (.not. p % alive) return - - ! Sample a scattering reaction and determine the secondary energy of the - ! exiting neutron - call scatter(p, energy_bin_avg) - - ! Play russian roulette if survival biasing is turned on - if (survival_biasing) then - call russian_roulette(p) - if (.not. p % alive) return - end if - - end subroutine sample_reaction - -end module physics_mg diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index bd7bcec11b..8651e29f18 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -8,14 +8,76 @@ #include "openmc/constants.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" +#include "openmc/material.h" #include "openmc/math_functions.h" #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" +#include "openmc/physics_common.h" #include "openmc/random_lcg.h" #include "openmc/settings.h" +#include "openmc/simulation.h" namespace openmc { +void +collision_mg(Particle* p, Bank* fission_bank, const int64_t fission_bank_size, + const double* energy_bin_avg, const MaterialMacroXS& material_xs) +{ + // Add to the collision counter for the particle + p->n_collision++; + + // Sample the reaction type + sample_reaction(p, fission_bank, fission_bank_size, energy_bin_avg, + material_xs); + + // Display information about collision + if ((settings::verbosity >= 10) || (openmc_trace)) { + std::stringstream msg; + msg << " Energy Group = " << p->g; + write_message(msg, 1); + } +} + +void +sample_reaction(Particle* p, Bank* fission_bank, + const int64_t fission_bank_size, const double* energy_bin_avg, + const MaterialMacroXS& material_xs) +{ + // Create fission bank sites. Note that while a fission reaction is sampled, + // it never actually "happens", i.e. the weight of the particle does not + // change when sampling fission sites. The following block handles all + // absorption (including fission) + + if (materials[p->material - 1]->fissionable) { + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + create_fission_sites(p, fission_bank, n_bank, fission_bank_size, + material_xs); + } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && + (settings::create_fission_neutrons)) { + create_fission_sites(p, p->secondary_bank, p->n_secondary, + MAX_SECONDARY, material_xs); + } + } + + // If survival biasing is being used, the following subroutine adjusts the + // weight of the particle. Otherwise, it checks to see if absorption occurs. + if (material_xs.absorption > 0.) { + absorption(p, material_xs); + } else { + p->absorb_wgt = 0.; + } + if (!p->alive) return; + + // Sample a scattering event to determine the energy of the exiting neutron + scatter(p, energy_bin_avg); + + // Play Russian roulette if survival biasing is turned on + if (settings::survival_biasing) { + russian_roulette(p); + if (!p->alive) return; + } +} + void scatter(Particle* p, const double* energy_bin_avg) { @@ -42,7 +104,8 @@ scatter(Particle* p, const double* energy_bin_avg) void create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, - int64_t& bank_array_size, MaterialMacroXS& material_xs) + const int64_t& bank_array_size, + const MaterialMacroXS& material_xs) { // TODO: Heat generation from fission @@ -146,7 +209,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, } void -absorption(Particle* p, MaterialMacroXS& material_xs) +absorption(Particle* p, const MaterialMacroXS& material_xs) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing diff --git a/src/tracking.F90 b/src/tracking.F90 index 8cf2fc05c7..02f10f7cc8 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -17,7 +17,7 @@ module tracking use nuclide_header use particle_header use physics, only: collision - use physics_mg, only: collision_mg + ! use physics_mg, only: collision_mg use random_lcg, only: prn, prn_set_stream use settings use simulation_header @@ -33,6 +33,19 @@ module tracking implicit none + interface + subroutine collision_mg(p, fission_bank, fission_bank_size, & + energy_bin_avg, material_xs) bind(C) + import Particle, Bank, C_INT64_T, C_DOUBLE, MaterialMacroXS + type(Particle), intent(inout) :: p + type(Bank), intent(inout) :: fission_bank(*) + integer(C_INT64_T), value, intent(in) :: fission_bank_size + real(C_DOUBLE), intent(in) :: energy_bin_avg(*) + type(MaterialMacroXS), intent(in) :: material_xs + end subroutine collision_mg + + end interface + contains !=============================================================================== @@ -226,7 +239,8 @@ contains if (run_CE) then call collision(p) else - call collision_mg(p) + call collision_mg(p, fission_bank, & + size(fission_bank, kind=C_INT64_T), energy_bin_avg, material_xs) end if ! Score collision estimator tallies -- this is done after a collision