mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 21:55:41 -04:00
Finish moving apply_derivative_to_score to C++
This commit is contained in:
Sterling Harper
parent
e29931fef7
commit
ddaed7311b
4 changed files with 275 additions and 409 deletions
|
|
@ -1,6 +1,8 @@
|
|||
#ifndef OPENMC_TALLIES_DERIVATIVE_H
|
||||
#define OPENMC_TALLIES_DERIVATIVE_H
|
||||
|
||||
#include "openmc/particle.h"
|
||||
|
||||
#include <unordered_map>
|
||||
#include <vector>
|
||||
|
||||
|
|
@ -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
|
||||
|
||||
//==============================================================================
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue