diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h index d9202433c3..da9fa60a6e 100644 --- a/include/openmc/tallies/derivative.h +++ b/include/openmc/tallies/derivative.h @@ -1,6 +1,8 @@ #ifndef OPENMC_TALLIES_DERIVATIVE_H #define OPENMC_TALLIES_DERIVATIVE_H +#include "openmc/particle.h" + #include #include @@ -23,6 +25,16 @@ extern "C" struct TallyDerivative { TallyDerivative(pugi::xml_node node); }; +//============================================================================== +// Non-method functions +//============================================================================== + +//! Scale the given score by its logarithmic derivative + +void +apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide, + double atom_density, int score_bin, double* score); + } // namespace openmc //============================================================================== diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index fd8d080ad4..b0ff83ea3a 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -106,8 +106,8 @@ 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, +void +apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide, double atom_density, int score_bin, double* score) { //TODO: off-by-one @@ -303,8 +303,268 @@ apply_derivative_to_score_c(Particle* p, int i_tally, int i_nuclide, break; //============================================================================ + // Temperature 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_T = (d_sigma_Mt_i / d_T) * N_i + // (1 / c) * (d_c / d_T) = (d_sigma_MT_i / d_T) / sigma_MT_i + // If the score is for the total material (i_nuclide = -1) + // (1 / c) * (d_c / d_T) = Sum_i((d_sigma_MT_i / d_T) * N_i) / Sigma_MT_i + // where i is the perturbed nuclide. The d_sigma_MT_i / d_T term is + // computed by multipole_deriv_eval. It only works for the resolved + // resonance range and requires multipole data. case DIFF_TEMPERATURE: + switch (tally.estimator_) { + + case ESTIMATOR_ANALOG: + { + // Find the index of the event nuclide. + int i; + for (i = 0; i < material.nuclide_.size(); ++i) + if (material.nuclide_[i] == p->event_nuclide-1) break; + + const auto& nuc {*data::nuclides[p->event_nuclide-1]}; + if (!multipole_in_range(&nuc, p->last_E)) { + *score *= flux_deriv; + break; + } + + switch (score_bin) { + + case SCORE_TOTAL: + if (simulation::micro_xs[p->event_nuclide-1].total) { + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + *score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i) + / simulation::material_xs.total; + } else { + *score *= flux_deriv; + } + break; + + case SCORE_SCATTER: + if (simulation::micro_xs[p->event_nuclide-1].total + - simulation::micro_xs[p->event_nuclide-1].absorption) { + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + *score *= flux_deriv + dsig_s * material.atom_density_(i) + / (simulation::material_xs.total + - simulation::material_xs.absorption); + } else { + *score *= flux_deriv; + } + break; + + case SCORE_ABSORPTION: + if (simulation::micro_xs[p->event_nuclide-1].absorption) { + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + *score *= flux_deriv + dsig_a * material.atom_density_(i) + / simulation::material_xs.absorption; + } else { + *score *= flux_deriv; + } + break; + + case SCORE_FISSION: + if (simulation::micro_xs[p->event_nuclide-1].fission) { + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + *score *= flux_deriv + dsig_f * material.atom_density_(i) + / simulation::material_xs.fission; + } else { + *score *= flux_deriv; + } + break; + + case SCORE_NU_FISSION: + if (simulation::micro_xs[p->event_nuclide-1].fission) { + double nu = simulation::micro_xs[p->event_nuclide-1].nu_fission + / simulation::micro_xs[p->event_nuclide-1].fission; + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + *score *= flux_deriv + nu * dsig_f * material.atom_density_(i) + / simulation::material_xs.nu_fission; + } else { + *score *= flux_deriv; + } + break; + + default: + fatal_error("Tally derivative not defined for a score on tally " + + std::to_string(tally.id_)); + } + } + break; + + case ESTIMATOR_COLLISION: + if (i_nuclide != -1) { + const auto& nuc {*data::nuclides[i_nuclide]}; + if (!multipole_in_range(&nuc, p->last_E)) { + *score *= flux_deriv; + return; + } + } + + switch (score_bin) { + + case SCORE_TOTAL: + if (i_nuclide == -1 && simulation::material_xs.total) { + double cum_dsig = 0; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto i_nuc = material.nuclide_[i]; + const auto& nuc {*data::nuclides[i_nuc]}; + if (multipole_in_range(&nuc, p->last_E) + && simulation::micro_xs[i_nuc].total) { + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + cum_dsig += (dsig_s + dsig_a) * material.atom_density_(i); + } + } + *score *= flux_deriv + cum_dsig / simulation::material_xs.total; + } else if (simulation::micro_xs[i_nuclide].total) { + const auto& nuc {*data::nuclides[i_nuclide]}; + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + *score *= flux_deriv + + (dsig_s + dsig_a) / simulation::micro_xs[i_nuclide].total; + } else { + *score *= flux_deriv; + } + break; + + case SCORE_SCATTER: + if (i_nuclide == -1 && (simulation::material_xs.total + - simulation::material_xs.absorption)) { + double cum_dsig = 0; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto i_nuc = material.nuclide_[i]; + const auto& nuc {*data::nuclides[i_nuc]}; + if (multipole_in_range(&nuc, p->last_E) + && (simulation::micro_xs[i_nuc].total + - simulation::micro_xs[i_nuc].absorption)) { + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + cum_dsig += dsig_s * material.atom_density_(i); + } + } + *score *= flux_deriv + cum_dsig / (simulation::material_xs.total + - simulation::material_xs.absorption); + } else if (simulation::micro_xs[i_nuclide].total + - simulation::micro_xs[i_nuclide].absorption) { + const auto& nuc {*data::nuclides[i_nuclide]}; + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + *score *= flux_deriv + dsig_s / (simulation::micro_xs[i_nuclide].total + - simulation::micro_xs[i_nuclide].absorption); + } else { + *score *= flux_deriv; + } + break; + + case SCORE_ABSORPTION: + if (i_nuclide == -1 && simulation::material_xs.absorption) { + double cum_dsig = 0; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto i_nuc = material.nuclide_[i]; + const auto& nuc {*data::nuclides[i_nuc]}; + if (multipole_in_range(&nuc, p->last_E) + && simulation::micro_xs[i_nuc].absorption) { + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + cum_dsig += dsig_a * material.atom_density_(i); + } + } + *score *= flux_deriv + cum_dsig / simulation::material_xs.absorption; + } else if (simulation::micro_xs[i_nuclide].absorption) { + const auto& nuc {*data::nuclides[i_nuclide]}; + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + *score *= flux_deriv + + dsig_a / simulation::micro_xs[i_nuclide].absorption; + } else { + *score *= flux_deriv; + } + break; + + case SCORE_FISSION: + if (i_nuclide == -1 && simulation::material_xs.fission) { + double cum_dsig = 0; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto i_nuc = material.nuclide_[i]; + const auto& nuc {*data::nuclides[i_nuc]}; + if (multipole_in_range(&nuc, p->last_E) + && simulation::micro_xs[i_nuc].fission) { + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + cum_dsig += dsig_f * material.atom_density_(i); + } + } + *score *= flux_deriv + cum_dsig / simulation::material_xs.fission; + } else if (simulation::micro_xs[i_nuclide].fission) { + const auto& nuc {*data::nuclides[i_nuclide]}; + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + *score *= flux_deriv + + dsig_f / simulation::micro_xs[i_nuclide].fission; + } else { + *score *= flux_deriv; + } + break; + + case SCORE_NU_FISSION: + if (i_nuclide == -1 && simulation::material_xs.nu_fission) { + double cum_dsig = 0; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto i_nuc = material.nuclide_[i]; + const auto& nuc {*data::nuclides[i_nuc]}; + if (multipole_in_range(&nuc, p->last_E) + && simulation::micro_xs[i_nuc].fission) { + double nu = simulation::micro_xs[i_nuc].nu_fission + / simulation::micro_xs[i_nuc].fission; + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + cum_dsig += nu * dsig_f * material.atom_density_(i); + } + } + *score *= flux_deriv + cum_dsig / simulation::material_xs.nu_fission; + } else if (simulation::micro_xs[i_nuclide].fission) { + const auto& nuc {*data::nuclides[i_nuclide]}; + double dsig_s, dsig_a, dsig_f; + std::tie(dsig_s, dsig_a, dsig_f) + = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + *score *= flux_deriv + + dsig_f / simulation::micro_xs[i_nuclide].fission; + } else { + *score *= flux_deriv; + } + break; + + default: + break; + } + break; + + default: + fatal_error("Differential tallies are only implemented for analog and " + "collision estimators."); + } break; } } @@ -344,7 +604,6 @@ score_track_derivative(Particle* p, double distance) 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 @@ -427,7 +686,7 @@ score_collision_derivative(Particle* p) 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). + // 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 diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index c5efaabf8e..d4fb204554 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -95,407 +95,6 @@ contains end if end subroutine init_tally_routines -!=============================================================================== -! APPLY_DERIVATIVE_TO_SCORE multiply the given score by its relative derivative -!=============================================================================== - - subroutine apply_derivative_to_score(p, i_tally, i_nuclide, atom_density, & - score_bin, score) bind(C) - 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 ! atom/b-cm - integer(C_INT), value, intent(in) :: score_bin - real(C_DOUBLE), intent(inout) :: score - - type(TallyDerivative), pointer :: deriv - integer :: l - integer :: i_nuc - logical :: scoring_diff_nuclide - 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 - - ! 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. - - !associate(deriv => tally_derivs(t % deriv())) - 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) - - case (DIFF_DENSITY) - return - - case (DIFF_NUCLIDE_DENSITY) - return - - !========================================================================= - ! Temperature 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_T = (d_sigma_Mt_i / d_T) * N_i - ! (1 / c) * (d_c / d_T) = (d_sigma_MT_i / d_T) / sigma_MT_i - ! If the score is for the total material (i_nuclide = -1) - ! (1 / c) * (d_c / d_T) = Sum_i((d_sigma_MT_i / d_T) * N_i) / Sigma_MT_i - ! where i is the perturbed nuclide. The d_sigma_MT_i / d_T term is - ! computed by multipole_deriv_eval. It only works for the resolved - ! resonance range and requires multipole data. - - case (DIFF_TEMPERATURE) - select case (t % estimator()) - - case (ESTIMATOR_ANALOG) - - select case (score_bin) - - case (SCORE_FLUX) - !score = score * flux_deriv - - case (SCORE_TOTAL) - if (material_id(p % material) == deriv % diff_material .and. & - micro_xs(p % event_nuclide) % total > ZERO) then - ! Search for the index of the perturbed nuclide. - do l = 1, material_nuclide_size(p % material) - if (material_nuclide(p % material, l) == p % event_nuclide) exit - end do - - dsig_s = ZERO - dsig_a = ZERO - associate (nuc => nuclides(p % event_nuclide)) - if (multipole_in_range(nuc % ptr, p % last_E)) then - call multipole_deriv_eval(nuc % ptr, p % last_E, & - p % sqrtkT, dsig_s, dsig_a, dsig_f) - end if - end associate - score = score * (flux_deriv & - + (dsig_s + dsig_a) * material_atom_density(p % material, l) & - / material_xs % total) - else - score = score * flux_deriv - end if - - case (SCORE_SCATTER) - if (material_id(p % material) == deriv % diff_material .and. & - (micro_xs(p % event_nuclide) % total & - - micro_xs(p % event_nuclide) % absorption) > ZERO) then - ! Search for the index of the perturbed nuclide. - do l = 1, material_nuclide_size(p % material) - if (material_nuclide(p % material, l) == p % event_nuclide) exit - end do - - dsig_s = ZERO - associate (nuc => nuclides(p % event_nuclide)) - if (multipole_in_range(nuc % ptr, p % last_E)) then - call multipole_deriv_eval(nuc % ptr, p % last_E, & - p % sqrtkT, dsig_s, dsig_a, dsig_f) - end if - end associate - score = score * (flux_deriv + dsig_s * material_atom_density(p % material, l) / & - (material_xs % total - material_xs % absorption)) - else - score = score * flux_deriv - end if - - case (SCORE_ABSORPTION) - if (material_id(p % material) == deriv % diff_material .and. & - micro_xs(p % event_nuclide) % absorption > ZERO) then - ! Search for the index of the perturbed nuclide. - do l = 1, material_nuclide_size(p % material) - if (material_nuclide(p % material, l) == p % event_nuclide) exit - end do - - dsig_a = ZERO - associate (nuc => nuclides(p % event_nuclide)) - if (multipole_in_range(nuc % ptr, p % last_E)) then - call multipole_deriv_eval(nuc % ptr, p % last_E, & - p % sqrtkT, dsig_s, dsig_a, dsig_f) - end if - end associate - score = score * (flux_deriv + dsig_a * material_atom_density(p % material, l) & - / material_xs % absorption) - else - score = score * flux_deriv - end if - - case (SCORE_FISSION) - if (material_id(p % material) == deriv % diff_material .and. & - micro_xs(p % event_nuclide) % fission > ZERO) then - ! Search for the index of the perturbed nuclide. - do l = 1, material_nuclide_size(p % material) - if (material_nuclide(p % material, l) == p % event_nuclide) exit - end do - - dsig_f = ZERO - associate (nuc => nuclides(p % event_nuclide)) - if (multipole_in_range(nuc % ptr, p % last_E)) then - call multipole_deriv_eval(nuc % ptr, p % last_E, & - p % sqrtkT, dsig_s, dsig_a, dsig_f) - end if - end associate - score = score * (flux_deriv & - + dsig_f * material_atom_density(p % material, l) / material_xs % fission) - else - score = score * flux_deriv - end if - - case (SCORE_NU_FISSION) - if (material_id(p % material) == deriv % diff_material .and. & - micro_xs(p % event_nuclide) % nu_fission > ZERO) then - ! Search for the index of the perturbed nuclide. - do l = 1, material_nuclide_size(p % material) - if (material_nuclide(p % material, l) == p % event_nuclide) exit - end do - - dsig_f = ZERO - associate (nuc => nuclides(p % event_nuclide)) - if (multipole_in_range(nuc % ptr, p % last_E)) then - call multipole_deriv_eval(nuc % ptr, p % last_E, & - p % sqrtkT, dsig_s, dsig_a, dsig_f) - end if - end associate - score = score * (flux_deriv & - + 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) - 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 - return - - case (SCORE_TOTAL) - if (i_nuclide == -1 .and. & - material_id(p % material) == deriv % diff_material .and. & - material_xs % total > ZERO) then - cum_dsig = ZERO - 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(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) & - * material_atom_density(p % material, l) - end if - end associate - end do - score = score * (flux_deriv & - + cum_dsig / material_xs % total) - else if (material_id(p % material) == deriv % diff_material & - .and. material_xs % total > ZERO) then - dsig_s = ZERO - dsig_a = ZERO - associate (nuc => nuclides(i_nuclide+1)) - if (multipole_in_range(nuc % ptr, p % last_E)) then - call multipole_deriv_eval(nuc % ptr, p % last_E, & - p % sqrtkT, dsig_s, dsig_a, dsig_f) - end if - end associate - score = score * (flux_deriv & - + (dsig_s + dsig_a) / micro_xs(i_nuclide+1) % total) - 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 - cum_dsig = ZERO - 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(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 * material_atom_density(p % material, l) - end if - end associate - end do - score = score * (flux_deriv + cum_dsig & - / (material_xs % total - material_xs % absorption)) - else if ( material_id(p % material) == deriv % diff_material & - .and. (material_xs % total - material_xs % absorption) > ZERO)& - then - dsig_s = ZERO - associate (nuc => nuclides(i_nuclide+1)) - if (multipole_in_range(nuc % ptr, p % last_E)) then - call multipole_deriv_eval(nuc % ptr, p % last_E, & - p % sqrtkT, dsig_s, dsig_a, dsig_f) - end if - end associate - score = score * (flux_deriv + dsig_s & - / (micro_xs(i_nuclide+1) % total & - - micro_xs(i_nuclide+1) % absorption)) - 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 - cum_dsig = ZERO - 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(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 * material_atom_density(p % material, l) - end if - end associate - end do - score = score * (flux_deriv & - + cum_dsig / material_xs % absorption) - else if (material_id(p % material) == deriv % diff_material & - .and. material_xs % absorption > ZERO) then - dsig_a = ZERO - associate (nuc => nuclides(i_nuclide+1)) - if (multipole_in_range(nuc % ptr, p % last_E)) then - call multipole_deriv_eval(nuc % ptr, p % last_E, & - p % sqrtkT, dsig_s, dsig_a, dsig_f) - end if - end associate - score = score * (flux_deriv & - + dsig_a / micro_xs(i_nuclide+1) % absorption) - 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 - cum_dsig = ZERO - 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(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 * material_atom_density(p % material, l) - end if - end associate - end do - score = score * (flux_deriv & - + cum_dsig / material_xs % fission) - else if (material_id(p % material) == deriv % diff_material & - .and. material_xs % fission > ZERO) then - dsig_f = ZERO - associate (nuc => nuclides(i_nuclide+1)) - if (multipole_in_range(nuc % ptr, p % last_E)) then - call multipole_deriv_eval(nuc % ptr, p % last_E, & - p % sqrtkT, dsig_s, dsig_a, dsig_f) - end if - end associate - score = score * (flux_deriv & - + dsig_f / micro_xs(i_nuclide+1) % fission) - 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 - cum_dsig = ZERO - 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(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 * material_atom_density(p % material, l) & - * micro_xs(i_nuc) % nu_fission & - / micro_xs(i_nuc) % fission - end if - end associate - end do - score = score * (flux_deriv & - + cum_dsig / material_xs % nu_fission) - else if (material_id(p % material) == deriv % diff_material & - .and. material_xs % nu_fission > ZERO) then - dsig_f = ZERO - associate (nuc => nuclides(i_nuclide+1)) - if (multipole_in_range(nuc % ptr, p % last_E)) then - call multipole_deriv_eval(nuc % ptr, p % last_E, & - p % sqrtkT, dsig_s, dsig_a, dsig_f) - end if - end associate - score = score * (flux_deriv & - + dsig_f / micro_xs(i_nuclide+1) % fission) - 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 - end select - end associate - end subroutine apply_derivative_to_score - !=============================================================================== ! ACCUMULATE_TALLIES accumulates the sum of the contributions from each history ! within the batch to a new random variable diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 3b94b70e3e..704bedaa50 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -38,10 +38,6 @@ namespace openmc { // Functions defined in Fortran //============================================================================== -extern "C" void -apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide, - double atom_density, int score_bin, double* score); - extern "C" int energy_filter_search(const EnergyFilter* filt, double val);