From cdb10ef2ccc3e5c28dd4956048c13c1edb1d8ec2 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 10 Feb 2019 18:43:48 -0500 Subject: [PATCH] Move score_*_derivative to C++ --- include/openmc/nuclide.h | 2 + src/nuclide.cpp | 2 +- src/tallies/derivative.cpp | 133 ++++++++++++++++++++++++++++++ src/tallies/tally.F90 | 161 +++---------------------------------- 4 files changed, 147 insertions(+), 151 deletions(-) diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index dd4dfc678..be81dd917 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -176,6 +176,8 @@ private: //! Checks for the right version of nuclear data within HDF5 files void check_data_version(hid_t file_id); +extern "C" bool multipole_in_range(const Nuclide* nuc, double E); + //============================================================================== // Global variables //============================================================================== diff --git a/src/nuclide.cpp b/src/nuclide.cpp index b507b9130..057b7c113 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -990,7 +990,7 @@ extern "C" void multipole_deriv_eval(Nuclide* nuc, double E, double sqrtkT, std::tie(*sig_s, *sig_a, *sig_f) = nuc->multipole_->evaluate_deriv(E, sqrtkT); } -extern "C" bool multipole_in_range(Nuclide* nuc, double E) +extern "C" bool multipole_in_range(const Nuclide* nuc, double E) { return nuc->multipole_ && E >= nuc->multipole_->E_min_&& E <= nuc->multipole_->E_max_; diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 5908ce338..9c9571bb2 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -1,9 +1,12 @@ #include "openmc/error.h" +#include "openmc/material.h" #include "openmc/nuclide.h" #include "openmc/settings.h" #include "openmc/tallies/derivative.h" #include "openmc/xml_interface.h" +#include + template class std::vector; namespace openmc { @@ -102,6 +105,136 @@ read_tally_derivatives(pugi::xml_node* node) fatal_error("Differential tallies not supported in multi-group mode"); } +//! Adjust diff tally flux derivatives for a particle tracking event. + +extern "C" void +score_track_derivative(Particle* p, double distance) +{ + // A void material cannot be perturbed so it will not affect flux derivatives. + if (p->material == MATERIAL_VOID) return; + //TODO: off-by-one + const Material& material {*model::materials[p->material-1]}; + + for (auto& deriv : model::tally_derivs) { + if (deriv.diff_material != material.id_) continue; + + switch (deriv.variable) { + + case DIFF_DENSITY: + // phi is proportional to e^(-Sigma_tot * dist) + // (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 -= distance * simulation::material_xs.total + / material.density_gpcc_; + break; + + case DIFF_NUCLIDE_DENSITY: + // phi is proportional to e^(-Sigma_tot * dist) + // (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist + // (1 / phi) * (d_phi / d_N) = - sigma_tot * dist + //TODO: off-by-one + deriv.flux_deriv -= distance + * simulation::micro_xs[deriv.diff_nuclide-1].total; + break; + + case DIFF_TEMPERATURE: + for (auto i = 0; i < material.nuclide_.size(); ++i) { + const auto& nuc {*data::nuclides[material.nuclide_[i]]}; + //TODO: off-by-one + if (multipole_in_range(&nuc, p->last_E)) { + // phi is proportional to e^(-Sigma_tot * dist) + // (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist + // (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->E, p->sqrtkT); + deriv.flux_deriv -= distance * (dsig_s + dsig_a) + * material.atom_density_(i); + } + } + break; + } + } +} + +//! Adjust diff tally flux derivatives for a particle scattering event. +// +//! Note that this subroutine will be called after absorption events in +//! addition to scattering events, but any flux derivatives scored after an +//! absorption will never be tallied. The paricle will be killed before any +//! further tallies are scored. + +extern "C" void +score_collision_derivative(Particle* p) +{ + // A void material cannot be perturbed so it will not affect flux derivatives. + if (p->material == MATERIAL_VOID) return; + //TODO: off-by-one + const Material& material {*model::materials[p->material-1]}; + + for (auto& deriv : model::tally_derivs) { + if (deriv.diff_material != material.id_) continue; + + switch (deriv.variable) { + + case DIFF_DENSITY: + // phi is proportional to Sigma_s + // (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s + // (1 / phi) * (d_phi / d_rho) = 1 / rho + deriv.flux_deriv += 1. / material.density_gpcc_; + break; + + case DIFF_NUCLIDE_DENSITY: + //TODO: off-by-one throughout on diff_nuclide + if (p->event_nuclide != deriv.diff_nuclide) continue; + // Find the index in this material for the diff_nuclide. + int i; + for (i = 0; i < material.nuclide_.size(); ++i) + if (material.nuclide_[i] == deriv.diff_nuclide - 1) break; + // Make sure we found the nuclide. + if (material.nuclide_[i] != deriv.diff_nuclide - 1) { + std::stringstream err_msg; + err_msg << "Could not find nuclide " + << data::nuclides[deriv.diff_nuclide-1]->name_ << " in material " + << material.id_ << " for tally derivative " << deriv.id; + fatal_error(err_msg); + } + // phi is proportional to Sigma_s + // (1 / phi) * (d_phi / d_N) = (d_Sigma_s / d_N) / Sigma_s + // (1 / phi) * (d_phi / d_N) = sigma_s / Sigma_s + // (1 / phi) * (d_phi / d_N) = 1 / N + deriv.flux_deriv += 1. / material.atom_density_(i); + break; + + case DIFF_TEMPERATURE: + // Loop over the material's nuclides until we find the event nuclide. + for (auto i_nuc : material.nuclide_) { + const auto& nuc {*data::nuclides[i_nuc]}; + //TODO: off-by-one + if (i_nuc == p->event_nuclide - 1 + && multipole_in_range(&nuc, p->last_E)) { + // phi is proportional to Sigma_s + // (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s + // (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s + const auto& micro_xs {simulation::micro_xs[i_nuc]}; + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + deriv.flux_deriv += dsig_s / (micro_xs.total - micro_xs.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 d_P(E'->E, uvw'->uvw) / d_T = 0 and only + // computing d_S(E') / d_T. Using this approximation in the vicinity + // of low-energy resonances causes errors (~2-5% for PWR pincell + // eigenvalue derivatives). + } + } + break; + } + } +} + //! Set the flux derivatives on differential tallies to zero. extern "C" void zero_flux_derivs() diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index a83a9c21b..854aab751 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -65,6 +65,17 @@ module tally type(Particle) :: p end subroutine + subroutine score_track_derivative(p, distance) bind(C) + import Particle, C_DOUBLE + type(Particle) :: p + real(C_DOUBLE), value :: distance + end subroutine + + subroutine score_collision_derivative(p) bind(C) + import Particle + type(Particle) :: p + end subroutine + subroutine zero_flux_derivs() bind(C) end subroutine end interface @@ -643,156 +654,6 @@ contains end associate end subroutine apply_derivative_to_score -!=============================================================================== -! SCORE_TRACK_DERIVATIVE Adjust flux derivatives on differential tallies to -! account for a neutron travelling through a perturbed material. -!=============================================================================== - - subroutine score_track_derivative(p, distance) - type(Particle), intent(in) :: p - real(8), intent(in) :: distance ! Neutron flight distance - - type(TallyDerivative), pointer :: deriv - integer :: i, l - real(8) :: dsig_s, dsig_a, dsig_f - - ! A void material cannot be perturbed so it will not affect flux derivatives - if (p % material == MATERIAL_VOID) return - - do i = 0, n_tally_derivs() - 1 - !associate(deriv => tally_derivs(i)) - deriv => tally_deriv_c(i) - select case (deriv % variable) - - case (DIFF_DENSITY) - if (material_id(p % material) == deriv % diff_material) then - ! phi is proportional to e^(-Sigma_tot * dist) - ! (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 / material_density_gpcc(p % material) - end if - - case (DIFF_NUCLIDE_DENSITY) - if (material_id(p % material) == deriv % diff_material) then - ! phi is proportional to e^(-Sigma_tot * dist) - ! (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist - ! (1 / phi) * (d_phi / d_N) = - sigma_tot * dist - deriv % flux_deriv = deriv % flux_deriv & - - distance * micro_xs(deriv % diff_nuclide) % total - end if - - case (DIFF_TEMPERATURE) - if (material_id(p % material) == deriv % diff_material) then - 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 - ! (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist - 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) * material_atom_density(p % material, l) - end if - end associate - end do - end if - end select - !end associate - end do - end subroutine score_track_derivative - -!=============================================================================== -! SCORE_COLLISION_DERIVATIVE Adjust flux derivatives on differential tallies to -! account for a neutron scattering in the perturbed material. Note that this -! subroutine will be called after absorption events in addition to scattering -! events, but any flux derivatives scored after an absorption will never be -! tallied. This is because the order of operations is -! 1. Particle is moved. -! 2. score_track_derivative is called. -! 3. Collision physics are computed, and the particle is labeled absorbed. -! 4. Analog- and collision-estimated tallies are scored. -! 5. This subroutine is called. -! 6. Particle is killed and no more tallies are scored. -! Hence, it is safe to assume that only derivative of the scattering cross -! section need to be computed here. -!=============================================================================== - - subroutine score_collision_derivative(p) - type(Particle), intent(in) :: p - - type(TallyDerivative), pointer :: deriv - 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 - if (p % material == MATERIAL_VOID) return - - do i = 0, n_tally_derivs() - 1 - !associate(deriv => tally_derivs(i)) - deriv => tally_deriv_c(i) - select case (deriv % variable) - - case (DIFF_DENSITY) - if (material_id(p % material) == deriv % diff_material) then - ! phi is proportional to Sigma_s - ! (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 / material_density_gpcc(p % material) - end if - - case (DIFF_NUCLIDE_DENSITY) - if (material_id(p % material) == deriv % diff_material & - .and. p % event_nuclide == deriv % diff_nuclide) then - ! Find the index in this material for the diff_nuclide. - 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 (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 - ! (1 / phi) * (d_phi / d_N) = (d_Sigma_s / d_N) / Sigma_s - ! (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 / material_atom_density(p % material, j) - end if - - case (DIFF_TEMPERATURE) - if (material_id(p % material) == deriv % diff_material) then - 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 - ! (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s - 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(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 - ! d_P(E'->E, uvw'->uvw) / d_T = 0 and only computing - ! d_S(E') / d_T. Using this approximation in the vicinity - ! of low-energy resonances causes errors (~2-5% for PWR - ! pincell eigenvalue derivatives). - end if - end associate - end do - end if - end select - !end associate - end do - end subroutine score_collision_derivative - !=============================================================================== ! ACCUMULATE_TALLIES accumulates the sum of the contributions from each history ! within the batch to a new random variable