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 diff --git a/src/enrichment_calculator.cc b/src/enrichment_calculator.cc index bcf5b35..45b98f2 100644 --- a/src/enrichment_calculator.cc +++ b/src/enrichment_calculator.cc @@ -13,38 +13,54 @@ 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( - 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) { + use_downblending(use_downblending), + 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. } - cyclus::compmath::Normalize(&this->feed_composition); + cyclus::compmath::Normalize(&feed_composition); + CalculateGammaAlphaStar_(); + + BuildMatchedAbundanceRatioCascade(); +} - IsotopesNucID(isotopes); +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +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()), + gamma_235(e.gamma_235) { CalculateGammaAlphaStar_(); BuildMatchedAbundanceRatioCascade(); } @@ -75,14 +91,8 @@ EnrichmentCalculator& EnrichmentCalculator::operator= ( use_downblending = e.use_downblending; - 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(); @@ -124,18 +134,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; @@ -185,12 +190,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; @@ -201,6 +216,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_(); @@ -279,27 +301,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 dc672dd..c8bad63 100644 --- a/src/enrichment_calculator.h +++ b/src/enrichment_calculator.h @@ -17,11 +17,9 @@ 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); void PPrint(); @@ -33,11 +31,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; } @@ -63,8 +62,7 @@ class EnrichmentCalculator { double swu = 0; // Separative work that has been performed // 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; @@ -82,7 +80,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 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..c1d60d3 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); } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -304,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; @@ -318,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); } @@ -448,22 +438,23 @@ 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::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, request_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 +474,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..da8ae54 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,14 +74,12 @@ 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); - 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 b958bd1..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. @@ -51,16 +50,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; @@ -81,11 +96,13 @@ 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; - } +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) { + i = (92*1000 + i) * 10000; + } + return isotopes; } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -102,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); @@ -181,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; @@ -208,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 2196e36..cc9c262 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,7 +24,7 @@ const double kEpsDouble = 1e-5; const double kEpsCompMap = 1e-5; const int kIterMax = 200; -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 f91be7e..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); @@ -43,6 +42,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); @@ -71,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;