Merge pull request #333 from rwcarlsen/fab-fix

fix fuel fab buffer composition assumption bug

src/fuel_fab.cc

@@ -225,7 +225,20 @@ void FuelFab::AcceptMatlTrades(const std::vector<
       throw cyclus::ValueError("cycamore::FuelFab was overmatched on requests");
     }
   }
+
   req_inventories_.clear();
+
+  // IMPORTANT - each buffer needs to be a single homogenous composition or
+  // the inventory mixing constraints for bids don't work
+  if (fill.count() > 1) {
+    fill.Push(cyclus::toolkit::Squash(fill.PopN(fill.count())));
+  }
+  if (fiss.count() > 1) {
+    fiss.Push(cyclus::toolkit::Squash(fiss.PopN(fiss.count())));
+  }
+  if (topup.count() > 1) {
+    topup.Push(cyclus::toolkit::Squash(topup.PopN(topup.count())));
+  }
 }
 
 std::set<cyclus::BidPortfolio<Material>::Ptr> FuelFab::GetMatlBids(
@@ -339,9 +352,10 @@ void FuelFab::GetMatlTrades(
         responses) {
   using cyclus::Trade;
 
+  // guard against cases where a buffer is empty - this is okay because some trades
+  // may not need that particular buffer.
   double w_fill = 0;
-  if (fill.count() >
-      0) {  // it's possible to only need fissile inventory for a trade
+  if (fill.count() > 0) {
     w_fill = CosiWeight(fill.Peek()->comp(), spectrum);
   }
   double w_topup = 0;
@@ -357,6 +371,7 @@ void FuelFab::GetMatlTrades(
   double tot = 0;
   for (int i = 0; i < trades.size(); i++) {
     Material::Ptr tgt = trades[i].request->target();
+
     double w_tgt = CosiWeight(tgt->comp(), spectrum);
     double qty = trades[i].amt;
     double wfiss = w_fiss;
@@ -371,10 +386,18 @@ void FuelFab::GetMatlTrades(
 
     if (fiss.count() == 0) {
       // use straight filler to satisfy this request
-      responses.push_back(std::make_pair(trades[i], fill.Pop(qty)));
+      double fillqty = qty;
+      if (std::abs(fillqty - fill.quantity()) < cyclus::eps()) {
+        fillqty = std::min(fill.quantity(), qty);
+      }
+      responses.push_back(std::make_pair(trades[i], fill.Pop(fillqty)));
     } else if (fill.count() == 0 && ValidWeights(w_fill, w_tgt, w_fiss)) {
       // use straight fissile to satisfy this request
-      responses.push_back(std::make_pair(trades[i], fiss.Pop(qty)));
+      double fissqty = qty;
+      if (std::abs(fissqty - fiss.quantity()) < cyclus::eps()) {
+        fissqty = std::min(fiss.quantity(), qty);
+      }
+      responses.push_back(std::make_pair(trades[i], fiss.Pop(fissqty)));
     } else if (ValidWeights(w_fill, w_tgt, w_fiss)) {
       double fiss_frac = HighFrac(w_fill, w_tgt, w_fiss);
       double fill_frac = LowFrac(w_fill, w_tgt, w_fiss);
@@ -383,10 +406,19 @@ void FuelFab::GetMatlTrades(
       fill_frac =
           AtomToMassFrac(fill_frac, fill.Peek()->comp(), fiss.Peek()->comp());
 
-      Material::Ptr m = fiss.Pop(fiss_frac * qty);
+      double fissqty = fiss_frac*qty;
+      if (std::abs(fissqty - fiss.quantity()) < cyclus::eps()) {
+        fissqty = std::min(fiss.quantity(), fiss_frac*qty);
+      }
+      double fillqty = fill_frac*qty;
+      if (std::abs(fillqty - fill.quantity()) < cyclus::eps()) {
+        fillqty = std::min(fill.quantity(), fill_frac*qty);
+      }
+
+      Material::Ptr m = fiss.Pop(fissqty);
       // this if block prevents zero qty ResBuf pop exceptions
       if (fill_frac > 0) {
-        m->Absorb(fill.Pop(fill_frac * qty));
+        m->Absorb(fill.Pop(fillqty));
       }
       responses.push_back(std::make_pair(trades[i], m));
     } else {
@@ -397,10 +429,19 @@ void FuelFab::GetMatlTrades(
       fiss_frac =
           AtomToMassFrac(fiss_frac, fiss.Peek()->comp(), topup.Peek()->comp());
 
-      Material::Ptr m = fiss.Pop(fiss_frac * qty);
+      double fissqty = fiss_frac*qty;
+      if (std::abs(fissqty - fiss.quantity()) < cyclus::eps()) {
+        fissqty = std::min(fiss.quantity(), fiss_frac*qty);
+      }
+      double topupqty = topup_frac*qty;
+      if (std::abs(topupqty - topup.quantity()) < cyclus::eps()) {
+        topupqty = std::min(topup.quantity(), topup_frac*qty);
+      }
+
+      Material::Ptr m = fiss.Pop(fissqty);
       // this if block prevents zero qty ResBuf pop exceptions
       if (topup_frac > 0) {
-        m->Absorb(topup.Pop(topup_frac * qty));
+        m->Absorb(topup.Pop(topupqty));
       }
       responses.push_back(std::make_pair(trades[i], m));
     }

src/fuel_fab_tests.cc

@@ -837,6 +837,52 @@ TEST(FuelFabTests, SwapTopup_FissConstrained) {
   EXPECT_NEAR(max_provide, m->quantity(), 1e-10) << "matched trade uses more fiss than available";
 }
 
+// Before this test and a fix, the fuel fab (partially) assumed each entire material
+// buffer had the same composition as the material on top of the buffer when
+// calculating stream mixing ratios for material to supply.  This problem was
+// compounded by the fact that material weights are computed on an atom basis
+// and mixing is done on a mass basis - corresponding conversions resulted in
+// the fab being matched for more than it could actually supply - due to
+// thinking it had an inventory of higher quality material than was actually
+// the case.  This test makes sure that doesn't happen again.
+TEST(FuelFabTests, HomogenousBuffers) {
+  std::string config = 
+     "<fill_commod>natu</fill_commod>"
+     "<fill_recipe>natu</fill_recipe>"
+     "<fill_size>40</fill_size>"
+     ""
+     "<fiss_commods> <val>stream1</val> </fiss_commods>"
+     "<fiss_size>4</fiss_size>"
+     "<fiss_recipe>spentuox</fiss_recipe>"
+     ""
+     "<outcommod>out</outcommod>"
+     "<spectrum>thermal</spectrum>"
+     "<throughput>1e10</throughput>"
+     ;
+
+  CompMap m;
+  m[id("u235")] = 7;
+  m[id("u238")] = 86;
+  // the zr90 is important to force the atom-mass basis conversion to push the
+  // dre to overmatch in the direction we want.
+  m[id("zr90")] = 7;
+  Composition::Ptr c = Composition::CreateFromMass(m);
+
+  int simdur = 5;
+  cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:FuelFab"), config, simdur);
+  sim.AddSource("stream1").start(0).lifetime(1).capacity(.01).recipe("special").Finalize();
+  sim.AddSource("stream1").start(1).lifetime(1).capacity(3.98).recipe("natu").Finalize();
+  sim.AddSource("natu").lifetime(1).Finalize();
+  sim.AddSink("out").start(2).capacity(4).lifetime(1).recipe("uox").Finalize();
+  sim.AddSink("out").start(2).capacity(4).lifetime(1).recipe("uox").Finalize();
+  sim.AddRecipe("uox", c_uox());
+  sim.AddRecipe("spentuox", c_pustream());
+  sim.AddRecipe("natu", c_natu());
+  sim.AddRecipe("special", c);
+  ASSERT_NO_THROW(sim.Run());
+}
+
 } // namespace fuelfabtests
 } // namespace cycamore
 
+