CYCAMORE
src/fuel_fab.cc
Go to the documentation of this file.
1 #include "fuel_fab.h"
2 
3 #include <sstream>
4 
5 using cyclus::Material;
6 using cyclus::Composition;
7 using pyne::simple_xs;
8 
9 #define SHOW(X) \
10  std::cout << std::setprecision(17) << __FILE__ << ":" << __LINE__ \
11  << ": " #X " = " << X << "\n"
12 
13 namespace cycamore {
14 
15 class FissConverter : public cyclus::Converter<cyclus::Material> {
16  public:
17  FissConverter(Composition::Ptr c_fill, Composition::Ptr c_fiss,
18  Composition::Ptr c_topup, std::string spectrum)
19  : c_fiss_(c_fiss), c_topup_(c_topup), c_fill_(c_fill), spec_(spectrum) {
20  w_fiss_ = CosiWeight(c_fiss, spectrum);
21  w_fill_ = CosiWeight(c_fill, spectrum);
22  w_topup_ = CosiWeight(c_topup, spectrum);
23  }
24 
25  virtual ~FissConverter() {}
26 
27  virtual double convert(
28  cyclus::Material::Ptr m, cyclus::Arc const* a = NULL,
29  cyclus::ExchangeTranslationContext<cyclus::Material> const* ctx =
30  NULL) const {
31  double w_tgt = CosiWeight(m->comp(), spec_);
32  if (ValidWeights(w_fill_, w_tgt, w_fiss_)) {
33  double frac = HighFrac(w_fill_, w_tgt, w_fiss_);
34  return AtomToMassFrac(frac, c_fiss_, c_fill_) * m->quantity();
35  } else if (ValidWeights(w_fiss_, w_tgt, w_topup_)) {
36  // use fiss inventory as filler, and topup as fissile
37  double frac = LowFrac(w_fiss_, w_tgt, w_topup_);
38  return AtomToMassFrac(frac, c_fiss_, c_topup_) * m->quantity();
39  } else {
40  // don't bid at all
41  return 1e200;
42  }
43  }
44 
45  private:
47  double w_fiss_;
48  double w_topup_;
49  double w_fill_;
50  Composition::Ptr c_fiss_;
51  Composition::Ptr c_fill_;
52  Composition::Ptr c_topup_;
53 };
54 
55 class FillConverter : public cyclus::Converter<cyclus::Material> {
56  public:
57  FillConverter(Composition::Ptr c_fill, Composition::Ptr c_fiss,
58  Composition::Ptr c_topup, std::string spectrum)
59  : c_fiss_(c_fiss), c_topup_(c_topup), c_fill_(c_fill), spec_(spectrum) {
60  w_fiss_ = CosiWeight(c_fiss, spectrum);
61  w_fill_ = CosiWeight(c_fill, spectrum);
62  w_topup_ = CosiWeight(c_topup, spectrum);
63  }
64 
65  virtual ~FillConverter() {}
66 
67  virtual double convert(
68  cyclus::Material::Ptr m, cyclus::Arc const* a = NULL,
69  cyclus::ExchangeTranslationContext<cyclus::Material> const* ctx =
70  NULL) const {
71  double w_tgt = CosiWeight(m->comp(), spec_);
72  if (ValidWeights(w_fill_, w_tgt, w_fiss_)) {
73  double frac = LowFrac(w_fill_, w_tgt, w_fiss_);
74  return AtomToMassFrac(frac, c_fill_, c_fiss_) * m->quantity();
75  } else if (ValidWeights(w_fiss_, w_tgt, w_topup_)) {
76  // switched fissile inventory to filler so don't need any filler inventory
77  return 0;
78  } else {
79  // don't bid at all
80  return 1e200;
81  }
82  }
83 
84  private:
86  double w_fiss_;
87  double w_topup_;
88  double w_fill_;
89  Composition::Ptr c_fiss_;
90  Composition::Ptr c_fill_;
91  Composition::Ptr c_topup_;
92 };
93 
94 class TopupConverter : public cyclus::Converter<cyclus::Material> {
95  public:
96  TopupConverter(Composition::Ptr c_fill, Composition::Ptr c_fiss,
97  Composition::Ptr c_topup, std::string spectrum)
98  : c_fiss_(c_fiss), c_topup_(c_topup), c_fill_(c_fill), spec_(spectrum) {
99  w_fiss_ = CosiWeight(c_fiss, spectrum);
100  w_fill_ = CosiWeight(c_fill, spectrum);
101  w_topup_ = CosiWeight(c_topup, spectrum);
102  }
103 
104  virtual ~TopupConverter() {}
105 
106  virtual double convert(
107  cyclus::Material::Ptr m, cyclus::Arc const* a = NULL,
108  cyclus::ExchangeTranslationContext<cyclus::Material> const* ctx =
109  NULL) const {
110  double w_tgt = CosiWeight(m->comp(), spec_);
111  if (ValidWeights(w_fill_, w_tgt, w_fiss_)) {
112  return 0;
113  } else if (ValidWeights(w_fiss_, w_tgt, w_topup_)) {
114  // switched fissile inventory to filler and topup as fissile
115  double frac = HighFrac(w_fiss_, w_tgt, w_topup_);
116  return AtomToMassFrac(frac, c_topup_, c_fiss_) * m->quantity();
117  } else {
118  // don't bid at all
119  return 1e200;
120  }
121  }
122 
123  private:
125  double w_fiss_;
126  double w_topup_;
127  double w_fill_;
128  Composition::Ptr c_fiss_;
129  Composition::Ptr c_fill_;
130  Composition::Ptr c_topup_;
131 };
132 
133 FuelFab::FuelFab(cyclus::Context* ctx)
134  : cyclus::Facility(ctx), fill_size(0), fiss_size(0), throughput(0) {}
135 
136 void FuelFab::EnterNotify() {
137  cyclus::Facility::EnterNotify();
138 
139  if (fiss_commod_prefs.empty()) {
140  for (int i = 0; i < fiss_commods.size(); i++) {
141  fiss_commod_prefs.push_back(cyclus::kDefaultPref);
142  }
143  } else if (fiss_commod_prefs.size() != fiss_commods.size()) {
144  std::stringstream ss;
145  ss << "prototype '" << prototype() << "' has " << fiss_commod_prefs.size()
146  << " fiss_commod_prefs vals, expected " << fiss_commods.size();
147  throw cyclus::ValidationError(ss.str());
148  }
149 
150  if (fill_commod_prefs.empty()) {
151  for (int i = 0; i < fill_commods.size(); i++) {
152  fill_commod_prefs.push_back(cyclus::kDefaultPref);
153  }
154  } else if (fill_commod_prefs.size() != fill_commods.size()) {
155  std::stringstream ss;
156  ss << "prototype '" << prototype() << "' has " << fill_commod_prefs.size()
157  << " fill_commod_prefs vals, expected " << fill_commods.size();
158  throw cyclus::ValidationError(ss.str());
159  }
160 }
161 
162 std::set<cyclus::RequestPortfolio<Material>::Ptr> FuelFab::GetMatlRequests() {
163  using cyclus::RequestPortfolio;
164 
165  std::set<RequestPortfolio<Material>::Ptr> ports;
166 
167  bool exclusive = false;
168 
169  if (fiss.space() > cyclus::eps_rsrc()) {
170  RequestPortfolio<Material>::Ptr port(new RequestPortfolio<Material>());
171 
172  Material::Ptr m = cyclus::NewBlankMaterial(fiss.space());
173  if (!fiss_recipe.empty()) {
174  Composition::Ptr c = context()->GetRecipe(fiss_recipe);
175  m = Material::CreateUntracked(fiss.space(), c);
176  }
177 
178  std::vector<cyclus::Request<Material>*> reqs;
179  for (int i = 0; i < fiss_commods.size(); i++) {
180  std::string commod = fiss_commods[i];
181  double pref = fiss_commod_prefs[i];
182  reqs.push_back(port->AddRequest(m, this, commod, pref, exclusive));
183  req_inventories_[reqs.back()] = "fiss";
184  }
185  port->AddMutualReqs(reqs);
186  ports.insert(port);
187  }
188 
189  if (fill.space() > cyclus::eps_rsrc()) {
190  RequestPortfolio<Material>::Ptr port(new RequestPortfolio<Material>());
191 
192  Material::Ptr m = cyclus::NewBlankMaterial(fill.space());
193  if (!fill_recipe.empty()) {
194  Composition::Ptr c = context()->GetRecipe(fill_recipe);
195  m = Material::CreateUntracked(fill.space(), c);
196  }
197 
198  std::vector<cyclus::Request<Material>*> reqs;
199  for (int i = 0; i < fill_commods.size(); i++) {
200  std::string commod = fill_commods[i];
201  double pref = fill_commod_prefs[i];
202  reqs.push_back(port->AddRequest(m, this, commod, pref, exclusive));
203  req_inventories_[reqs.back()] = "fill";
204  }
205  port->AddMutualReqs(reqs);
206  ports.insert(port);
207  }
208 
209  if (topup.space() > cyclus::eps_rsrc()) {
210  RequestPortfolio<Material>::Ptr port(new RequestPortfolio<Material>());
211 
212  Material::Ptr m = cyclus::NewBlankMaterial(topup.space());
213  if (!topup_recipe.empty()) {
214  Composition::Ptr c = context()->GetRecipe(topup_recipe);
215  m = Material::CreateUntracked(topup.space(), c);
216  }
217  cyclus::Request<Material>* r =
218  port->AddRequest(m, this, topup_commod, topup_pref, exclusive);
219  req_inventories_[r] = "topup";
220  ports.insert(port);
221  }
222 
223  return ports;
224 }
225 
226 bool Contains(std::vector<std::string> vec, std::string s) {
227  for (int i = 0; i < vec.size(); i++) {
228  if (vec[i] == s) {
229  return true;
230  }
231  }
232  return false;
233 }
234 
236  const std::vector<std::pair<cyclus::Trade<Material>, Material::Ptr> >&
237  responses) {
238  std::vector<std::pair<cyclus::Trade<cyclus::Material>,
239  cyclus::Material::Ptr> >::const_iterator trade;
240 
241  for (trade = responses.begin(); trade != responses.end(); ++trade) {
242  std::string commod = trade->first.request->commodity();
243  double req_qty = trade->first.request->target()->quantity();
244  cyclus::Request<Material>* req = trade->first.request;
245  Material::Ptr m = trade->second;
246  if (req_inventories_[req] == "fill") {
247  fill.Push(m);
248  } else if (req_inventories_[req] == "topup") {
249  topup.Push(m);
250  } else if (req_inventories_[req] == "fiss") {
251  fiss.Push(m);
252  } else {
253  throw cyclus::ValueError("cycamore::FuelFab was overmatched on requests");
254  }
255  }
256 
257  req_inventories_.clear();
258 
259  // IMPORTANT - each buffer needs to be a single homogenous composition or
260  // the inventory mixing constraints for bids don't work
261  if (fill.count() > 1) {
262  fill.Push(cyclus::toolkit::Squash(fill.PopN(fill.count())));
263  }
264  if (fiss.count() > 1) {
265  fiss.Push(cyclus::toolkit::Squash(fiss.PopN(fiss.count())));
266  }
267  if (topup.count() > 1) {
268  topup.Push(cyclus::toolkit::Squash(topup.PopN(topup.count())));
269  }
270 }
271 
272 std::set<cyclus::BidPortfolio<Material>::Ptr> FuelFab::GetMatlBids(
273  cyclus::CommodMap<Material>::type& commod_requests) {
274  using cyclus::BidPortfolio;
275 
276  std::set<BidPortfolio<Material>::Ptr> ports;
277  std::vector<cyclus::Request<Material>*>& reqs = commod_requests[outcommod];
278 
279  if (throughput == 0) {
280  return ports;
281  } else if (reqs.size() == 0) {
282  return ports;
283  }
284 
285  double w_fill = 0;
286  Composition::Ptr
287  c_fill; // no default needed - this is non-optional parameter
288  if (fill.count() > 0) {
289  c_fill = fill.Peek()->comp();
290  w_fill = CosiWeight(c_fill, spectrum);
291  } else {
292  c_fill = context()->GetRecipe(fill_recipe);
293  w_fill = CosiWeight(c_fill, spectrum);
294  }
295 
296  double w_topup = 0;
297  Composition::Ptr c_topup = c_fill;
298  if (topup.count() > 0) {
299  c_topup = topup.Peek()->comp();
300  w_topup = CosiWeight(c_topup, spectrum);
301  } else if (!topup_recipe.empty()) {
302  c_topup = context()->GetRecipe(topup_recipe);
303  w_topup = CosiWeight(c_topup, spectrum);
304  }
305 
306  double w_fiss =
307  w_fill; // this allows trading just fill with no fiss inventory
308  Composition::Ptr c_fiss = c_fill;
309  if (fiss.count() > 0) {
310  c_fiss = fiss.Peek()->comp();
311  w_fiss = CosiWeight(c_fiss, spectrum);
312  } else if (!fiss_recipe.empty()) {
313  c_fiss = context()->GetRecipe(fiss_recipe);
314  w_fiss = CosiWeight(c_fiss, spectrum);
315  }
316 
317  BidPortfolio<Material>::Ptr port(new BidPortfolio<Material>());
318  for (int j = 0; j < reqs.size(); j++) {
319  cyclus::Request<Material>* req = reqs[j];
320 
321  Composition::Ptr tgt = req->target()->comp();
322  double w_tgt = CosiWeight(tgt, spectrum);
323  double tgt_qty = req->target()->quantity();
324  if (ValidWeights(w_fill, w_tgt, w_fiss)) {
325  double fiss_frac = HighFrac(w_fill, w_tgt, w_fiss);
326  double fill_frac = 1 - fiss_frac;
327  fiss_frac = AtomToMassFrac(fiss_frac, c_fiss, c_fill);
328  fill_frac = AtomToMassFrac(fill_frac, c_fill, c_fiss);
329  Material::Ptr m1 = Material::CreateUntracked(fiss_frac * tgt_qty, c_fiss);
330  Material::Ptr m2 = Material::CreateUntracked(fill_frac * tgt_qty, c_fill);
331  m1->Absorb(m2);
332 
333  bool exclusive = false;
334  port->AddBid(req, m1, this, exclusive);
335  } else if (topup.count() > 0 && ValidWeights(w_fiss, w_tgt, w_topup)) {
336  // only bid with topup if we have filler - otherwise we might be able to
337  // meet target with filler when we get it. we should only use topup
338  // when the fissile has too poor neutronics.
339  double topup_frac = HighFrac(w_fiss, w_tgt, w_topup);
340  double fiss_frac = 1 - topup_frac;
341  fiss_frac = AtomToMassFrac(fiss_frac, c_fiss, c_topup);
342  topup_frac = AtomToMassFrac(topup_frac, c_topup, c_fiss);
343  Material::Ptr m1 =
344  Material::CreateUntracked(topup_frac * tgt_qty, c_topup);
345  Material::Ptr m2 = Material::CreateUntracked(fiss_frac * tgt_qty, c_fiss);
346  m1->Absorb(m2);
347 
348  bool exclusive = false;
349  port->AddBid(req, m1, this, exclusive);
350  } else if (fiss.count() > 0 && fill.count() > 0 ||
351  fiss.count() > 0 && topup.count() > 0) {
352  // else can't meet the target weight - don't bid. Just a plain else
353  // doesn't work because we set w_fiss = w_fill if we don't have any fiss
354  // or fill inventory.
355  std::stringstream ss;
356  ss << "prototype '" << prototype()
357  << "': Input stream weights/reactivity do not span "
358  "the requested material weight.";
359  cyclus::Warn<cyclus::VALUE_WARNING>(ss.str());
360  }
361  }
362 
363  cyclus::Converter<Material>::Ptr fissconv(
364  new FissConverter(c_fill, c_fiss, c_topup, spectrum));
365  cyclus::Converter<Material>::Ptr fillconv(
366  new FillConverter(c_fill, c_fiss, c_topup, spectrum));
367  cyclus::Converter<Material>::Ptr topupconv(
368  new TopupConverter(c_fill, c_fiss, c_topup, spectrum));
369  // important! - the std::max calls prevent CapacityConstraint throwing a zero
370  // cap exception
371  cyclus::CapacityConstraint<Material> fissc(std::max(fiss.quantity(), 1e-10),
372  fissconv);
373  cyclus::CapacityConstraint<Material> fillc(std::max(fill.quantity(), 1e-10),
374  fillconv);
375  cyclus::CapacityConstraint<Material> topupc(std::max(topup.quantity(), 1e-10),
376  topupconv);
377  port->AddConstraint(fillc);
378  port->AddConstraint(fissc);
379  port->AddConstraint(topupc);
380 
381  cyclus::CapacityConstraint<Material> cc(throughput);
382  port->AddConstraint(cc);
383  ports.insert(port);
384  return ports;
385 }
386 
388  const std::vector<cyclus::Trade<Material> >& trades,
389  std::vector<std::pair<cyclus::Trade<Material>, Material::Ptr> >&
390  responses) {
391  using cyclus::Trade;
392 
393  // guard against cases where a buffer is empty - this is okay because some
394  // trades may not need that particular buffer.
395  double w_fill = 0;
396  if (fill.count() > 0) {
397  w_fill = CosiWeight(fill.Peek()->comp(), spectrum);
398  }
399  double w_topup = 0;
400  if (topup.count() > 0) {
401  w_topup = CosiWeight(topup.Peek()->comp(), spectrum);
402  }
403  double w_fiss = 0;
404  if (fiss.count() > 0) {
405  w_fiss = CosiWeight(fiss.Peek()->comp(), spectrum);
406  }
407 
408  std::vector<cyclus::Trade<cyclus::Material> >::const_iterator it;
409  double tot = 0;
410  for (int i = 0; i < trades.size(); i++) {
411  Material::Ptr tgt = trades[i].request->target();
412 
413  double w_tgt = CosiWeight(tgt->comp(), spectrum);
414  double qty = trades[i].amt;
415  double wfiss = w_fiss;
416 
417  tot += qty;
418  if (tot > throughput + cyclus::eps_rsrc()) {
419  std::stringstream ss;
420  ss << "FuelFab was matched above throughput limit: " << tot << " > "
421  << throughput;
422  throw cyclus::ValueError(ss.str());
423  }
424 
425  if (fiss.count() == 0) {
426  // use straight filler to satisfy this request
427  double fillqty = qty;
428  if (std::abs(fillqty - fill.quantity()) < cyclus::eps_rsrc()) {
429  fillqty = std::min(fill.quantity(), qty);
430  }
431  responses.push_back(
432  std::make_pair(trades[i], fill.Pop(fillqty, cyclus::eps_rsrc())));
433  } else if (fill.count() == 0 && ValidWeights(w_fill, w_tgt, w_fiss)) {
434  // use straight fissile to satisfy this request
435  double fissqty = qty;
436  if (std::abs(fissqty - fiss.quantity()) < cyclus::eps_rsrc()) {
437  fissqty = std::min(fiss.quantity(), qty);
438  }
439  responses.push_back(
440  std::make_pair(trades[i], fiss.Pop(fissqty, cyclus::eps_rsrc())));
441  } else if (ValidWeights(w_fill, w_tgt, w_fiss)) {
442  double fiss_frac = HighFrac(w_fill, w_tgt, w_fiss);
443  double fill_frac = LowFrac(w_fill, w_tgt, w_fiss);
444  fiss_frac =
445  AtomToMassFrac(fiss_frac, fiss.Peek()->comp(), fill.Peek()->comp());
446  fill_frac =
447  AtomToMassFrac(fill_frac, fill.Peek()->comp(), fiss.Peek()->comp());
448 
449  double fissqty = fiss_frac * qty;
450  if (std::abs(fissqty - fiss.quantity()) < cyclus::eps_rsrc()) {
451  fissqty = std::min(fiss.quantity(), fiss_frac * qty);
452  }
453  double fillqty = fill_frac * qty;
454  if (std::abs(fillqty - fill.quantity()) < cyclus::eps_rsrc()) {
455  fillqty = std::min(fill.quantity(), fill_frac * qty);
456  }
457 
458  Material::Ptr m = fiss.Pop(fissqty, cyclus::eps_rsrc());
459  // this if block prevents zero qty ResBuf pop exceptions
460  if (fill_frac > 0) {
461  m->Absorb(fill.Pop(fillqty, cyclus::eps_rsrc()));
462  }
463  responses.push_back(std::make_pair(trades[i], m));
464  } else {
465  double topup_frac = HighFrac(w_fiss, w_tgt, w_topup);
466  double fiss_frac = 1 - topup_frac;
467  topup_frac =
468  AtomToMassFrac(topup_frac, topup.Peek()->comp(), fiss.Peek()->comp());
469  fiss_frac =
470  AtomToMassFrac(fiss_frac, fiss.Peek()->comp(), topup.Peek()->comp());
471 
472  double fissqty = fiss_frac * qty;
473  if (std::abs(fissqty - fiss.quantity()) < cyclus::eps_rsrc()) {
474  fissqty = std::min(fiss.quantity(), fiss_frac * qty);
475  }
476  double topupqty = topup_frac * qty;
477  if (std::abs(topupqty - topup.quantity()) < cyclus::eps_rsrc()) {
478  topupqty = std::min(topup.quantity(), topup_frac * qty);
479  }
480 
481  Material::Ptr m = fiss.Pop(fissqty, cyclus::eps_rsrc());
482  // this if block prevents zero qty ResBuf pop exceptions
483  if (topup_frac > 0) {
484  m->Absorb(topup.Pop(topupqty, cyclus::eps_rsrc()));
485  }
486  responses.push_back(std::make_pair(trades[i], m));
487  }
488  }
489 }
490 
491 extern "C" cyclus::Agent* ConstructFuelFab(cyclus::Context* ctx) {
492  return new FuelFab(ctx);
493 }
494 
495 // Returns the weight of c using 1 group cross sections of type spectrum
496 // which must be one of:
497 //
498 // * thermal
499 // * thermal_maxwell_ave
500 // * fission_spectrum_ave
501 // * resonance_integral
502 // * fourteen_MeV
503 //
504 // The weight is calculated as "(nu*sigma_f - sigma_a) * N". Since weights
505 // are computed based on nuclide atom fractions, corresponding computed
506 // material/mixing fractions will also be atom-based naturally and will need
507 // to be converted to mass-based for actual material object mixing.
508 double CosiWeight(cyclus::Composition::Ptr c, const std::string& spectrum) {
509  cyclus::CompMap cm = c->atom();
510  cyclus::compmath::Normalize(&cm);
511 
512  if (spectrum == "thermal") {
513  double nu_pu239 = 2.85;
514  double nu_u233 = 2.5;
515  double nu_u235 = 2.43;
516  double nu_u238 = 0;
517  double nu_pu241 = nu_pu239;
518 
519  static std::map<int, double> absorb_xs;
520  static std::map<int, double> fiss_xs;
521  static double p_u238 = 0;
522  static double p_pu239 = 0;
523  if (p_u238 == 0) {
524  double fiss_u238 = simple_xs(922380000, "fission", "thermal");
525  double absorb_u238 = simple_xs(922380000, "absorption", "thermal");
526  p_u238 = nu_u238 * fiss_u238 - absorb_u238;
527 
528  double fiss_pu239 = simple_xs(942390000, "fission", "thermal");
529  double absorb_pu239 = simple_xs(942390000, "absorption", "thermal");
530  p_pu239 = nu_pu239 * fiss_pu239 - absorb_pu239;
531  }
532 
533  cyclus::CompMap::iterator it;
534  double w = 0;
535  for (it = cm.begin(); it != cm.end(); ++it) {
536  cyclus::Nuc nuc = it->first;
537  double nu = 0;
538  if (nuc == 922350000) {
539  nu = nu_u235;
540  } else if (nuc == 922330000) {
541  nu = nu_u233;
542  } else if (nuc == 942390000) {
543  nu = nu_pu239;
544  } else if (nuc == 942410000) {
545  nu = nu_pu241;
546  }
547 
548  double fiss = 0;
549  double absorb = 0;
550  if (absorb_xs.count(nuc) == 0) {
551  try {
552  fiss = simple_xs(nuc, "fission", "thermal");
553  absorb = simple_xs(nuc, "absorption", "thermal");
554  absorb_xs[nuc] = absorb;
555  fiss_xs[nuc] = fiss;
556  } catch (pyne::InvalidSimpleXS err) {
557  fiss = 0;
558  absorb = 0;
559  }
560  } else {
561  fiss = fiss_xs[nuc];
562  absorb = absorb_xs[nuc];
563  }
564 
565  double p = nu * fiss - absorb;
566  w += it->second * (p - p_u238) / (p_pu239 - p_u238);
567  }
568  return w;
569  } else if (spectrum == "fission_spectrum_ave") {
570  double nu_pu239 = 3.1;
571  double nu_u233 = 2.63;
572  double nu_u235 = 2.58;
573  double nu_u238 = 0;
574  double nu_pu241 = nu_pu239;
575 
576  static std::map<int, double> absorb_xs;
577  static std::map<int, double> fiss_xs;
578  static double p_u238 = 0;
579  static double p_pu239 = 0;
580  if (p_u238 == 0) {
581  double fiss_u238 =
582  simple_xs(922380000, "fission", "fission_spectrum_ave");
583  double absorb_u238 =
584  simple_xs(922380000, "absorption", "fission_spectrum_ave");
585  p_u238 = nu_u238 * fiss_u238 - absorb_u238;
586 
587  double fiss_pu239 =
588  simple_xs(942390000, "fission", "fission_spectrum_ave");
589  double absorb_pu239 =
590  simple_xs(942390000, "absorption", "fission_spectrum_ave");
591  p_pu239 = nu_pu239 * fiss_pu239 - absorb_pu239;
592  }
593 
594  cyclus::CompMap::iterator it;
595  double w = 0;
596  for (it = cm.begin(); it != cm.end(); ++it) {
597  cyclus::Nuc nuc = it->first;
598  double nu = 0;
599  if (nuc == 922350000) {
600  nu = nu_u235;
601  } else if (nuc == 922330000) {
602  nu = nu_u233;
603  } else if (nuc == 942390000) {
604  nu = nu_pu239;
605  } else if (nuc == 942410000) {
606  nu = nu_pu241;
607  }
608 
609  double fiss = 0;
610  double absorb = 0;
611  if (absorb_xs.count(nuc) == 0) {
612  try {
613  fiss = simple_xs(nuc, "fission", "fission_spectrum_ave");
614  absorb = simple_xs(nuc, "absorption", "fission_spectrum_ave");
615  absorb_xs[nuc] = absorb;
616  fiss_xs[nuc] = fiss;
617  } catch (pyne::InvalidSimpleXS err) {
618  fiss = 0;
619  absorb = 0;
620  }
621  } else {
622  fiss = fiss_xs[nuc];
623  absorb = absorb_xs[nuc];
624  }
625 
626  double p = nu * fiss - absorb;
627  w += it->second * (p - p_u238) / (p_pu239 - p_u238);
628  }
629  return w;
630  } else {
631  double nu_pu239 = 3.1;
632  double nu_u233 = 2.63;
633  double nu_u235 = 2.58;
634  double nu_u238 = 0;
635  double nu_pu241 = nu_pu239;
636 
637  double fiss_u238 = simple_xs(922380000, "fission", spectrum);
638  double absorb_u238 = simple_xs(922380000, "absorption", spectrum);
639  double p_u238 = nu_u238 * fiss_u238 - absorb_u238;
640 
641  double fiss_pu239 = simple_xs(942390000, "fission", spectrum);
642  double absorb_pu239 = simple_xs(942390000, "absorption", spectrum);
643  double p_pu239 = nu_pu239 * fiss_pu239 - absorb_pu239;
644 
645  cyclus::CompMap::iterator it;
646  double w = 0;
647  for (it = cm.begin(); it != cm.end(); ++it) {
648  cyclus::Nuc nuc = it->first;
649  double nu = 0;
650  if (nuc == 922350000) {
651  nu = nu_u235;
652  } else if (nuc == 922330000) {
653  nu = nu_u233;
654  } else if (nuc == 942390000) {
655  nu = nu_pu239;
656  } else if (nuc == 942410000) {
657  nu = nu_pu241;
658  }
659 
660  double fiss = 0;
661  double absorb = 0;
662  try {
663  fiss = simple_xs(nuc, "fission", spectrum);
664  absorb = simple_xs(nuc, "absorption", spectrum);
665  } catch (pyne::InvalidSimpleXS err) {
666  fiss = 0;
667  absorb = 0;
668  }
669 
670  double p = nu * fiss - absorb;
671  w += it->second * (p - p_u238) / (p_pu239 - p_u238);
672  }
673  return w;
674  }
675 }
676 
677 // Convert an atom frac (n1/(n1+n2) to a mass frac (m1/(m1+m2) given
678 // corresponding compositions c1 and c2.
679 double AtomToMassFrac(double atomfrac, Composition::Ptr c1,
680  Composition::Ptr c2) {
681  cyclus::CompMap n1 = c1->atom();
682  cyclus::CompMap n2 = c2->atom();
683  cyclus::compmath::Normalize(&n1, atomfrac);
684  cyclus::compmath::Normalize(&n2, 1 - atomfrac);
685 
686  cyclus::CompMap::iterator it;
687 
688  double mass1 = 0;
689  for (it = n1.begin(); it != n1.end(); ++it) {
690  mass1 += it->second * pyne::atomic_mass(it->first);
691  }
692 
693  double mass2 = 0;
694  for (it = n2.begin(); it != n2.end(); ++it) {
695  mass2 += it->second * pyne::atomic_mass(it->first);
696  }
697 
698  return mass1 / (mass1 + mass2);
699 }
700 
701 double HighFrac(double w_low, double w_target, double w_high, double eps) {
702  if (!ValidWeights(w_low, w_target, w_high)) {
703  throw cyclus::ValueError("low and high weights cannot meet target");
704  } else if (w_low == w_high && w_target == w_low) {
705  return 1;
706  }
707  double f = std::abs((w_target - w_low) / (w_high - w_low));
708  if (1 - f < eps) {
709  return 1;
710  } else if (f < eps) {
711  return 0;
712  }
713  return f;
714 }
715 
716 double LowFrac(double w_low, double w_target, double w_high, double eps) {
717  return 1 - HighFrac(w_low, w_target, w_high, eps);
718 }
719 
720 // Returns true if the given weights can be used to linearly interpolate valid
721 // mixing fractions of the filler and fissile streams to hit the target.
722 bool ValidWeights(double w_low, double w_target, double w_high) {
723  // w_tgt must be in between w_fill and w_fiss for the equivalence
724  // interpolation to work.
725  return w_low <= w_target && w_target <= w_high;
726 }
727 
728 } // namespace cycamore
FissConverter(Composition::Ptr c_fill, Composition::Ptr c_fiss, Composition::Ptr c_topup, std::string spectrum)
Definition: src/fuel_fab.cc:17
double HighFrac(double w_low, double w_target, double w_high, double eps)
bool Contains(std::vector< std::string > vec, std::string s)
double AtomToMassFrac(double atomfrac, Composition::Ptr c1, Composition::Ptr c2)
double CosiWeight(cyclus::Composition::Ptr c, const std::string &spectrum)
FuelFab takes in 2 streams of material and mixes them in ratios in order to supply material that matc...
home mouginot work app cyclus cycamore _tmp_build cycamore reactor cc
cycamore::GrowthRegion string
cyclus::Agent * ConstructFuelFab(cyclus::Context *ctx)
virtual void AcceptMatlTrades(const std::vector< std::pair< cyclus::Trade< cyclus::Material >, cyclus::Material::Ptr > > &responses)
TopupConverter(Composition::Ptr c_fill, Composition::Ptr c_fiss, Composition::Ptr c_topup, std::string spectrum)
Definition: src/fuel_fab.cc:96
FuelFab(cyclus::Context *ctx)
FillConverter(Composition::Ptr c_fill, Composition::Ptr c_fiss, Composition::Ptr c_topup, std::string spectrum)
Definition: src/fuel_fab.cc:57
virtual void GetMatlTrades(const std::vector< cyclus::Trade< cyclus::Material > > &trades, std::vector< std::pair< cyclus::Trade< cyclus::Material >, cyclus::Material::Ptr > > &responses)
virtual double convert(cyclus::Material::Ptr m, cyclus::Arc const *a=NULL, cyclus::ExchangeTranslationContext< cyclus::Material > const *ctx=NULL) const
bool ValidWeights(double w_low, double w_target, double w_high)
virtual double convert(cyclus::Material::Ptr m, cyclus::Arc const *a=NULL, cyclus::ExchangeTranslationContext< cyclus::Material > const *ctx=NULL) const
Definition: src/fuel_fab.cc:67
virtual double convert(cyclus::Material::Ptr m, cyclus::Arc const *a=NULL, cyclus::ExchangeTranslationContext< cyclus::Material > const *ctx=NULL) const
Definition: src/fuel_fab.cc:27
virtual std::set< cyclus::RequestPortfolio< cyclus::Material >::Ptr > GetMatlRequests()
double LowFrac(double w_low, double w_target, double w_high, double eps)
virtual std::set< cyclus::BidPortfolio< cyclus::Material >::Ptr > GetMatlBids(cyclus::CommodMap< cyclus::Material >::type &commod_requests)