CYCLUS
material.cc
Go to the documentation of this file.
1 #include "material.h"
2 
3 #include <math.h>
4 
5 #include "comp_math.h"
6 #include "context.h"
7 #include "decayer.h"
8 #include "error.h"
9 #include "logger.h"
10 
11 namespace cyclus {
12 
13 const ResourceType Material::kType = "Material";
14 
16 
18  Composition::Ptr c) {
19  Material::Ptr m(new Material(creator->context(), quantity, c));
20  m->tracker_.Create(creator);
21  return m;
22 }
23 
25  Composition::Ptr c) {
26  Material::Ptr m(new Material(NULL, quantity, c));
27  return m;
28 }
29 
30 int Material::qual_id() const {
31  return comp_->id();
32 }
33 
35  return Material::kType;
36 }
37 
39  Material* m = new Material(*this);
40  Resource::Ptr c(m);
41  m->tracker_.DontTrack();
42  return c;
43 }
44 
45 void Material::Record(Context* ctx) const {
46  // Note that no time field is needed because the resource ID changes
47  // every time the resource changes - state_id by itself is already unique.
48  ctx_->NewDatum("MaterialInfo")
49  ->AddVal("ResourceId", state_id())
50  ->AddVal("PrevDecayTime", prev_decay_time_)
51  ->Record();
52 
53  comp_->Record(ctx);
54 }
55 
56 std::string Material::units() const {
57  return "kg";
58 }
59 
60 double Material::quantity() const {
61  return qty_;
62 }
63 
65  return boost::static_pointer_cast<Resource>(ExtractQty(qty));
66 }
67 
69  return ExtractComp(qty, comp_);
70 }
71 
73  double threshold) {
74  if (qty_ < qty) {
75  throw ValueError("mass extraction causes negative quantity");
76  }
77 
78  // TODO: decide if ExtractComp should force lazy-decay by calling comp()
79  if (comp_ != c) {
80  CompMap v(comp_->mass());
81  compmath::Normalize(&v, qty_);
82  CompMap otherv(c->mass());
83  compmath::Normalize(&otherv, qty);
84  CompMap newv = compmath::Sub(v, otherv);
85  compmath::ApplyThreshold(&newv, threshold);
86  comp_ = Composition::CreateFromMass(newv);
87  }
88 
89  qty_ -= qty;
90 
91  Material::Ptr other(new Material(ctx_, qty, c));
92 
93  // Decay called on the extracted material should have the same dt as for
94  // this material regardless of composition.
95  other->prev_decay_time_ = prev_decay_time_;
96 
97  tracker_.Extract(&other->tracker_);
98 
99  return other;
100 }
101 
103  // these calls force lazy evaluation if in lazy decay mode
104  Composition::Ptr c0 = comp();
105  Composition::Ptr c1 = mat->comp();
106 
107  if (c0 != c1) {
108  CompMap v(c0->mass());
109  compmath::Normalize(&v, qty_);
110  CompMap otherv(c1->mass());
111  compmath::Normalize(&otherv, mat->qty_);
112  comp_ = Composition::CreateFromMass(compmath::Add(v, otherv));
113  }
114 
115  // Set the decay time to the value of the material that had the larger
116  // quantity. This helps avoid inheriting erroneous prev decay times if, for
117  // example, you absorb a material into a zero-quantity material that had a
118  // prev decay time prior to the current simulation time step.
119  if (qty_ < mat->qty_) {
120  prev_decay_time_ = mat->prev_decay_time_;
121  }
122 
123  qty_ += mat->qty_;
124  mat->qty_ = 0;
125  tracker_.Absorb(&mat->tracker_);
126 }
127 
129  comp_ = c;
130  tracker_.Modify();
131 
132  // Presumably the user has chosen the new composition to be accurate for
133  // the current simulation time. The next call to decay should not include
134  // accumulated decay delta t from a composition that no longer exists in
135  // this material. This ---------+ condition allows testing to work.
136  // |
137  // |
138  // V
139  if (ctx_ != NULL && ctx_->time() > prev_decay_time_) {
140  prev_decay_time_ = ctx_->time();
141  }
142 }
143 
144 void Material::Decay(int curr_time) {
145  if (ctx_ != NULL && ctx_->sim_info().decay == "never") {
146  return;
147  } else if (curr_time < 0 && ctx_ == NULL) {
148  throw ValueError("decay cannot use default time with NULL context");
149  }
150 
151  if (curr_time < 0) {
152  curr_time = ctx_->time();
153  }
154 
155  int dt = curr_time - prev_decay_time_;
156  if (dt == 0) {
157  return;
158  }
159 
160  double eps = 1e-3;
161  const CompMap c = comp_->atom();
162 
163  // If composition has too many nuclides (i.e. > 100), it is cheaper to
164  // just do the decay rather than check all the decay constants.
165  bool decay = c.size() > 100;
166 
167  uint64_t secs_per_timestep = kDefaultTimeStepDur;
168  if (ctx_ != NULL) {
169  secs_per_timestep = ctx_->sim_info().dt;
170  }
171 
172  if (!decay) {
173  // Only do the decay calc if one of the nuclides would change in number
174  // density more than fraction eps.
175  // i.e. decay if (1 - eps) > exp(-lambda*dt)
176  CompMap::const_reverse_iterator it;
177  for (it = c.rbegin(); it != c.rend(); ++it) {
178  int nuc = it->first;
179  double lambda_timesteps = pyne::decay_const(nuc) * static_cast<double>(secs_per_timestep);
180  double change = 1.0 - std::exp(-lambda_timesteps * static_cast<double>(dt));
181  if (change >= eps) {
182  decay = true;
183  break;
184  }
185  }
186  if (!decay) {
187  return;
188  }
189  }
190 
191  prev_decay_time_ = curr_time; // this must go before Transmute call
192  Composition::Ptr decayed = comp_->Decay(dt, secs_per_timestep);
193  Transmute(decayed);
194 }
195 
197  double decay_heat = 0.;
198  // Pyne decay heat operates with grams, cyclus generally in kilograms.
199  pyne::Material p_map = pyne::Material(comp_->mass(), qty_ * 1000);
200  std::map<int, double> dec_heat = p_map.decay_heat();
201  for (auto nuc : dec_heat) {
202  decay_heat += nuc.second;
203  }
204  return decay_heat;
205 }
206 
208  throw Error("comp() const is deprecated - use non-const comp() function."
209  " Recompilation should fix the problem.");
210 }
211 
213  if (ctx_ != NULL && ctx_->sim_info().decay == "lazy") {
214  Decay(-1);
215  }
216  return comp_;
217 }
218 
220  : qty_(quantity),
221  comp_(c),
222  tracker_(ctx, this),
223  ctx_(ctx),
224  prev_decay_time_(0) {
225  if (ctx != NULL) {
226  prev_decay_time_ = ctx->time();
227  } else {
228  tracker_.DontTrack();
229  }
230 }
231 
234  return Material::CreateUntracked(quantity, comp);
235 }
236 
237 } // namespace cyclus
void Absorb(ResTracker *absorbed)
Should be called when a resource is combined with another.
Definition: res_tracker.cc:57
void Transmute(Composition::Ptr c)
Changes the material&#39;s composition to c without changing its mass.
Definition: material.cc:128
boost::shared_ptr< Composition > Ptr
Definition: composition.h:43
virtual int qual_id() const
Returns the id of the material&#39;s internal nuclide composition.
Definition: material.cc:30
CompMap Sub(const CompMap &v1, const CompMap &v2)
Does component-wise subtraction of the nuclide quantities of v1 and v2 and returns the result...
Definition: comp_math.cc:27
A generic mechanism to manually manage exceptions.
Definition: error.h:12
void Decay(int curr_time)
Updates the material&#39;s composition by performing a decay calculation.
Definition: material.cc:144
Composition::Ptr comp()
Returns the nuclide composition of this material.
Definition: material.cc:212
For values that are too big, too small, etc.
Definition: error.h:41
boost::shared_ptr< Material > Ptr
Definition: material.h:75
static const ResourceType kType
Definition: material.h:76
std::string ResourceType
Definition: resource.h:12
Ptr ExtractComp(double qty, Composition::Ptr c, double threshold=eps_rsrc())
Creates a new material by extracting from this one.
Definition: material.cc:72
virtual double quantity() const
Returns the mass of this material in kg.
Definition: material.cc:60
std::map< Nuc, double > CompMap
a raw definition of nuclides and corresponding (dimensionless quantities).
Definition: composition.h:17
virtual void Record(Context *ctx) const
Records the internal nuclide composition of this resource.
Definition: material.cc:45
static Ptr CreateUntracked(double quantity, Composition::Ptr c)
Creates a new material resource that does not actually exist as part of the simulation and is untrack...
Definition: material.cc:24
Ptr ExtractQty(double qty)
Same as ExtractComp with c = this->comp().
Definition: material.cc:68
Resource defines an abstract interface implemented by types that are offered, requested, and transferred between simulation agents.
Definition: resource.h:19
virtual const ResourceType type() const
Returns Material::kType.
Definition: material.cc:34
uint64_t dt
Duration in seconds of a single time step in the simulation.
Definition: context.h:99
virtual int time()
Returns the current simulation timestep.
Definition: context.cc:225
static Ptr Create(Agent *creator, double quantity, Composition::Ptr c)
Creates a new material resource that is "live" and tracked.
Definition: material.cc:17
CompMap Add(const CompMap &v1, const CompMap &v2)
Does component-wise subtraction of the nuclide quantities of v1 and v2 and returns the result...
Definition: comp_math.cc:18
const uint64_t kDefaultTimeStepDur
Definition: context.h:21
comp_map decay_heat()
Calculates the decay heat of a material based on the composition and each nuclide&#39;s mass...
Definition: pyne.cc:18136
virtual Resource::Ptr ExtractRes(double qty)
Splits the resource and returns the extracted portion as a new resource object.
Definition: material.cc:64
double DecayHeat()
Returns a double with the decay heat of the material in units of W/kg.
Definition: material.cc:196
struct pyne::decay decay
a struct matching the &#39;/decay/decays&#39; table in nuc_data.h5.
void DontTrack()
Prevent a resource&#39;s heritage from being tracked and recorded.
Definition: res_tracker.cc:14
Datum * AddVal(const char *field, boost::spirit::hold_any val, std::vector< int > *shape=NULL)
Add an arbitrary field-value pair to the datum.
Definition: datum.cc:12
std::string decay
"manual" if use of the decay function is allowed, "never" otherwise
Definition: context.h:77
virtual Resource::Ptr Clone() const
Creates an untracked copy of this material object.
Definition: material.cc:38
void Modify()
Should be called when the state of a resource changes (e.g.
Definition: res_tracker.cc:32
Material::Ptr NewBlankMaterial(double quantity)
Creates and returns a new material with the specified quantity and a default, meaningless composition...
Definition: material.cc:232
void Extract(ResTracker *removed)
Should be called when a resource has some quantity removed from it (e.g.
Definition: res_tracker.cc:42
void Absorb(Ptr mat)
Combines material mat with this one. mat&#39;s quantity becomes zero.
Definition: material.cc:102
Code providing rudimentary logging capability for the Cyclus core.
A simulation context provides access to necessary simulation-global functions and state...
Definition: context.h:128
Material(Context *ctx, double quantity, Composition::Ptr c)
Definition: material.cc:219
boost::shared_ptr< Resource > Ptr
Definition: resource.h:24
The abstract base class used by all types of agents that live and interact in a simulation.
Definition: agent.h:51
taken directly from OsiSolverInterface.cpp on 2/17/14 from https://projects.coin-or.org/Osi/browser/trunk.
Definition: agent.cc:14
The material class is primarily responsible for enabling basic material manipulation while helping en...
Definition: material.h:71
Context * context() const
Returns this agent&#39;s simulation context.
Definition: agent.h:369
double decay_const(int nuc)
Returns the decay constant for a nuclide nuc.
Definition: pyne.cc:11914
void Record()
Record this datum to its Recorder.
Definition: datum.cc:24
SimInfo sim_info() const
Return static simulation info.
Definition: context.h:238
Material composed of nuclides.
Definition: pyne.h:4800
Datum * NewDatum(std::string title)
See Recorder::NewDatum documentation.
Definition: context.cc:237
virtual std::string units() const
Returns "kg".
Definition: material.cc:56
const int state_id() const
Returns the unique id corresponding to this resource and its current state.
Definition: resource.h:39
virtual ~Material()
Definition: material.cc:15
void Normalize(CompMap *v, double val)
The sum of quantities of all nuclides of v is normalized to val.
Definition: comp_math.cc:63
double eps()
a generic epsilon value
Definition: cyc_limits.h:12
static Ptr CreateFromMass(CompMap v)
Creates a new composition from v with its components having appropriate mass-based ratios...
Definition: composition.cc:26
void ApplyThreshold(CompMap *v, double threshold)
All nuclides with quantities below threshold will have their quantity set to zero.
Definition: comp_math.cc:45