CYCAMORE
src/enrichment.cc
Go to the documentation of this file.
1 // Implements the Enrichment class
2 #include "enrichment.h"
3 
4 #include <algorithm>
5 #include <cmath>
6 #include <limits>
7 #include <sstream>
8 #include <vector>
9 
10 #include <boost/lexical_cast.hpp>
11 
12 namespace cycamore {
13 
14 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
15 Enrichment::Enrichment(cyclus::Context* ctx)
16  : cyclus::Facility(ctx),
17  tails_assay(0),
18  swu_capacity(0),
19  max_enrich(1),
20  initial_feed(0),
21  feed_commod(""),
22  feed_recipe(""),
23  product_commod(""),
24  tails_commod(""),
25  order_prefs(true) {}
26 
27 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29 
30 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
32  std::stringstream ss;
33  ss << cyclus::Facility::str() << " with enrichment facility parameters:"
34  << " * SWU capacity: " << SwuCapacity()
35  << " * Tails assay: " << tails_assay << " * Feed assay: " << FeedAssay()
36  << " * Input cyclus::Commodity: " << feed_commod
37  << " * Output cyclus::Commodity: " << product_commod
38  << " * Tails cyclus::Commodity: " << tails_commod;
39  return ss.str();
40 }
41 
42 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43 void Enrichment::Build(cyclus::Agent* parent) {
44  using cyclus::Material;
45 
46  Facility::Build(parent);
47  if (initial_feed > 0) {
48  inventory.Push(Material::Create(this, initial_feed,
49  context()->GetRecipe(feed_recipe)));
50  }
51 
52  LOG(cyclus::LEV_DEBUG2, "EnrFac") << "Enrichment "
53  << " entering the simuluation: ";
54  LOG(cyclus::LEV_DEBUG2, "EnrFac") << str();
55 }
56 
57 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59 
60 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61 void Enrichment::Tock() {
62  using cyclus::toolkit::RecordTimeSeries;
63  LOG(cyclus::LEV_INFO4, "EnrFac") << prototype() << " used "
64  << intra_timestep_swu_ << " SWU";
65  RecordTimeSeries<cyclus::toolkit::ENRICH_SWU>(this, intra_timestep_swu_);
66  LOG(cyclus::LEV_INFO4, "EnrFac") << prototype() << " used "
67  << intra_timestep_feed_ << " feed";
68  RecordTimeSeries<cyclus::toolkit::ENRICH_FEED>(this, intra_timestep_feed_);
69 }
70 
71 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72 std::set<cyclus::RequestPortfolio<cyclus::Material>::Ptr>
74  using cyclus::Material;
75  using cyclus::RequestPortfolio;
76  using cyclus::Request;
77 
78  std::set<RequestPortfolio<Material>::Ptr> ports;
79  RequestPortfolio<Material>::Ptr port(new RequestPortfolio<Material>());
80  Material::Ptr mat = Request_();
81  double amt = mat->quantity();
82 
83  if (amt > cyclus::eps_rsrc()) {
84  port->AddRequest(mat, this, feed_commod);
85  ports.insert(port);
86  }
87 
88  return ports;
89 }
90 
91 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92 bool SortBids(cyclus::Bid<cyclus::Material>* i,
93  cyclus::Bid<cyclus::Material>* j) {
94  cyclus::Material::Ptr mat_i = i->offer();
95  cyclus::Material::Ptr mat_j = j->offer();
96 
97  cyclus::toolkit::MatQuery mq_i(mat_i);
98  cyclus::toolkit::MatQuery mq_j(mat_j);
99 
100  return ((mq_i.mass(922350000) / mq_i.qty()) <=
101  (mq_j.mass(922350000) / mq_j.qty()));
102 }
103 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104 // Sort offers of input material to have higher preference for more
105 // U-235 content
107  cyclus::PrefMap<cyclus::Material>::type& prefs) {
108  using cyclus::Bid;
109  using cyclus::Material;
110  using cyclus::Request;
111 
112  if (order_prefs == false) {
113  return;
114  }
115 
116  cyclus::PrefMap<cyclus::Material>::type::iterator reqit;
117 
118  // Loop over all requests
119  for (reqit = prefs.begin(); reqit != prefs.end(); ++reqit) {
120  std::vector<Bid<Material>*> bids_vector;
121  std::map<Bid<Material>*, double>::iterator mit;
122  for (mit = reqit->second.begin(); mit != reqit->second.end(); ++mit) {
123  Bid<Material>* bid = mit->first;
124  bids_vector.push_back(bid);
125  }
126  std::sort(bids_vector.begin(), bids_vector.end(), SortBids);
127 
128  // Assign preferences to the sorted vector
129  double n_bids = bids_vector.size();
130  bool u235_mass = 0;
131 
132  for (int bidit = 0; bidit < bids_vector.size(); bidit++) {
133  int new_pref = bidit + 1;
134 
135  // For any bids with U-235 qty=0, set pref to zero.
136  if (!u235_mass) {
137  cyclus::Material::Ptr mat = bids_vector[bidit]->offer();
138  cyclus::toolkit::MatQuery mq(mat);
139  if (mq.mass(922350000) == 0) {
140  new_pref = -1;
141  } else {
142  u235_mass = true;
143  }
144  }
145  (reqit->second)[bids_vector[bidit]] = new_pref;
146  } // each bid
147  } // each Material Request
148 }
149 
150 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152  const std::vector<std::pair<cyclus::Trade<cyclus::Material>,
153  cyclus::Material::Ptr> >& responses) {
154  // see
155  // http://stackoverflow.com/questions/5181183/boostshared-ptr-and-inheritance
156  std::vector<std::pair<cyclus::Trade<cyclus::Material>,
157  cyclus::Material::Ptr> >::const_iterator it;
158  for (it = responses.begin(); it != responses.end(); ++it) {
159  AddMat_(it->second);
160  }
161 }
162 
163 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164 std::set<cyclus::BidPortfolio<cyclus::Material>::Ptr> Enrichment::GetMatlBids(
165  cyclus::CommodMap<cyclus::Material>::type& out_requests) {
166  using cyclus::Bid;
167  using cyclus::BidPortfolio;
168  using cyclus::CapacityConstraint;
169  using cyclus::Converter;
170  using cyclus::Material;
171  using cyclus::Request;
172  using cyclus::toolkit::MatVec;
173 
174  std::set<BidPortfolio<Material>::Ptr> ports;
175 
176  if ((out_requests.count(tails_commod) > 0) && (tails.quantity() > 0)) {
177  BidPortfolio<Material>::Ptr tails_port(new BidPortfolio<Material>());
178 
179  std::vector<Request<Material>*>& tails_requests =
180  out_requests[tails_commod];
181  std::vector<Request<Material>*>::iterator it;
182  for (it = tails_requests.begin(); it != tails_requests.end(); ++it) {
183  // offer bids for all tails material, keeping discrete quantities
184  // to preserve possible variation in composition
185  MatVec mats = tails.PopN(tails.count());
186  tails.Push(mats);
187  for (int k = 0; k < mats.size(); k++) {
188  Material::Ptr m = mats[k];
189  Request<Material>* req = *it;
190  tails_port->AddBid(req, m, this);
191  }
192  }
193  // overbidding (bidding on every offer)
194  // add an overall capacity constraint
195  CapacityConstraint<Material> tails_constraint(tails.quantity());
196  tails_port->AddConstraint(tails_constraint);
197  LOG(cyclus::LEV_INFO5, "EnrFac") << prototype()
198  << " adding tails capacity constraint of "
199  << tails.capacity();
200  ports.insert(tails_port);
201  }
202 
203  if ((out_requests.count(product_commod) > 0) && (inventory.quantity() > 0)) {
204  BidPortfolio<Material>::Ptr commod_port(new BidPortfolio<Material>());
205 
206  std::vector<Request<Material>*>& commod_requests =
207  out_requests[product_commod];
208  std::vector<Request<Material>*>::iterator it;
209  for (it = commod_requests.begin(); it != commod_requests.end(); ++it) {
210  Request<Material>* req = *it;
211  Material::Ptr mat = req->target();
212  double request_enrich = cyclus::toolkit::UraniumAssay(mat);
213  if (ValidReq(req->target()) &&
214  ((request_enrich < max_enrich) ||
215  (cyclus::AlmostEq(request_enrich, max_enrich)))) {
216  Material::Ptr offer = Offer_(req->target());
217  commod_port->AddBid(req, offer, this);
218  }
219  }
220 
221  Converter<Material>::Ptr sc(new SWUConverter(FeedAssay(), tails_assay));
222  Converter<Material>::Ptr nc(new NatUConverter(FeedAssay(), tails_assay));
223  CapacityConstraint<Material> swu(swu_capacity, sc);
224  CapacityConstraint<Material> natu(inventory.quantity(), nc);
225  commod_port->AddConstraint(swu);
226  commod_port->AddConstraint(natu);
227 
228  LOG(cyclus::LEV_INFO5, "EnrFac")
229  << prototype() << " adding a swu constraint of " << swu.capacity();
230  LOG(cyclus::LEV_INFO5, "EnrFac")
231  << prototype() << " adding a natu constraint of " << natu.capacity();
232  ports.insert(commod_port);
233  }
234  return ports;
235 }
236 
237 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238 bool Enrichment::ValidReq(const cyclus::Material::Ptr mat) {
239  cyclus::toolkit::MatQuery q(mat);
240  double u235 = q.atom_frac(922350000);
241  double u238 = q.atom_frac(922380000);
242  return (u238 > 0 && u235 / (u235 + u238) > tails_assay);
243 }
244 
245 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
247  const std::vector<cyclus::Trade<cyclus::Material> >& trades,
248  std::vector<std::pair<cyclus::Trade<cyclus::Material>,
249  cyclus::Material::Ptr> >& responses) {
250  using cyclus::Material;
251  using cyclus::Trade;
252 
255 
256  std::vector<Trade<Material> >::const_iterator it;
257  for (it = trades.begin(); it != trades.end(); ++it) {
258  double qty = it->amt;
259  std::string commod_type = it->bid->request()->commodity();
260  Material::Ptr response;
261 
262  // Figure out whether material is tails or enriched,
263  // if tails then make transfer of material
264  if (commod_type == tails_commod) {
265  LOG(cyclus::LEV_INFO5, "EnrFac")
266  << prototype() << " just received an order"
267  << " for " << it->amt << " of " << tails_commod;
268  double pop_qty = std::min(qty, tails.quantity());
269  response = tails.Pop(pop_qty, cyclus::eps_rsrc());
270  } else {
271  LOG(cyclus::LEV_INFO5, "EnrFac")
272  << prototype() << " just received an order"
273  << " for " << it->amt << " of " << product_commod;
274  response = Enrich_(it->bid->offer(), qty);
275  }
276  responses.push_back(std::make_pair(*it, response));
277  }
278 
279  if (cyclus::IsNegative(tails.quantity())) {
280  std::stringstream ss;
281  ss << "is being asked to provide more than its current inventory.";
282  throw cyclus::ValueError(Agent::InformErrorMsg(ss.str()));
283  }
284  if (cyclus::IsNegative(current_swu_capacity)) {
285  throw cyclus::ValueError("EnrFac " + prototype() +
286  " is being asked to provide more than" +
287  " its SWU capacity.");
288  }
289 }
290 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
291 void Enrichment::AddMat_(cyclus::Material::Ptr mat) {
292  // Elements and isotopes other than U-235, U-238 are sent directly to tails
293  cyclus::CompMap cm = mat->comp()->atom();
294  bool extra_u = false;
295  bool other_elem = false;
296  for (cyclus::CompMap::const_iterator it = cm.begin(); it != cm.end(); ++it) {
297  if (pyne::nucname::znum(it->first) == 92) {
298  if (pyne::nucname::anum(it->first) != 235 &&
299  pyne::nucname::anum(it->first) != 238 && it->second > 0) {
300  extra_u = true;
301  }
302  } else if (it->second > 0) {
303  other_elem = true;
304  }
305  }
306  if (extra_u) {
307  cyclus::Warn<cyclus::VALUE_WARNING>(
308  "More than 2 isotopes of U. "
309  "Istopes other than U-235, U-238 are sent directly to tails.");
310  }
311  if (other_elem) {
312  cyclus::Warn<cyclus::VALUE_WARNING>(
313  "Non-uranium elements are "
314  "sent directly to tails.");
315  }
316 
317  LOG(cyclus::LEV_INFO5, "EnrFac") << prototype() << " is initially holding "
318  << inventory.quantity() << " total.";
319 
320  try {
321  inventory.Push(mat);
322  } catch (cyclus::Error& e) {
323  e.msg(Agent::InformErrorMsg(e.msg()));
324  throw e;
325  }
326 
327  LOG(cyclus::LEV_INFO5, "EnrFac")
328  << prototype() << " added " << mat->quantity() << " of " << feed_commod
329  << " to its inventory, which is holding " << inventory.quantity()
330  << " total.";
331 }
332 
333 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
334 cyclus::Material::Ptr Enrichment::Request_() {
335  double qty = std::max(0.0, inventory.capacity() - inventory.quantity());
336  return cyclus::Material::CreateUntracked(qty,
337  context()->GetRecipe(feed_recipe));
338 }
339 
340 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
341 cyclus::Material::Ptr Enrichment::Offer_(cyclus::Material::Ptr mat) {
342  cyclus::toolkit::MatQuery q(mat);
343  cyclus::CompMap comp;
344  comp[922350000] = q.atom_frac(922350000);
345  comp[922380000] = q.atom_frac(922380000);
346  return cyclus::Material::CreateUntracked(
347  mat->quantity(), cyclus::Composition::CreateFromAtom(comp));
348 }
349 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
350 cyclus::Material::Ptr Enrichment::Enrich_(cyclus::Material::Ptr mat,
351  double qty) {
352  using cyclus::Material;
353  using cyclus::ResCast;
354  using cyclus::toolkit::Assays;
355  using cyclus::toolkit::UraniumAssay;
356  using cyclus::toolkit::SwuRequired;
357  using cyclus::toolkit::FeedQty;
358  using cyclus::toolkit::TailsQty;
359 
360  // get enrichment parameters
361  Assays assays(FeedAssay(), UraniumAssay(mat), tails_assay);
362  double swu_req = SwuRequired(qty, assays);
363  double natu_req = FeedQty(qty, assays);
364 
365  // Determine the composition of the natural uranium
366  // (ie. U-235+U-238/TotalMass)
367  double pop_qty = inventory.quantity();
368  Material::Ptr natu_matl = inventory.Pop(pop_qty, cyclus::eps_rsrc());
369  inventory.Push(natu_matl);
370 
371  cyclus::toolkit::MatQuery mq(natu_matl);
372  std::set<cyclus::Nuc> nucs;
373  nucs.insert(922350000);
374  nucs.insert(922380000);
375  double natu_frac = mq.mass_frac(nucs);
376  double feed_req = natu_req / natu_frac;
377 
378  // pop amount from inventory and blob it into one material
379  Material::Ptr r;
380  try {
381  // required so popping doesn't take out too much
382  if (cyclus::AlmostEq(feed_req, inventory.quantity())) {
383  r = cyclus::toolkit::Squash(inventory.PopN(inventory.count()));
384  } else {
385  r = inventory.Pop(feed_req, cyclus::eps_rsrc());
386  }
387  } catch (cyclus::Error& e) {
388  NatUConverter nc(FeedAssay(), tails_assay);
389  std::stringstream ss;
390  ss << " tried to remove " << feed_req << " from its inventory of size "
391  << inventory.quantity()
392  << " and the conversion of the material into natu is "
393  << nc.convert(mat);
394  throw cyclus::ValueError(Agent::InformErrorMsg(ss.str()));
395  }
396 
397  // "enrich" it, but pull out the composition and quantity we require from the
398  // blob
399  cyclus::Composition::Ptr comp = mat->comp();
400  Material::Ptr response = r->ExtractComp(qty, comp);
401  tails.Push(r);
402 
403  current_swu_capacity -= swu_req;
404 
405  intra_timestep_swu_ += swu_req;
406  intra_timestep_feed_ += feed_req;
407  RecordEnrichment_(feed_req, swu_req);
408 
409  LOG(cyclus::LEV_INFO5, "EnrFac") << prototype()
410  << " has performed an enrichment: ";
411  LOG(cyclus::LEV_INFO5, "EnrFac") << " * Feed Qty: " << feed_req;
412  LOG(cyclus::LEV_INFO5, "EnrFac") << " * Feed Assay: "
413  << assays.Feed() * 100;
414  LOG(cyclus::LEV_INFO5, "EnrFac") << " * Product Qty: " << qty;
415  LOG(cyclus::LEV_INFO5, "EnrFac") << " * Product Assay: "
416  << assays.Product() * 100;
417  LOG(cyclus::LEV_INFO5, "EnrFac") << " * Tails Qty: "
418  << TailsQty(qty, assays);
419  LOG(cyclus::LEV_INFO5, "EnrFac") << " * Tails Assay: "
420  << assays.Tails() * 100;
421  LOG(cyclus::LEV_INFO5, "EnrFac") << " * SWU: " << swu_req;
422  LOG(cyclus::LEV_INFO5, "EnrFac") << " * Current SWU capacity: "
424 
425  return response;
426 }
427 
428 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
429 void Enrichment::RecordEnrichment_(double natural_u, double swu) {
430  using cyclus::Context;
431  using cyclus::Agent;
432 
433  LOG(cyclus::LEV_DEBUG1, "EnrFac") << prototype()
434  << " has enriched a material:";
435  LOG(cyclus::LEV_DEBUG1, "EnrFac") << " * Amount: " << natural_u;
436  LOG(cyclus::LEV_DEBUG1, "EnrFac") << " * SWU: " << swu;
437 
438  Context* ctx = Agent::context();
439  ctx->NewDatum("Enrichments")
440  ->AddVal("ID", id())
441  ->AddVal("Time", ctx->time())
442  ->AddVal("Natural_Uranium", natural_u)
443  ->AddVal("SWU", swu)
444  ->Record();
445 }
446 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
447 double Enrichment::FeedAssay() {
448  using cyclus::Material;
449 
450  if (inventory.empty()) {
451  return 0;
452  }
453  double pop_qty = inventory.quantity();
454  cyclus::Material::Ptr fission_matl =
455  inventory.Pop(pop_qty, cyclus::eps_rsrc());
456  inventory.Push(fission_matl);
457  return cyclus::toolkit::UraniumAssay(fission_matl);
458 }
459 
460 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
461 extern "C" cyclus::Agent* ConstructEnrichment(cyclus::Context* ctx) {
462  return new Enrichment(ctx);
463 }
464 
465 } // namespace cycamore
double FeedAssay()
calculates the feed assay based on the unenriched inventory
cyclus::Material::Ptr Offer_(cyclus::Material::Ptr req)
Generates a material offer for a given request.
virtual ~Enrichment()
Destructor for the Enrichment class.
virtual void Build(cyclus::Agent *parent)
perform module-specific tasks when entering the simulation
virtual void Tick()
Each facility is prompted to do its beginning-of-time-step stuff at the tick of the timer...
bool SortBids(cyclus::Bid< cyclus::Material > *i, cyclus::Bid< cyclus::Material > *j)
cyclus::Agent * ConstructEnrichment(cyclus::Context *ctx)
virtual void AdjustMatlPrefs(cyclus::PrefMap< cyclus::Material >::type &prefs)
The Enrichment adjusts preferences for offers of natural uranium it has received to maximize U-235 co...
virtual void AcceptMatlTrades(const std::vector< std::pair< cyclus::Trade< cyclus::Material >, cyclus::Material::Ptr > > &responses)
The Enrichment place accepted trade Materials in their Inventory.
bool ValidReq(const cyclus::Material::Ptr mat)
Determines if a particular material is a valid request to respond to.
cyclus::Material::Ptr Request_()
generates a request for this facility given its current state.
cyclus::Material::Ptr Enrich_(cyclus::Material::Ptr mat, double qty)
cycamore::GrowthRegion string
virtual std::set< cyclus::BidPortfolio< cyclus::Material >::Ptr > GetMatlBids(cyclus::CommodMap< cyclus::Material >::type &commod_requests)
Responds to each request for this facility&#39;s commodity.
virtual std::string str()
Print information about this agent.
Enrichment(cyclus::Context *ctx)
Constructor for the Enrichment class.
void AddMat_(cyclus::Material::Ptr mat)
adds a material into the natural uranium inventory
virtual void Tock()
Each facility is prompted to its end-of-time-step stuff on the tock of the timer. ...
virtual std::set< cyclus::RequestPortfolio< cyclus::Material >::Ptr > GetMatlRequests()
The Enrichment request Materials of its given commodity.
virtual void GetMatlTrades(const std::vector< cyclus::Trade< cyclus::Material > > &trades, std::vector< std::pair< cyclus::Trade< cyclus::Material >, cyclus::Material::Ptr > > &responses)
respond to each trade with a material enriched to the appropriate level given this facility&#39;s invento...
cyclus::toolkit::ResBuf< cyclus::Material > tails
void RecordEnrichment_(double natural_u, double swu)
records and enrichment with the cyclus::Recorder
cyclus::toolkit::ResBuf< cyclus::Material > inventory