Skip to content

Loading…

Removed decay data, now with pyne! #762

Merged
merged 10 commits into from

3 participants

@scopatz
Cyclus member

This PR replaces the decayInfo.dat file with nuclear data coming from pyne. A custom cyclus_nuc_data.h5 has been created and uploaded to our webspace. The install process will automatically download and install this file. This data file may be generated by the following command:

scopatz@ares ~/pyne staging $ nuc_data_make -o cyclus_nuc_data.h5 
    -m atomic_mass,scattering_lengths,decay,simple_xs,materials,eaf,wimsd_fpy,nds_fpy

This PR also exposed what I would consider some fundamental problems with how the current methodology works. These will be addressed in separate issues as the purpose here is only to integrate with pyne data. As much as I wanted to throw my machine out of the 4th floor window of ERB during this process, the unite tests now pass :).

@scopatz scopatz added this to the v0.5 milestone
@gidden gidden commented on the diff
cli/cyclus.cc
@@ -167,6 +167,9 @@ int main(int argc, char* argv[]) {
std::string path = Env::PathBase(argv[0]);
Env::SetCyclusRelPath(path);
+ // tell pyne about the path to nuc data
+ pyne::NUC_DATA_PATH = Env::GetBuildPath() + "/share/cyclus_nuc_data.h5";
@gidden Cyclus member
gidden added a note

should likely use a boost file system path for os robustness

@scopatz Cyclus member
scopatz added a note

There is a lot of this kind of code in cyclus right now. Maybe we could make an issue and tackle it later.

@gidden Cyclus member
gidden added a note

I see thats the case with finding mass.sqlite as well. Would you like to make the issue? We could also make an issue about using pyne nuclide masses as well and get rid of mass.sqlite.

@scopatz Cyclus member
scopatz added a note

I see thats the case with finding mass.sqlite as well. Would you like to make the issue?

I'd rather if you did. Thanks!

We could also make an issue about using pyne nuclide masses as well and get rid of mass.sqlite.

I am going to do this tomorrow :)

@gidden Cyclus member
gidden added a note

done

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@gidden gidden commented on the diff
src/CMakeLists.txt
@@ -8,6 +8,16 @@ CONFIGURE_FILE(suffix.h.in "${CMAKE_CURRENT_SOURCE_DIR}/suffix.h" @ONLY)
EXECUTE_PROCESS(COMMAND git describe --tags OUTPUT_VARIABLE core_version OUTPUT_STRIP_TRAILING_WHITESPACE)
CONFIGURE_FILE(version.h.in "${CMAKE_CURRENT_SOURCE_DIR}/version.h" @ONLY)
+if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/cyclus_nuc_data.h5")
+ message("-- Found '${CMAKE_CURRENT_SOURCE_DIR}/cyclus_nuc_data.h5'")
+else()
+ message("-- Downloading http://regtests.fuelcycle.org/cyclus_nuc_data.h5 to "
+ "'${CMAKE_CURRENT_SOURCE_DIR}/cyclus_nuc_data.h5'...")
+ file(DOWNLOAD "http://regtests.fuelcycle.org/cyclus_nuc_data.h5"
@gidden Cyclus member
gidden added a note

TIL. Neat!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@gidden gidden commented on an outdated diff
src/composition.cc
@@ -105,11 +105,13 @@ void Composition::Record(Context* ctx) {
Composition::Ptr Composition::NewDecay(int delta) {
int tot_decay = prev_decay_ + delta;
- double months_per_year = 12;
- double years = double(delta) / months_per_year;
+
+ if (atom_.size() == 0)
+ // FIXME this is only here for testing
+ return (Composition::Ptr) new Composition(tot_decay, decay_line_);
@gidden Cyclus member
gidden added a note
@rwcarlsen Cyclus member

Also - need curlies on multi-line if. Also - the casting syntax is strange and non-obvious. You should use the constructor syntax - which is what is really going on:

return Composition::Ptr(new Composition(tot_decay, decay_line_));
@rwcarlsen Cyclus member

@scopatz I just noticed a bug here - could you throw in this line at the top of the NewDecay function:

  atom(); // force evaluation of atom-composition if not calculated already
@scopatz Cyclus member
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@gidden gidden commented on the diff
src/decayer.cc
((105 lines not shown))
- nuc));
- }
- }
- temp[i] = std::make_pair(nuc, branchRatio);
- }
-
- daughters_[jcol] = temp;
- ++jcol; // set next column
- decayInfo >> nuc; // get next parent
- nuc *= 10000; // put in id form
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void Decayer::AddNucToMaps(int nuc) {
+ int i;
+ int col;
+ int daughter;
+ std::set<int> daughters;
@gidden Cyclus member
gidden added a note

I believe in a PyNE PR, someone suggested either from_nuc/to_nuc, parent/child, or mother/daughter, in that preference order =)

