From 4d80d54ffe806b1342a1322ea1f2008cc9be03d5 Mon Sep 17 00:00:00 2001 From: maxschalz <58850267+maxschalz@users.noreply.github.com> Date: Wed, 24 Feb 2021 14:54:19 +0100 Subject: [PATCH 1/6] Fix typos in example input files --- input/facility.py | 2 +- input/main.py | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/input/facility.py b/input/facility.py index ec539ca..0558f03 100644 --- a/input/facility.py +++ b/input/facility.py @@ -18,7 +18,7 @@ def facility(): { "name": "DepletedUSink", "config": {"Sink": { - "in_commods": {"val": ["DepeletedU"]}, + "in_commods": {"val": ["DepletedU"]}, "recipe_name": "DepletedURecipe" }} }, diff --git a/input/main.py b/input/main.py index e9433cd..f127fda 100644 --- a/input/main.py +++ b/input/main.py @@ -1,3 +1,7 @@ +import sys +sys.path.append('./') +sys.path.append('./input') + import archetypes import commodity import control From ae6904dffc850e93e2118b158d06dc0f327c215e Mon Sep 17 00:00:00 2001 From: maxschalz <58850267+maxschalz@users.noreply.github.com> Date: Fri, 26 Feb 2021 10:40:49 +0100 Subject: [PATCH 2/6] Add copy constructor, prevent uninitialised variables This commit adds a proper copy constructor in addition to the copy assignment operator. Moreover, a lot of variables were uninitialised which has been fixed. Finally, the `isotopes` vector in `EnrichmentCalculator` now is a `const` variable. --- src/enrichment_calculator.cc | 44 ++++++++++++++++++++++-------------- src/enrichment_calculator.h | 3 ++- src/miso_helper.cc | 30 ++++++++++++++++++++++-- src/miso_helper.h | 2 ++ src/miso_helper_tests.cc | 1 + 5 files changed, 60 insertions(+), 20 deletions(-) diff --git a/src/enrichment_calculator.cc b/src/enrichment_calculator.cc index bcf5b35..17e14e0 100644 --- a/src/enrichment_calculator.cc +++ b/src/enrichment_calculator.cc @@ -13,16 +13,16 @@ namespace misoenrichment { // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -EnrichmentCalculator::EnrichmentCalculator() { - IsotopesNucID(isotopes); -} +// Constructor delegation only possible from C++11 onwards, CMake checks if +// C++11 is supported. +EnrichmentCalculator::EnrichmentCalculator() : EnrichmentCalculator(1.4) {} // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -EnrichmentCalculator::EnrichmentCalculator(double gamma_235) : - gamma_235(gamma_235) { - IsotopesNucID(isotopes); - CalculateGammaAlphaStar_(); -} +// Constructor delegation only possible from C++11 onwards, CMake checks if +// C++11 is supported. +EnrichmentCalculator::EnrichmentCalculator(double gamma_235) : + EnrichmentCalculator(misotest::comp_reprocessedU(), 0.05, 0.003, + gamma_235, 1, 1, 1e299, true) {} // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - EnrichmentCalculator::EnrichmentCalculator( @@ -37,14 +37,29 @@ EnrichmentCalculator::EnrichmentCalculator( target_feed_qty(feed_qty), target_product_qty(product_qty), max_swu(max_swu), - use_downblending(use_downblending) { + use_downblending(use_downblending), + isotopes(IsotopesNucID_vector()) { if (feed_qty==1e299 && product_qty==1e299 && max_swu==1e299) { // TODO think about whether one or two of these variables have to be // defined. Additionally, add an exception that should be thrown. } cyclus::compmath::Normalize(&this->feed_composition); + CalculateGammaAlphaStar_(); - IsotopesNucID(isotopes); + BuildMatchedAbundanceRatioCascade(); +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +EnrichmentCalculator::EnrichmentCalculator(const EnrichmentCalculator& e) : + feed_composition(e.feed_composition), + product_composition(e.product_composition), + tails_composition(e.tails_composition), + target_product_assay(e.target_product_assay), + target_tails_assay(e.target_tails_assay), + target_feed_qty(e.target_feed_qty), + target_product_qty(e.target_product_qty), max_swu(e.max_swu), + use_downblending(e.use_downblending), isotopes(IsotopesNucID_vector()), + gamma_235(e.gamma_235) { CalculateGammaAlphaStar_(); BuildMatchedAbundanceRatioCascade(); } @@ -75,14 +90,9 @@ EnrichmentCalculator& EnrichmentCalculator::operator= ( use_downblending = e.use_downblending; - IsotopesNucID(isotopes); + //IsotopesNucID(isotopes); gamma_235 = e.gamma_235; - separation_factors = CalculateSeparationFactor(gamma_235); - for (int i : isotopes) { - // E. von Halle Eq. (15) - alpha_star[i] = separation_factors[i] - / std::sqrt(separation_factors[IsotopeToNucID(235)]); - } + CalculateGammaAlphaStar_(); // TODO Check why the recalculated variables are not copied BuildMatchedAbundanceRatioCascade(); diff --git a/src/enrichment_calculator.h b/src/enrichment_calculator.h index dc672dd..e09e974 100644 --- a/src/enrichment_calculator.h +++ b/src/enrichment_calculator.h @@ -22,6 +22,7 @@ class EnrichmentCalculator { // TODO in the above constructor it might not make sense to keep the // default arguments for feed_qty and for product_qty. This will be // determined in later steps of the implementation. + EnrichmentCalculator(const EnrichmentCalculator& e); EnrichmentCalculator& operator= (const EnrichmentCalculator& e); void PPrint(); @@ -64,7 +65,7 @@ class EnrichmentCalculator { // in kg SWU timestep^-1 // TODO declare vector as const? - std::vector isotopes; + const std::vector isotopes; std::map separation_factors; std::map alpha_star; diff --git a/src/miso_helper.cc b/src/miso_helper.cc index b958bd1..7933764 100644 --- a/src/miso_helper.cc +++ b/src/miso_helper.cc @@ -51,16 +51,32 @@ cyclus::Composition::Ptr comp_depletedU() { // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - cyclus::Composition::Ptr comp_natU() { - // TODO check where this is used and if the hydrogen could be removed cyclus::CompMap comp; comp[922340000] = 5.5e-3; comp[922350000] = 0.711; comp[922380000] = 99.2835; - comp[10010000] = 10; // insert hydrogen to check if it is filtered out return cyclus::Composition::CreateFromMass(comp); }; +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +cyclus::Composition::Ptr comp_reprocessedU() { + // Composition taken from Exploring Uranium Resource Constraints on + // Fissile Material Production in Pakistan. Zia Mian, A. H. Nayyar and + // R. Rajaraman. Science and Global Security 17, 2009: p.87. + // DOI: 10.1080/08929880902975834 + + cyclus::CompMap comp; + comp[922320000] = 1.013e-10; + comp[922330000] = 2.550e-9; + comp[922340000] = 0.005; + comp[922350000] = 0.616; + comp[922360000] = 0.015; + comp[922380000] = 99.364; + + return cyclus::Composition::CreateFromMass(comp); +} + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - cyclus::Composition::Ptr comp_weapongradeU() { cyclus::CompMap comp; @@ -88,6 +104,16 @@ void IsotopesNucID(std::vector& isotopes) { } } +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +std::vector IsotopesNucID_vector() { + int iso[6] = {232, 233, 234, 235, 236, 238}; + std::vector isotopes(iso, iso + sizeof(iso)/sizeof(int)); + for (int& i : isotopes) { + i = (92*1000 + i) * 10000; + } + return isotopes; +} + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - int IsotopeToNucID(int isotope) { std::vector isotopes = {232, 233, 234, 235, 236, 238}; diff --git a/src/miso_helper.h b/src/miso_helper.h index 2196e36..40b20fc 100644 --- a/src/miso_helper.h +++ b/src/miso_helper.h @@ -14,6 +14,7 @@ namespace misotest { bool CompareCompMap(cyclus::CompMap cm1, cyclus::CompMap cm2); cyclus::Composition::Ptr comp_depletedU(); cyclus::Composition::Ptr comp_natU(); +cyclus::Composition::Ptr comp_reprocessedU(); cyclus::Composition::Ptr comp_weapongradeU(); cyclus::Material::Ptr mat_natU(); @@ -23,6 +24,7 @@ const double kEpsDouble = 1e-5; const double kEpsCompMap = 1e-5; const int kIterMax = 200; +std::vector IsotopesNucID_vector(); void IsotopesNucID(std::vector& isotopes); int IsotopeToNucID(int isotope); int NucIDToIsotope(int nuc_id); diff --git a/src/miso_helper_tests.cc b/src/miso_helper_tests.cc index f91be7e..0060aad 100644 --- a/src/miso_helper_tests.cc +++ b/src/miso_helper_tests.cc @@ -43,6 +43,7 @@ TEST(MIsoHelperTest, CheckFractionsComposition) { + 0.992835/pyne::atomic_mass(922380000)); cyclus::Composition::Ptr comp = misotest::comp_natU(); + std::cout << "TODO Add composition with non-uranium elements\n"; EXPECT_DOUBLE_EQ(MIsoAtomAssay(comp), expected_atom235); EXPECT_DOUBLE_EQ(MIsoAtomFrac(comp, 922350000), expected_atom235); EXPECT_DOUBLE_EQ(MIsoMassAssay(comp), expected_mass235); From 390e275f6d35d681df14e700cb6b65a8e21d9ca8 Mon Sep 17 00:00:00 2001 From: maxschalz <58850267+maxschalz@users.noreply.github.com> Date: Fri, 26 Feb 2021 11:08:49 +0100 Subject: [PATCH 3/6] Add pass by reference to EnrichmentCalculator::ValueFunction --- src/enrichment_calculator.cc | 36 +++++++++++++++++++++++++++--------- src/enrichment_calculator.h | 2 +- 2 files changed, 28 insertions(+), 10 deletions(-) diff --git a/src/enrichment_calculator.cc b/src/enrichment_calculator.cc index 17e14e0..ad79ae3 100644 --- a/src/enrichment_calculator.cc +++ b/src/enrichment_calculator.cc @@ -289,27 +289,45 @@ void EnrichmentCalculator::CalculateSwu_() { // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - double EnrichmentCalculator::ValueFunction_( - cyclus::CompMap composition) { + const cyclus::CompMap& composition) { const int NUCID_235 = IsotopeToNucID(235); const int NUCID_238 = IsotopeToNucID(238); - double value = 0.; - - std::vector::iterator it; - for (it = isotopes.begin(); it != isotopes.end(); it++) { - double k = (separation_factors[*it]-1) + + try { + composition.at(NUCID_235); + composition.at(NUCID_238); + } catch (const std::out_of_range& err) { + if (composition.size()==0) { + // This case can happen, e.g., during initalisation and is not a bug. + return value; + } + // Else, throw an error. + std::stringstream msg; + msg << "No U-235 or U-238 present in composition passed to " + << "'EnrichmentCalculator::ValueFunction_'."; + throw cyclus::KeyError(msg.str()); + } + + for (int i : isotopes) { + try { + composition.at(i); + } catch (const std::out_of_range& err) { + continue; + } + double k = (separation_factors[i]-1) / (separation_factors[NUCID_235]-1); if (cyclus::AlmostEq(k, 0.5)) { // This formula is not included in de la Garza 1963, it is taken // from the preceding article, see Eq. (26) in: // A. de la Garza et al., 'Multicomponent isotope separation in // cascades'. Chemical Engineering Science 15, pp. 188-209 (1961). - value += std::log(composition[*it] / composition[NUCID_238]); + value += std::log(composition.at(i) / composition.at(NUCID_238)); } else { - value += composition[*it] / (2*k - 1); + value += composition.at(i) / (2*k - 1); } } - value *= std::log(composition[NUCID_235] / composition[NUCID_238]); + value *= std::log(composition.at(NUCID_235) / composition.at(NUCID_238)); return value; } diff --git a/src/enrichment_calculator.h b/src/enrichment_calculator.h index e09e974..0f246da 100644 --- a/src/enrichment_calculator.h +++ b/src/enrichment_calculator.h @@ -83,7 +83,7 @@ class EnrichmentCalculator { void Downblend_(); void CalculateSums(double& sum_e, double& sum_s); - double ValueFunction_(cyclus::CompMap composition); + double ValueFunction_(const cyclus::CompMap& composition); }; } // namespace misoenrichment From 7c943a7fd51fd90f162d0631ba3cf9170e149631 Mon Sep 17 00:00:00 2001 From: maxschalz <58850267+maxschalz@users.noreply.github.com> Date: Tue, 9 Mar 2021 09:25:44 +0100 Subject: [PATCH 4/6] Fix memory errors Segmentation faults kept appearing. To fix this, uninitialised values have been set properly, the copying and passing of values to-and-fro `EnrichmentCalculator` and `MIsoEnrich` has been adapted and checked and minor functions such as `EnrichmentCalculator::ProductOutput` have been added. In this commit, it has also been ensured that `EnrichmentCalculator` only in- and outputs compositions under the form of `Composition::Ptr` instead of `CompMap`s. Thus, all handling of the compositions (e.g., as mass fractions or atom fractions) is done inside the calculator, independently of the other classes. This commit fixes the memory errors. Also, one error due to Cyclus is fixed (inside this module, not for the whole of Cyclus). See from line 197 onwards in `src/EnrichmentCalculator.cc`, taking into account the 'double free' error caused by `cyclus::Material::ExtractComp` and `cyclus::compmath::ApplyThreshold`. See also Cyclus issue #1524: https://github.com/cyclus/cyclus/issues/1524 --- src/enrichment_calculator.cc | 39 ++++++++++++++++++--------- src/enrichment_calculator.h | 11 ++++---- src/enrichment_calculator_tests.cc | 24 ++++++++++------- src/enrichment_calculator_tests.h | 3 ++- src/miso_enrich.cc | 43 +++++++++++++----------------- src/miso_enrich.h | 10 +++---- 6 files changed, 71 insertions(+), 59 deletions(-) diff --git a/src/enrichment_calculator.cc b/src/enrichment_calculator.cc index ad79ae3..81efa78 100644 --- a/src/enrichment_calculator.cc +++ b/src/enrichment_calculator.cc @@ -26,16 +26,17 @@ EnrichmentCalculator::EnrichmentCalculator(double gamma_235) : // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - EnrichmentCalculator::EnrichmentCalculator( - cyclus::Composition::Ptr feed_composition, + cyclus::Composition::Ptr feed_comp, double target_product_assay, double target_tails_assay, double gamma_235, double feed_qty, double product_qty, double max_swu, bool use_downblending) : - feed_composition(feed_composition->atom()), + feed_composition(feed_comp->atom()), target_product_assay(target_product_assay), target_tails_assay(target_tails_assay), gamma_235(gamma_235), target_feed_qty(feed_qty), target_product_qty(product_qty), + feed_qty(0.), product_qty(0.), max_swu(max_swu), use_downblending(use_downblending), isotopes(IsotopesNucID_vector()) { @@ -43,7 +44,7 @@ EnrichmentCalculator::EnrichmentCalculator( // TODO think about whether one or two of these variables have to be // defined. Additionally, add an exception that should be thrown. } - cyclus::compmath::Normalize(&this->feed_composition); + cyclus::compmath::Normalize(&feed_composition); CalculateGammaAlphaStar_(); BuildMatchedAbundanceRatioCascade(); @@ -134,18 +135,13 @@ void EnrichmentCalculator::SetInput( double new_feed_qty, double new_product_qty, double new_max_swu, double new_gamma_235, bool new_use_downblending) { - // This temporary variable is needed because feed_comp->atom() returns - // a const cyclus::CompMap object and hence it cannot be normalised. - // However, the normalisation is needed to be able to compare it to the - // current feed_composition - cyclus::CompMap new_compmap = new_feed_composition->atom(); - cyclus::compmath::Normalize(&new_compmap); + feed_composition = new_feed_composition->atom(); + cyclus::compmath::Normalize(&feed_composition); if (new_gamma_235 != gamma_235) { gamma_235 = new_gamma_235; CalculateGammaAlphaStar_(); } - feed_composition = new_compmap; target_product_assay = new_target_product_assay; target_tails_assay = new_target_tails_assay; @@ -195,12 +191,22 @@ void EnrichmentCalculator::SetInput( // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - void EnrichmentCalculator::EnrichmentOutput( - cyclus::CompMap& product_comp, cyclus::CompMap& tails_comp, + cyclus::Composition::Ptr& product_comp, cyclus::Composition::Ptr& tails_comp, double& feed_used, double& swu_used, double& product_produced, double& tails_produced, int& n_enrich, int& n_strip) { - product_comp = product_composition; - tails_comp = tails_composition; + // This step is needed to prevent a 'double free' error. It is caused for + // unknown reasons by cyclus::Material::ExtractComp and + // cyclus::compmath::ApplyThreshold. See also Cyclus issue #1524: + // https://github.com/cyclus/cyclus/issues/1524 + cyclus::CompMap cm; + for (const auto& x : product_composition) { + if (x.second > 0) { + cm[x.first] = x.second; + } + } + product_comp = cyclus::Composition::CreateFromAtom(cm); + tails_comp = cyclus::Composition::CreateFromAtom(tails_composition); feed_used = feed_qty; swu_used = swu; @@ -211,6 +217,13 @@ void EnrichmentCalculator::EnrichmentOutput( n_strip = n_stripping; } +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +void EnrichmentCalculator::ProductOutput( + cyclus::Composition::Ptr& old_product_comp, double& old_product_qty) { + old_product_comp = cyclus::Composition::CreateFromAtom(product_composition); + old_product_qty = product_qty; +} + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - void EnrichmentCalculator::BuildMatchedAbundanceRatioCascade() { CalculateNStages_(); diff --git a/src/enrichment_calculator.h b/src/enrichment_calculator.h index 0f246da..0edd8a3 100644 --- a/src/enrichment_calculator.h +++ b/src/enrichment_calculator.h @@ -34,11 +34,12 @@ class EnrichmentCalculator { double new_feed_qty, double new_product_qty, double new_max_swu, double gamma_235, bool use_downblending); - void EnrichmentOutput( - cyclus::CompMap& product_comp, cyclus::CompMap& tails_comp, - double& feed_used, double& swu_used, double& product_produced, - double& tails_produced, int& n_enrich, int& n_strip); - + void EnrichmentOutput(cyclus::Composition::Ptr& product_comp, + cyclus::Composition::Ptr& tails_comp, double& feed_used, + double& swu_used, double& product_produced, + double& tails_produced, int& n_enrich, int& n_strip); + void ProductOutput(cyclus::Composition::Ptr&, double&); + inline double FeedUsed() { return feed_qty; } inline double SwuUsed() { return swu; } diff --git a/src/enrichment_calculator_tests.cc b/src/enrichment_calculator_tests.cc index 4ccd1a5..61d760e 100644 --- a/src/enrichment_calculator_tests.cc +++ b/src/enrichment_calculator_tests.cc @@ -63,6 +63,8 @@ EnrichmentCalculatorTest::EnrichmentCalculatorTest() : use_downblending); e.EnrichmentOutput(product_comp, tails_comp, feed_qty, swu_used, product_qty, tails_qty, n_enriching, n_stripping); + product_cm = product_comp->atom(); + tails_cm = tails_comp->atom(); } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -70,10 +72,10 @@ EnrichmentCalculatorTest::EnrichmentCalculatorTest() : TEST_F(EnrichmentCalculatorTest, AssignmentOperator) { EnrichmentCalculator e2( cyclus::Composition::CreateFromAtom(weapons_grade_U()), 0.95, 0.1, 1.1, - 1., true); + 1., 1e299, 1e299, true); e2 = e; - cyclus::CompMap product_comp2, tails_comp2; + cyclus::Composition::Ptr product_comp2, tails_comp2; double feed_qty2, product_qty2, tails_qty2, swu_used2; int n_enriching2, n_stripping2; @@ -84,8 +86,10 @@ TEST_F(EnrichmentCalculatorTest, AssignmentOperator) { // assignment operator, but it does check the results. It is expected // that if the assignment should not have worked correctly then the // results would be wrong as well. - EXPECT_TRUE(misotest::CompareCompMap(product_comp2, product_comp)); - EXPECT_TRUE(misotest::CompareCompMap(tails_comp2, tails_comp)); + EXPECT_TRUE(misotest::CompareCompMap(product_comp2->atom(), + product_comp->atom())); + EXPECT_TRUE(misotest::CompareCompMap(tails_comp2->atom(), + tails_comp->atom())); EXPECT_DOUBLE_EQ(feed_qty2, feed_qty); EXPECT_DOUBLE_EQ(product_qty2, product_qty); EXPECT_DOUBLE_EQ(tails_qty2, tails_qty); @@ -97,8 +101,8 @@ TEST_F(EnrichmentCalculatorTest, AssignmentOperator) { // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - TEST_F(EnrichmentCalculatorTest, Concentrations) { EXPECT_TRUE(misotest::CompareCompMap(expect_product_comp, - product_comp)); - EXPECT_TRUE(misotest::CompareCompMap(expect_tails_comp, tails_comp)); + product_cm)); + EXPECT_TRUE(misotest::CompareCompMap(expect_tails_comp, tails_cm)); } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -123,8 +127,8 @@ TEST_F(EnrichmentCalculatorTest, NumberStages) { TEST_F(EnrichmentCalculatorTest, Downblending) { double target_product_assay = MIsoAssay(weapons_grade_U()) - 0.001; - cyclus::CompMap bl_product_comp, bl_tails_comp; - cyclus::CompMap bl_product_comp2; + cyclus::Composition::Ptr bl_product_comp, bl_tails_comp; + cyclus::Composition::Ptr bl_product_comp2; double bl_feed_qty, bl_product_qty, bl_tails_qty, bl_swu_used; double bl_feed_qty2, bl_product_qty2; int dummy_int; @@ -136,7 +140,7 @@ TEST_F(EnrichmentCalculatorTest, Downblending) { bl_swu_used, bl_product_qty, bl_tails_qty, dummy_int, dummy_int); - EXPECT_DOUBLE_EQ(target_product_assay, MIsoAssay(bl_product_comp)); + EXPECT_DOUBLE_EQ(target_product_assay, MIsoAtomAssay(bl_product_comp)); EXPECT_DOUBLE_EQ(expect_feed_qty, bl_feed_qty); EXPECT_TRUE(bl_product_qty > expect_product_qty); @@ -150,7 +154,7 @@ TEST_F(EnrichmentCalculatorTest, Downblending) { bl_swu_used, bl_product_qty2, bl_tails_qty, dummy_int, dummy_int); - EXPECT_DOUBLE_EQ(target_product_assay, MIsoAssay(bl_product_comp2)); + EXPECT_DOUBLE_EQ(target_product_assay, MIsoAtomAssay(bl_product_comp2)); EXPECT_DOUBLE_EQ(bl_feed_qty, bl_feed_qty2); EXPECT_DOUBLE_EQ(bl_product_qty, bl_product_qty2); } diff --git a/src/enrichment_calculator_tests.h b/src/enrichment_calculator_tests.h index db4af8e..12785c7 100644 --- a/src/enrichment_calculator_tests.h +++ b/src/enrichment_calculator_tests.h @@ -23,7 +23,8 @@ class EnrichmentCalculatorTest : public ::testing::Test { EnrichmentCalculator e; // Values to be calculated by EnrichmentCalculator - cyclus::CompMap product_comp, tails_comp; + cyclus::Composition::Ptr product_comp, tails_comp; + cyclus::CompMap product_cm, tails_cm; double feed_qty, product_qty, tails_qty, swu_used; int n_enriching, n_stripping; diff --git a/src/miso_enrich.cc b/src/miso_enrich.cc index ba02e22..5106fcc 100644 --- a/src/miso_enrich.cc +++ b/src/miso_enrich.cc @@ -31,10 +31,7 @@ MIsoEnrich::MIsoEnrich(cyclus::Context* ctx) latitude(0.0), longitude(0.0), coordinates(latitude, longitude), - use_downblending(true) { - - enrichment_calc = EnrichmentCalculator(gamma_235); -} + use_downblending(true) {} // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - MIsoEnrich::~MIsoEnrich() {} @@ -283,20 +280,17 @@ std::set::Ptr> // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - cyclus::Material::Ptr MIsoEnrich::Offer_( cyclus::Material::Ptr mat) { - cyclus::CompMap product_comp, dummy_comp; - double dummy_double; + cyclus::Composition::Ptr product_comp; double feed_qty = feed_inv[feed_idx].quantity(); + double product_assay = MIsoAtomAssay(mat); double product_qty = mat->quantity(); - int dummy_int; - - enrichment_calc.SetInput(feed_inv_comp[feed_idx], MIsoAtomAssay(mat), - tails_assay, feed_qty, product_qty, - swu_capacity, gamma_235, use_downblending); - enrichment_calc.EnrichmentOutput(product_comp, dummy_comp, dummy_double, - dummy_double, product_qty, dummy_double, - dummy_int, dummy_int); - return cyclus::Material::CreateUntracked( - product_qty, cyclus::Composition::CreateFromAtom(product_comp)); + + EnrichmentCalculator e(feed_inv_comp[feed_idx], product_assay, + tails_assay, gamma_235, feed_qty, product_qty, + swu_capacity, use_downblending); + e.ProductOutput(product_comp, product_qty); + + return cyclus::Material::CreateUntracked(product_qty, product_comp); } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -450,20 +444,21 @@ void MIsoEnrich::AddMat_(cyclus::Material::Ptr mat) { cyclus::Material::Ptr MIsoEnrich::Enrich_( cyclus::Material::Ptr mat, double qty) { - cyclus::CompMap product_comp, tails_comp; + cyclus::Composition::Ptr product_comp; + cyclus::Composition::Ptr tails_comp; double feed_required, swu_required, product_qty, tails_qty; int n_enriching, n_stripping; double feed_assay = MIsoAtomAssay(feed_inv_comp[feed_idx]); + double feed_qty = feed_inv[feed_idx].quantity(); double product_assay = MIsoAtomAssay(mat); // In the following line, the enrichment is calculated but it is not yet // performed! - enrichment_calc.SetInput(feed_inv_comp[feed_idx], product_assay, - tails_assay, feed_inv[feed_idx].quantity(), - qty, current_swu_capacity, gamma_235, - use_downblending); - enrichment_calc.EnrichmentOutput(product_comp, tails_comp, feed_required, + EnrichmentCalculator e(feed_inv_comp[feed_idx], product_assay, + tails_assay, gamma_235, feed_qty, qty, + swu_capacity, use_downblending); + e.EnrichmentOutput(product_comp, tails_comp, feed_required, swu_required, product_qty, tails_qty, n_enriching, n_stripping); // Now, perform the enrichment by popping the feed and converting it to @@ -483,8 +478,8 @@ cyclus::Material::Ptr MIsoEnrich::Enrich_( << feed_inv[feed_idx].quantity(); throw cyclus::ValueError(cyclus::Agent::InformErrorMsg(ss.str())); } - cyclus::Material::Ptr response = pop_mat->ExtractComp( - product_qty, cyclus::Composition::CreateFromAtom(product_comp)); + cyclus::Material::Ptr response = pop_mat->ExtractComp(product_qty, + product_comp); tails_inv.Push(pop_mat); current_swu_capacity -= swu_required; diff --git a/src/miso_enrich.h b/src/miso_enrich.h index 0ee8ce2..e0542fa 100644 --- a/src/miso_enrich.h +++ b/src/miso_enrich.h @@ -27,12 +27,11 @@ class SwuConverter : public cyclus::Converter { cyclus::Material::Ptr m, cyclus::Arc const * a = NULL, cyclus::ExchangeTranslationContext const * ctx = NULL) const { - EnrichmentCalculator e; double product_qty = m->quantity(); double product_assay = MIsoAtomAssay(m); - e.SetInput(feed_comp_, product_assay, tails_assay_, 1e299, product_qty, - 1e299, gamma_235_, use_downblending); + EnrichmentCalculator e(feed_comp_, product_assay, tails_assay_, gamma_235_, + 1e299, product_qty, 1e299, use_downblending); double swu_used = e.SwuUsed(); return swu_used; @@ -75,9 +74,8 @@ class FeedConverter : public cyclus::Converter { double product_qty = m->quantity(); double product_assay = MIsoAtomAssay(m); - EnrichmentCalculator e; - e.SetInput(feed_comp_, product_assay, tails_assay_, 1e299, product_qty, - 1e299, gamma_235_, use_downblending); + EnrichmentCalculator e(feed_comp_, product_assay, tails_assay_, gamma_235_, + 1e299, product_qty, 1e299, use_downblending); double feed_used = e.FeedUsed(); cyclus::toolkit::MatQuery mq(m); From 2b69debff83906322783f01271190b5e884b113d Mon Sep 17 00:00:00 2001 From: maxschalz <58850267+maxschalz@users.noreply.github.com> Date: Tue, 9 Mar 2021 09:35:31 +0100 Subject: [PATCH 5/6] Adapt IsotopesNucID function to return a vector Instead of filling an already existing vector passed by reference to `IsotopesNucID`, the function now returns a new (filled) vector. This seems neater and I do not know why I have not done it like this from the start. --- src/enrichment_calculator.cc | 5 ++--- src/miso_enrich.h | 3 +-- src/miso_helper.cc | 24 ++++++------------------ src/miso_helper.h | 3 +-- src/miso_helper_tests.cc | 6 ++---- 5 files changed, 12 insertions(+), 29 deletions(-) diff --git a/src/enrichment_calculator.cc b/src/enrichment_calculator.cc index 81efa78..45b98f2 100644 --- a/src/enrichment_calculator.cc +++ b/src/enrichment_calculator.cc @@ -39,7 +39,7 @@ EnrichmentCalculator::EnrichmentCalculator( feed_qty(0.), product_qty(0.), max_swu(max_swu), use_downblending(use_downblending), - isotopes(IsotopesNucID_vector()) { + isotopes(IsotopesNucID()) { if (feed_qty==1e299 && product_qty==1e299 && max_swu==1e299) { // TODO think about whether one or two of these variables have to be // defined. Additionally, add an exception that should be thrown. @@ -59,7 +59,7 @@ EnrichmentCalculator::EnrichmentCalculator(const EnrichmentCalculator& e) : target_tails_assay(e.target_tails_assay), target_feed_qty(e.target_feed_qty), target_product_qty(e.target_product_qty), max_swu(e.max_swu), - use_downblending(e.use_downblending), isotopes(IsotopesNucID_vector()), + use_downblending(e.use_downblending), isotopes(IsotopesNucID()), gamma_235(e.gamma_235) { CalculateGammaAlphaStar_(); BuildMatchedAbundanceRatioCascade(); @@ -91,7 +91,6 @@ EnrichmentCalculator& EnrichmentCalculator::operator= ( use_downblending = e.use_downblending; - //IsotopesNucID(isotopes); gamma_235 = e.gamma_235; CalculateGammaAlphaStar_(); diff --git a/src/miso_enrich.h b/src/miso_enrich.h index e0542fa..da8ae54 100644 --- a/src/miso_enrich.h +++ b/src/miso_enrich.h @@ -79,8 +79,7 @@ class FeedConverter : public cyclus::Converter { double feed_used = e.FeedUsed(); cyclus::toolkit::MatQuery mq(m); - std::vector isotopes; - IsotopesNucID(isotopes); + std::vector isotopes(IsotopesNucID()); std::set nucs(isotopes.begin(), isotopes.end()); double feed_uranium_frac = mq.atom_frac(nucs); diff --git a/src/miso_helper.cc b/src/miso_helper.cc index 7933764..412a7a2 100644 --- a/src/miso_helper.cc +++ b/src/miso_helper.cc @@ -13,8 +13,7 @@ namespace misotest { // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - bool CompareCompMap(cyclus::CompMap cm1, cyclus::CompMap cm2) { - std::vector isotopes; - IsotopesNucID(isotopes); + std::vector isotopes(IsotopesNucID()); std::vector::iterator it; // The following for-loop has been added to ensure that the all of the // uranium keys are present in both compmaps, else the comparison fails. @@ -97,15 +96,7 @@ cyclus::Material::Ptr mat_natU() { } // namespace misotest // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -void IsotopesNucID(std::vector& isotopes) { - isotopes = {232, 233, 234, 235, 236, 238}; - for (int i = 0; i < isotopes.size(); i++) { - isotopes[i] = (92*1000 + isotopes[i]) * 10000; - } -} - -// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -std::vector IsotopesNucID_vector() { +const std::vector IsotopesNucID() { int iso[6] = {232, 233, 234, 235, 236, 238}; std::vector isotopes(iso, iso + sizeof(iso)/sizeof(int)); for (int& i : isotopes) { @@ -128,8 +119,7 @@ int IsotopeToNucID(int isotope) { // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - int NucIDToIsotope(int nuc_id) { - std::vector isotopes; - IsotopesNucID(isotopes); + std::vector isotopes(IsotopesNucID()); std::vector::iterator it; it = std::find(isotopes.begin(), isotopes.end(), nuc_id); @@ -207,8 +197,7 @@ double MIsoAssay(cyclus::CompMap compmap) { // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - double MIsoFrac(cyclus::CompMap compmap, int isotope) { - std::vector isotopes; - IsotopesNucID(isotopes); + std::vector isotopes(IsotopesNucID()); double isotope_assay = 0; double uranium_atom_frac = 0; @@ -234,10 +223,9 @@ double MIsoFrac(cyclus::CompMap compmap, int isotope) { // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - std::map CalculateSeparationFactor(double gamma_235) { - std::vector isotopes; - IsotopesNucID(isotopes); + std::vector isotopes(IsotopesNucID()); std::map separation_factors; - + // We consider U-238 to be the key component hence the mass differences // are calculated with respect to this isotope. for (int i : isotopes) { diff --git a/src/miso_helper.h b/src/miso_helper.h index 40b20fc..cc9c262 100644 --- a/src/miso_helper.h +++ b/src/miso_helper.h @@ -24,8 +24,7 @@ const double kEpsDouble = 1e-5; const double kEpsCompMap = 1e-5; const int kIterMax = 200; -std::vector IsotopesNucID_vector(); -void IsotopesNucID(std::vector& isotopes); +const std::vector IsotopesNucID(); int IsotopeToNucID(int isotope); int NucIDToIsotope(int nuc_id); int ResBufIdx( diff --git a/src/miso_helper_tests.cc b/src/miso_helper_tests.cc index 0060aad..4d5ae69 100644 --- a/src/miso_helper_tests.cc +++ b/src/miso_helper_tests.cc @@ -25,8 +25,7 @@ TEST(MIsoHelperTest, ChooseCorrectResBuf) { // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - TEST(MIsoHelperTest, NucIDConversion) { - std::vector isotopes; - IsotopesNucID(isotopes); + std::vector isotopes(IsotopesNucID()); for (int i : isotopes) { int isotope = NucIDToIsotope(i); @@ -72,8 +71,7 @@ TEST(MIsoHelperTest, SeparationFactor) { std::map separation_factor = CalculateSeparationFactor( gamma_235); - std::vector isotopes; - IsotopesNucID(isotopes); + std::vector isotopes(IsotopesNucID()); std::map expected; expected[922320000] = 2.2; expected[922330000] = 2.0; From 21cb43e7e73b1f79b50f6ac51625011a7c2cab1c Mon Sep 17 00:00:00 2001 From: maxschalz <58850267+maxschalz@users.noreply.github.com> Date: Tue, 9 Mar 2021 09:38:36 +0100 Subject: [PATCH 6/6] Remove unnecessary comments, add minor tweaks Nothing crazy going on, adding and removing some comments and changing some variable names to improve clarity and readability of the code. --- src/enrichment_calculator.h | 8 ++------ src/miso_enrich.cc | 8 ++------ 2 files changed, 4 insertions(+), 12 deletions(-) diff --git a/src/enrichment_calculator.h b/src/enrichment_calculator.h index 0edd8a3..c8bad63 100644 --- a/src/enrichment_calculator.h +++ b/src/enrichment_calculator.h @@ -17,11 +17,8 @@ class EnrichmentCalculator { EnrichmentCalculator(cyclus::Composition::Ptr feed_comp, double target_product_assay, double target_tails_assay, double gamma, - double feed_qty=1e299, double product_qty=1e299, - double max_swu=1e299, bool use_downblending=true); - // TODO in the above constructor it might not make sense to keep the - // default arguments for feed_qty and for product_qty. This will be - // determined in later steps of the implementation. + double feed_qty, double product_qty, + double max_swu, bool use_downblending=true); EnrichmentCalculator(const EnrichmentCalculator& e); EnrichmentCalculator& operator= (const EnrichmentCalculator& e); @@ -65,7 +62,6 @@ class EnrichmentCalculator { double swu = 0; // Separative work that has been performed // in kg SWU timestep^-1 - // TODO declare vector as const? const std::vector isotopes; std::map separation_factors; std::map alpha_star; diff --git a/src/miso_enrich.cc b/src/miso_enrich.cc index 5106fcc..c1d60d3 100644 --- a/src/miso_enrich.cc +++ b/src/miso_enrich.cc @@ -298,7 +298,6 @@ bool MIsoEnrich::ValidReq_(const cyclus::Material::Ptr& req_mat) { double u_235 = MIsoAtomAssay(req_mat); double u_238 = MIsoAtomFrac(req_mat, IsotopeToNucID(238)); -// bool u_238_present = u_238 > 0 && !cyclus::AlmostEq(u_238, 0); bool u_238_present = u_238 > 0; bool not_depleted = u_235 > tails_assay; bool possible_enrichment = u_235 < max_enrich; @@ -312,9 +311,6 @@ bool SortBids(cyclus::Bid* i, cyclus::Material::Ptr mat_i = i->offer(); cyclus::Material::Ptr mat_j = j->offer(); - // TODO cycamore uses mass(U235) compared to total mass. This would also - // include possible non-U elements that are sent directly to tails. - // Because of this, they should not be considered here IMO. return MIsoAtomAssay(mat_i) <= MIsoAtomAssay(mat_j); } @@ -442,7 +438,7 @@ void MIsoEnrich::AddMat_(cyclus::Material::Ptr mat) { // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - cyclus::Material::Ptr MIsoEnrich::Enrich_( - cyclus::Material::Ptr mat, double qty) { + cyclus::Material::Ptr mat, double request_qty) { cyclus::Composition::Ptr product_comp; cyclus::Composition::Ptr tails_comp; @@ -456,7 +452,7 @@ cyclus::Material::Ptr MIsoEnrich::Enrich_( // In the following line, the enrichment is calculated but it is not yet // performed! EnrichmentCalculator e(feed_inv_comp[feed_idx], product_assay, - tails_assay, gamma_235, feed_qty, qty, + tails_assay, gamma_235, feed_qty, request_qty, swu_capacity, use_downblending); e.EnrichmentOutput(product_comp, tails_comp, feed_required, swu_required, product_qty, tails_qty,