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