diff --git a/.gitignore b/.gitignore index d5193fa..b362b7a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +.vscode build/ *~ src/tmp/* \ No newline at end of file diff --git a/src/.ipynb_checkpoints/CascadeConfig-checkpoint.cc b/src/.ipynb_checkpoints/CascadeConfig-checkpoint.cc new file mode 100644 index 0000000..4b5f955 --- /dev/null +++ b/src/.ipynb_checkpoints/CascadeConfig-checkpoint.cc @@ -0,0 +1,426 @@ +// Implements the CascadeEnrich class +#include "CascadeConfig.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mbmore { + +CascadeConfig::CascadeConfig() + : n_enrich(0), + n_strip(0), + n_machines(0), + feed_flow(0), + feed_assay(0), + design_product_assay(0), + design_tail_assay(0) {} + +CascadeConfig::CascadeConfig(CentrifugeConfig centrifuge_, double f_assay, + double p_assay, double t_assay, + double max_feed_flow, int max_centrifuge, + double precision) { + centrifuge = centrifuge_; + feed_assay = f_assay; + design_product_assay = p_assay; + design_tail_assay = t_assay; + + feed_flow = max_feed_flow; + n_machines = max_centrifuge; + BuildIdealCascade(f_assay, p_assay, t_assay, precision); + ScaleCascade(max_feed_flow, max_centrifuge); +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Calculate steady-state flow rates into each stage +// Linear system of equations in form AX = B, where A is nxn square matrix +// of linear equations for the flow rates of each stage and B are the external +// feeds for the stage. External feed is zero for all stages accept cascade +// feed stage (F_0) stages start with last strip stage [-2, -1, 0, 1, 2] +// http://www.physics.utah.edu/~detar/phys6720/handouts/lapack.html +// +void CascadeConfig::CalcFeedFlows() { + + int n_stages = n_enrich + n_strip; + int max_stages = n_stages; + + // Build Array with pointers + double flow_eqns[max_stages][max_stages]; + double flows[1][max_stages]; + + // build matrix of equations in this pattern + // [[ -1, 1-cut, 0, 0, 0] [[0] + // [cut, -1, 1-cut, 0, 0] [0] + // [ 0, cut, -1, 1-cut, 0] * X = [-1*feed] + // [ 0, 0, cut, -1, 1-cut] [0] + // [ 0, 0, 0, cut, -1]] [0]] + // + for (int row_idx = 0; row_idx < max_stages; row_idx++) { + // fill the array with zeros, then update individual elements as nonzero + flows[0][row_idx] = 0; + for (int fill_idx = 0; fill_idx < max_stages; fill_idx++) { + flow_eqns[fill_idx][row_idx] = 0; + } + // Required to do the artificial 'Max Stages' defn. Only calculate + // non-zero matrix elements where stages really exist. + if (row_idx < n_stages) { + int stg_i = row_idx - n_strip; + int col_idx = n_strip + stg_i; + flow_eqns[col_idx][row_idx] = -1.; + if (row_idx != 0) { + std::map::iterator it = stgs_config.find(stg_i - 1); + if (it != stgs_config.end()) { + flow_eqns[col_idx - 1][row_idx] = it->second.cut(); + } + } + if (row_idx != n_stages - 1) { + std::map::iterator it = stgs_config.find(stg_i + 1); + if (it != stgs_config.end()) { + flow_eqns[col_idx + 1][row_idx] = (1 - it->second.cut()); + } + } + // Add the external feed for the cascade + if (stg_i == 0) { + flows[0][row_idx] = -1. * feed_flow; + } + } + } + + // LAPACK solver variables + int nrhs = 1; // 1 column solution + int lda = max_stages; // must be >= MAX(1,N) + int ldb = max_stages; // must be >= MAX(1,N) + int ipiv[max_stages]; + int info; + + // Solve the linear system + dgesv_(&n_stages, &nrhs, &flow_eqns[0][0], &lda, ipiv, &flows[0][0], &ldb, + &info); + + // Check for success + if (info != 0) { + std::cerr << "LAPACK linear solver dgesv returned error " << info << "\n"; + } + + int i = 0; + std::map::iterator it; + for (it = stgs_config.begin(); it != stgs_config.end(); it++) { + it->second.feed_flow(flows[0][i]); + i++; + } +} + +void CascadeConfig::BuildIdealCascade(double f_assay, double product_assay, + double waste_assay, double precision) { + std::map ideal_stgs; + int ideal_enrich_stage = 0; + int ideal_strip_stage = 0; + + // Initialisation of Feeding stage (I == 0) + StageConfig stg(centrifuge, f_assay, precision); + int stg_i = 0; + ideal_stgs[stg_i] = stg; + double ref_alpha = ideal_stgs[0].alpha(); + double ref_du = ideal_stgs[0].DU(); + // Calculate number of enriching stages + while (stg.product_assay() < product_assay) { + stg.feed_assay(stg.product_assay()); + stg.BuildIdealStg(); + stg_i++; + ideal_stgs.insert(std::make_pair(stg_i, stg)); + } + n_enrich = stg_i + 1; + // reset + stg_i = 0; + stg = ideal_stgs[stg_i]; + // Calculate number of stripping stages + while (stg.tail_assay() > waste_assay) { + stg.feed_assay(stg.tail_assay()); + stg.BuildIdealStg(); + stg_i--; + ideal_stgs.insert(std::make_pair(stg_i, stg)); + } + n_strip = -stg_i; + stgs_config = ideal_stgs; +} + +void CascadeConfig::PopulateStages() { + int n_stages = n_enrich + n_strip; + + for (int i = 0; i < n_stages; i++) { + int curr_stage = i - n_strip; + std::map::iterator it = stgs_config.find(curr_stage); + if (it == stgs_config.end()) { + std::cout << "Bad Stage number" << std::endl; + exit(1); + } + it->second.MachinesNeededPerStage(); + } +} + +int CascadeConfig::FindTotalMachines() { + int machines = 0; + std::map::iterator it; + for (it = stgs_config.begin(); it != stgs_config.end(); it++) { + if (it->second.n_machines() == -1) it->second.MachinesNeededPerStage(); + machines += it->second.n_machines(); + } + return machines; +} + +void CascadeConfig::ScaleCascade(double max_feed, int max_centrifuges) { + // Determine the ideal steady-state feed flows for this cascade design given + // the maximum potential design feed rate + feed_flow = max_feed; + CalcFeedFlows(); + PopulateStages(); + + int machines_needed = FindTotalMachines(); + + if (max_centrifuges == -1) { + n_machines = machines_needed; + return; + } + + // Do design parameters require more centrifuges than what is available? + double max_feed_from_machine = max_feed; + while (machines_needed > max_centrifuges) { + double scaling_ratio = (double)machines_needed / (double)max_centrifuges; + max_feed_from_machine = max_feed_from_machine / scaling_ratio; + feed_flow = max_feed_from_machine; + CalcFeedFlows(); + PopulateStages(); + machines_needed = FindTotalMachines(); + } + n_machines = machines_needed; +} + +CascadeConfig CascadeConfig::ModelMisuseCascade(double f_assay, + int modeling_opt, + double precision) { + CascadeConfig misuse_cascade = (*this); + misuse_cascade.PropagateAssay(f_assay); + + // Default: alpha & cut = constant; beta = varying + // Case 1: alpha = beta, cut = varying + // Case 2: gamma & cut = constant, alpha & veta = varying + switch (modeling_opt) { + default: + misuse_cascade.ComputeAssayByAlpha(f_assay, precision); + break; + + case 1: + misuse_cascade.UpdateCut(); + misuse_cascade.UpdateFlow(); + break; + + case 2: + misuse_cascade.ComputeAssayByGamma(f_assay, precision); + break; + } + return misuse_cascade; +} + +void CascadeConfig::UpdateFlow() { + // recompute the flow according to the new cuts + (*this).CalcFeedFlows(); + + double ratio = 1; + std::map::iterator it; + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + std::map::iterator it_real = + (*this).stgs_config.find(it->first); + double max_stg_flow = + it_real->second.n_machines() * it_real->second.centrifuge.feed; + double stg_flow_ratio = it->second.feed_flow() / max_stg_flow; + if (ratio > stg_flow_ratio) { + ratio = stg_flow_ratio; + } + } + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + it->second.feed_flow(it->second.feed_flow() * ratio); + } + (*this).feed_flow *= ratio; +} + +void CascadeConfig::UpdateCut() { + std::map::iterator it; + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + it->second.CutByAlphaBeta(); + } +} + +void CascadeConfig::PropagateAssay(double f_assay) { + // Initialise Feeding stage + std::map::iterator it = (*this).stgs_config.find(0); + it->second.feed_assay(f_assay); + it->second.ProductAssay(); + it->second.TailAssay(); + + // Propagate initialisation to all stages + for (int i = 0; i < (*this).n_enrich; i++) { + it = (*this).stgs_config.find(i); + std::map::iterator it_feed; + + // Enrich stage -> feed = Product from Previous Stage + it_feed = (*this).stgs_config.find(it->first - 1); + if (it->first > 0) { + if (it_feed != (*this).stgs_config.end()) { + it->second.feed_assay(it_feed->second.product_assay()); + } + // Update Product and Tail assay from feed assay + it->second.ProductAssay(); + it->second.TailAssay(); + } + } + for (int i = 1; i <= (*this).n_strip; i++) { + it = (*this).stgs_config.find(-i); + std::map::iterator it_feed; + + // Striping stage -> feed = tails from Next Stage + it_feed = (*this).stgs_config.find(it->first + 1); + if (it->first < 0) { + if (it_feed != (*this).stgs_config.end()) { + it->second.feed_assay(it_feed->second.tail_assay()); + } + // Update Product and Tail assay from feed assay + it->second.ProductAssay(); + it->second.TailAssay(); + } + } +} + +void CascadeConfig::ComputeAssayByAlpha(double f_assay, double precision) { + CascadeConfig previous_cascade; + while (DeltaEnrichment((*this), previous_cascade, precision) > precision) { + previous_cascade = (*this); + (*this).IterateEnrichment(f_assay); + (*this).UpdateByAlpha(); + } +} + +void CascadeConfig::ComputeAssayByGamma(double f_assay, double precision) { + CascadeConfig previous_cascade; + while (DeltaEnrichment((*this), previous_cascade, precision) > precision) { + previous_cascade = (*this); + (*this).IterateEnrichment(f_assay); + (*this).UpdateByGamma(); + } +} + +double CascadeConfig::DeltaEnrichment(CascadeConfig a_enrichments, + CascadeConfig p_enrichments, + double precision) { + if (p_enrichments.n_enrich == 0) { + return 100. * std::abs(precision); + } + double square_feed_diff = 0; + double square_product_diff = 0; + double square_waste_diff = 0; + std::map::iterator it; + for (it = a_enrichments.stgs_config.begin(); + it != a_enrichments.stgs_config.end(); it++) { + int i = it->first; + std::map::iterator it2 = + p_enrichments.stgs_config.find(it->first); + if (it2 != p_enrichments.stgs_config.end()) { + square_feed_diff += + pow(it->second.feed_assay() - it2->second.feed_assay(), 2); + square_product_diff += + pow(it->second.product_assay() - it2->second.product_assay(), 2); + square_waste_diff += + pow(it->second.tail_assay() - it2->second.tail_assay(), 2); + } + } + return square_feed_diff + square_product_diff + square_waste_diff; +} + +void CascadeConfig::IterateEnrichment(double f_assay) { + CascadeConfig updated_enrichment = (*this); + + // mixing variables + double down_assay = 0; + double up_assay = 0; + double down_flow = 0; + double up_flow = 0; + double stg_feed_flow = 0; + std::map::iterator it; + + // Get the Flow and Assay quantity + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + int i = it->first; + std::map::iterator it_up = + (*this).stgs_config.find(i + 1); + std::map::iterator it_down = + (*this).stgs_config.find(i - 1); + down_assay = 0; + up_assay = 0; + down_flow = 0; + up_flow = 0; + + if (it_down != (*this).stgs_config.end()) { + down_assay = it_down->second.product_assay(); + down_flow = it_down->second.feed_flow() * it_down->second.cut(); + } + if (it_up != (*this).stgs_config.end()) { + up_assay = it_up->second.tail_assay(); + up_flow = it_up->second.feed_flow() * (1 - it_up->second.cut()); + } + + // Mix the Product and the Tail to have the correct Feed Assay + double stg_f_assay = + (down_assay * down_flow + up_assay * up_flow) / (down_flow + up_flow); + if (i == 0) { // add Feed flow in the entry stage + stg_f_assay = (down_assay * down_flow + up_assay * up_flow + + f_assay * (*this).feed_flow) / + (down_flow + up_flow + (*this).feed_flow); + stg_feed_flow = down_flow + up_flow + (*this).feed_flow; + } + + std::map::iterator it_new = + updated_enrichment.stgs_config.find(i); + + // Update Stage feed assay + it_new->second.feed_assay(stg_f_assay); + } + + (*this).stgs_config = updated_enrichment.stgs_config; +} + +void CascadeConfig::UpdateByAlpha(){ + std::map::iterator it; + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + // Update Beta values (from feed) -- Alpha & Cut are cte + it->second.BetaByAlphaAndCut(); + // Recompute Product Assay and Tail Assay + it->second.ProductAssay(); + it->second.TailAssay(); + } +} + +void CascadeConfig::UpdateByGamma(){ + std::map::iterator it; + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + double gamma = it->second.alpha() * it->second.beta(); + // Recompute Product Assay + it->second.ProductAssayByGamma(gamma); + // Alpha by Product assay + it->second.AlphaByProductAssay(); + // Beta and tail assay... + it->second.BetaByAlphaAndCut(); + it->second.TailAssay(); + } +} + +} // namespace mbmore diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 23ea67f..29c1b60 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -3,7 +3,12 @@ USE_CYCLUS("mbmore" "mytest") USE_CYCLUS("mbmore" "behavior_functions") USE_CYCLUS("mbmore" "RandomEnrich") + USE_CYCLUS("mbmore" "CentrifugeConfig") +USE_CYCLUS("mbmore" "StageConfig") +USE_CYCLUS("mbmore" "CascadeConfig") +USE_CYCLUS("mbmore" "CascadeEnrich") + USE_CYCLUS("mbmore" "RandomSink") USE_CYCLUS("mbmore" "StateInst") USE_CYCLUS("mbmore" "InteractRegion") diff --git a/src/CascadeConfig.cc b/src/CascadeConfig.cc new file mode 100644 index 0000000..4b5f955 --- /dev/null +++ b/src/CascadeConfig.cc @@ -0,0 +1,426 @@ +// Implements the CascadeEnrich class +#include "CascadeConfig.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mbmore { + +CascadeConfig::CascadeConfig() + : n_enrich(0), + n_strip(0), + n_machines(0), + feed_flow(0), + feed_assay(0), + design_product_assay(0), + design_tail_assay(0) {} + +CascadeConfig::CascadeConfig(CentrifugeConfig centrifuge_, double f_assay, + double p_assay, double t_assay, + double max_feed_flow, int max_centrifuge, + double precision) { + centrifuge = centrifuge_; + feed_assay = f_assay; + design_product_assay = p_assay; + design_tail_assay = t_assay; + + feed_flow = max_feed_flow; + n_machines = max_centrifuge; + BuildIdealCascade(f_assay, p_assay, t_assay, precision); + ScaleCascade(max_feed_flow, max_centrifuge); +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Calculate steady-state flow rates into each stage +// Linear system of equations in form AX = B, where A is nxn square matrix +// of linear equations for the flow rates of each stage and B are the external +// feeds for the stage. External feed is zero for all stages accept cascade +// feed stage (F_0) stages start with last strip stage [-2, -1, 0, 1, 2] +// http://www.physics.utah.edu/~detar/phys6720/handouts/lapack.html +// +void CascadeConfig::CalcFeedFlows() { + + int n_stages = n_enrich + n_strip; + int max_stages = n_stages; + + // Build Array with pointers + double flow_eqns[max_stages][max_stages]; + double flows[1][max_stages]; + + // build matrix of equations in this pattern + // [[ -1, 1-cut, 0, 0, 0] [[0] + // [cut, -1, 1-cut, 0, 0] [0] + // [ 0, cut, -1, 1-cut, 0] * X = [-1*feed] + // [ 0, 0, cut, -1, 1-cut] [0] + // [ 0, 0, 0, cut, -1]] [0]] + // + for (int row_idx = 0; row_idx < max_stages; row_idx++) { + // fill the array with zeros, then update individual elements as nonzero + flows[0][row_idx] = 0; + for (int fill_idx = 0; fill_idx < max_stages; fill_idx++) { + flow_eqns[fill_idx][row_idx] = 0; + } + // Required to do the artificial 'Max Stages' defn. Only calculate + // non-zero matrix elements where stages really exist. + if (row_idx < n_stages) { + int stg_i = row_idx - n_strip; + int col_idx = n_strip + stg_i; + flow_eqns[col_idx][row_idx] = -1.; + if (row_idx != 0) { + std::map::iterator it = stgs_config.find(stg_i - 1); + if (it != stgs_config.end()) { + flow_eqns[col_idx - 1][row_idx] = it->second.cut(); + } + } + if (row_idx != n_stages - 1) { + std::map::iterator it = stgs_config.find(stg_i + 1); + if (it != stgs_config.end()) { + flow_eqns[col_idx + 1][row_idx] = (1 - it->second.cut()); + } + } + // Add the external feed for the cascade + if (stg_i == 0) { + flows[0][row_idx] = -1. * feed_flow; + } + } + } + + // LAPACK solver variables + int nrhs = 1; // 1 column solution + int lda = max_stages; // must be >= MAX(1,N) + int ldb = max_stages; // must be >= MAX(1,N) + int ipiv[max_stages]; + int info; + + // Solve the linear system + dgesv_(&n_stages, &nrhs, &flow_eqns[0][0], &lda, ipiv, &flows[0][0], &ldb, + &info); + + // Check for success + if (info != 0) { + std::cerr << "LAPACK linear solver dgesv returned error " << info << "\n"; + } + + int i = 0; + std::map::iterator it; + for (it = stgs_config.begin(); it != stgs_config.end(); it++) { + it->second.feed_flow(flows[0][i]); + i++; + } +} + +void CascadeConfig::BuildIdealCascade(double f_assay, double product_assay, + double waste_assay, double precision) { + std::map ideal_stgs; + int ideal_enrich_stage = 0; + int ideal_strip_stage = 0; + + // Initialisation of Feeding stage (I == 0) + StageConfig stg(centrifuge, f_assay, precision); + int stg_i = 0; + ideal_stgs[stg_i] = stg; + double ref_alpha = ideal_stgs[0].alpha(); + double ref_du = ideal_stgs[0].DU(); + // Calculate number of enriching stages + while (stg.product_assay() < product_assay) { + stg.feed_assay(stg.product_assay()); + stg.BuildIdealStg(); + stg_i++; + ideal_stgs.insert(std::make_pair(stg_i, stg)); + } + n_enrich = stg_i + 1; + // reset + stg_i = 0; + stg = ideal_stgs[stg_i]; + // Calculate number of stripping stages + while (stg.tail_assay() > waste_assay) { + stg.feed_assay(stg.tail_assay()); + stg.BuildIdealStg(); + stg_i--; + ideal_stgs.insert(std::make_pair(stg_i, stg)); + } + n_strip = -stg_i; + stgs_config = ideal_stgs; +} + +void CascadeConfig::PopulateStages() { + int n_stages = n_enrich + n_strip; + + for (int i = 0; i < n_stages; i++) { + int curr_stage = i - n_strip; + std::map::iterator it = stgs_config.find(curr_stage); + if (it == stgs_config.end()) { + std::cout << "Bad Stage number" << std::endl; + exit(1); + } + it->second.MachinesNeededPerStage(); + } +} + +int CascadeConfig::FindTotalMachines() { + int machines = 0; + std::map::iterator it; + for (it = stgs_config.begin(); it != stgs_config.end(); it++) { + if (it->second.n_machines() == -1) it->second.MachinesNeededPerStage(); + machines += it->second.n_machines(); + } + return machines; +} + +void CascadeConfig::ScaleCascade(double max_feed, int max_centrifuges) { + // Determine the ideal steady-state feed flows for this cascade design given + // the maximum potential design feed rate + feed_flow = max_feed; + CalcFeedFlows(); + PopulateStages(); + + int machines_needed = FindTotalMachines(); + + if (max_centrifuges == -1) { + n_machines = machines_needed; + return; + } + + // Do design parameters require more centrifuges than what is available? + double max_feed_from_machine = max_feed; + while (machines_needed > max_centrifuges) { + double scaling_ratio = (double)machines_needed / (double)max_centrifuges; + max_feed_from_machine = max_feed_from_machine / scaling_ratio; + feed_flow = max_feed_from_machine; + CalcFeedFlows(); + PopulateStages(); + machines_needed = FindTotalMachines(); + } + n_machines = machines_needed; +} + +CascadeConfig CascadeConfig::ModelMisuseCascade(double f_assay, + int modeling_opt, + double precision) { + CascadeConfig misuse_cascade = (*this); + misuse_cascade.PropagateAssay(f_assay); + + // Default: alpha & cut = constant; beta = varying + // Case 1: alpha = beta, cut = varying + // Case 2: gamma & cut = constant, alpha & veta = varying + switch (modeling_opt) { + default: + misuse_cascade.ComputeAssayByAlpha(f_assay, precision); + break; + + case 1: + misuse_cascade.UpdateCut(); + misuse_cascade.UpdateFlow(); + break; + + case 2: + misuse_cascade.ComputeAssayByGamma(f_assay, precision); + break; + } + return misuse_cascade; +} + +void CascadeConfig::UpdateFlow() { + // recompute the flow according to the new cuts + (*this).CalcFeedFlows(); + + double ratio = 1; + std::map::iterator it; + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + std::map::iterator it_real = + (*this).stgs_config.find(it->first); + double max_stg_flow = + it_real->second.n_machines() * it_real->second.centrifuge.feed; + double stg_flow_ratio = it->second.feed_flow() / max_stg_flow; + if (ratio > stg_flow_ratio) { + ratio = stg_flow_ratio; + } + } + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + it->second.feed_flow(it->second.feed_flow() * ratio); + } + (*this).feed_flow *= ratio; +} + +void CascadeConfig::UpdateCut() { + std::map::iterator it; + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + it->second.CutByAlphaBeta(); + } +} + +void CascadeConfig::PropagateAssay(double f_assay) { + // Initialise Feeding stage + std::map::iterator it = (*this).stgs_config.find(0); + it->second.feed_assay(f_assay); + it->second.ProductAssay(); + it->second.TailAssay(); + + // Propagate initialisation to all stages + for (int i = 0; i < (*this).n_enrich; i++) { + it = (*this).stgs_config.find(i); + std::map::iterator it_feed; + + // Enrich stage -> feed = Product from Previous Stage + it_feed = (*this).stgs_config.find(it->first - 1); + if (it->first > 0) { + if (it_feed != (*this).stgs_config.end()) { + it->second.feed_assay(it_feed->second.product_assay()); + } + // Update Product and Tail assay from feed assay + it->second.ProductAssay(); + it->second.TailAssay(); + } + } + for (int i = 1; i <= (*this).n_strip; i++) { + it = (*this).stgs_config.find(-i); + std::map::iterator it_feed; + + // Striping stage -> feed = tails from Next Stage + it_feed = (*this).stgs_config.find(it->first + 1); + if (it->first < 0) { + if (it_feed != (*this).stgs_config.end()) { + it->second.feed_assay(it_feed->second.tail_assay()); + } + // Update Product and Tail assay from feed assay + it->second.ProductAssay(); + it->second.TailAssay(); + } + } +} + +void CascadeConfig::ComputeAssayByAlpha(double f_assay, double precision) { + CascadeConfig previous_cascade; + while (DeltaEnrichment((*this), previous_cascade, precision) > precision) { + previous_cascade = (*this); + (*this).IterateEnrichment(f_assay); + (*this).UpdateByAlpha(); + } +} + +void CascadeConfig::ComputeAssayByGamma(double f_assay, double precision) { + CascadeConfig previous_cascade; + while (DeltaEnrichment((*this), previous_cascade, precision) > precision) { + previous_cascade = (*this); + (*this).IterateEnrichment(f_assay); + (*this).UpdateByGamma(); + } +} + +double CascadeConfig::DeltaEnrichment(CascadeConfig a_enrichments, + CascadeConfig p_enrichments, + double precision) { + if (p_enrichments.n_enrich == 0) { + return 100. * std::abs(precision); + } + double square_feed_diff = 0; + double square_product_diff = 0; + double square_waste_diff = 0; + std::map::iterator it; + for (it = a_enrichments.stgs_config.begin(); + it != a_enrichments.stgs_config.end(); it++) { + int i = it->first; + std::map::iterator it2 = + p_enrichments.stgs_config.find(it->first); + if (it2 != p_enrichments.stgs_config.end()) { + square_feed_diff += + pow(it->second.feed_assay() - it2->second.feed_assay(), 2); + square_product_diff += + pow(it->second.product_assay() - it2->second.product_assay(), 2); + square_waste_diff += + pow(it->second.tail_assay() - it2->second.tail_assay(), 2); + } + } + return square_feed_diff + square_product_diff + square_waste_diff; +} + +void CascadeConfig::IterateEnrichment(double f_assay) { + CascadeConfig updated_enrichment = (*this); + + // mixing variables + double down_assay = 0; + double up_assay = 0; + double down_flow = 0; + double up_flow = 0; + double stg_feed_flow = 0; + std::map::iterator it; + + // Get the Flow and Assay quantity + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + int i = it->first; + std::map::iterator it_up = + (*this).stgs_config.find(i + 1); + std::map::iterator it_down = + (*this).stgs_config.find(i - 1); + down_assay = 0; + up_assay = 0; + down_flow = 0; + up_flow = 0; + + if (it_down != (*this).stgs_config.end()) { + down_assay = it_down->second.product_assay(); + down_flow = it_down->second.feed_flow() * it_down->second.cut(); + } + if (it_up != (*this).stgs_config.end()) { + up_assay = it_up->second.tail_assay(); + up_flow = it_up->second.feed_flow() * (1 - it_up->second.cut()); + } + + // Mix the Product and the Tail to have the correct Feed Assay + double stg_f_assay = + (down_assay * down_flow + up_assay * up_flow) / (down_flow + up_flow); + if (i == 0) { // add Feed flow in the entry stage + stg_f_assay = (down_assay * down_flow + up_assay * up_flow + + f_assay * (*this).feed_flow) / + (down_flow + up_flow + (*this).feed_flow); + stg_feed_flow = down_flow + up_flow + (*this).feed_flow; + } + + std::map::iterator it_new = + updated_enrichment.stgs_config.find(i); + + // Update Stage feed assay + it_new->second.feed_assay(stg_f_assay); + } + + (*this).stgs_config = updated_enrichment.stgs_config; +} + +void CascadeConfig::UpdateByAlpha(){ + std::map::iterator it; + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + // Update Beta values (from feed) -- Alpha & Cut are cte + it->second.BetaByAlphaAndCut(); + // Recompute Product Assay and Tail Assay + it->second.ProductAssay(); + it->second.TailAssay(); + } +} + +void CascadeConfig::UpdateByGamma(){ + std::map::iterator it; + for (it = (*this).stgs_config.begin(); it != (*this).stgs_config.end(); + it++) { + double gamma = it->second.alpha() * it->second.beta(); + // Recompute Product Assay + it->second.ProductAssayByGamma(gamma); + // Alpha by Product assay + it->second.AlphaByProductAssay(); + // Beta and tail assay... + it->second.BetaByAlphaAndCut(); + it->second.TailAssay(); + } +} + +} // namespace mbmore diff --git a/src/CascadeConfig.h b/src/CascadeConfig.h new file mode 100644 index 0000000..9116512 --- /dev/null +++ b/src/CascadeConfig.h @@ -0,0 +1,94 @@ +#ifndef MBMORE_SRC_CASCADE_H_ +#define MBMORE_SRC_CASCADE_H_ + +#include +#include + +#include "CentrifugeConfig.h" +#include "StageConfig.h" + +namespace mbmore { +// LAPACK solver for system of linear equations +extern "C" { +void dgesv_(int *n, int *nrhs, double *a, int *lda, int *ipivot, double *b, + int *ldb, int *info); +} + +class CascadeConfig { + public: + CascadeConfig() ; + CascadeConfig(CentrifugeConfig centrifuge, double f_assay, double p_assay, + double t_assay, double max_feed_flow, int max_centrifuge, + double precision = 1e-8); + // Build a full cascade such as all stage follow alpha = beta = const. Get + // alpha/beta value from feeding stage. From the design feed/product/assay + void BuildIdealCascade(double f_assay, double p_assay, double w_assay, + double precision = 1e-8); + // Get the total number of machine in the Cascade + int FindTotalMachines(); + + // Solve the flow matrix from the stages cuts + void CalcFeedFlows(); + // Establish machines needed in each stage + void PopulateStages(); + + // Scale the Cascade to meet the limitation in max feed or max centrifuges + void ScaleCascade(double max_feed, int max_centrifuges); + + CascadeConfig ModelMisuseCascade(double f_assay, int modeling_opt = 0, double precision = 1e-8); + + // Compute the response of the cascade to a non ideal feed assay + void PropagateAssay(double f_assay); + + + // Propagate iterratively the assay assuming + // Alpha constant + void ComputeAssayByAlpha(double f_assay, double precision = 1e-8); + // Gamma constant + void ComputeAssayByGamma(double f_assay, double precision = 1e-8); + + + void UpdateCut(); + void UpdateFlow(); + + double FeedFlow() { return feed_flow; } + // Configuration of the centrifuges in the stages + CentrifugeConfig centrifuge; + // Map of all the stage configuration + std::map stgs_config; + // number of enrich stages + int n_enrich; + // number of stripping stages + int n_strip; + + private: + // total number of machine in the Cascade + int n_machines; + // real feed flow (constrained by the cascade design/total number of + // machine/max feed flow) + double feed_flow; + + //design feed assay + double feed_assay; + //design product assay + double design_product_assay; + //design tail assay + double design_tail_assay; + + // Method to check the assays different between 2 cascades + double DeltaEnrichment(CascadeConfig actual_enrichments, + CascadeConfig previous_enrichment, + double precision); + + // method computing one iteration, of the algorithm used to get the response + // to non ideal feed assay + void IterateEnrichment(double feed_assay); + // recompute product assay, beta and tails assay from alpha, cut and feed assay + void UpdateByAlpha(); + //recompute alpha, product assay, beta and tails assay from gamma, cut and feed assay + void UpdateByGamma(); +}; + +} // namespace mbmore + +#endif // MBMORE_SRC_CASCADE_H_ diff --git a/src/CascadeConfig_tests.cc b/src/CascadeConfig_tests.cc new file mode 100644 index 0000000..ca7a93d --- /dev/null +++ b/src/CascadeConfig_tests.cc @@ -0,0 +1,202 @@ +#include + +#include "CascadeConfig.h" + +#include "agent_tests.h" +#include "context.h" +#include "facility_tests.h" + +namespace mbmore { + +// Benchmarked using a regression test (expected values calculated manually) +namespace cascadeconfig_tests { +// Fixed for a cascade separating out U235 from U238 in UF6 gas +const double M = 0.352; // kg/mol UF6 +const double dM = 0.003; // kg/mol U238 - U235 +const double x = 1000; // Pressure ratio (Glaser) + +// General cascade assumptions +const double flow_ratio = 2.0; +const double cut = 0.5; + +// Centrifuge parameters based on Glaser SGS 2009 paper +const double v_a = 485; // m/s +const double height = 0.5; // meters +const double diameter = 0.15; // meters +const double feed_m = 15 * 60 * 60 / ((1e3) * 60 * 60 * 1000.0); // kg/sec +const double temp = 320.0; // Kelvin + +// Cascade params used in calculating expected values +const double feed_assay = 0.0071; +const double product_assay = 0.035; +const double waste_assay = 0.001; +const double feed_c = 739 / (30.4 * 24 * 60 * 60); // kg/month -> kg/sec +const double product_c = 77 / (30.4 * 24 * 60 * 60); // kg/month -> kg/sec +CentrifugeConfig centrifuge(v_a, height, diameter, feed_m, temp, M, dM, x, + flow_ratio); + +// del U=8.638345e-08 alpha=1.130517 +double delU = centrifuge.ComputeDeltaU(cut); + +const double tol_assay = 1e-5; +const double tol_qty = 1e-6; +const double tol_num = 1e-2; + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Ideal cascade design, and then using away from ideal design +// +TEST(CascadeStage_Test, TestCascade) { + // These expected values are found by regression test + // and are working as intended as of committing this line. + CascadeConfig cascade; + cascade.centrifuge = centrifuge; + cascade.BuildIdealCascade(feed_assay, product_assay, waste_assay, 1e-8); + int expected_n_enrich_stage = 13; + int expected_n_strip_stage = 15; + // integer + int n_stage_enrich = cascade.n_enrich; + int n_stage_waste = cascade.n_strip; + + EXPECT_EQ(n_stage_enrich, expected_n_enrich_stage); + EXPECT_EQ(n_stage_waste, expected_n_strip_stage); + + // Now test assays when cascade is modified away from ideal design + // (cascade optimized for natural uranium feed, now use 20% enriched) + double feed_assay_mod = 0.20; + cascade.ScaleCascade(feed_c, 1000000); + CascadeConfig cascade_non_ideal = + cascade.ModelMisuseCascade(feed_assay_mod, 0, 1e-31); + + double mod_product_assay = + cascade_non_ideal.stgs_config[n_stage_enrich - 1].product_assay(); + double mod_waste_assay = + cascade_non_ideal.stgs_config[-n_stage_waste].product_assay(); + + double expected_mod_product_assay = 0.791461; + EXPECT_NEAR(mod_product_assay, expected_mod_product_assay, tol_assay); + + double expected_mod_waste_assay = 0.097997; + EXPECT_NEAR(mod_waste_assay, expected_mod_waste_assay, tol_assay); +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Tests the steady state flow rates for a cascade +// +TEST(CascadeStage_Test, TestCascadeDesign) { + double fa = 0.10; + double pa = 0.20; + double wa = 0.05; + + std::vector expected_flows = { + 0.00030693, 0.00061387, 0.0009208, 0.00122774, 0.00153467, + 0.00127889, 0.00102311, 0.00076734, 0.00051156, 0.00025578}; + + std::vector expected_machines = {25, 46, 66, 84, 100, 115, + 94, 74, 57, 41}; + + CascadeConfig cascade(centrifuge, fa, pa, wa, feed_c, 1000000); + + for (int i = 0; i < expected_flows.size(); i++) { + EXPECT_NEAR(cascade.stgs_config[i - cascade.n_strip].feed_flow(), + expected_flows[i], tol_num); + int nmach = cascade.stgs_config[i - cascade.n_strip].n_machines(); + EXPECT_EQ(nmach, expected_machines[i]); + } + + // not enough machines + int max_centrifuges = 40; + cascade.ScaleCascade(feed_c, max_centrifuges); + int expected_tot_mach = 40; + double expected_opt_feed = 1.261064e-05; + + EXPECT_EQ(expected_tot_mach, cascade.FindTotalMachines()); + EXPECT_NEAR(expected_opt_feed, cascade.FeedFlow(), tol_qty); + + // more machines than requested capacity + max_centrifuges = 800; + cascade.ScaleCascade(feed_c, max_centrifuges); + expected_tot_mach = 741; + expected_opt_feed = 0.0002814; + + EXPECT_EQ(expected_tot_mach, cascade.FindTotalMachines()); + EXPECT_NEAR(expected_opt_feed, cascade.FeedFlow(), tol_qty); +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Tests feed flows in misue model A, where alpha & cut are constant +// +TEST(CascadeStage_Test, TestUpdateAssay) { + double fa = 0.10; + double pa = 0.20; + double wa = 0.05; + + CascadeConfig cascade(centrifuge, fa, pa, wa, feed_c, 100); + double product_assay = + cascade.stgs_config[cascade.n_enrich - 1].product_assay(); + double tail_assay = cascade.stgs_config[-cascade.n_strip].tail_assay(); + double product_flow = cascade.stgs_config[cascade.n_enrich - 1].feed_flow() * + cascade.stgs_config[cascade.n_enrich - 1].cut(); + double tail_flow = cascade.stgs_config[-cascade.n_strip].feed_flow() * + (1 - cascade.stgs_config[-cascade.n_strip].cut()); + + double feed_from_assay = + product_flow * (product_assay - tail_assay) / (fa - tail_assay); + double tail_from_assay = + product_flow * (product_assay - fa) / (fa - tail_assay); + + EXPECT_NEAR(cascade.FeedFlow(), feed_from_assay, 1e-3); + EXPECT_NEAR(tail_flow, tail_from_assay, 1e-3); + + fa = 0.2; + cascade = cascade.ModelMisuseCascade(fa, 0, 1e-17); + product_assay = cascade.stgs_config[cascade.n_enrich - 1].product_assay(); + tail_assay = cascade.stgs_config[-cascade.n_strip].tail_assay(); + product_flow = cascade.stgs_config[cascade.n_enrich - 1].feed_flow() * + cascade.stgs_config[cascade.n_enrich - 1].cut(); + tail_flow = cascade.stgs_config[-cascade.n_strip].feed_flow() * + (1 - cascade.stgs_config[-cascade.n_strip].cut()); + feed_from_assay = + product_flow * (product_assay - tail_assay) / (fa - tail_assay); + tail_from_assay = product_flow * (product_assay - fa) / (fa - tail_assay); + + EXPECT_NEAR(cascade.FeedFlow(), feed_from_assay, 1e-3); + EXPECT_NEAR(tail_flow, tail_from_assay, 1e-3); +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Test feed flows in misuse model B, where alpha = beta +// +TEST(CascadeStage_Test, TestUpdateAlphaBetaFix) { + double fa = 0.10; + double pa = 0.20; + double wa = 0.05; + + CascadeConfig cascade(centrifuge, fa, pa, wa, feed_c, 100); + double product_assay = + cascade.stgs_config[cascade.n_enrich - 1].product_assay(); + double tail_assay = cascade.stgs_config[-cascade.n_strip].tail_assay(); + double product_flow = cascade.stgs_config[cascade.n_enrich - 1].feed_flow() * + cascade.stgs_config[cascade.n_enrich - 1].cut(); + double tail_flow = cascade.stgs_config[-cascade.n_strip].feed_flow() * + (1 - cascade.stgs_config[-cascade.n_strip].cut()); + + double feed_from_assay = + product_flow * (product_assay - tail_assay) / (fa - tail_assay); + double tail_from_assay = + product_flow * (product_assay - fa) / (fa - tail_assay); + + EXPECT_NEAR(cascade.FeedFlow(), feed_from_assay, 1e-3); + EXPECT_NEAR(tail_flow, tail_from_assay, 1e-3); + + fa = 0.2; + cascade = cascade.ModelMisuseCascade(fa, 1, 1e-17); + double alpha_ref = cascade.stgs_config[0].alpha(); + std::map::iterator it; + for (it = cascade.stgs_config.begin(); it != cascade.stgs_config.end(); it++){ + EXPECT_EQ(alpha_ref, it->second.alpha()); + EXPECT_EQ(alpha_ref, it->second.beta()); + } +} + +} // namespace cascadeconfig_tests +} // namespace mbmore diff --git a/src/CascadeEnrich.cc b/src/CascadeEnrich.cc index 64dc116..99d7c8c 100644 --- a/src/CascadeEnrich.cc +++ b/src/CascadeEnrich.cc @@ -1,48 +1,48 @@ // Implements the CascadeEnrich class #include "CascadeEnrich.h" #include "behavior_functions.h" -#include "enrich_functions.h" #include "sim_init.h" #include +#include #include #include #include #include -#include - namespace mbmore { // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CascadeEnrich::CascadeEnrich(cyclus::Context* ctx) +CascadeEnrich::CascadeEnrich(cyclus::Context* ctx) : cyclus::Facility(ctx), - feed_recipe(""), - max_centrifuges(), - design_feed_assay(), - design_product_assay(), - design_tails_assay(), - centrifuge_velocity(485.0), - temp(320.0), - height(0.5), - diameter(0.15), - machine_feed(15), - max_enrich(1), - design_feed_flow(0), - feed_commod(""), - product_commod(""), - tails_commod(""), - order_prefs(true) {} + feed_recipe(""), + max_centrifuges(), + design_feed_assay(), + design_product_assay(), + design_tails_assay(), + centrifuge_velocity(485.0), + temp(320.0), + height(0.5), + diameter(0.15), + machine_feed(15), + max_enrich(1), + design_feed_flow(100), + L_over_F(2), + feed_commod(""), + product_commod(""), + tails_commod(""), + misuse_model(0), + order_prefs(true) { + secpertimestep = (*this).context()->sim_info().dt; + } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CascadeEnrich::~CascadeEnrich() {} // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - std::string CascadeEnrich::str() { std::stringstream ss; - ss << cyclus::Facility::str() - << " with enrichment facility parameters:" - << " * SWU capacity: " << SwuCapacity() - << " * Tails assay: " << tails_assay << " * Feed assay: " << FeedAssay() + ss << cyclus::Facility::str() << " with enrichment facility parameters:" + << " * Tails assay: " << design_tails_assay << " * Feed assay: " << design_feed_assay << " * Input cyclus::Commodity: " << feed_commod << " * Output cyclus::Commodity: " << product_commod << " * Tails cyclus::Commodity: " << tails_commod; @@ -50,72 +50,60 @@ std::string CascadeEnrich::str() { } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -void CascadeEnrich::Build(cyclus::Agent* parent) { +void CascadeEnrich::EnterNotify() { using cyclus::Material; - - tails_assay = design_tails_assay; - - // Calculate ideal machine performance - double design_delU = CalcDelU(centrifuge_velocity, height, diameter, - Mg2kgPerSec(machine_feed), temp, - cut, eff, M, dM, x, flow_internal); - double design_alpha = AlphaBySwu(design_delU, Mg2kgPerSec(machine_feed), - cut, M); - - // Design ideal cascade based on target feed flow and product assay - std::pair n_stages = - FindNStages(design_alpha, design_feed_assay, design_product_assay, - design_tails_assay); - - // TODO DELETE THIS, STAGES ARE ALREADY INTS - // set as internal state variables - // int truncates but we need # of stages to assure target values, - // so if the number is 5.1 we need 6. - // n_enrich_stages = int(n_stages.first) + 1; - // n_strip_stages = int(n_stages.second) + 1; - n_enrich_stages = n_stages.first; - n_strip_stages = n_stages.second; - - std::pair cascade_info = DesignCascade(FlowPerSec(design_feed_flow), - design_alpha, - design_delU, - cut, max_centrifuges, - n_stages); - - max_feed_inventory = FlowPerMon(cascade_info.second); - // Number of machines times swu per machine - SwuCapacity(cascade_info.first * FlowPerMon(design_delU)); - - Facility::Build(parent); + cyclus::Facility::EnterNotify(); + centrifuge = CentrifugeConfig(); + // Update Centrifuge parameter from the user input: + centrifuge.v_a = centrifuge_velocity; + centrifuge.height = height; + centrifuge.diameter = diameter; + centrifuge.feed = machine_feed / 1000 / 1000; + centrifuge.temp = temp; + centrifuge.flow_ratio = L_over_F; + secpertimestep = (*this).context()->sim_info().dt; + + cascade = CascadeConfig(centrifuge, design_feed_assay, design_product_assay, + design_tails_assay, design_feed_flow, + max_centrifuges, precision); + + std::map::iterator it; + for (it = cascade.stgs_config.begin(); it != cascade.stgs_config.end(); + it++) { + std::cout << "stg " << it->first; + std::cout << " FA: " << it->second.feed_assay(); + std::cout << " PA: " << it->second.product_assay(); + std::cout << " TA: " << it->second.tail_assay(); + std::cout << " feed_flow: " << it->second.feed_flow(); + std::cout << " cut: " << it->second.cut(); + std::cout << " alpha: " << it->second.alpha(); + std::cout << " beta: " << it->second.beta(); + std::cout << " machine: " << it->second.n_machines(); + std::cout << std::endl; + } + std::cout << "Dsign Feed Flow " << FlowPerDt(cascade.FeedFlow()) << std::endl; + if (max_feed_inventory > 0) { + inventory.capacity(max_feed_inventory); + } if (initial_feed > 0) { - inventory.Push( - Material::Create( - this, initial_feed, context()->GetRecipe(feed_recipe))); + inventory.Push(Material::Create(this, initial_feed, + context()->GetRecipe(feed_recipe))); } - LOG(cyclus::LEV_DEBUG2, "EnrFac") << "CascadeEnrich " - << " entering the simuluation: "; + << " entering the simuluation: "; LOG(cyclus::LEV_DEBUG2, "EnrFac") << str(); } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -void CascadeEnrich::Tick() { - - current_swu_capacity = SwuCapacity(); - - } +void CascadeEnrich::Tick() {} // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - void CascadeEnrich::Tock() { using cyclus::toolkit::RecordTimeSeries; - LOG(cyclus::LEV_INFO4, "EnrFac") << prototype() << " used " - << intra_timestep_swu_ << " SWU"; - RecordTimeSeries(this, intra_timestep_swu_); LOG(cyclus::LEV_INFO4, "EnrFac") << prototype() << " used " << intra_timestep_feed_ << " feed"; RecordTimeSeries(this, intra_timestep_feed_); - } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -137,9 +125,10 @@ CascadeEnrich::GetMatlRequests() { return ports; } + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - // Sort offers of input material to have higher preference for more -// U-235 content +// U-235 content void CascadeEnrich::AdjustMatlPrefs( cyclus::PrefMap::type& prefs) { using cyclus::Bid; @@ -149,7 +138,6 @@ void CascadeEnrich::AdjustMatlPrefs( if (order_prefs == false) { return; } - cyclus::PrefMap::type::iterator reqit; // Loop over all requests @@ -168,7 +156,7 @@ void CascadeEnrich::AdjustMatlPrefs( bool u235_mass = 0; for (int bidit = 0; bidit < bids_vector.size(); bidit++) { - int new_pref = bidit + 1; + int new_pref = bidit + 10; // For any bids with U-235 qty=0, set pref to zero. if (!u235_mass) { @@ -188,15 +176,16 @@ void CascadeEnrich::AdjustMatlPrefs( // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - void CascadeEnrich::AcceptMatlTrades( const std::vector, - cyclus::Material::Ptr> >& responses) { + cyclus::Material::Ptr>>& responses) { // see // http://stackoverflow.com/questions/5181183/boostshared-ptr-and-inheritance std::vector, - cyclus::Material::Ptr> >::const_iterator it; + cyclus::Material::Ptr>>::const_iterator it; for (it = responses.begin(); it != responses.end(); ++it) { AddMat_(it->second); } } + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - void CascadeEnrich::AddMat_(cyclus::Material::Ptr mat) { // Elements and isotopes other than U-235, U-238 are sent directly to tails @@ -239,9 +228,11 @@ void CascadeEnrich::AddMat_(cyclus::Material::Ptr mat) { << " to its inventory, which is holding " << inventory.quantity() << " total."; } + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - std::set::Ptr> -CascadeEnrich::GetMatlBids(cyclus::CommodMap::type& out_requests) { +CascadeEnrich::GetMatlBids( + cyclus::CommodMap::type& out_requests) { using cyclus::Bid; using cyclus::BidPortfolio; using cyclus::CapacityConstraint; @@ -287,27 +278,34 @@ CascadeEnrich::GetMatlBids(cyclus::CommodMap::type& out_reques std::vector*>::iterator it; for (it = commod_requests.begin(); it != commod_requests.end(); ++it) { Request* req = *it; - Material::Ptr mat = req->target(); - double request_enrich = cyclus::toolkit::UraniumAssay(mat); - if (ValidReq(req->target()) && - ((request_enrich < max_enrich) || - (cyclus::AlmostEq(request_enrich, max_enrich)))) { - Material::Ptr offer = Offer_(req->target()); - commod_port->AddBid(req, offer, this); - } + Material::Ptr offer = Offer_(req->target()); + // The offer might not match the required enrichment ! It just produces + // what it can according to the cascade configuration and the feed assays + commod_port->AddBid(req, offer, this); } - Converter::Ptr sc(new SWUConverter(FeedAssay(), tails_assay)); - Converter::Ptr nc(new NatUConverter(FeedAssay(), tails_assay)); - CapacityConstraint swu(swu_capacity, sc); - CapacityConstraint natu(inventory.quantity(), nc); - commod_port->AddConstraint(swu); - commod_port->AddConstraint(natu); - - LOG(cyclus::LEV_INFO5, "EnrFac") - << prototype() << " adding a swu constraint of " << swu.capacity(); + // overbidding (bidding on every offer) + // add an overall production capacity constraint + + // correct the actual inventory quantity by the amount of Uranium in it... + double feed_qty = inventory.quantity(); + Material::Ptr natu_matl = inventory.Pop(feed_qty, cyclus::eps_rsrc()); + inventory.Push(natu_matl); + cyclus::toolkit::MatQuery mq(natu_matl); + std::set nucs; + nucs.insert(922350000); + nucs.insert(922380000); + double u_frac = mq.mass_frac(nucs); + double cor_feed_qty = feed_qty * u_frac; + double production_capacity = + ProductFlow(std::min(cor_feed_qty, MaxFeedFlow(FeedAssay(feed_qty)))); + cyclus::CapacityConstraint production_contraint( + production_capacity); + commod_port->AddConstraint(production_contraint); LOG(cyclus::LEV_INFO5, "EnrFac") - << prototype() << " adding a natu constraint of " << natu.capacity(); + << prototype() << " adding production capacity constraint of " + << production_capacity; + ports.insert(commod_port); } return ports; @@ -315,21 +313,18 @@ CascadeEnrich::GetMatlBids(cyclus::CommodMap::type& out_reques // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - void CascadeEnrich::GetMatlTrades( - const std::vector >& trades, + const std::vector>& trades, std::vector, - cyclus::Material::Ptr> >& responses) { + cyclus::Material::Ptr>>& responses) { using cyclus::Material; using cyclus::Trade; - intra_timestep_swu_ = 0; intra_timestep_feed_ = 0; - - std::vector >::const_iterator it; + std::vector>::const_iterator it; for (it = trades.begin(); it != trades.end(); ++it) { double qty = it->amt; std::string commod_type = it->bid->request()->commodity(); Material::Ptr response; - // Figure out whether material is tails or enriched, // if tails then make transfer of material if (commod_type == tails_commod) { @@ -352,28 +347,24 @@ void CascadeEnrich::GetMatlTrades( ss << "is being asked to provide more than its current inventory."; throw cyclus::ValueError(Agent::InformErrorMsg(ss.str())); } - if (cyclus::IsNegative(current_swu_capacity)) { - throw cyclus::ValueError("EnrFac " + prototype() + - " is being asked to provide more than" + - " its SWU capacity."); - } } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - cyclus::Material::Ptr CascadeEnrich::Enrich_(cyclus::Material::Ptr mat, - double qty) { + double qty) { using cyclus::Material; using cyclus::ResCast; using cyclus::toolkit::Assays; using cyclus::toolkit::UraniumAssay; - using cyclus::toolkit::SwuRequired; using cyclus::toolkit::FeedQty; using cyclus::toolkit::TailsQty; // get enrichment parameters - Assays assays(FeedAssay(), UraniumAssay(mat), tails_assay); - double swu_req = SwuRequired(qty, assays); - double natu_req = FeedQty(qty, assays); + double feed_qty = FeedRequired(qty); + double feed_assay = FeedAssay(feed_qty); + double product_assay = ProductAssay(feed_assay); + double tails_assay = TailsAssay(FeedAssay(feed_qty)); + double tails_mass = TailsFlow(feed_qty); // Determine the composition of the natural uranium // (ie. U-235+U-238/TotalMass) @@ -386,8 +377,7 @@ cyclus::Material::Ptr CascadeEnrich::Enrich_(cyclus::Material::Ptr mat, nucs.insert(922350000); nucs.insert(922380000); double natu_frac = mq.mass_frac(nucs); - double feed_req = natu_req / natu_frac; - + double feed_req = feed_qty / natu_frac; // pop amount from inventory and blob it into one material Material::Ptr r; try { @@ -398,12 +388,9 @@ cyclus::Material::Ptr CascadeEnrich::Enrich_(cyclus::Material::Ptr mat, r = inventory.Pop(feed_req, cyclus::eps_rsrc()); } } catch (cyclus::Error& e) { - NatUConverter nc(FeedAssay(), tails_assay); std::stringstream ss; ss << " tried to remove " << feed_req << " from its inventory of size " - << inventory.quantity() - << " and the conversion of the material into natu is " - << nc.convert(mat); + << inventory.quantity(); throw cyclus::ValueError(Agent::InformErrorMsg(ss.str())); } @@ -413,46 +400,35 @@ cyclus::Material::Ptr CascadeEnrich::Enrich_(cyclus::Material::Ptr mat, Material::Ptr response = r->ExtractComp(qty, comp); tails.Push(r); - current_swu_capacity -= swu_req; - - intra_timestep_swu_ += swu_req; - intra_timestep_feed_ += feed_req; - RecordEnrichment_(feed_req, swu_req); + RecordEnrichment_(feed_req); LOG(cyclus::LEV_INFO5, "EnrFac") << prototype() << " has performed an enrichment: "; LOG(cyclus::LEV_INFO5, "EnrFac") << " * Feed Qty: " << feed_req; - LOG(cyclus::LEV_INFO5, "EnrFac") << " * Feed Assay: " - << assays.Feed() * 100; + LOG(cyclus::LEV_INFO5, "EnrFac") << " * Feed Assay: " << feed_assay * 100; LOG(cyclus::LEV_INFO5, "EnrFac") << " * Product Qty: " << qty; LOG(cyclus::LEV_INFO5, "EnrFac") << " * Product Assay: " - << assays.Product() * 100; - LOG(cyclus::LEV_INFO5, "EnrFac") << " * Tails Qty: " - << TailsQty(qty, assays); - LOG(cyclus::LEV_INFO5, "EnrFac") << " * Tails Assay: " - << assays.Tails() * 100; - LOG(cyclus::LEV_INFO5, "EnrFac") << " * SWU: " << swu_req; - LOG(cyclus::LEV_INFO5, "EnrFac") << " * Current SWU capacity: " - << current_swu_capacity; + << product_assay * 100; + LOG(cyclus::LEV_INFO5, "EnrFac") << " * Tails Qty: " << tails_mass; + LOG(cyclus::LEV_INFO5, "EnrFac") << " * Tails Assay: " << tails_assay * 100; return response; } - // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -void CascadeEnrich::RecordEnrichment_(double natural_u, double swu) { + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +void CascadeEnrich::RecordEnrichment_(double natural_u) { using cyclus::Context; using cyclus::Agent; LOG(cyclus::LEV_DEBUG1, "EnrFac") << prototype() << " has enriched a material:"; LOG(cyclus::LEV_DEBUG1, "EnrFac") << " * Amount: " << natural_u; - LOG(cyclus::LEV_DEBUG1, "EnrFac") << " * SWU: " << swu; Context* ctx = Agent::context(); ctx->NewDatum("Enrichments") ->AddVal("ID", id()) ->AddVal("Time", ctx->time()) ->AddVal("Natural_Uranium", natural_u) - ->AddVal("SWU", swu) ->Record(); } @@ -462,15 +438,20 @@ cyclus::Material::Ptr CascadeEnrich::Request_() { return cyclus::Material::CreateUntracked(qty, context()->GetRecipe(feed_recipe)); } + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - cyclus::Material::Ptr CascadeEnrich::Offer_(cyclus::Material::Ptr mat) { - cyclus::toolkit::MatQuery q(mat); + double feed_assay = FeedAssay(mat->quantity()); + double product_assay = ProductAssay(feed_assay); + cyclus::CompMap comp; - comp[922350000] = q.atom_frac(922350000); - comp[922380000] = q.atom_frac(922380000); + comp[922350000] = product_assay; + comp[922380000] = 1 - product_assay; + return cyclus::Material::CreateUntracked( mat->quantity(), cyclus::Composition::CreateFromAtom(comp)); } + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - bool CascadeEnrich::ValidReq(const cyclus::Material::Ptr mat) { cyclus::toolkit::MatQuery q(mat); @@ -478,23 +459,79 @@ bool CascadeEnrich::ValidReq(const cyclus::Material::Ptr mat) { double u238 = q.atom_frac(922380000); return (u238 > 0 && u235 / (u235 + u238) > tails_assay); } + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - double CascadeEnrich::FeedAssay() { +double CascadeEnrich::FeedAssay(double quantity) { using cyclus::Material; if (inventory.empty()) { return 0; } - double pop_qty = inventory.quantity(); + double pop_qty = std::min(inventory.quantity(), quantity); cyclus::Material::Ptr fission_matl = inventory.Pop(pop_qty, cyclus::eps_rsrc()); inventory.Push(fission_matl); return cyclus::toolkit::UraniumAssay(fission_matl); } +double CascadeEnrich::ProductAssay(double feed_assay) { + CascadeConfig cascade_tmp = cascade.ModelMisuseCascade(feed_assay, misuse_model, precision); + return cascade_tmp.stgs_config.rbegin()->second.product_assay(); +} + +double CascadeEnrich::TailsAssay(double feed_assay) { + CascadeConfig cascade_tmp = cascade.ModelMisuseCascade(feed_assay, misuse_model, precision); + return cascade_tmp.stgs_config.begin()->second.tail_assay(); +} + +double CascadeEnrich::MaxFeedFlow(double feed_assay){ + CascadeConfig cascade_tmp = cascade.ModelMisuseCascade(feed_assay, misuse_model, precision); + + return FlowPerDt(cascade_tmp.FeedFlow()); + +} + +double CascadeEnrich::FeedRequired(double prod_qty) { + double max_feed_flow = MaxFeedFlow(FeedAssay(inventory.quantity())); + double max_product_flow = ProductFlow(max_feed_flow); + double feed_required = max_feed_flow / max_product_flow * prod_qty; + + max_feed_flow = MaxFeedFlow(FeedAssay(feed_required)); + max_product_flow = ProductFlow(max_feed_flow); + double corrected_feed_required = max_feed_flow / max_product_flow * prod_qty; + double diff_feed = std::abs(feed_required - corrected_feed_required); + + while (diff_feed > precision) { + // reset feed_required + feed_required = corrected_feed_required; + + max_feed_flow = MaxFeedFlow(FeedAssay(feed_required)); + max_product_flow = ProductFlow(max_feed_flow); + corrected_feed_required = max_feed_flow / max_product_flow * prod_qty; + diff_feed = std::abs(feed_required - corrected_feed_required); + } + + return corrected_feed_required; +} + +double CascadeEnrich::ProductFlow(double feed_flow) { + double feed_assay = FeedAssay(feed_flow); + double feed_ratio = feed_flow / MaxFeedFlow(feed_assay); + CascadeConfig cascade_tmp = cascade.ModelMisuseCascade(feed_assay, misuse_model, precision); + + StageConfig last_stg = cascade_tmp.stgs_config.rbegin()->second; + double product_flow = last_stg.feed_flow() * last_stg.cut(); + return feed_ratio * FlowPerDt(product_flow); +} + +double CascadeEnrich::TailsFlow(double feed_flow) { + // this assume mass flow conservation + return feed_flow - ProductFlow(feed_flow); +} + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - extern "C" cyclus::Agent* ConstructCascadeEnrich(cyclus::Context* ctx) { return new CascadeEnrich(ctx); } - + } // namespace mbmore diff --git a/src/CascadeEnrich.h b/src/CascadeEnrich.h index eceb46a..5549dbb 100644 --- a/src/CascadeEnrich.h +++ b/src/CascadeEnrich.h @@ -5,6 +5,8 @@ #include "cyclus.h" #include "sim_init.h" +#include "CascadeConfig.h" +#include "CentrifugeConfig.h" /* Working with cycamore Develop build: 3ada148442de636d @@ -37,138 +39,69 @@ B) Change cascade feed assay */ namespace mbmore { -/// @class SWUConverter -/// -/// @brief The SWUConverter is a simple Converter class for material to -/// determine the amount of SWU required for their proposed enrichment -class SWUConverter : public cyclus::Converter { - public: - SWUConverter(double feed_commod, double tails) - : feed_(feed_commod), tails_(tails) {} - virtual ~SWUConverter() {} - - /// @brief provides a conversion for the SWU required - virtual double convert( - cyclus::Material::Ptr m, cyclus::Arc const* a = NULL, - cyclus::ExchangeTranslationContext const* ctx = - NULL) const { - cyclus::toolkit::Assays assays(feed_, cyclus::toolkit::UraniumAssay(m), - tails_); - return cyclus::toolkit::SwuRequired(m->quantity(), assays); - } - - /// @returns true if Converter is a SWUConverter and feed and tails equal - virtual bool operator==(Converter& other) const { - SWUConverter* cast = dynamic_cast(&other); - return cast != NULL && feed_ == cast->feed_ && tails_ == cast->tails_; - } - - private: - double feed_, tails_; -}; - - -/// @class NatUConverter -/// -/// @brief The NatUConverter is a simple Converter class for material to -/// determine the amount of natural uranium required for their proposed -/// enrichment -class NatUConverter : public cyclus::Converter { - public: - NatUConverter(double feed_commod, double tails) - : feed_(feed_commod), tails_(tails) {} - virtual ~NatUConverter() {} - - // virtual std::string version() { return CYCAMORE_VERSION; } - - /// @brief provides a conversion for the amount of natural Uranium required - virtual double convert( - cyclus::Material::Ptr m, cyclus::Arc const* a = NULL, - cyclus::ExchangeTranslationContext const* ctx = - NULL) const { - cyclus::toolkit::Assays assays(feed_, cyclus::toolkit::UraniumAssay(m), - tails_); - cyclus::toolkit::MatQuery mq(m); - std::set nucs; - nucs.insert(922350000); - nucs.insert(922380000); - - double natu_frac = mq.mass_frac(nucs); - double natu_req = cyclus::toolkit::FeedQty(m->quantity(), assays); - return natu_req / natu_frac; - } - - /// @returns true if Converter is a NatUConverter and feed and tails equal - virtual bool operator==(Converter& other) const { - NatUConverter* cast = dynamic_cast(&other); - return cast != NULL && feed_ == cast->feed_ && tails_ == cast->tails_; - } - - private: - double feed_, tails_; -}; - class CascadeEnrich : public cyclus::Facility { #pragma cyclus note { \ "niche": "enrichment facility", \ "doc": \ "The CascadeEnrich facility based on the Cycamore Enrich facility. " \ - "timesteps (see README for full implementation documentation ",\ + "timesteps (see README for full implementation documentation ",\ } public: // --- Module Members --- - /// Constructor for the CascadeEnrich class - /// @param ctx the cyclus context for access to simulation-wide parameters + /// Constructor for the CascadeEnrich class + /// @param ctx the cyclus context for access to simulation-wide parameters CascadeEnrich(cyclus::Context* ctx); - /// Destructor for the CascadeEnrich class + /// Destructor for the CascadeEnrich class virtual ~CascadeEnrich(); #pragma cyclus - /// Print information about this agent + /// Print information about this agent virtual std::string str(); // --- - // --- Facility Members --- - /// perform module-specific tasks when entering the simulation - virtual void Build(cyclus::Agent* parent); - // --- + /// --- Facility Members --- + /// Perform module-specific tasks when entering the simulation + /// virtual void Build(cyclus::Agent* parent); + /// --- // --- Agent Members --- - /// Each facility is prompted to do its beginning-of-time-step - /// stuff at the tick of the timer. - - /// @param time is the time to perform the tick + /// Each facility is prompted to do its beginning-of-time-step + /// stuff at the tick of the timer. + virtual void EnterNotify(); + /// @param time is the time to perform the tick virtual void Tick(); - /// Each facility is prompted to its end-of-time-step - /// stuff on the tock of the timer. + /// Each facility is prompted to its end-of-time-step + /// stuff on the tock of the timer. - /// @param time is the time to perform the tock + /// @param time is the time to perform the tock virtual void Tock(); - /// @brief The Enrichment request Materials of its given - /// commodity. + /// @brief The Enrichment request Materials of its given commodity. virtual std::set::Ptr> GetMatlRequests(); /// @brief The Enrichment adjusts preferences for offers of - /// natural uranium it has received to maximize U-235 content - /// Any offers that have zero U-235 content are not accepted + /// natural uranium it has received to maximize U-235 content. + /// Any offers that have zero U-235 content are not accepted. + /// + /// @param prefs is the preference map with all requests. virtual void AdjustMatlPrefs(cyclus::PrefMap::type& prefs); - /// @brief The Enrichment place accepted trade Materials in their - /// Inventory + /// @brief The Enrichment place accepted trade Materials in their Inventory. virtual void AcceptMatlTrades( const std::vector, cyclus::Material::Ptr> >& responses); /// @brief Responds to each request for this facility's commodity. If a given /// request is more than this facility's inventory or SWU capacity, it will - /// offer its minimum of its capacities. + /// offer the minimum of its capacities. + /// + /// @param commod_requests is a cyclus map of requests for this facility virtual std::set::Ptr> GetMatlBids( cyclus::CommodMap::type& commod_requests); @@ -188,33 +121,28 @@ class CascadeEnrich : public cyclus::Facility { inventory.capacity(size); } - inline void SwuCapacity(double capacity) { - swu_capacity = capacity; - current_swu_capacity = swu_capacity; - } - - inline double SwuCapacity() const { return swu_capacity; } - // TODO: MAKE THESE CONVERSIONS TOOLKIT FNS and have them explicitly check // timestep duration - // Convert input file flows kg/mon to SI units + /// --- Convert input file flows kg/mon to SI units --- inline double FlowPerSec(double flow_per_mon) { - return flow_per_mon / secpermonth; + return flow_per_mon / secpertimestep; } - inline double FlowPerMon(double flow_per_sec) { - return flow_per_sec * secpermonth; + inline double FlowPerDt(double flow_per_sec) { + return flow_per_sec * secpertimestep; } inline double Mg2kgPerSec(double feed_mg_per_sec) { - return feed_mg_per_sec / (1e6); + return feed_mg_per_sec / (1e6); } - - /// @brief Determines if a particular material is a valid request to respond - /// to. Valid requests must contain U235 and U238 and must have a relative - /// U235-to-U238 ratio less than this facility's tails_assay(). - /// @return true if the above description is met by the material + /// --- + + /// @brief Determines if a particular material is a valid request to respond + /// to. Valid requests must contain U235 and U238 and must have a relative + /// U235-to-U238 ratio less than this facility's tails_assay(). + /// + /// @return true if the above description is met by the material bool ValidReq(const cyclus::Material::Ptr mat); inline const cyclus::toolkit::ResBuf& Tails() const { @@ -222,142 +150,185 @@ class CascadeEnrich : public cyclus::Facility { } private: - /// @brief calculates the feed assay based on the unenriched inventory - double FeedAssay(); + /// @brief calculates the feed assay based on the unenriched inventory + double FeedAssay(double quantity); - /// @brief adds a material into the natural uranium inventory - /// @throws if the material is not the same composition as the feed_recipe + /// @brief adds a material into the natural uranium inventory + /// @throws if the material is not the same composition as the feed_recipe void AddMat_(cyclus::Material::Ptr mat); - /// @brief generates a request for this facility given its current state. - /// Quantity of the material will be equal to remaining inventory size. + /// @brief generates a request for this facility given its current state. + /// Quantity of the material will be equal to remaining inventory size. cyclus::Material::Ptr Request_(); - /// @brief Generates a material offer for a given request. The response - /// composition will be comprised only of U235 and U238 at their relative - /// ratio in the requested material. The response quantity will be the - /// same as the requested commodity. + /// @brief Generates a material offer for a given request. The response + /// composition will be comprised only of U235 and U238 at their relative + /// ratio in the requested material. The response quantity will be the + /// same as the requested commodity. /// - /// @param req the requested material being responded to + /// @param req the requested material being responded to cyclus::Material::Ptr Offer_(cyclus::Material::Ptr req); cyclus::Material::Ptr Enrich_(cyclus::Material::Ptr mat, double qty); - /// @brief records and enrichment with the cyclus::Recorder - void RecordEnrichment_(double natural_u, double swu); + /// @brief records and enrichment with the cyclus::Recorder + void RecordEnrichment_(double natural_u); + + + // Not physical constants but fixed assumptions for a cascade separating out + // U235 from U238 in UF6 gas group all the characteristic of a centrifuges. + CentrifugeConfig centrifuge; + CascadeConfig cascade; + double precision = 1e-15; + + + double secpertimestep = 60.*60.*24.*(365.25/12.); // Set to design_tails at beginning of simulation. Gets reset if // facility is used off-design - double tails_assay; + double tails_assay; - // These state variables are constrained by the design input params at - // the start of the simulation: - // Set by max feed for an individual machine - double design_delU; - double design_alpha; - // Set by design assays (feed, product, tails) int n_enrich_stages; int n_strip_stages; // Set by maximum allowable centrifuges - double max_feed_inventory; - double swu_capacity; - - // Not physical constants but fixed assumptions for a cascade separating - // out U235 from U238 in UF6 gas - const double M = 0.352; // kg/mol UF6 - const double dM = 0.003; // kg/mol U238 - U235 - const double x = 1000; // Pressure ratio (Glaser) + double ProductAssay(double feed_assay); + double ProductFlow(double feed_flow); + double TailsAssay(double feed_assay); + double TailsFlow(double feed_flow); + double FeedRequired(double prod_qty); + double MaxFeedFlow(double feed_assay); - const double flow_internal = 2.0; // can vary from 2-4 - const double eff = 1.0; // typical efficiencies <0.6 - const double cut = 0.5; // target for ideal cascade - - const double secpermonth = 60*60*24*(365.25/12); - - private: #pragma cyclus var { \ "tooltip" : "feed recipe", \ "doc" : "recipe for enrichment facility feed commodity", \ - "uilabel" : "Feed Recipe", "uitype" : "recipe" } + "uilabel" : "Feed Recipe", \ + "uitype" : "recipe" } std::string feed_recipe; #pragma cyclus var { \ - "default": 0, "tooltip": "initial uranium reserves (kg)", \ + "default": 0, \ + "tooltip": "initial uranium reserves (kg)", \ "uilabel": "Initial Feed Inventory", \ "doc": "amount of natural uranium stored at the enrichment " \ "facility at the beginning of the simulation (kg)" } double initial_feed; + #pragma cyclus var { \ + "default": 1e299, "tooltip": "max inventory of feed material (kg)", \ + "uilabel": "Maximum Feed Inventory", \ + "uitype": "range", \ + "range": [0.0, 1e299], \ + "doc": "maximum total inventory of natural uranium in " \ + "the enrichment facility (kg)" \ + } + double max_feed_inventory; + #pragma cyclus var { \ - "default": 0, "tooltip": "design feed flow (kg/mon)", \ + "default": 0.000015, \ + "tooltip": "design feed flow (kg/s)", \ "uilabel": "Design Feed Flow", \ - "doc": "Target amount of feed material to be processed by the" \ - " facility (kg/mon). Either this or max_centrifuges is used to constrain" \ - " the cascade design" } + "doc": "Target amount of feed material to be processed by the " \ + "facility (kg/s). Either this or max_centrifuges is used to constrain " \ + "the cascade design" } double design_feed_flow; #pragma cyclus var { \ - "default" : 0, "tooltip" : "number of centrifuges available ", \ + "default" : -1, \ + "tooltip" : "number of centrifuges available ", \ "uilabel" : "Number of Centrifuges", \ "doc" : "number of centrifuges available to make the cascade" } int max_centrifuges; // TODO: USE FEED RECIPE TO DETERMINE FEED ASSAY!!! #pragma cyclus var { \ - "default": 0.0071, "tooltip": "initial uranium reserves (kg)", \ + "default": 0.0071, \ + "tooltip": "initial uranium reserves (kg)", \ "uilabel": "Initial feed assay", \ "doc": "desired fraction of U235 in feed material (should be consistent "\ "with feed recipe" } double design_feed_assay; #pragma cyclus var { \ - "default" : 0.035, "tooltip" : "Initial target product assay", \ + "default" : 0.035, \ + "tooltip" : "Initial target product assay", \ "uilabel" : "Target product assay", \ "doc" : "desired fraction of U235 in product" } double design_product_assay; #pragma cyclus var { \ - "default" : 0.003, "tooltip" : "Initial target tails assay", \ + "default" : 0.003, \ + "tooltip" : "Initial target tails assay", \ "uilabel" : "Target tails assay", \ "doc" : "desired fraction of U235 in tails" } double design_tails_assay; #pragma cyclus var { \ - "default" : 320.0, "tooltip" : "Centrifuge temperature (Kelvin)", \ + "default" : 320.0, \ + "tooltip" : "Centrifuge temperature (Kelvin)", \ "uilabel" : "Centrifuge temperature (Kelvin)", \ "doc" : "temperature at which centrifuges are operated (Kelvin)" } double temp; #pragma cyclus var { \ - "default" : 485.0, "tooltip" : "Centrifuge velocity (m/s)", \ + "default" : 485.0, \ + "tooltip" : "Centrifuge velocity (m/s)", \ "uilabel" : "Centrifuge velocity (m/s)", \ - "doc" : "operational centrifuge velocity (m/s) at the outer radius (a)"} + "doc" : "operational centrifuge velocity (m/s) at the outer radius (a)" } double centrifuge_velocity; #pragma cyclus var { \ - "default" : 0.5, "tooltip" : "Centrifuge height (m)", \ + "default" : 0.5, \ + "tooltip" : "Centrifuge height (m)", \ "uilabel" : "Centrifuge height (m)", \ - "doc" : "height of centrifuge (m)"} + "doc" : "height of centrifuge (m)" } double height; #pragma cyclus var { \ - "default" : 0.15, "tooltip" : "Centrifuge diameter (m)", \ + "default" : 0.15, \ + "tooltip" : "Centrifuge diameter (m)", \ "uilabel" : "Centrifuge diameter (m)", \ - "doc" : "diameter of centrifuge (m)"} + "doc" : "diameter of centrifuge (m)" } double diameter; #pragma cyclus var { \ - "default" : 15.0, "tooltip" : "Centrifuge feed rate (mg/sec)", \ + "default" : 2, \ + "tooltip" : "Centrifuge L/F* ", \ + "uilabel" : "Centrifuge Countercurrent to feed ratio", \ + "doc" : "Countercurrent to feed ratio" } + double L_over_F; + +#pragma cyclus var { \ + "default" : 15.0, \ + "tooltip" : "Centrifuge feed rate (mg/sec)", \ "uilabel" : "Max feed rate for single centrifuge (mg/sec)", \ - "doc" : "maximum feed rate for a single centrifuge (mg/sec)"} + "doc" : "maximum feed rate for a single centrifuge (mg/sec)" } double machine_feed; +#pragma cyclus var { \ + "default" : 1000., \ + "tooltip" : "Internal Pressure Ratio", \ + "uilabel" : "Centrifuge internal pressure ratio", \ + "doc" : "Pressure ratio parameter for centrifuge design"\ + "that drives the r1 / r2 ratio."} + double x; + + #pragma cyclus var { \ + "default": 0, \ + "userlevel": 10, \ + "tooltip": "Modeling option for misuse calculation", \ + "uilabel": "Misuse Modeling", \ + "doc": "Misuse modeling option: " \ + " - 0: alpha-theta fix -- beta varies, " \ + " - 1: alpha=beta fix -- theta varies, " \ + " - 2: alpha*beta fix " } + int misuse_model; + + - // Input params from cycamore::Enrichment #pragma cyclus var { \ "default": 1, \ @@ -386,34 +357,34 @@ class CascadeEnrich : public cyclus::Facility { #pragma cyclus var { \ "tooltip" : "feed commodity", \ "doc" : "feed commodity that the enrichment facility accepts", \ - "uilabel" : "Feed Commodity", "uitype" : "incommodity" } + "uilabel" : "Feed Commodity", \ + "uitype" : "incommodity" } std::string feed_commod; #pragma cyclus var { \ "tooltip" : "product commodity", \ "doc" : "product commodity that the enrichment facility generates", \ - "uilabel" : "Product Commodity", "uitype" : "outcommodity" } + "uilabel" : "Product Commodity", \ + "uitype" : "outcommodity" } std::string product_commod; #pragma cyclus var { \ "tooltip" : "tails commodity", \ "doc" : "tails commodity supplied by enrichment facility", \ - "uilabel" : "Tails Commodity", "uitype" : "outcommodity" } + "uilabel" : "Tails Commodity", \ + "uitype" : "outcommodity" } std::string tails_commod; - double current_swu_capacity; #pragma cyclus var {} cyclus::toolkit::ResBuf tails; // depleted u // used to total intra-timestep swu and natu usage for meeting requests - // these help enable time series generation. - double intra_timestep_swu_; double intra_timestep_feed_; // END LEGACY -#pragma cyclus var { 'capacity' : 'max_feed_inventory' } cyclus::toolkit::ResBuf inventory; // natural u friend class CascadeEnrichTest; diff --git a/src/CascadeEnrich_tests.cc b/src/CascadeEnrich_tests.cc index cbcf745..41cfab5 100644 --- a/src/CascadeEnrich_tests.cc +++ b/src/CascadeEnrich_tests.cc @@ -20,9 +20,9 @@ using cyclus::Material; namespace mbmore { -namespace cascadenrichtest{ +namespace cascadenrichtest { - Composition::Ptr c_nou235() { +Composition::Ptr c_nou235() { cyclus::CompMap m; m[922380000] = 1.0; return Composition::CreateFromMass(m); @@ -62,10 +62,10 @@ TEST_F(CascadeEnrichTest, RequestQty) { " natu1 " " enr_u " " tails " - " 1.0 " - " 0.003 "; + " 0.003 " + " 1000 "; - int simdur = 10; + int simdur = 1; cyclus::MockSim sim(cyclus::AgentSpec(":mbmore:CascadeEnrich"), config, simdur); sim.AddRecipe("natu1", cascadenrichtest::c_natu1()); @@ -80,14 +80,14 @@ TEST_F(CascadeEnrichTest, RequestQty) { Material::Ptr m = sim.GetMaterial(qr.GetVal("ResourceId")); // Should be only one transaction into the EF, - // and it should be exactly 1kg of natu - EXPECT_EQ(1.0, qr.rows.size()); - EXPECT_NEAR(1.0, m->quantity(), 1e-10) + // and it should be exactly 207.8928179411368kg of natu + EXPECT_EQ(1, qr.rows.size()); + EXPECT_NEAR(1000, m->quantity(), 1e-3) << "matched trade provides the wrong quantity of material"; } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -TEST_F(CascadeEnrichTest, CheckSWUConstraint) { +TEST_F(CascadeEnrichTest, CheckFlowConstraint) { // Tests that request for enrichment that exceeds the SWU constraint // fulfilled only up to the available SWU. // Also confirms that initial_feed flag works. @@ -99,8 +99,7 @@ TEST_F(CascadeEnrichTest, CheckSWUConstraint) { " enr_u " " tails " " 0.003 " - " 1000 " - " 195 "; + " 1000 "; int simdur = 1; @@ -110,7 +109,7 @@ TEST_F(CascadeEnrichTest, CheckSWUConstraint) { sim.AddRecipe("natu1", cascadenrichtest::c_natu1()); sim.AddRecipe("heu", cascadenrichtest::c_heu()); - sim.AddSink("enr_u").recipe("heu").capacity(10).Finalize(); + sim.AddSink("enr_u").recipe("heu").capacity(4).Finalize(); int id = sim.Run(); @@ -120,8 +119,8 @@ TEST_F(CascadeEnrichTest, CheckSWUConstraint) { Material::Ptr m = sim.GetMaterial(qr.GetVal("ResourceId")); EXPECT_EQ(1.0, qr.rows.size()); - EXPECT_NEAR(5.0, m->quantity(), 0.1) - << "traded quantity exceeds SWU constraint"; + EXPECT_NEAR(4, m->quantity(), 0.1) + << "traded quantity differ from flow contraints"; } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -135,7 +134,8 @@ TEST_F(CascadeEnrichTest, CheckCapConstraint) { " enr_u " " tails " " 0.003 " - " 243 "; + " 200 " + " 50 "; int simdur = 1; @@ -145,7 +145,7 @@ TEST_F(CascadeEnrichTest, CheckCapConstraint) { sim.AddRecipe("natu1", cascadenrichtest::c_natu1()); sim.AddRecipe("heu", cascadenrichtest::c_heu()); - sim.AddSink("enr_u").recipe("heu").capacity(10).Finalize(); + sim.AddSink("enr_u").recipe("heu").capacity(100).Finalize(); int id = sim.Run(); @@ -155,14 +155,14 @@ TEST_F(CascadeEnrichTest, CheckCapConstraint) { Material::Ptr m = sim.GetMaterial(qr.GetVal("ResourceId")); EXPECT_EQ(1.0, qr.rows.size()); - EXPECT_NEAR(5.0, m->quantity(), 0.01) + EXPECT_NEAR(24.691, m->quantity(), 0.01) << "traded quantity exceeds capacity constraint"; } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -TEST_F(CascadeEnrichTest, RequestEnrich) { - // this tests verifies that requests for output material exceeding - // the maximum allowed enrichment are not fulfilled. +TEST_F(CascadeEnrichTest, FeedFlowTime) { + // Tests that simulations using different timesteps reach the + // same equilibrium value, all other parameters held the same. std::string config = " natu " @@ -170,17 +170,22 @@ TEST_F(CascadeEnrichTest, RequestEnrich) { " enr_u " " tails " " 0.003 " - " 0.20 "; + " 1e09 " + " 300"; + + int simdur = 1; - int simdur = 2; + //-------------------------------------------------------------------- + // BEGIN: Daily timestep sim cyclus::MockSim sim(cyclus::AgentSpec(":mbmore:CascadeEnrich"), config, - simdur); + simdur); + cyclus::SimInfo SI(1, 0, 1, "", "never"); + SI.dt = 3600*24; + cyclus::Context* ctx = sim.context(); + ctx->InitSim(SI); + sim.AddRecipe("natu1", cascadenrichtest::c_natu1()); - sim.AddRecipe("leu", cascadenrichtest::c_leu()); sim.AddRecipe("heu", cascadenrichtest::c_heu()); - - sim.AddSource("natu").recipe("natu1").Finalize(); - sim.AddSink("enr_u").recipe("leu").capacity(1.0).Finalize(); sim.AddSink("enr_u").recipe("heu").Finalize(); int id = sim.Run(); @@ -190,22 +195,34 @@ TEST_F(CascadeEnrichTest, RequestEnrich) { QueryResult qr = sim.db().Query("Transactions", &conds); Material::Ptr m = sim.GetMaterial(qr.GetVal("ResourceId")); - // Should be only one transaction out of the EF, - // and it should be 1kg of LEU - EXPECT_EQ(1.0, qr.rows.size()); - EXPECT_NEAR(1.0, m->quantity(), 0.01) - << "Not providing the requested quantity"; + double day_q = m->quantity(); - CompMap got = m->comp()->mass(); - CompMap want = cascadenrichtest::c_leu()->mass(); - cyclus::compmath::Normalize(&got); - cyclus::compmath::Normalize(&want); + // END: Daily timestep sim + //-------------------------------------------------------------------- + // BEGIN: Yearly timestep sim + cyclus::MockSim sim2(cyclus::AgentSpec(":mbmore:CascadeEnrich"), config, + simdur); + cyclus::SimInfo SI2(1, 0, 1, "", "never"); + SI2.dt = 3600*24*365.25; + cyclus::Context* ctx2 = sim2.context(); + ctx2->InitSim(SI2); - CompMap::iterator it; - for (it = want.begin(); it != want.end(); ++it) { - EXPECT_DOUBLE_EQ(it->second, got[it->first]) - << "nuclide qty off: " << pyne::nucname::name(it->first); - } + sim2.AddRecipe("natu1", cascadenrichtest::c_natu1()); + sim2.AddRecipe("heu", cascadenrichtest::c_heu()); + sim2.AddSink("enr_u").recipe("heu").Finalize(); + + int id2 = sim2.Run(); + + std::vector conds2; + conds2.push_back(Cond("Commodity", "==", std::string("enr_u"))); + QueryResult qr2 = sim2.db().Query("Transactions", &conds2); + Material::Ptr m2 = sim2.GetMaterial(qr2.GetVal("ResourceId")); + + double year_q = m2->quantity(); + // END: Yearly timestep sim + + // Check if, scaling for timestep, the same quantity is produced + EXPECT_NEAR(day_q*365.25,year_q, 0.01); } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -226,7 +243,7 @@ TEST_F(CascadeEnrichTest, TradeTails) { sim.AddRecipe("natu1", cascadenrichtest::c_natu1()); sim.AddRecipe("leu", cascadenrichtest::c_leu()); - sim.AddSource("natu").recipe("natu1").Finalize(); + sim.AddSource("natu").recipe("natu1").capacity(200).Finalize(); sim.AddSink("enr_u").recipe("leu").Finalize(); sim.AddSink("tails").Finalize(); @@ -258,7 +275,7 @@ TEST_F(CascadeEnrichTest, TailsQty) { sim.AddRecipe("natu1", cascadenrichtest::c_natu1()); sim.AddRecipe("leu", cascadenrichtest::c_leu()); - sim.AddSource("natu").recipe("natu1").Finalize(); + sim.AddSource("natu").recipe("natu1").capacity(200).Finalize(); sim.AddSink("enr_u").recipe("leu").capacity(0.5).Finalize(); sim.AddSink("enr_u").recipe("leu").capacity(0.5).Finalize(); sim.AddSink("tails").Finalize(); @@ -270,7 +287,8 @@ TEST_F(CascadeEnrichTest, TailsQty) { QueryResult qr = sim.db().Query("Transactions", &conds); Material::Ptr m = sim.GetMaterial(qr.GetVal("ResourceId")); - // Should be 2 tails transactions, one from each LEU sink, each 4.084kg. + // Should be 2 tails transactions, one from each LEU sink, each 3.787453kg. + // (1* (0.0403193-0.00271456) / (0.0071-0.00271456) -1)/2. = 3.787453 EXPECT_EQ(2, qr.rows.size()); cyclus::SqlStatement::Ptr stmt = sim.db().db().Prepare( @@ -280,7 +298,7 @@ TEST_F(CascadeEnrichTest, TailsQty) { stmt->BindText(1, "tails"); stmt->Step(); - EXPECT_NEAR(8.168, stmt->GetDouble(0), 0.01) + EXPECT_NEAR(7.100, stmt->GetDouble(0), 0.01) << "Not providing the requested quantity"; } @@ -291,7 +309,7 @@ TEST_F(CascadeEnrichTest, BidPrefs) { std::string config = " natu " - " natu1 " + " natu2 " " enr_u " " tails " " 0.003 " @@ -303,9 +321,10 @@ TEST_F(CascadeEnrichTest, BidPrefs) { sim.AddRecipe("natu1", cascadenrichtest::c_natu1()); sim.AddRecipe("natu2", cascadenrichtest::c_natu2()); + sim.AddSource("natu").recipe("natu1").capacity(200).Finalize(); sim.AddSource("natu").recipe("natu1").capacity(1).Finalize(); - sim.AddSource("natu").recipe("natu2").capacity(1).Finalize(); + sim.AddSource("natu").recipe("natu2").capacity(2).Finalize(); int id = sim.Run(); @@ -422,7 +441,7 @@ void CascadeEnrichTest::InitParameters() { tails_assay = 0.002; swu_capacity = 100; //** - inv_size = 5; + inv_size = 105.5; reserves = 105.5; } @@ -436,7 +455,6 @@ void CascadeEnrichTest::SetUpSource() { src_facility->tails_assay = tails_assay; src_facility->max_enrich = max_enrich; src_facility->SetMaxInventorySize(inv_size); - src_facility->SwuCapacity(swu_capacity); src_facility->initial_feed = reserves; } @@ -521,158 +539,6 @@ TEST_F(CascadeEnrichTest, ValidReq) { } // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -TEST_F(CascadeEnrichTest, ConstraintConverters) { - // Tests the SWU and NatU converters to make sure that amount of - // feed and SWU required are correct to fulfill the enrichment request. - using cyclus::CompMap; - using cyclus::Material; - using cyclus::toolkit::MatQuery; - using cyclus::Composition; - using cyclus::toolkit::Assays; - using cyclus::toolkit::UraniumAssay; - using cyclus::toolkit::SwuRequired; - using cyclus::toolkit::FeedQty; - using cyclus::toolkit::MatQuery; - using mbmore::SWUConverter; - cyclus::Env::SetNucDataPath(); - - double qty = 5; // 5 kg - double product_assay = 0.05; // of 5 w/o enriched U - CompMap v; - v[922350000] = product_assay; - v[922380000] = 1 - product_assay; - v[94239] = 0.5; // 94239 shouldn't be taken into account - Material::Ptr target = - Material::CreateUntracked(qty, Composition::CreateFromMass(v)); - - std::set nucs; - nucs.insert(922350000); - nucs.insert(922380000); - - MatQuery mq(target); - double mass_frac = mq.mass_frac(nucs); - - SWUConverter swuc(feed_assay, tails_assay); - NatUConverter natuc(feed_assay, tails_assay); - - Material::Ptr offer = DoOffer(target); - - EXPECT_NEAR(swuc.convert(target), swuc.convert(offer), 0.001); - EXPECT_NEAR(natuc.convert(target) * mass_frac, natuc.convert(offer), 0.001); -} - -// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -TEST_F(CascadeEnrichTest, Enrich) { - // this test asks the facility to enrich a material that results in an amount - // of natural uranium required that is exactly its inventory level. that - // inventory will be comprised of two materials to test the manifest/absorb - // strategy employed in Enrich_. - using cyclus::CompMap; - using cyclus::Material; - using cyclus::toolkit::MatQuery; - using cyclus::Composition; - using cyclus::toolkit::Assays; - using cyclus::toolkit::UraniumAssay; - using cyclus::toolkit::SwuRequired; - using cyclus::toolkit::FeedQty; - - double qty = 5; // kg - double product_assay = 0.05; // of 5 w/o enriched U - cyclus::CompMap v; - v[922350000] = product_assay; - v[922380000] = 1 - product_assay; - // target qty need not be = to request qty - Material::Ptr target = cyclus::Material::CreateUntracked( - qty + 10, cyclus::Composition::CreateFromMass(v)); - - Assays assays(feed_assay, UraniumAssay(target), tails_assay); - double swu_req = SwuRequired(qty, assays); - double natu_req = FeedQty(qty, assays); - double tails_qty = TailsQty(qty, assays); - - double swu_cap = swu_req * 5; - src_facility->SwuCapacity(swu_cap); - src_facility->SetMaxInventorySize(natu_req); - DoAddMat(GetMat(natu_req / 2)); - DoAddMat(GetMat(natu_req / 2)); - - Material::Ptr response; - EXPECT_NO_THROW(response = DoEnrich(target, qty)); - EXPECT_DOUBLE_EQ(src_facility->Tails().quantity(), tails_qty); - - MatQuery q(response); - EXPECT_EQ(response->quantity(), qty); - EXPECT_EQ(q.mass_frac(922350000), product_assay); - EXPECT_EQ(q.mass_frac(922380000), 1 - product_assay); - - // test too much natu request - DoAddMat(GetMat(natu_req - 1)); - EXPECT_THROW(response = DoEnrich(target, qty), cyclus::Error); -} - -// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -TEST_F(CascadeEnrichTest, Response) { - // this test asks the facility to respond to multiple requests for enriched - // uranium. two requests are provided, whose total equals the swu capacity of - // the facility while not exceeding its inventory capacity (that's taken care - // of in the Enrich tests). - // - // note that response quantity and quality need not be tested, because they - // are covered by the Enrich and RequestEnrich tests - using cyclus::Bid; - using cyclus::CompMap; - using cyclus::Material; - using cyclus::Request; - using cyclus::Trade; - using cyclus::toolkit::Assays; - using cyclus::toolkit::FeedQty; - using cyclus::toolkit::SwuRequired; - using cyclus::toolkit::UraniumAssay; - - // problem set up - std::vector > trades; - std::vector< - std::pair, cyclus::Material::Ptr> > - responses; - - double qty = 5; // kg - double trade_qty = qty / 3; - double product_assay = 0.05; // of 5 w/o enriched U - - cyclus::CompMap v; - v[922350000] = product_assay; - v[922380000] = 1 - product_assay; - // target qty need not be = to request qty - Material::Ptr target = cyclus::Material::CreateUntracked( - qty + 10, cyclus::Composition::CreateFromMass(v)); - - Assays assays(feed_assay, UraniumAssay(target), tails_assay); - double swu_req = SwuRequired(qty, assays); - double natu_req = FeedQty(qty, assays); - - src_facility->SetMaxInventorySize(natu_req * 4); // not capacitated by nat - src_facility->SwuCapacity(swu_req); // swu capacitated - - src_facility->GetMatlTrades(trades, responses); - - // set up state - DoAddMat(GetMat(natu_req * 2)); - - src_facility->GetMatlTrades(trades, responses); - - Request* req = - Request::Create(target, trader, product_commod); - Bid* bid = Bid::Create(req, target, src_facility); - Trade trade(req, bid, trade_qty); - trades.push_back(trade); - - // 2 trades, SWU = SWU cap - ASSERT_GT(src_facility->SwuCapacity() - 2 * swu_req / 3, -1 * cyclus::eps()); - trades.push_back(trade); - responses.clear(); - EXPECT_NO_THROW(src_facility->GetMatlTrades(trades, responses)); - EXPECT_EQ(responses.size(), 2); -} } // namespace cycamore diff --git a/src/CascadeEnrich_tests.h b/src/CascadeEnrich_tests.h index 22a650d..025911b 100644 --- a/src/CascadeEnrich_tests.h +++ b/src/CascadeEnrich_tests.h @@ -11,6 +11,7 @@ #include "material.h" #include "CascadeEnrich.h" +#include "CascadeConfig.h" namespace mbmore { diff --git a/src/CentrifugeConfig.cc b/src/CentrifugeConfig.cc index 6ed994f..43fcb00 100644 --- a/src/CentrifugeConfig.cc +++ b/src/CentrifugeConfig.cc @@ -8,123 +8,127 @@ #include #include -/* +/* @article{alexander_glaser_characteristics_2008, - title = {Characteristics of the {Gas} {Centrifuge} for {Uranium} {Enrichment} and {Their} {Relevance} for {Nuclear} {Weapon} {Proliferation}}, - issn = {1547-7800}, - doi = {10.1080/08929880802335998}, - number = {16}, - journal = {Science and Global Security}, - author = {Alexander Glaser}, - year = {2008}, - pages = {1--25}, - annote = {From Tamara}, - file = {2008aglaser_sgsvol16.pdf:/Users/mouginot/Zotero/storage/UAGA49FS/2008aglaser_sgsvol16.pdf:application/pdf} + title = {Characteristics of the {Gas} {Centrifuge} for {Uranium} +{Enrichment} and {Their} {Relevance} for {Nuclear} {Weapon} {Proliferation}}, + issn = {1547-7800}, + doi = {10.1080/08929880802335998}, + number = {16}, + journal = {Science and Global Security}, + author = {Alexander Glaser}, + year = {2008}, + pages = {1--25}, + annote = {From Tamara}, + file = +{2008aglaser_sgsvol16.pdf:/Users/mouginot/Zotero/storage/UAGA49FS/2008aglaser_sgsvol16.pdf:application/pdf} } */ namespace mbmore { CentrifugeConfig::CentrifugeConfig() { -// Standard Values for a P1-type centrifuges + // Standard Values for a P1-type centrifuges v_a = 320; height = 1.8; diameter = 0.10; feed = 15. / 1000. / 1000.; temp = 320; - eff = 1.0; x = 1000; - flow_internal = 2.0; + flow_ratio = 2.0; + design_cut = 0.5; // default UF6 0.352 kg/mol M = 0.352; - + // default 0.003kg/mol (u235/U238) dM = 0.003; // SEE GLASER P2 D_rho = 2.2e-5; // kg/m/s gas_const = 8.314; // J/K/mol M_238 = 0.238; // kg/mol - secpermonth = 60. * 60. * 24. * (365.25 / 12.); } -CentrifugeConfig::CentrifugeConfig(double v_a_, double h_, double d_, double feed_, double T_, - double eff_, double M_, double dM_, double x_, double flow_){ +CentrifugeConfig::CentrifugeConfig(double v_a_, double h_, double d_, + double feed_, double T_, double M_, + double dM_, double x_, double flow_) { v_a = v_a_; height = h_; diameter = d_; feed = feed_; temp = T_; - eff = eff_; x = x_; - flow_internal = flow_; + flow_ratio = flow_; + design_cut = 0.5; M = M_; dM = dM_; D_rho = 2.2e-5; // kg/m/s gas_const = 8.314; // J/K/mol M_238 = 0.238; // kg/mol - secpermonth = 60. * 60. * 24. * (365.25 / 12.); } double CentrifugeConfig::ComputeDeltaU(double cut) { - // Inputs that are effectively constants: double a = diameter / 2.0; // outer radius - // withdrawl radius for heavy isotope + // withdrawal radius for heavy isotope // Glaser 2009 says operationally it ranges from 0.96-0.99 double r2 = 0.975 * a; // Glaser: if v_a < 380 then r1_over_r2 = cte = 0.534 double r1_over_r2 = 0; - if (v_a > 380){ - r1_over_r2 = std::sqrt( - 1.0 - (2.0 * gas_const * temp * log(x)) / M / (v_a*v_a) ); + if (v_a > 380) { + r1_over_r2 = + std::sqrt(1.0 - (2.0 * gas_const * temp * log(x)) / M / (v_a * v_a)); } else { r1_over_r2 = 0.534; } - double r1 = r2 * r1_over_r2; // withdrawl radius for ligher isotope + double r1 = r2 * r1_over_r2; // withdrawal radius for lighter isotope // Glaser eqn 12 // Vertical location of feed - double Z_p = height * (1.0 - cut) * (1.0 + flow_internal) / - (1.0 - cut + flow_internal); + double Z_p = + height * (1.0 - design_cut) * (1.0 + flow_ratio) / (1.0 - design_cut + flow_ratio); // Glaser eqn 3 // Converting from molecular mass to atomic mass (assuming U238) // Warning: This assumption is only valid at low enrichment // TODO: EXPLICITLY CALCULATE ATOMIC MASS GIVEN FEED ASSAY - double C1 = (2.0 * M_PI * (D_rho * M_238 / M) / (log(r2 / r1))); - // double C1 = (2.0 * M_PI * D_rho / (log(r2 / r1))); - double A_p = C1 * (1.0 / feed) * - (cut / ((1.0 + flow_internal) * (1.0 - cut + flow_internal))); - double A_w = C1 * (1.0 / feed) * - ((1.0 - cut) / (flow_internal * (1.0 - cut + flow_internal))); + + // converting feed from UF6 to U + double ratio_UF6_U = M_238 / M; + double feed_U = feed * ratio_UF6_U; + + double C1 = (2.0 * M_PI * D_rho * ratio_UF6_U / (log(r2 / r1))); + double A_p = C1 * (1.0 / feed_U) * + (cut / ((1.0 + flow_ratio) * (1.0 - cut + flow_ratio))); + double A_w = C1 * (1.0 / feed_U) * + ((1.0 - cut) / (flow_ratio * (1.0 - cut + flow_ratio))); double C_therm = CalcCTherm(v_a, temp, dM); // defining terms in the Ratz equation - double r1_over_r2_sq = r1_over_r2*r1_over_r2; - double C_scale = pow((r2 / a), 4) * (1 - r1_over_r2_sq)*(1 - r1_over_r2_sq); - double bracket1 = (1 + flow_internal) / cut; + double r1_over_r2_sq = pow(r1_over_r2, 2); + double C_scale = pow((r2 / a), 4) * pow(1 - r1_over_r2_sq, 2); + double bracket1 = (1 + flow_ratio) / cut; double exp1 = exp(-1.0 * A_p * Z_p); - double bracket2 = flow_internal / (1 - cut); + double bracket2 = flow_ratio / (1 - cut); double exp2 = exp(-1.0 * A_w * (height - Z_p)); // Glaser eqn 10 (Ratz Equation) double major_term = - 0.5 * cut * (1.0 - cut) * (C_therm*C_therm) * C_scale * + 0.5 * cut * (1.0 - cut) * (C_therm * C_therm) * C_scale * pow(((bracket1 * (1 - exp1)) + (bracket2 * (1 - exp2))), 2); // kg/s - double del_U = feed * major_term * eff; // kg/s + double del_U = feed_U * major_term; // kg/s return del_U; } double CentrifugeConfig::CalcCTherm(double v_a, double temp, double dM) { - return dM * (v_a*v_a) / (2.0 * gas_const * temp); + return dM * (v_a * v_a) / (2.0 * gas_const * temp); } double CentrifugeConfig::CalcV(double assay) { diff --git a/src/CentrifugeConfig.h b/src/CentrifugeConfig.h index f04f0e2..ea107d9 100644 --- a/src/CentrifugeConfig.h +++ b/src/CentrifugeConfig.h @@ -12,27 +12,27 @@ class CentrifugeConfig { public: CentrifugeConfig(); CentrifugeConfig(double v_a, double h, double d, double feed, double T, - double eff, double M, double dM, double x, double i_flow); + double M, double dM, double x, double i_flow); double ComputeDeltaU(double cut); // compute the solution of the Raetz equation using all the paramters valus and the provided cut value double dM; // molar mass diff between the 2 components of the gas default 0.003kg/mol (u235/U238) double M; // gas molar mass, default UF6 0.352 kg/mol double x; // pressure ration -> drive r1/r2 ratio - double flow_internal; // countercurrent-to-feed ratios k range between 2 and 4 + double flow_ratio; // countercurrent-to-feed ratios k range between 2 and 4 double v_a; // peripheral velocities double height; // centrifugre height double diameter; // centrifuge diameter double feed; // feed flow rate double temp; // average gas temperature - double eff; // efficiency private: double D_rho; // kg/m/s double gas_const; // J/K/mol double M_238; // kg/mol - double secpermonth; // obvious ?! + + double design_cut; // used for calculating the rectifier length // two method to compute some part of the solution to the Raetz equation double CalcCTherm(double v_a, double temp, double dM); diff --git a/src/CentrifugeConfig_tests.cc b/src/CentrifugeConfig_tests.cc index b22421a..ee571e6 100644 --- a/src/CentrifugeConfig_tests.cc +++ b/src/CentrifugeConfig_tests.cc @@ -18,11 +18,11 @@ TEST(CentrifugeConfig_Test, TestSWU) { double x = 1000; // Pressure ratio (Glaser) // General cascade assumptions - double flow_internal = 2.0; - double eff = 1.0; + double flow_ratio = 2.0; double cut = 0.5; // Centrifuge params used in Python test code + // found at CNERG/enrich_calc:mbmore_unit_tests // (based on Glaser SGS 2009 paper) double v_a = 485; // m/s double height = 0.5; // meters @@ -30,12 +30,12 @@ TEST(CentrifugeConfig_Test, TestSWU) { double feed_m = 15 * 60 * 60 / ((1e3) * 60 * 60 * 1000.0); // kg/sec double temp = 320.0; // Kelvin - CentrifugeConfig centrifuge(v_a, height, diameter, feed_m, temp, eff, M, dM, - x, flow_internal); - // del U=7.0323281e-08 alpha=1.16321 + CentrifugeConfig centrifuge(v_a, height, diameter, feed_m, temp, M, dM, + x, flow_ratio); + // dU=8.638345e-08 alpha=1.130517 double delU = centrifuge.ComputeDeltaU(cut); - double pycode_U = 7.03232816847e-08; + double pycode_U = 8.638345e-08; double tol = 1e-9; EXPECT_NEAR(delU, pycode_U, tol); diff --git a/src/StageConfig.cc b/src/StageConfig.cc new file mode 100644 index 0000000..b7d5b44 --- /dev/null +++ b/src/StageConfig.cc @@ -0,0 +1,188 @@ +// Implements the CascadeEnrich class +#include "StageConfig.h" +#include "cyclus.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace mbmore { + +StageConfig::StageConfig(CentrifugeConfig cent, double f_assay, + double precision__, double feed__) + : centrifuge(cent), + feed_assay_(f_assay), + feed_flow_(feed__), + precision_(precision__) { + BuildIdealStg(); +} +StageConfig::StageConfig(double f_assay, double feed__, double cut__, + double DU__, double alpha__, double precision__) + : centrifuge(), + feed_assay_(f_assay), + feed_flow_(feed__), + precision_(precision__), + cut_(cut__), + DU_(DU__), + alpha_(alpha__) { + // if alpha is not provided, compute it from the dU + if (alpha_ == -1) { + AlphaByDU(); + } + + BetaByAlphaAndCut(); + ProductAssay(); + TailAssay(); +} + +// Search for cut where alpha = beta by starting at cut = 0.1, 0.9 +void StageConfig::CutForIdealStg() { + // Calculating high and low parameters + double low_cut = cut_ = 0.1; + DU_ = centrifuge.ComputeDeltaU(cut_); + AlphaByDU(); + BetaByAlphaAndCut(); + double low_alpha = alpha_; + double low_beta = beta_; + double high_cut = cut_ = 0.9; + DU_ = centrifuge.ComputeDeltaU(cut_); + AlphaByDU(); + BetaByAlphaAndCut(); + double high_alpha = alpha_; + double high_beta = beta_; + + double p_alpha, p_beta, p_cut; + // Set initial guess to closer cut values + if (std::abs(low_alpha - low_beta) < std::abs(high_alpha - high_beta)) { + alpha_ = low_alpha; + beta_ = low_beta; + cut_ = low_cut; + + p_alpha = high_alpha; + p_beta = high_beta; + p_cut = high_cut; + } else { + alpha_ = high_alpha; + beta_ = high_beta; + cut_ = high_cut; + + p_alpha = low_alpha; + p_beta = low_beta; + p_cut = low_cut; + } + + double p_alpha_minus_beta = p_alpha - p_beta; + while (std::abs(alpha_ - beta_) > precision_) { + // a*cut +b =y + double alpha_minus_beta = alpha_ - beta_; + double a = (p_alpha_minus_beta - alpha_minus_beta) / (p_cut - cut_); + double b = alpha_minus_beta - cut_ * a; + // old = new + p_alpha_minus_beta = alpha_minus_beta; + p_cut = cut_; + // targeting alpha_minus_beta = 0 + cut_ = -b / a; + DU_ = centrifuge.ComputeDeltaU(cut_); + AlphaByDU(); + BetaByAlphaAndCut(); + } +} + +void StageConfig::ProductAssay() { + double ratio = alpha_ * feed_assay_ / (1. - feed_assay_); + product_assay_ = ratio / (1. + ratio); +} + +void StageConfig::TailAssay() { + double A = (feed_assay_ / (1. - feed_assay_)) / beta_; + tail_assay_ = A / (1. + A); +} + +void StageConfig::AlphaByDU() { + double M = centrifuge.M; // UF6 kg/mol + double M_238 = 0.238; // U kg/mol + double ratio_UF6_U = M_238 / M; + + // converting feed from UF6 to U + double feed_U = centrifuge.feed * ratio_UF6_U; + // "Uranium Enrichment By Gas Centrifuge" D.G. Avery & E. Davies pg. 18 + alpha_ = 1. + std::sqrt((2. * DU_ * (1. - cut_) / (cut_ * feed_U))); +} + +void StageConfig::AlphaByProductAssay() { + alpha_ = + (1 - feed_assay_) / feed_assay_ * product_assay_ / (1 - product_assay_); +} + +void StageConfig::BetaByAlphaAndCut() { + ProductAssay(); + double waste_assay = (feed_assay_ - cut_ * product_assay_) / (1. - cut_); + + beta_ = feed_assay_ / (1. - feed_assay_) * (1. - waste_assay) / waste_assay; +} + +void StageConfig::ProductAssayByGamma(double gamma) { + product_assay_ = (-sqrt(pow((feed_assay_ - cut_) * gamma, 2) + + gamma * (2 * feed_assay_ + 2 * cut_ - + 2 * pow(feed_assay_, 2) - 2 * pow(cut_, 2)) + + pow(feed_assay_ + cut_ - 1, 2)) + + gamma * (feed_assay_ + cut_) - feed_assay_ - cut_ + 1) / + (2 * cut_ * (gamma - 1)); +} + +void StageConfig::CutByAlphaBeta() { + ProductAssay(); + TailAssay(); + + cut_ = (feed_assay_ - tail_assay_) / (product_assay_ - tail_assay_); +} + +void StageConfig::BuildIdealStg() { + if (DU_ == -1 || alpha_ == -1) { + CutForIdealStg(); + DU_ = centrifuge.ComputeDeltaU(cut_); + AlphaByDU(); + } + + beta_ = alpha_; + CutByAlphaBeta(); + ProductAssay(); + TailAssay(); +} + +void StageConfig::MachinesNeededPerStage(double tolerance) { + // n_machines: the denominator should be equal to the + // centrifuge feed flow (centrifuge.feed). + + double M = centrifuge.M; // UF6 kg/mol + double M_238 = 0.238; // U kg/mol + double ratio_UF6_U = M_238 / M; + + // "Uranium Enrichment By Gas Centrifuge" D.G. Avery & E. Davies pg. 18 + double centrifuge_feed_flow = + 2 * DU_ * ((1 - cut_) / cut_) / pow((alpha_ - 1.), 2.) / ratio_UF6_U; + double n_exact = (*this).feed_flow_ / (centrifuge_feed_flow); + + // Adds a machine if fractional amount is needed + n_machines_ = int(n_exact); + if (std::abs(n_exact - n_machines_) > tolerance) { + n_machines_ = int(n_exact) + 1; + } +} +double StageConfig::SWU() { + using cyclus::toolkit::Assays; + using cyclus::toolkit::SwuRequired; + if (product_assay_ == -1 || feed_assay_ == -1 || tail_assay_ == 1 || + feed_flow_ == -1) + return -1; + + return SwuRequired(feed_flow_, + Assays(feed_assay_, product_assay_, tail_assay_)); +} + +} // namespace mbmore diff --git a/src/StageConfig.h b/src/StageConfig.h new file mode 100644 index 0000000..59f10b8 --- /dev/null +++ b/src/StageConfig.h @@ -0,0 +1,110 @@ +#ifndef MBMORE_SRC_STAGE_H_ +#define MBMORE_SRC_STAGE_H_ + +#include +#include "CentrifugeConfig.h" + +namespace mbmore { +class CascadeConfig; + +class StageConfig { + + private: + + // Compute the cut to ensure alpha = beta (from dU) + void CutForIdealStg(); + + // Precision used for the cut calculation defautl 1e-8 + double precision_ = 1e-8; + + // cut value of the stage + double cut_; + // dU value of the stage (calculated form the centrifuges config with the cut) + double DU_ = -1; + // Feed to Product enrichment ratio + double alpha_ = -1; + // Feed to Tail enrichment ratio + double beta_; + // Feed flow (g/s) + double feed_flow_; + + // number of centrifuges in the stage + int n_machines_ = -1; + + // Feed assay + double feed_assay_; + // Product assay + double product_assay_; + // Tail assay + double tail_assay_; + + public: + // Setup a empty stage + StageConfig() { ; } + // Design an ideal stage for a specific feed assay and feed flow + StageConfig(CentrifugeConfig cent, double f_assay, + double precision = 1e-8, double feed_flow = -1); + // Design an ideal stage for a specific feed assay and feed flow + StageConfig(double f_assay, double feed_flow, double cut_, double DU_, + double alpha_ = -1, double precision = 1e-8); + + // Build a stage assuming alpha = beta + // (if cut is not defined, compute the cut to make it so) + void BuildIdealStg(); + + // calculate Alpha value using the dU, Cut and the centrifuge feed flow value + void AlphaByDU(); + + // calculate Alpha from the feed assay and the product assay; + void AlphaByProductAssay(); + + // Compute Beta value from alpha value, cut and feed assay + void BetaByAlphaAndCut(); + // recompute Cut value assuming Alpha and Beta fixed + void CutByAlphaBeta(); + + // return the SWU of the stage + double SWU(); + + // Compute Product from gamma + void ProductAssayByGamma(double gamma); + + // Compute Product assay from feed assay and alpha + void ProductAssay(); + // Compute Waste assay from feed assay and beta + void TailAssay(); + + // Return the minimum number of centrifuges required to meet the feed flow + void MachinesNeededPerStage(double tolerance = 0.01); + + // Configuration of all the centrifuges in the stage + // Default centrifuge initialized + CentrifugeConfig centrifuge; + + // Setter methods + void precision(double p) {precision_ = p;} + void cut(double c) {cut_ = c;} + void DU(double du) {DU_ = du;} + void alpha(double a) {alpha_ = a;} + void beta(double b) {beta_ = b;} + void feed_flow(double f) {feed_flow_ = f;} + void feed_assay(double fa) {feed_assay_ = fa;} + double product_assay(double pa) {product_assay_ = pa;} + double tail_assay(double ta) {tail_assay_ = ta;} + + // Getter methods + double precision() {return precision_;} + double cut() {return cut_;} + double DU() {return DU_;} + double alpha() {return alpha_;} + double beta() {return beta_;} + double feed_flow() {return feed_flow_;} + int n_machines() {return n_machines_;} + double feed_assay() {return feed_assay_;} + double product_assay() {return product_assay_;} + double tail_assay() {return tail_assay_;} +}; + +} // namespace mbmore + +#endif // MBMORE_SRC_STAGE_H_ diff --git a/src/StageConfig_tests.cc b/src/StageConfig_tests.cc new file mode 100644 index 0000000..eeae9e0 --- /dev/null +++ b/src/StageConfig_tests.cc @@ -0,0 +1,207 @@ +#include + +#include "StageConfig.h" + +#include "agent_tests.h" +#include "context.h" +#include "facility_tests.h" + +namespace mbmore { + +// Benchmarked using a regression test (expected values calculated manually) +// Calculations can be found at CNERG/enrich_calc:mbmore_unit_tests +namespace stageconfig_test { +// Fixed for a cascade separating out U235 from U238 in UF6 gas +double M = 0.352; // kg/mol UF6 +double dM = 0.003; // kg/mol U238 - U235 +double x = 1000; // Pressure ratio (Glaser) + +// General cascade assumptions +double flow_ratio = 2.0; +double cut = 0.5; + +// Centrifgue parameters based on Glaser SGS 2009 paper +double v_a = 485; // m/s +double height = 0.5; // meters +double diameter = 0.15; // meters +double feed_m = 15 * 60 * 60 / ((1e3) * 60 * 60 * 1000.0); // kg/sec +double temp = 320.0; // Kelvin + +// Cascade params used in in calculating expected values +const double feed_assay = 0.0071; +const double prod_assay = 0.035; +const double waste_assay = 0.001; +const double feed_c = 739 / (30.4 * 24 * 60 * 60); // kg/month -> kg/sec +const double product_c = 77 / (30.4 * 24 * 60 * 60); // kg/month -> kg/sec +CentrifugeConfig centrifuge(v_a, height, diameter, feed_m, temp, M, dM, x, + flow_ratio); + +// del U=8.638345e-08 alpha=1.130517 +double delU = centrifuge.ComputeDeltaU(cut); + +const double tol_assay = 1e-5; +const double tol_qty = 1e-6; +const double tol_num = 1e-2; + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Find product assay from separation factor alpha +TEST(StageConfig_Test, TestAssays) { + double cur_alpha = 1.4; + double cur_f_assay = 0.007; + + StageConfig stage(cur_f_assay, feed_m, cut, delU, cur_alpha, 1e-16); + stage.ProductAssay(); + + // N_prime = alpha*R / ( 1+alpha*R) + double target_prod_assay = 0.009773; + double tol = 1e-6; + + EXPECT_NEAR(stage.product_assay(), target_prod_assay, tol); + + double n_stages = 5; + double target_w_assay = 0.004227; + stage.BetaByAlphaAndCut(); + stage.TailAssay(); + + EXPECT_NEAR(stage.tail_assay(), target_w_assay, tol); +} +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Calculate ideal SWU params of single machinefeed_assay +// (separation potential delU and separation factor alpha) +TEST(StageConfig_Test, TestSWU) { + double expected_U = 8.638345e-08; + double tol = 1e-9; + + StageConfig stage(feed_assay, feed_m, cut, delU, -1, 1e-16); + + double expected_alpha = 1.130517; + double tol_alpha = 1e-2; + EXPECT_NEAR(stage.alpha(), expected_alpha, tol_alpha); +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Calculate the product assay for an ideal stage configuration. +TEST(StageConfig_Test, TestIdeal) { + StageConfig stage_ideal(feed_assay, feed_m, cut, -1, -1, 1e-16); + + // Only setting precision for building ideal stage + stage_ideal.precision(1e-3); + stage_ideal.BuildIdealStg(); + stage_ideal.precision(1e-16); + + double expected_alpha = 1.130517; + double tol_alpha = 1e-2; + + // These expected values are found by regression test + // and are working as intended as of committing this line. + double expected_cut = 0.4685952; + double tol_cut = 1e-3; + double expected_U = 8.281007e-08; + double tol_DU = 1e-9; + + EXPECT_NEAR(stage_ideal.alpha(), expected_alpha, tol_alpha); + EXPECT_NEAR(stage_ideal.cut(), expected_cut, tol_cut); + EXPECT_NEAR(stage_ideal.DU(), expected_U, tol_DU); +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +TEST(StageConfig_Test, ProductAssayByGamma) { + double gamma = 1.3798316056650026; + double target_product_assay = 0.00822; + double theta_ = 0.46040372309; + double feed_assay_ = 0.007; + StageConfig stage(feed_assay_, feed_c, theta_, delU, -1, 1e-16); + stage.ProductAssayByGamma(gamma); + + EXPECT_NEAR(target_product_assay,stage.product_assay(), 1e-5); +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +TEST(StageConfig_Test, AlphaByProductAssay) { + double gamma = 1.3798316056650026; + double target_product_assay = 0.00821; + double theta_ = 0.46040372309; + double feed_assay_ = 0.007; + StageConfig stage(feed_assay_, feed_c, theta_, delU, -1, 1e-16); + stage.feed_assay(0.1); + stage.product_assay(0.3); + double alpha_ = 0.3/(1-0.3)*(1-0.1)/0.1; + stage.AlphaByProductAssay(); + + EXPECT_NEAR(alpha_,stage.alpha(), 1e-5); +} + + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Determine the output of the first enrich/strip stage of a cascade +// based on the design params for the cascade +TEST(StageConfig_Test, TestStages) { + StageConfig stage(feed_assay, feed_c, cut, delU, -1, 1e-16); + + stage.ProductAssay(); + stage.MachinesNeededPerStage(); + + double expected_product_assay_s = 0.0080192; + + // Calculated using equations from 2009 Glaser paper + int expected_n_mach_e = ceil(feed_c/feed_m); + + EXPECT_NEAR(stage.n_machines(), expected_n_mach_e, tol_num); + EXPECT_NEAR(stage.product_assay(), expected_product_assay_s, tol_assay); +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Verify that, for a P1-type centrifuge, results are in agreement +// with those presented by Glaser's paper. +TEST(StageConfig_Test, AlphaBeta) { + // P1-type parameters + v_a = 320; // m/s + x = 1000.; + flow_ratio = 2.0; + temp = 320; // K + M = 0.352; // kg/mol UF6 + dM = 0.003; // kg/mol U238 - U235 + cut = 0.5; + double Z = 1.8; // m + double d = 0.10; // m + double feed = 12.6 / 1000. / 1000.; // mg/s -> kg/s + double f_assay = 0.00720; + + CentrifugeConfig cent(v_a, Z, d, feed, temp, M, dM, x, flow_ratio); + + double delU = cent.ComputeDeltaU(cut); + + // Glaser's paper value is delU_e = 2.5 / (365.25 * 24 * 60 * 60) + // The answer below is apprx. 6% different. + // WIP to reconcile the difference. + double delU_e = 8.40508e-8; + EXPECT_NEAR(delU, delU_e, 1e-8); + + StageConfig stg(f_assay, feed, cut, delU); + + double ab_e = 1.29; + double tol_ab = 0.01; + EXPECT_NEAR(stg.alpha() * stg.beta(), ab_e, tol_ab); + +} + +// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +// Verify that the centrifuge feed flow can be reconstructed if alpha +// is not provided. +TEST(StageConfig_Test, FeedFlow) { + StageConfig stage(feed_assay, feed_c, cut, delU, -1, 1e-16); + + double M = centrifuge.M; // UF6 kg/mol + double M_238 = 0.238; // U kg/mol + double ratio_UF6_U = M_238 / M; + + // "Uranium Enrichment By Gas Centrifuge" D.G. Avery & E. Davies pg. 18 + double centrifuge_feed_flow = + 2 * stage.DU() * ((1 - stage.cut()) / stage.cut()) + / pow((stage.alpha() - 1.), 2.) / ratio_UF6_U; + + EXPECT_NEAR(centrifuge_feed_flow,stage.centrifuge.feed,1e-6); +} + +} // namespace enrichfunctiontests +} // namespace mbmore diff --git a/src/enrich_functions.cc b/src/enrich_functions.cc index 29a5471..c0c593b 100644 --- a/src/enrich_functions.cc +++ b/src/enrich_functions.cc @@ -15,25 +15,6 @@ double gas_const = 8.314; // J/K/mol double M_238 = 0.238; // kg/mol -double AlphaBySwu(double del_U, double feed, double cut, double M) { - double alpha = 1 + std::sqrt((2 * (del_U / M) * (1 - cut) / (cut * feed))); - return alpha; -} - -// per machine -double ProductAssayByAlpha(double alpha, double feed_assay) { - // Possibly incorrect is commented out ? - // double ratio = (1.0 - feed_assay) / (alpha * feed_assay); - // return 1.0 / (ratio + 1.0); - double ratio = alpha * feed_assay / (1.0 - feed_assay); - return ratio / (1 + ratio); -} - -double WasteAssayByAlpha(double alpha, double feed_assay) { - double A = (feed_assay / (1 - feed_assay)) / alpha; - return A / (1 + A); -} - // This equation can only be used in the limit where the separation factor // (alpha) is very close to one, which is not true for modern gas centrifuges // DO NOT USE THIS EQUATION!!! @@ -63,38 +44,6 @@ std::pair StagesPerCascade(double alpha, double feed_assay, // Determine number of stages required to reach ideal cascade product assay // (requires integer number of stages, so output may exceed target assay) -std::pair FindNStages(double alpha, double feed_assay, - double product_assay, double waste_assay) { - using std::pair; - - int ideal_enrich_stage = 0; - int ideal_strip_stage = 0; - double stage_feed_assay = feed_assay; - double stage_product_assay = feed_assay; - double stage_waste_assay = feed_assay; // start w/waste of 1st enrich stage - - // Calculate number of enriching stages - while (stage_product_assay < product_assay) { - stage_product_assay = ProductAssayByAlpha(alpha, stage_feed_assay); - if (ideal_enrich_stage == 0) { - stage_waste_assay = WasteAssayByAlpha(alpha, stage_feed_assay); - } - ideal_enrich_stage += 1; - stage_feed_assay = stage_product_assay; - } - // Calculate number of stripping stages - stage_feed_assay = stage_waste_assay; - while (stage_waste_assay > waste_assay) { - stage_waste_assay = WasteAssayByAlpha(alpha, stage_feed_assay); - ideal_strip_stage += 1; - stage_feed_assay = stage_waste_assay; - } - - std::pair stages = - std::make_pair(ideal_enrich_stage, ideal_strip_stage); - return stages; -} - double ProductAssayFromNStages(double alpha, double feed_assay, double enrich_stages) { double A = @@ -154,86 +103,6 @@ double DelUByCascadeConfig(double product_assay, double waste_assay, return U_cascade / feed_assay; } -// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -// Calculate steady-state flow rates into each cascade stage -// Linear system of equations in form AX = B, where A is nxn square matrix -// of linear equations for the flow rates of each stage and B are the external -// feeds for the stage. External feed is zero for all stages accept cascade -// feed stage (F_0) stages start with last strip stage [-2, -1, 0, 1, 2] -// http://www.physics.utah.edu/~detar/phys6720/handouts/lapack.html -// -std::vector CalcFeedFlows(std::pair n_st, double cascade_feed, - double cut) { - // This is the Max # of stages in cascade. It cannot be passed in due to - // how memory is allocated and so must be hardcoded. It's been chosen - // to be much larger than it should ever need to be - int max_stages = 100; - - int n_enrich = n_st.first; - int n_strip = n_st.second; - int n_stages = n_st.first + n_st.second; - - // LAPACK takes the external flow feeds as B, and then returns a modified - // version of the same array now representing the solution flow rates. - - // Build Array with pointers - double flow_eqns[max_stages][max_stages]; - double flows[1][max_stages]; - - // build matrix of equations in this pattern - // [[ -1, 1-cut, 0, 0, 0] [[0] - // [cut, -1, 1-cut, 0, 0] [0] - // [ 0, cut, -1, 1-cut, 0] * X = [-1*cascade_feed] - // [ 0, 0, cut, -1, 1-cut] [0] - // [ 0, 0, 0, cut, -1]] [0]] - // - for (int row_idx = 0; row_idx < max_stages; row_idx++) { - // fill the array with zeros, then update individual elements as nonzero - flows[0][row_idx] = 0; - for (int fill_idx = 0; fill_idx < max_stages; fill_idx++) { - flow_eqns[fill_idx][row_idx] = 0; - } - // Required do to the artificial 'Max Stages' defn. Only calculate - // non-zero matrix elements where stages really exist. - if (row_idx < n_stages) { - int i = row_idx - n_strip; - int col_idx = n_strip + i; - flow_eqns[col_idx][row_idx] = -1; - if (col_idx != 0) { - flow_eqns[col_idx - 1][row_idx] = cut; - } - if (col_idx != n_stages - 1) { - flow_eqns[col_idx + 1][row_idx] = (1 - cut); - } - // Add the external feed for the cascade - if (i == 0) { - flows[0][row_idx] = -1 * cascade_feed; - } - } - } - - // LAPACK solver variables - int nrhs = 1; // 1 column solution - int lda = max_stages; // must be >= MAX(1,N) - int ldb = max_stages; // must be >= MAX(1,N) - int ipiv[max_stages]; - int info; - - // Solve the linear system - dgesv_(&n_stages, &nrhs, &flow_eqns[0][0], &lda, ipiv, &flows[0][0], &ldb, - &info); - - // Check for success - if (info != 0) { - std::cerr << "LAPACK linear solver dgesv returned error " << info << "\n"; - } - - std::vector final_flows; - for (int i = 0; i < n_stages; i++) { - final_flows.push_back(flows[0][i]); - } - return final_flows; -} // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - // Determine number of machines in each stage of the cascade, and total // output flow from each stage @@ -298,16 +167,6 @@ std::vector> CalcStageFeatures( return stage_info; } -// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -// Determine total number of machines in the cascade from machines per stage -int FindTotalMachines(std::vector> stage_info) { - int machines_needed = 0; - std::vector>::const_iterator it; - for (it = stage_info.begin(); it != stage_info.end(); it++) { - machines_needed += it->first; - } - return machines_needed; -} // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - std::pair DesignCascade(double design_feed, diff --git a/src/enrich_functions.h b/src/enrich_functions.h index db6e707..cba6f99 100644 --- a/src/enrich_functions.h +++ b/src/enrich_functions.h @@ -15,26 +15,6 @@ extern "C" { } - // Calculates the separations factor given the ideal separation energy of a - // single machine (del_U has units of moles/sec) - // Avery p 18 - double AlphaBySwu(double del_U, double feed, double cut, double M); - - // Calculates the assay of the product given the assay - // of the feed and the theoretical separation factor of the machine - // Avery p 57 - double ProductAssayByAlpha(double alpha, double feed_assay); - - // Calculates the assay of the waste given the assay - // of the feed and the theoretical separation factor of the machine - // Avery p 59 (per machine) - double WasteAssayByAlpha(double alpha, double feed_assay); - - // Calculates the number of stages needed in a cascade given the separation - // potential of a single centrifuge and the material assays - std::pair - FindNStages(double alpha, double feed_assay, double product_assay, - double Nwc); // Calculates the product assay after N enriching stages double ProductAssayFromNStages(double alpha, double feed_assay, @@ -83,10 +63,6 @@ extern "C" { double product_flow, double waste_flow, double feed_assay); - // Solves system of linear eqns to determine steady state flow rates - // in each stage of cascade - std::vector CalcFeedFlows(std::pair n_st, - double cascade_feed, double cut); // Determines the number of machines and product in each stage based // on the steady-state flows defined for the cascade. @@ -96,8 +72,6 @@ extern "C" { std::pair n_st, std::vector feed_flow); - // Determine total number of machines in the cascade from machines per stage - int FindTotalMachines(std::vector> stage_info); std::pair DesignCascade( double design_feed, double design_alpha, double design_delU, double cut, diff --git a/src/enrich_functions_tests.cc b/src/enrich_functions_tests.cc deleted file mode 100644 index 6421450..0000000 --- a/src/enrich_functions_tests.cc +++ /dev/null @@ -1,208 +0,0 @@ -#include - -#include "enrich_functions.h" - -#include "agent_tests.h" -#include "context.h" -#include "facility_tests.h" - -namespace mbmore { - - // Benchmarked against mbmore_enrich_compare.ipynb - // https://github.com/mbmcgarry/data_analysis/tree/master/enrich_calcs - namespace enrichfunctiontests { - // Fixed for a cascade separating out U235 from U238 in UF6 gas - const double M = 0.352; // kg/mol UF6 - const double dM = 0.003; // kg/mol U238 - U235 - const double x = 1000; // Pressure ratio (Glaser) - - // General cascade assumptions - const double flow_internal = 2.0 ; - const double eff = 1.0; - const double cut = 0.5; - - // Centrifuge params used in Python test code - // (based on Glaser SGS 2009 paper) - const double v_a = 485; // m/s - const double height = 0.5; // meters - const double diameter = 0.15; // meters - const double feed_m = 15 * 60 * 60 / ((1e3) * 60 * 60 * 1000.0); // kg/sec - const double temp = 320.0 ; //Kelvin - - // Cascade params used in Python test code (Enrichment_Calculations.ipynb) - const double feed_assay = 0.0071; - const double product_assay = 0.035; - const double waste_assay = 0.001; - const double feed_c = 739 / (30.4 * 24 * 60 * 60); // kg/month -> kg/sec - const double product_c = 77 / (30.4 * 24 * 60 * 60); // kg/month -> kg/sec - - //del U=7.0323281e-08 alpha=1.16321 - double delU = CalcDelU(v_a, height, diameter, feed_m, temp, cut, eff, - M, dM, x, flow_internal); - - double alpha = AlphaBySwu(delU, feed_m, cut, M); - const double tol_assay = 1e-5; - const double tol_qty = 1e-6; - const double tol_num = 1e-2; - -// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -// Find product assay from separation factor alpha -TEST(Enrich_Functions_Test, TestAssays) { - - double cur_alpha = 1.4; - double cur_f_assay = 0.007; - - double cpp_assay = ProductAssayByAlpha(cur_alpha, cur_f_assay); - - double pycode_assay = 0.009772636; - double tol = 1e-6; - - EXPECT_NEAR(cpp_assay, pycode_assay, tol); - - double n_stages = 5; - double pycode_w_assay = 0.00095311 ; - - double cpp_w_assay = WasteAssayFromNStages(cur_alpha, cur_f_assay, n_stages); - - EXPECT_NEAR(cpp_w_assay, pycode_w_assay, tol); - -} - -// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -// Ideal cascade design, and then using away from ideal design - -TEST(Enrich_Functions_Test, TestCascade) { - double n_machines = MachinesPerCascade(delU, product_assay, - waste_assay, feed_c, product_c); - - double pycode_n_mach = 25990.392; - EXPECT_NEAR(n_machines, pycode_n_mach, tol_num); - - std::pair n_stages = FindNStages(alpha, feed_assay, - product_assay, - waste_assay); - int pycode_n_enrich_stage = 11; - int pycode_n_strip_stage = 13; - - - // int n_stage_enrich = (int) n_stages.first + 1; // Round up to next integer - // int n_stage_waste = (int) n_stages.second + 1; // Round up to next integer - int n_stage_enrich = n_stages.first; - int n_stage_waste = n_stages.second; - - EXPECT_EQ(n_stage_enrich, pycode_n_enrich_stage); - EXPECT_EQ(n_stage_waste, pycode_n_strip_stage); - - // Now test assays when cascade is modified away from ideal design - // (cascade optimized for natural uranium feed, now use 20% enriched - double feed_assay_mod = 0.20; - - double mod_product_assay = ProductAssayFromNStages(alpha, feed_assay_mod, - n_stage_enrich); - double mod_waste_assay = WasteAssayFromNStages(alpha, feed_assay_mod, - n_stage_waste); - - std::cout << "alpha " << alpha << " feed " << feed_assay_mod << " nstage " << n_stage_enrich << " unrounded stages " << n_stages.first << std::endl; - double pycode_mod_product_assay = 0.60085; - double pycode_mod_waste_assay = 0.0290846; - EXPECT_NEAR(mod_product_assay, pycode_mod_product_assay, tol_assay); - EXPECT_NEAR(mod_waste_assay, pycode_mod_waste_assay, tol_assay); - } - - // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - // Determine the output of the first enrich/strip stage of a cascade - // based on the design params for the cascade - TEST(Enrich_Functions_Test, TestStages) { - double product_assay_s = ProductAssayByAlpha(alpha, feed_assay); - double n_mach_e = MachinesPerStage(alpha, delU, feed_c); - double product_s = ProductPerEnrStage(alpha, feed_assay, - product_assay_s, feed_c); - - double enrich_waste = feed_c - product_s; - double enrich_waste_assay = (feed_c * feed_assay - product_s * - product_assay_s)/enrich_waste; - - double pycode_product_assay_s = 0.0082492; - double pycode_n_mach_e = 53.287; - double pycode_product_s = 0.0001408; - - EXPECT_NEAR(n_mach_e, pycode_n_mach_e, tol_num); - EXPECT_NEAR(product_assay_s, pycode_product_assay_s, tol_assay); - EXPECT_NEAR(product_s, pycode_product_s, tol_qty); - - double n_mach_w = MachinesPerStage(alpha, delU, enrich_waste); - double strip_waste_assay = WasteAssayByAlpha(alpha, enrich_waste_assay); - - // This AVERY EQN doesn't work for some reason - // double strip_waste = WastePerStripStage(alpha, enrich_waste_assay, - // strip_waste_assay, enrich_waste); - - double pycode_n_mach_w = 26.6127; - double pycode_waste_assay_s = 0.005117; - // double pycode_waste_s = 8.60660553717e-05; - - EXPECT_NEAR(n_mach_w, pycode_n_mach_w, tol_num); - EXPECT_NEAR(strip_waste_assay, pycode_waste_assay_s, tol_assay); - // EXPECT_NEAR(strip_waste, pycode_waste_s, tol_qty); - - } - -// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -// tests the steady state flow rates for a cascade -// -TEST(Enrich_Functions_Test, TestCascadeDesign) { - double fa = 0.10; - double pa = 0.20; - double wa = 0.05; - - std::vector pycode_flows = {0.00030693, 0.00061387, 0.0009208 , - 0.00122774, 0.00153467, 0.00127889, - 0.00102311, 0.00076734, 0.00051156, - 0.00025578}; - - - std::vector pycode_machines={59, 117, 175, 233, 291, 243, 194, - 146, 97, 49}; - - std::pair n_stages = FindNStages(alpha, fa, pa, wa); - std::vector flows = CalcFeedFlows(n_stages, feed_c, cut); - - // if # Machines for the stage is within tol_num of an integer - // then round down. Otherwise round up to the next integer machine to - // preserve steady-state flow calculations. - std::vector> stage_info = - CalcStageFeatures(fa, alpha, delU, cut, n_stages, flows); - - for (int i = 0; i < pycode_flows.size(); i++){ - EXPECT_NEAR(flows[i], pycode_flows[i], tol_num); - int nmach = stage_info[i].first; - EXPECT_EQ(nmach, pycode_machines[i]); - } - - // not enough machines - int max_centrifuges = 80; - std::pair design_params = DesignCascade(feed_c, alpha, delU, - cut, max_centrifuges, - n_stages); - int py_tot_mach = 79; - double py_opt_feed = 1.30116169899e-05; - - EXPECT_EQ(py_tot_mach, design_params.first); - EXPECT_NEAR(py_opt_feed, design_params.second, tol_qty); - - // more machines than requested capacity - max_centrifuges = 1000; - design_params = DesignCascade(feed_c, alpha, delU, - cut, max_centrifuges, - n_stages); - py_tot_mach = 986; - py_opt_feed = 0.000172728; - - EXPECT_EQ(py_tot_mach, design_params.first); - EXPECT_NEAR(py_opt_feed, design_params.second, tol_qty); - -} - - - } // namespace enrichfunctiontests -} // namespace mbmore