mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Start moving apply_derivative_to_score to C++
This commit is contained in:
Sterling Harper
parent
cdb10ef2cc
commit
e29931fef7
2 changed files with 237 additions and 191 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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. &
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue