From 38fa90b3e42d525db6f1c5da9ec70f27acd34fe5 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Thu, 12 Dec 2019 16:02:06 +0000 Subject: [PATCH] Moved flux derivatives to the particle. Diff tallies are still not getting right answer, but they look closer... --- include/openmc/particle.h | 2 ++ include/openmc/tallies/derivative.h | 12 ++------ src/particle.cpp | 4 ++- src/simulation.cpp | 44 ++++++++++++--------------- src/tallies/derivative.cpp | 46 +++++++++++++++-------------- src/tallies/tally.cpp | 5 +--- 6 files changed, 52 insertions(+), 61 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index f818299cd..6d1e05a49 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -313,6 +313,8 @@ public: std::vector secondary_bank_; int64_t current_work_; // current work index + + std::vector flux_derivs_; // Derivatives of the current particle's weight }; } // namespace openmc diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h index d937116b1..a0db07da2 100644 --- a/include/openmc/tallies/derivative.h +++ b/include/openmc/tallies/derivative.h @@ -19,7 +19,6 @@ struct TallyDerivative { int variable; //!< Independent variable (like temperature) int diff_material; //!< Material this derivative is applied to int diff_nuclide; //!< Nuclide this material is applied to - double flux_deriv; //!< Derivative of the current particle's weight TallyDerivative() {} explicit TallyDerivative(pugi::xml_node node); @@ -46,16 +45,16 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, //! further tallies are scored. // //! \param p The particle being tracked -void score_collision_derivative(const Particle* p); +void score_collision_derivative(Particle* p); //! Adjust diff tally flux derivatives for a particle tracking event. // //! \param p The particle being tracked //! \param distance The distance in [cm] traveled by the particle -void score_track_derivative(const Particle* p, double distance); +void score_track_derivative(Particle* p, double distance); //! Set the flux derivatives on differential tallies to zero. -void zero_flux_derivs(); +void zero_flux_derivs(std::vector v); } // namespace openmc @@ -63,15 +62,10 @@ void zero_flux_derivs(); // Global variables //============================================================================== -// Explicit vector template specialization declaration of threadprivate variable -// outside of the openmc namespace for the picky Intel compiler. -//extern template class std::vector; - namespace openmc { namespace model { extern std::vector tally_derivs; -//#pragma omp threadprivate(tally_derivs) extern std::unordered_map tally_deriv_map; } // namespace model diff --git a/src/particle.cpp b/src/particle.cpp index d6867da5c..2fd9425df 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -109,6 +109,7 @@ Particle::from_source(const Bank* src) material_ = C_NONE; n_collision_ = 0; fission_ = false; + std::fill(flux_derivs_.begin(), flux_derivs_.end(), 0.0); // Copy attributes from source bank site type_ = src->particle; @@ -154,7 +155,8 @@ Particle::transport() if (write_track_) add_particle_track(); // Every particle starts with no accumulated flux derivative. - if (!model::active_tallies.empty()) zero_flux_derivs(); + if (!model::active_tallies.empty()) + flux_derivs_.resize(model::tally_derivs.size(), 0.0); while (true) { // Set the random number stream diff --git a/src/simulation.cpp b/src/simulation.cpp index dfb664ef9..1b65231de 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1280,35 +1280,29 @@ void initialize_history(Particle* p, int64_t index_source) } } -// Display message if high verbosity or trace is on - if (settings::verbosity >= 9 || simulation::trace) { - write_message("Simulating Particle " + std::to_string(p->id_)); - } + // Display message if high verbosity or trace is on + if (settings::verbosity >= 9 || simulation::trace) { + write_message("Simulating Particle " + std::to_string(p->id_)); + } - // // Initialize number of events to zero - // int n_event = 0; + // Add paricle's starting weight to count for normalizing tallies later + #pragma omp atomic + simulation::total_weight += p->wgt_; - // Add paricle's starting weight to count for normalizing tallies later - #pragma omp atomic - simulation::total_weight += p->wgt_; + // Force calculation of cross-sections by setting last energy to zero + if (settings::run_CE) { + for (auto& micro : p->neutron_xs_) micro.last_E = 0.0; + } - // Force calculation of cross-sections by setting last energy to zero - if (settings::run_CE) { - for (auto& micro : p->neutron_xs_) micro.last_E = 0.0; - } + // Prepare to write out particle track. + if (p->write_track_) add_particle_track(); - // Prepare to write out particle track. - if (p->write_track_) add_particle_track(); - - // Every particle starts with no accumulated flux derivative. - if (!model::active_tallies.empty()) zero_flux_derivs(); - - /* - std::cout << "Initialized particle " << particle_seed << " with E = " << p->E_ << " and Position {" << - p->r().x << ", " << - p->r().y << ", " << - p->r().z << "}" << std::endl; - */ + // Every particle starts with no accumulated flux derivative. + if (!model::active_tallies.empty()) + { + p->flux_derivs_.resize(model::tally_derivs.size(), 0.0); + std::fill(p->flux_derivs_.begin(), p->flux_derivs_.end(), 0.0); + } } int overall_generation() diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index c2e70ab0d..8cf571afd 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -81,13 +81,9 @@ TallyDerivative::TallyDerivative(pugi::xml_node node) void read_tally_derivatives(pugi::xml_node node) { - // Populate the derivatives array. This must be done in parallel because - // the derivatives are threadprivate. - //#pragma omp parallel - { - for (auto deriv_node : node.children("derivative")) - model::tally_derivs.emplace_back(deriv_node); - } + // Populate the derivatives array. + for (auto deriv_node : node.children("derivative")) + model::tally_derivs.emplace_back(deriv_node); // Fill the derivative map. for (auto i = 0; i < model::tally_derivs.size(); ++i) { @@ -119,8 +115,8 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, // 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; + const auto& deriv {model::tally_derivs[tally.deriv_]}; + const double& flux_deriv {p->flux_derivs_[tally.deriv_]}; // Handle special cases where we know that d_c/d_p must be zero. if (score_bin == SCORE_FLUX) { @@ -567,13 +563,15 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, } void -score_track_derivative(const Particle* p, double distance) +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; const Material& material {*model::materials[p->material_]}; - - for (auto& deriv : model::tally_derivs) { + + for (int idx = 0; idx < model::tally_derivs.size(); idx++) { + const auto& deriv {model::tally_derivs[idx]}; + double& flux_deriv = p->flux_derivs_[idx]; if (deriv.diff_material != material.id_) continue; switch (deriv.variable) { @@ -582,7 +580,7 @@ score_track_derivative(const Particle* p, double distance) // 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 * p->macro_xs_.total + flux_deriv -= distance * p->macro_xs_.total / material.density_gpcc_; break; @@ -590,7 +588,7 @@ score_track_derivative(const Particle* p, double distance) // 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 -= distance + flux_deriv -= distance * p->neutron_xs_[deriv.diff_nuclide].total; break; @@ -604,7 +602,7 @@ score_track_derivative(const Particle* p, double distance) 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) + flux_deriv -= distance * (dsig_s + dsig_a) * material.atom_density_(i); } } @@ -613,14 +611,18 @@ score_track_derivative(const Particle* p, double distance) } } -void score_collision_derivative(const Particle* p) +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; const Material& material {*model::materials[p->material_]}; - for (auto& deriv : model::tally_derivs) { + //for (auto& deriv : model::tally_derivs) { + for (int idx = 0; idx < model::tally_derivs.size(); idx++) { + const auto& deriv = model::tally_derivs[idx]; + double& flux_deriv = p->flux_derivs_[idx]; + if (deriv.diff_material != material.id_) continue; switch (deriv.variable) { @@ -629,7 +631,7 @@ void score_collision_derivative(const Particle* p) // 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_; + flux_deriv += 1. / material.density_gpcc_; break; case DIFF_NUCLIDE_DENSITY: @@ -650,7 +652,7 @@ void score_collision_derivative(const Particle* p) // (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); + flux_deriv += 1. / material.atom_density_(i); break; case DIFF_TEMPERATURE: @@ -665,7 +667,7 @@ void score_collision_derivative(const Particle* p) double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); - deriv.flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption); + 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, u'->u) = Sigma_s(E') * P(E'->E, u'->u). @@ -680,9 +682,9 @@ void score_collision_derivative(const Particle* p) } } -void zero_flux_derivs() +void zero_flux_derivs(std::vector v) { - for (auto& deriv : model::tally_derivs) deriv.flux_deriv = 0.; + std::fill(v.begin(), v.end(), 0); } }// namespace openmc diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 9b438ab34..4c1563377 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -1042,10 +1042,7 @@ setup_active_tallies() void free_memory_tally() { - //#pragma omp parallel - { - model::tally_derivs.clear(); - } + model::tally_derivs.clear(); model::tally_deriv_map.clear(); model::tally_filters.clear();