@scopatz Cyclus member
scopatz added a note

I feel it safe to call this legacy code that matches the names around it :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@gidden gidden commented on an outdated diff
src/decayer.cc
@@ -162,7 +110,7 @@ void Decayer::GetResult(CompMap& comp) {
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void Decayer::BuildDecayMatrix() {
- double decayConst = 0; // decay constant, in inverse years
+ double decayConst = 0; // decay constant, in inverse secs
@gidden Cyclus member
gidden added a note
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@gidden gidden commented on an outdated diff
src/decayer.cc
@@ -173,6 +121,9 @@ void Decayer::BuildDecayMatrix() {
while (parent_iter != parent_.end()) {
jcol = parent_iter->second.first; // determines column index
decayConst = parent_iter->second.second;
+ // Gross heuristic for mostly stable nuclides 2903040000 sec / 100 years
+ if ((long double) exp(-2903040000 * decayConst) == 0.0)
@gidden Cyclus member
gidden added a note

c++ style casts

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@gidden gidden commented on the diff
src/material.cc
@@ -128,7 +127,8 @@ void Material::Decay(int curr_time) {
CompMap::const_iterator it;
for (it = c.end(); it != c.begin(); --it) {
int nuc = it->first;
- double lambda_months = Decayer::DecayConstant(nuc) / 12;
+ // 2419200 == secs / month
+ double lambda_months = pyne::decay_const(nuc) * 2419200;
@gidden Cyclus member
gidden added a note

shouldn't [months] = [secs] / [secs / month]?

@gidden Cyclus member
gidden added a note

errr, [months] = 1 / ([secs / month] * [1 / secs])

@scopatz Cyclus member
scopatz added a note

lambda: [1/months] = [1/sec] * [sec/month]

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@gidden
Cyclus member

we definitely need decayer unit tests..

@rwcarlsen
Cyclus member

+1 for tests. I also would like to award anthony a code destroyer medal. Code is the enemy.

@gidden - who is in your new profile pic?

@gidden
Cyclus member
@scopatz
Cyclus member

+1 for tests.

I completely agree, but I think that is a separate issue.

I also would like to award anthony a code destroyer medal. Code is the enemy.

Thanks! Hopefully I'll get rid of even more tomorrow.

I think this is ready for another look.

@rwcarlsen rwcarlsen commented on an outdated diff
src/composition.cc
@@ -105,11 +105,14 @@ void Composition::Record(Context* ctx) {
Composition::Ptr Composition::NewDecay(int delta) {
int tot_decay = prev_decay_ + delta;
- double months_per_year = 12;
- double years = double(delta) / months_per_year;
+ atom(); // force evaluation of atom-composition if not calculated already
+
+ // FIXME this is only here for testing
+ if (atom_.size() == 0)
+ return static_cast<Composition::Ptr>(new Composition(tot_decay, decay_line_));
@rwcarlsen Cyclus member

This should be:

return Composition::Ptr(new Composition(tot_decay, decay_line_));

Composition::Ptr is a smart pointer - so we want to call its constructor. A static cast is wrong.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@scopatz
Cyclus member

ping.

@gidden
Cyclus member

Looks good. Thanks, @scopatz!

@gidden gidden merged commit 825f02b into cyclus:develop

1 check passed

Details default build and test completed successfully
@scopatz scopatz deleted the scopatz:rm-decay branch
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Mar 4, 2014
  1. @scopatz

    part way through

    scopatz committed
  2. @scopatz

    so sleepy

    scopatz committed
  3. @scopatz

    first cute

    scopatz committed
  4. @scopatz

    compiles and downloads file

    scopatz committed
Commits on Mar 5, 2014
  1. @scopatz

    at least cyclus runs

    scopatz committed
  2. @scopatz

    yay! unit tests pass!

    scopatz committed
  3. @scopatz
  4. @scopatz

    some feedback fixes

    scopatz committed
  5. @scopatz
  6. @scopatz

    fixed comp

    scopatz committed
Showing with 113 additions and 266 deletions.
  1. +1 −0 .gitignore
  2. +3 −0 cli/cyclus.cc
  3. +13 −3 src/CMakeLists.txt
  4. +6 −3 src/composition.cc
  5. +0 −111 src/decayInfo.dat
  6. +57 −106 src/decayer.cc
  7. +21 −39 src/decayer.h
  8. +2 −2 src/material.cc
  9. +3 −1 src/uniform_taylor.cc
  10. +2 −0 src/uniform_taylor.h
  11. +5 −1 tests/composition_tests.cc
View
1 .gitignore
@@ -21,3 +21,4 @@ suffix.h
*.html
*.tex
src/Core/Utility/Env.cpp
+src/cyclus_nuc_data.h5
View
3 cli/cyclus.cc
@@ -167,6 +167,9 @@ int main(int argc, char* argv[]) {
std::string path = Env::PathBase(argv[0]);
Env::SetCyclusRelPath(path);
+ // tell pyne about the path to nuc data
+ pyne::NUC_DATA_PATH = Env::GetBuildPath() + "/share/cyclus_nuc_data.h5";
@gidden Cyclus member
gidden added a note

should likely use a boost file system path for os robustness

@scopatz Cyclus member
scopatz added a note

There is a lot of this kind of code in cyclus right now. Maybe we could make an issue and tackle it later.

@gidden Cyclus member
gidden added a note

I see thats the case with finding mass.sqlite as well. Would you like to make the issue? We could also make an issue about using pyne nuclide masses as well and get rid of mass.sqlite.

@scopatz Cyclus member
scopatz added a note

I see thats the case with finding mass.sqlite as well. Would you like to make the issue?

I'd rather if you did. Thanks!

We could also make an issue about using pyne nuclide masses as well and get rid of mass.sqlite.

I am going to do this tomorrow :)

@gidden Cyclus member
gidden added a note

done

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+
// read input file and setup simulation
std::string inputFile = vm["input-file"].as<std::string>();
XMLFileLoader* loader;
View
16 src/CMakeLists.txt
@@ -8,6 +8,16 @@ CONFIGURE_FILE(suffix.h.in "${CMAKE_CURRENT_SOURCE_DIR}/suffix.h" @ONLY)
EXECUTE_PROCESS(COMMAND git describe --tags OUTPUT_VARIABLE core_version OUTPUT_STRIP_TRAILING_WHITESPACE)
CONFIGURE_FILE(version.h.in "${CMAKE_CURRENT_SOURCE_DIR}/version.h" @ONLY)
+if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/cyclus_nuc_data.h5")
+ message("-- Found '${CMAKE_CURRENT_SOURCE_DIR}/cyclus_nuc_data.h5'")
+else()
+ message("-- Downloading http://regtests.fuelcycle.org/cyclus_nuc_data.h5 to "
+ "'${CMAKE_CURRENT_SOURCE_DIR}/cyclus_nuc_data.h5'...")
+ file(DOWNLOAD "http://regtests.fuelcycle.org/cyclus_nuc_data.h5"
@gidden Cyclus member
gidden added a note

TIL. Neat!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ "${CMAKE_CURRENT_SOURCE_DIR}/cyclus_nuc_data.h5")
+ message("-- ...download complete!")
+endif()
+
CONFIGURE_FILE(
"${CMAKE_CURRENT_SOURCE_DIR}/mass.sqlite"
"${CYCLUS_BINARY_DIR}/share/mass.sqlite"
@@ -15,8 +25,8 @@ CONFIGURE_FILE(
)
CONFIGURE_FILE(
- "${CMAKE_CURRENT_SOURCE_DIR}/decayInfo.dat"
- "${CYCLUS_BINARY_DIR}/share/decayInfo.dat"
+ "${CMAKE_CURRENT_SOURCE_DIR}/cyclus_nuc_data.h5"
+ "${CYCLUS_BINARY_DIR}/share/cyclus_nuc_data.h5"
COPYONLY
)
@@ -32,7 +42,7 @@ INSTALL(FILES
"${CMAKE_CURRENT_SOURCE_DIR}/mass.sqlite"
"${CMAKE_CURRENT_SOURCE_DIR}/cyclus.rng.in"
"${CMAKE_CURRENT_SOURCE_DIR}/cyclus-flat.rng.in"
- "${CMAKE_CURRENT_SOURCE_DIR}/decayInfo.dat"
+ "${CMAKE_CURRENT_SOURCE_DIR}/cyclus_nuc_data.h5"
DESTINATION share
COMPONENT core
)
View
9 src/composition.cc
@@ -105,11 +105,14 @@ Composition::Composition(int prev_decay, ChainPtr decay_line)
Composition::Ptr Composition::NewDecay(int delta) {
int tot_decay = prev_decay_ + delta;
- double months_per_year = 12;
- double years = double(delta) / months_per_year;
+ atom(); // force evaluation of atom-composition if not calculated already
+
+ // FIXME this is only here for testing, see issue #761
+ if (atom_.size() == 0)
+ return Composition::Ptr(new Composition(tot_decay, decay_line_));
Decayer handler(atom_);
- handler.Decay(years);
+ handler.Decay(2419200.0 * delta); // 2419200 == secs / month
// the new composition is a part of this decay chain and so is created with a
// pointer to the exact same decay_line_.
View
111 src/decayInfo.dat
@@ -1,111 +0,0 @@
-2004 0 0
-88226 4.332170E-04 1
- 82210 1
-88228 1.205473E-01 1
- 90228 1
-82206 0 0
-82207 0 0
-82208 0 0
-82210 3.113869E-02 1
- 82206 1
-90228 3.623352E-01 1
- 82208 1
-90229 9.443422E-05 1
- 83209 1
-90230 9.001911E-06 1
- 88226 1
-90232 4.933432E-11 1
- 88228 1
-83209 0 0
-89227 3.183956E-02 1
- 82207 1
-91231 2.310491E-05 1
- 89227 1
-92232 9.627044E-03 1
- 90228 1
-92233 4.332170E-06 1
- 90229 1
-92234 2.888113E-06 1
- 90230 1
-92235 9.848639E-10 1
- 91231 1
-92236 2.960270E-08 1
- 90232 1
-92238 1.540327E-10 1
- 92234 1
-93237 3.300701E-07 1
- 92233 1
-94238 7.876673E-03 1
- 92234 1
-94239 2.872435E-05 1
- 92235 1
-94240 1.055179E-04 1
- 92236 1
-94241 4.813522E-02 1
- 95241 1
-94242 1.824072E-06 1
- 92238 1
-94244 8.391612E-09 1
- 94240 1
-95241 1.604507E-03 1
- 93237 1
-952421 4.560179E-03 3
- 96242 0.822865
- 94242 0.172135
- 94238 0.005
-95243 9.392238E-05 1
- 94239 1
-96242 1.553202E+00 1
- 94238 1
-96243 2.432095E-02 2
- 94239 0.998
- 95243 0.002
-96244 3.827428E-02 1
- 94240 1
-96245 8.154673E-05 1
- 94241 1
-96246 1.459257E-04 1
- 94242 1
-96247 4.443251E-08 1
- 95243 1
-96248 2.044682E-06 1
- 94244 1
-96250 1.004561E-04 2
- 96246 0.86
- 98250 0.14
-98249 1.977031E-03 1
- 96245 1
-98250 5.331901E-02 1
- 96246 1
-98251 7.701635E-04 1
- 96247 1
-98252 2.626552E-01 1
- 96248 1
-1003 5.626195E-02 1
- 8881159 1
-6014 1.212856E-04 1
- 8881159 1
-6129 0 0
-36081 3.013683E-06 1
- 8881159 1
-36085 6.441888E-02 1
- 8881159 1
-36849 0 0
-38090 2.408434E-02 1
- 8881159 1
-38889 0 0
-43099 3.254212E-06 1
- 8881159 1
-43989 0 0
-53129 4.414950E-08 1
- 8881159 1
-531279 0 0
-55134 3.356645E-01 1
- 8881159 1
-55135 3.013683E-07 1
- 8881159 1
-55137 2.305112E-02 1
- 8881159 1
-551339 0 0
-999249 0 0
-8881159 0 0
View
163 src/decayer.cc
@@ -11,7 +11,6 @@
namespace cyclus {
-bool Decayer::decay_info_loaded_ = false;
ParentMap Decayer::parent_ = ParentMap();
DaughtersMap Decayer::daughters_ = DaughtersMap();
Matrix Decayer::decay_matrix_ = Matrix();
@@ -19,120 +18,69 @@ NucList Decayer::nuclides_tracked_ = NucList();
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Decayer::Decayer(const CompMap& comp) {
- if (!decay_info_loaded_) {
- Decayer::LoadDecayInfo();
- decay_info_loaded_ = true;
- }
+ int nuc;
+ int col;
+ long double atom_count;
+ bool needs_build = false;
- pre_vect_ = Vector(parent_.size(), 1);
std::map<int, double>::const_iterator comp_iter = comp.begin();
for (comp_iter = comp.begin(); comp_iter != comp.end(); ++comp_iter) {
- int nuc = comp_iter->first;
- long double atom_count = comp_iter->second;
-
- // if the nuclide is tracked in the decay matrix
- if (parent_.count(nuc) > 0) {
- int col = parent_[nuc].first; // get Vector position
- pre_vect_(col, 1) = atom_count;
- // if it is not in the decay matrix, then it is added as a stable nuclide
- } else {
- double decayConst = 0;
- int col = parent_.size() + 1;
- parent_[nuc] = std::make_pair(col, decayConst); // add nuclide to parent map
-
- int nDaughters = 0;
- std::vector< std::pair<int, double> > temp(nDaughters);
- daughters_[col] = temp; // add nuclide to daughters map
-
- std::vector<long double> row(1, atom_count);
- pre_vect_.AddRow(row); // add nuclide to the end of the Vector
+ nuc = comp_iter->first;
+ atom_count = comp_iter->second;
+ if (!IsNucTracked(nuc)) {
+ needs_build = true;
+ AddNucToMaps(nuc);
}
}
-}
-// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-void Decayer::LoadDecayInfo() {
- std::string path = Env::GetBuildPath() + "/share/decayInfo.dat";
- std::ifstream decayInfo(path.c_str());
-
- if (!decayInfo.is_open()) {
- throw IOError("Could not find file 'decayInfo.dat'.");
- }
+ if (needs_build)
+ BuildDecayMatrix();
- int jcol = 1;
- int nuc = 0;
- int nDaughters = 0;
- double decayConst = 0; // decay constant, in inverse years
- double branchRatio = 0;
-
- decayInfo >> nuc; // get first parent
- nuc *= 10000; // put in id form
-
- // checks to see if there are nuclides in 'decayInfo.dat'
- if (decayInfo.eof()) {
- std::string err_msg = "There are no nuclides in the 'decayInfo.dat' file";
- throw ValidationError(err_msg);
+ pre_vect_ = Vector(parent_.size(), 1);
+ for (comp_iter = comp.begin(); comp_iter != comp.end(); ++comp_iter) {
+ nuc = comp_iter->first;
+ atom_count = comp_iter->second;
+ col = parent_[nuc].first;
+ pre_vect_(col, 1) = atom_count;
}
+}
- // processes 'decayInfo.dat'
- while (!decayInfo.eof()) {
- if (parent_.find(nuc) != parent_.end()) {
- std::string err_msg;
- err_msg = "A duplicate parent nuclide was found in 'decayInfo.dat'";
- throw ValidationError(err_msg);
- }
-
- // make parent
- decayInfo >> decayConst;
- decayInfo >> nDaughters;
- AddNucToList(nuc);
-
- parent_[nuc] = std::make_pair(jcol, decayConst);
-
- // make daughters
- std::vector< std::pair<int, double> > temp(nDaughters);
- for (int i = 0; i < nDaughters; ++i) {
- decayInfo >> nuc;
- nuc *= 10000; // put in id form
- decayInfo >> branchRatio;
- AddNucToList(nuc);
-
- // checks for duplicate daughter nuclides
- for (int j = 0; j < nDaughters; ++j) {
- if (temp[j].first == nuc) {
- throw ValidationError(
- std::string("A duplicate daughter nuclide, %i , was found in decayInfo.dat",
- nuc));
- }
- }
- temp[i] = std::make_pair(nuc, branchRatio);
- }
-
- daughters_[jcol] = temp;
- ++jcol; // set next column
- decayInfo >> nuc; // get next parent
- nuc *= 10000; // put in id form
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void Decayer::AddNucToMaps(int nuc) {
+ int i;
+ int col;
+ int daughter;
+ std::set<int> daughters;
@gidden Cyclus member
gidden added a note

I believe in a PyNE PR, someone suggested either from_nuc/to_nuc, parent/child, or mother/daughter, in that preference order =)

@scopatz Cyclus member
scopatz added a note

I feel it safe to call this legacy code that matches the names around it :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+ std::set<int>::iterator d;
+
+ if (IsNucTracked(nuc))
+ return;
+
+ col = parent_.size() + 1;
+ parent_[nuc] = std::make_pair(col, pyne::decay_const(nuc));
+ AddNucToList(nuc);
+
+ i = 0;
+ daughters = pyne::decay_children(nuc);
+ std::vector< std::pair<int, double> > dvec(daughters.size());
+ for (d = daughters.begin(); d != daughters.end(); ++d) {
+ daughter = *d;
+ AddNucToMaps(daughter);
+ dvec[i] = std::make_pair<int, double>(daughter, pyne::branch_ratio(nuc, daughter));
+ i++;
}
- BuildDecayMatrix();
+ daughters_[col] = dvec;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-double Decayer::DecayConstant(int nuc) {
- if (!decay_info_loaded_) {
- Decayer::LoadDecayInfo();
- decay_info_loaded_ = true;
- }
- if (parent_.count(nuc) > 0) {
- return parent_[nuc].second;
- }
- return 0;
+bool Decayer::IsNucTracked(int nuc) {
+ return (find(nuclides_tracked_.begin(), nuclides_tracked_.end(), nuc)
+ != nuclides_tracked_.end());
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void Decayer::AddNucToList(int nuc) {
- bool exists = (find(nuclides_tracked_.begin(), nuclides_tracked_.end(), nuc)
- != nuclides_tracked_.end());
- if (!exists) {
+ if (!IsNucTracked(nuc)) {
nuclides_tracked_.push_back(nuc);
}
}
@@ -162,7 +110,7 @@ void Decayer::GetResult(CompMap& comp) {
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
void Decayer::BuildDecayMatrix() {
- double decayConst = 0; // decay constant, in inverse years
+ double decay_const = 0; // decay constant, in inverse secs
int jcol = 1;
int n = parent_.size();
decay_matrix_ = Matrix(n, n);
@@ -172,13 +120,16 @@ void Decayer::BuildDecayMatrix() {
// populates the decay matrix column by column
while (parent_iter != parent_.end()) {
jcol = parent_iter->second.first; // determines column index
- decayConst = parent_iter->second.second;
- decay_matrix_(jcol, jcol) = -1 * decayConst; // sets A(i,i) value
+ decay_const = parent_iter->second.second;
+ // Gross heuristic for mostly stable nuclides 2903040000 sec / 100 years
+ if (static_cast<long double>(exp(-2903040000 * decay_const)) == 0.0)
+ decay_const = 0.0;
+ decay_matrix_(jcol, jcol) = -1 * decay_const; // sets A(i,i) value
// processes the vector in the daughters map if it is not empty
if (!daughters_.find(jcol)->second.empty()) {
// an iterator that points to 1st daughter in the vector
- // pair<nuclide,branchratio>
+ // pair<nuclide,branch_ratio>
std::vector< std::pair<int, double> >::const_iterator
nuc_iter = daughters_.find(jcol)->second.begin();
@@ -186,8 +137,8 @@ void Decayer::BuildDecayMatrix() {
while (nuc_iter != daughters_.find(jcol)->second.end()) {
int nuc = nuc_iter->first;
int irow = parent_.find(nuc)->second.first; // determines row index
- double branchRatio = nuc_iter->second;
- decay_matrix_(irow, jcol) = branchRatio * decayConst; // sets A(i,j) value
+ double branch_ratio = nuc_iter->second;
+ decay_matrix_(irow, jcol) = branch_ratio * decay_const; // sets A(i,j) value
++nuc_iter; // get next daughter
}
@@ -197,9 +148,9 @@ void Decayer::BuildDecayMatrix() {
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-void Decayer::Decay(double years) {
+void Decayer::Decay(double secs) {
// solves the decay equation for the final composition
- post_vect_ = UniformTaylor::MatrixExpSolver(decay_matrix_, pre_vect_, years);
+ post_vect_ = UniformTaylor::MatrixExpSolver(decay_matrix_, pre_vect_, secs);
}
} // namespace cyclus
View
60 src/decayer.h
@@ -3,39 +3,31 @@
#define CYCLUS_SRC_DECAYER_H_
#include <map>
+#include <set>
#include "composition.h"
+#include "pyne.h"
#include "use_matrix_lib.h"
namespace cyclus {
-/**
- A map type to represent all of the parent nuclides tracked. The key
- for this map type is the parent's Nuc number, and the value is a pair
- that contains the corresponding decay matrix column and decay
- constant associated with that parent.
- */
+/// A map type to represent all of the parent nuclides tracked. The key
+/// for this map type is the parent's Nuc number, and the value is a pair
+/// that contains the corresponding decay matrix column and decay
+/// constant associated with that parent.
typedef std::map< int, std::pair<int, double> > ParentMap;
-/**
- A map type to represent all of the daughter nuclides tracked. The
- key for this map type is the decay matrix column associated with the
- parent, and the value is a vector of pairs of all the daughters for
- that parent. Each of the daughters are represented by a pair that
- contains the daughter's Nuc number and its branching ratio.
- */
+/// A map type to represent all of the daughter nuclides tracked. The
+/// key for this map type is the decay matrix column associated with the
+/// parent, and the value is a vector of pairs of all the daughters for
+/// that parent. Each of the daughters are represented by a pair that
+/// contains the daughter's Nuc number and its branching ratio.
typedef std::map<int, std::vector<std::pair<int, double> > > DaughtersMap;
typedef std::vector<int> NucList;
class Decayer {
public:
- /// Returns the decay constant for specified nuclide in inverse years.
- static double DecayConstant(int nuc);
-
- /**
- default constructor
- */
Decayer(const CompMap& comp);
/**
@@ -43,11 +35,9 @@ class Decayer {
*/
void GetResult(CompMap& comp);
- /**
- decay the material
- @param years the number of years to decay
- */
- void Decay(double years);
+ /// decay the material
+ /// @param secs the number of seconds to decay
+ void Decay(double secs);
/**
the number of tracked nuclides
@@ -72,13 +62,6 @@ class Decayer {
static void BuildDecayMatrix();
/**
- Reads the decay information found in the 'decayInfo.dat' file
- into the parent and daughters maps. Uses these maps to create the
- decay matrix.
- */
- static void LoadDecayInfo();
-
- /**
The CompMap's parent
*/
static ParentMap parent_;
@@ -100,19 +83,18 @@ class Decayer {
Vector post_vect_;
/**
- whether the decay information is loaded
- */
- static bool decay_info_loaded_;
-
- /**
the list of tracked nuclides
*/
static NucList nuclides_tracked_;
- /**
- Add the nuclide to our list of tracked nuclides IFF it is not in the list.
- */
+ /// Add the nuclide to the parent/daughter maps IFF it is not in the tracked list.
+ static void AddNucToMaps(int nuc);
+
+ /// Add the nuclide to our list of tracked nuclides IFF it is not in the list.
static void AddNucToList(int nuc);
+
+ /// Checks if the nuclide is tracked
+ static bool IsNucTracked(int nuc);
};
} // namespace cyclus
View
4 src/material.cc
@@ -77,7 +77,6 @@ Material::Ptr Material::ExtractComp(double qty, Composition::Ptr c,
if (qty_ < qty) {
throw ValueError("mass extraction causes negative quantity");
}
-
if (comp_ != c) {
CompMap v(comp_->mass());
compmath::Normalize(&v, qty_);
@@ -128,7 +127,8 @@ void Material::Decay(int curr_time) {
CompMap::const_iterator it;
for (it = c.end(); it != c.begin(); --it) {
int nuc = it->first;
- double lambda_months = Decayer::DecayConstant(nuc) / 12;
+ // 2419200 == secs / month
+ double lambda_months = pyne::decay_const(nuc) * 2419200;
@gidden Cyclus member
gidden added a note

shouldn't [months] = [secs] / [secs / month]?

@gidden Cyclus member
gidden added a note

errr, [months] = 1 / ([secs / month] * [1 / secs])

@scopatz Cyclus member
scopatz added a note

lambda: [1/months] = [1/sec] * [sec/month]

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
if (eps <= 1 - std::exp(-lambda_months * dt)) {
decay = true;
View
4 src/uniform_taylor.cc
@@ -17,7 +17,9 @@ Vector UniformTaylor::MatrixExpSolver(const Matrix& A, const Vector& x_o,
// checks if the dimensions of A and x_o are compatible for matrix-vector
// computations
if (x_o.NumRows() != n) {
- std::string error = "Error: Matrix-Vector dimensions are not compatible.";
+ std::string error = "Error: Matrix-Vector dimensions are not compatible: " + \
+ boost::lexical_cast<std::string>( x_o.NumRows()) + \
+ " rows vs " + boost::lexical_cast<std::string>(n) + " nuclides.";
throw ValueError(error);
}
View
2 src/uniform_taylor.h
@@ -2,6 +2,8 @@
#ifndef CYCLUS_SRC_UNIFORM_TAYLOR_H_
#define CYCLUS_SRC_UNIFORM_TAYLOR_H_
+#include <boost/lexical_cast.hpp>
+
#include "use_matrix_lib.h"
namespace cyclus {
View
6 tests/composition_tests.cc
@@ -3,7 +3,9 @@
#include <gtest/gtest.h>
#include "composition.h"
+#include "env.h"
#include "mass_table.h"
+#include "pyne.h"
class TestComp : public cyclus::Composition {
public:
@@ -46,6 +48,8 @@ TEST(CompositionTests, create_mass) {
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
TEST(CompositionTests, lineage) {
using cyclus::Composition;
+ // tell pyne about the path to nuc data
+ pyne::NUC_DATA_PATH = cyclus::Env::GetBuildPath() + "/share/cyclus_nuc_data.h5";
TestComp c;
@@ -65,4 +69,4 @@ TEST(CompositionTests, lineage) {
EXPECT_EQ(chain[3 * dt], dec4);
EXPECT_EQ(dec4, dec5);
}
-
+
Something went wrong with that request. Please try again.