diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 9c9571bb2..fd8d080ad 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -2,6 +2,7 @@ #include "openmc/material.h" #include "openmc/nuclide.h" #include "openmc/settings.h" +#include "openmc/tallies/tally.h" #include "openmc/tallies/derivative.h" #include "openmc/xml_interface.h" @@ -105,6 +106,209 @@ read_tally_derivatives(pugi::xml_node* node) fatal_error("Differential tallies not supported in multi-group mode"); } +extern "C" void +apply_derivative_to_score_c(Particle* p, int i_tally, int i_nuclide, + double atom_density, int score_bin, double* score) +{ + //TODO: off-by-one + const Tally& tally {*model::tallies[i_tally-1]}; + + if (*score == 0) return; + + // If our score was previously c then the new score is + // c * (1/f * d_f/d_p + 1/c * d_c/d_p) + // where (1/f * d_f/d_p) is the (logarithmic) flux derivative and p is the + // perturbated variable. + + const auto& deriv {model::tally_derivs[tally.deriv_]}; + auto flux_deriv = deriv.flux_deriv; + + // Handle special cases where we know that d_c/d_p must be zero. + if (score_bin == SCORE_FLUX) { + *score *= flux_deriv; + return; + } else if (p->material == MATERIAL_VOID) { + *score *= flux_deriv; + return; + } + //TODO: off-by-one + const Material& material {*model::materials[p->material-1]}; + if (material.id_ != deriv.diff_material) { + *score *= flux_deriv; + return; + } + + switch (deriv.variable) { + + //============================================================================ + // Density derivative: + // c = Sigma_MT + // c = sigma_MT * N + // c = sigma_MT * rho * const + // d_c / d_rho = sigma_MT * const + // (1 / c) * (d_c / d_rho) = 1 / rho + + case DIFF_DENSITY: + switch (tally.estimator_) { + + case ESTIMATOR_ANALOG: + case ESTIMATOR_COLLISION: + switch (score_bin) { + + case SCORE_TOTAL: + case SCORE_SCATTER: + case SCORE_ABSORPTION: + case SCORE_FISSION: + case SCORE_NU_FISSION: + *score *= flux_deriv + 1. / material.density_gpcc_; + break; + + default: + fatal_error("Tally derivative not defined for a score on tally " + + std::to_string(tally.id_)); + } + break; + + default: + fatal_error("Differential tallies are only implemented for analog and " + "collision estimators."); + } + break; + + //============================================================================ + // Nuclide density derivative: + // If we are scoring a reaction rate for a single nuclide then + // c = Sigma_MT_i + // c = sigma_MT_i * N_i + // d_c / d_N_i = sigma_MT_i + // (1 / c) * (d_c / d_N_i) = 1 / N_i + // If the score is for the total material (i_nuclide = -1) + // c = Sum_i(Sigma_MT_i) + // d_c / d_N_i = sigma_MT_i + // (1 / c) * (d_c / d_N) = sigma_MT_i / Sigma_MT + // where i is the perturbed nuclide. + + case DIFF_NUCLIDE_DENSITY: + //TODO: off-by-one throughout on diff_nuclide + switch (tally.estimator_) { + + case ESTIMATOR_ANALOG: + if (p->event_nuclide != deriv.diff_nuclide) { + *score *= flux_deriv; + return; + } + + switch (score_bin) { + + case SCORE_TOTAL: + case SCORE_SCATTER: + case SCORE_ABSORPTION: + case SCORE_FISSION: + case SCORE_NU_FISSION: + { + // Find the index of the perturbed nuclide. + int i; + for (i = 0; i < material.nuclide_.size(); ++i) + if (material.nuclide_[i] == deriv.diff_nuclide - 1) break; + *score *= flux_deriv + 1. / material.atom_density_(i); + } + break; + + default: + fatal_error("Tally derivative not defined for a score on tally " + + std::to_string(tally.id_)); + } + break; + + case ESTIMATOR_COLLISION: + switch (score_bin) { + + case SCORE_TOTAL: + if (i_nuclide == -1 && simulation::material_xs.total) { + *score *= flux_deriv + + simulation::micro_xs[deriv.diff_nuclide-1].total + / simulation::material_xs.total; + } else if (i_nuclide == deriv.diff_nuclide-1 + && simulation::micro_xs[i_nuclide].total) { + *score *= flux_deriv + 1. / atom_density; + } else { + *score *= flux_deriv; + } + break; + + case SCORE_SCATTER: + if (i_nuclide == -1 && (simulation::material_xs.total + - simulation::material_xs.absorption)) { + *score *= flux_deriv + + (simulation::micro_xs[deriv.diff_nuclide-1].total + - simulation::micro_xs[deriv.diff_nuclide-1].absorption) + / (simulation::material_xs.total + - simulation::material_xs.absorption); + } else if (i_nuclide == deriv.diff_nuclide-1) { + *score *= flux_deriv + 1. / atom_density; + } else { + *score *= flux_deriv; + } + break; + + case SCORE_ABSORPTION: + if (i_nuclide == -1 && simulation::material_xs.absorption) { + *score *= flux_deriv + + simulation::micro_xs[deriv.diff_nuclide-1].absorption + / simulation::material_xs.absorption; + } else if (i_nuclide == deriv.diff_nuclide-1 + && simulation::micro_xs[i_nuclide].absorption) { + *score *= flux_deriv + 1. / atom_density; + } else { + *score *= flux_deriv; + } + break; + + case SCORE_FISSION: + if (i_nuclide == -1 && simulation::material_xs.fission) { + *score *= flux_deriv + + simulation::micro_xs[deriv.diff_nuclide-1].fission + / simulation::material_xs.fission; + } else if (i_nuclide == deriv.diff_nuclide-1 + && simulation::micro_xs[i_nuclide].fission) { + *score *= flux_deriv + 1. / atom_density; + } else { + *score *= flux_deriv; + } + break; + + case SCORE_NU_FISSION: + if (i_nuclide == -1 && simulation::material_xs.nu_fission) { + *score *= flux_deriv + + simulation::micro_xs[deriv.diff_nuclide-1].nu_fission + / simulation::material_xs.nu_fission; + } else if (i_nuclide == deriv.diff_nuclide-1 + && simulation::micro_xs[i_nuclide].nu_fission) { + *score *= flux_deriv + 1. / atom_density; + } else { + *score *= flux_deriv; + } + break; + + default: + fatal_error("Tally derivative not defined for a score on tally " + + std::to_string(tally.id_)); + } + break; + + default: + fatal_error("Differential tallies are only implemented for analog and " + "collision estimators."); + } + break; + + //============================================================================ + + case DIFF_TEMPERATURE: + break; + } +} + //! Adjust diff tally flux derivatives for a particle tracking event. extern "C" void diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 854aab751..c5efaabf8 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -115,6 +115,22 @@ contains real(8) :: flux_deriv real(8) :: dsig_s, dsig_a, dsig_f, cum_dsig + interface + subroutine apply_derivative_to_score_c(p, i_tally, i_nuclide, & + atom_density, score_bin, score) bind(C) + import Particle, C_INT, C_DOUBLE + type(Particle), intent(in) :: p + integer(C_INT), value, intent(in) :: i_tally + integer(C_INT), value, intent(in) :: i_nuclide + real(C_DOUBLE), value, intent(in) :: atom_density + integer(C_INT), value, intent(in) :: score_bin + real(C_DOUBLE), intent(inout) :: score + end subroutine + end interface + + call apply_derivative_to_score_c(p, i_tally, i_nuclide, atom_density, & + score_bin, score) + associate (t => tallies(i_tally) % obj) if (score == ZERO) return @@ -128,201 +144,26 @@ contains deriv => tally_deriv_c(t % deriv()) flux_deriv = deriv % flux_deriv + if (p % material == MATERIAL_VOID) then + return + else if (material_id(p % material) /= deriv % diff_material) then + return + end if + + if (deriv % variable == DIFF_NUCLIDE_DENSITY) then + if (t % estimator() == ESTIMATOR_ANALOG) then + if (p % event_nuclide /= deriv % diff_nuclide) return + end if + end if + !select case (tally_derivs(t % deriv()) % variable) select case (deriv % variable) - !========================================================================= - ! Density derivative: - ! c = Sigma_MT - ! c = sigma_MT * N - ! c = sigma_MT * rho * const - ! d_c / d_rho = sigma_MT * const - ! (1 / c) * (d_c / d_rho) = 1 / rho - case (DIFF_DENSITY) - select case (t % estimator()) - - case (ESTIMATOR_ANALOG) - - select case (score_bin) - - case (SCORE_FLUX) - score = score * flux_deriv - - case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & - SCORE_NU_FISSION) - if (material_id(p % material) == deriv % diff_material) then - score = score * (flux_deriv + ONE & - / material_density_gpcc(p % material)) - else - score = score * flux_deriv - end if - - case default - call fatal_error('Tally derivative not defined for a score on & - &tally ' // trim(to_str(t % id()))) - end select - - case (ESTIMATOR_COLLISION) - - select case (score_bin) - - case (SCORE_FLUX) - score = score * flux_deriv - - case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & - SCORE_NU_FISSION) - if (material_id(p % material) == deriv % diff_material) then - score = score * (flux_deriv + ONE & - / material_density_gpcc(p % material)) - else - score = score * flux_deriv - end if - - case default - call fatal_error('Tally derivative not defined for a score on & - &tally ' // trim(to_str(t % id()))) - end select - - case default - call fatal_error("Differential tallies are only implemented for & - &analog and collision estimators.") - end select - - !========================================================================= - ! Nuclide density derivative: - ! If we are scoring a reaction rate for a single nuclide then - ! c = Sigma_MT_i - ! c = sigma_MT_i * N_i - ! d_c / d_N_i = sigma_MT_i - ! (1 / c) * (d_c / d_N_i) = 1 / N_i - ! If the score is for the total material (i_nuclide = -1) - ! c = Sum_i(Sigma_MT_i) - ! d_c / d_N_i = sigma_MT_i - ! (1 / c) * (d_c / d_N) = sigma_MT_i / Sigma_MT - ! where i is the perturbed nuclide. + return case (DIFF_NUCLIDE_DENSITY) - select case (t % estimator()) - - case (ESTIMATOR_ANALOG) - - select case (score_bin) - - case (SCORE_FLUX) - score = score * flux_deriv - - case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & - SCORE_NU_FISSION) - if (material_id(p % material) == deriv % diff_material & - .and. p % event_nuclide == deriv % diff_nuclide) then - ! Search for the index of the perturbed nuclide. - 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 / material_atom_density(p % material, l)) - else - score = score * flux_deriv - end if - - case default - call fatal_error('Tally derivative not defined for a score on & - &tally ' // trim(to_str(t % id()))) - end select - - case (ESTIMATOR_COLLISION) - scoring_diff_nuclide = & - (material_id(p % material) == deriv % diff_material) & - .and. (i_nuclide+1 == deriv % diff_nuclide) - - select case (score_bin) - - case (SCORE_FLUX) - score = score * flux_deriv - - case (SCORE_TOTAL) - if (i_nuclide == -1 .and. & - material_id(p % material) == deriv % diff_material .and. & - material_xs % total /= ZERO) then - score = score * (flux_deriv & - + micro_xs(deriv % diff_nuclide) % total & - / material_xs % total) - else if (scoring_diff_nuclide .and. & - micro_xs(deriv % diff_nuclide) % total /= ZERO) then - score = score * (flux_deriv + ONE / atom_density) - else - score = score * flux_deriv - end if - - case (SCORE_SCATTER) - if (i_nuclide == -1 .and. & - material_id(p % material) == deriv % diff_material .and. & - material_xs % total - material_xs % absorption /= ZERO) then - score = score * (flux_deriv & - + (micro_xs(deriv % diff_nuclide) % total & - - micro_xs(deriv % diff_nuclide) % absorption) & - / (material_xs % total - material_xs % absorption)) - else if (scoring_diff_nuclide .and. & - (micro_xs(deriv % diff_nuclide) % total & - - micro_xs(deriv % diff_nuclide) % absorption) /= ZERO) then - score = score * (flux_deriv + ONE / atom_density) - else - score = score * flux_deriv - end if - - case (SCORE_ABSORPTION) - if (i_nuclide == -1 .and. & - material_id(p % material) == deriv % diff_material .and. & - material_xs % absorption /= ZERO) then - score = score * (flux_deriv & - + micro_xs(deriv % diff_nuclide) % absorption & - / material_xs % absorption ) - else if (scoring_diff_nuclide .and. & - micro_xs(deriv % diff_nuclide) % absorption /= ZERO) then - score = score * (flux_deriv + ONE / atom_density) - else - score = score * flux_deriv - end if - - case (SCORE_FISSION) - if (i_nuclide == -1 .and. & - material_id(p % material) == deriv % diff_material .and. & - material_xs % fission /= ZERO) then - score = score * (flux_deriv & - + micro_xs(deriv % diff_nuclide) % fission & - / material_xs % fission) - else if (scoring_diff_nuclide .and. & - micro_xs(deriv % diff_nuclide) % fission /= ZERO) then - score = score * (flux_deriv + ONE / atom_density) - else - score = score * flux_deriv - end if - - case (SCORE_NU_FISSION) - if (i_nuclide == -1 .and. & - material_id(p % material) == deriv % diff_material .and. & - material_xs % nu_fission /= ZERO) then - score = score * (flux_deriv & - + micro_xs(deriv % diff_nuclide) % nu_fission & - / material_xs % nu_fission) - else if (scoring_diff_nuclide .and. & - micro_xs(deriv % diff_nuclide) % nu_fission /= ZERO) then - score = score * (flux_deriv + ONE / atom_density) - else - score = score * flux_deriv - end if - - case default - call fatal_error('Tally derivative not defined for a score on & - &tally ' // trim(to_str(t % id()))) - end select - - case default - call fatal_error("Differential tallies are only implemented for & - &analog and collision estimators.") - end select + return !========================================================================= ! Temperature derivative: @@ -345,7 +186,7 @@ contains select case (score_bin) case (SCORE_FLUX) - score = score * flux_deriv + !score = score * flux_deriv case (SCORE_TOTAL) if (material_id(p % material) == deriv % diff_material .and. & @@ -467,7 +308,8 @@ contains select case (score_bin) case (SCORE_FLUX) - score = score * flux_deriv + !score = score * flux_deriv + return case (SCORE_TOTAL) if (i_nuclide == -1 .and. &