tstTransitionSystemAdapter_AmpxN.cpp
./Manager/libld/tests/tstTransitionSystemAdapter_AmpxN.cpp
#include <iomanip>
#include "AmpxLoader/AmpxLoader.h"
#include "Nemesis/gtest/nemesis_gtest.hh"
#include "Origen/Core/TestPaths.h"
#include "Origen/Core/dc/FakeFactory.h"
#include "Origen/Core/fn/print.h"
#include "Origen/Core/io/DecayResourceIO.h"
#include "Origen/Core/io/LibraryIO.h"
#include "Origen/Core/io/ReactionResourceIO.h"
#include "Origen/Core/io/YieldResourceIO.h"
#include "Origen/Core/re/ReactionResource.h"
#include "Origen/Core/ts/Channel.h"
#include "Origen/Core/ts/ChannelType.h"
#include "Origen/Core/ts/Nuclide.h"
#include "Origen/Core/ts/Transition.h"
#include "Origen/Core/ts/TransitionSystem_Gen.h"
#include "Origen/Manager/libld/LibraryBuilder.h"
#include "ScaleData/Core/Utils.h"
#include "ScaleUtils/IO/DB.h"
// enable extra printing
std::shared_ptr<AmpxLibrary> loadAmpx( std::string ampxlib )
{
AmpxLoader loader;
std::shared_ptr<AmpxLibrary> idnxs;
bool success = loader.load( QString( ampxlib.c_str() ) );
if( success ) idnxs = std::shared_ptr<AmpxLibrary>( loader.getLibrary() );
return idnxs;
}
{
protected:
{
Origen::YieldResourceIO yr_io;
Origen::DecayResourceIO dr_io;
rr = std::make_shared<Origen::ReactionResource>();
// load library directly to idnxs for testing
TIMER_START( total_time );
}
{
TIMER_STOP( total_time );
TIMER_RECORD( "tst.total", total_time );
Origen::printTimingReport( std::cout );
nemesis::Timing_Diagnostics::delete_timers();
}
bool success;
nemesis::Timer total_time;
ScaleUtils::IO::DB opts;
std::shared_ptr<AmpxLibrary> idnxs;
};
const int search_mt,
const AmpxLibrary& nxs,
const std::vector<float>& flux )
{
SCOPED_TIMER( "tst.explicit_calc_og_xs" );
double xs1 = 0;
double sum_flux = 0;
for( size_t g = 0; g < flux.size(); ++g )
{
sum_flux += flux[g];
}
int num_ampx_nuclides = nxs.getNumberNuclides();
for( int j = 0; j < num_ampx_nuclides; j++ )
{
LibraryNuclide* ampx_nuclide = nxs.getNuclideAt( j );
int iza = ampx_nuclide->getId();
if( iza != search_iza ) continue;
int num_1dxs = ampx_nuclide->getNumNeutron1dData();
for( int k = 0; k < num_1dxs; k++ )
{
CrossSection1d* xs = ampx_nuclide->getNeutron1dDataAt( k );
int mt = xs->getMt();
bool partial;
int final_isomeric_state, initial_isomeric_state;
ScaleData::Utils::transform_mt(
mt, partial, final_isomeric_state, initial_isomeric_state );
if( mt == search_mt )
{
for( size_t g = 0; g < flux.size(); ++g )
{
xs1 += xs->getValues()[g] * flux[g];
}
jDebugTaggedLine(
"found iza=" << iza << " mt=" << mt << " xs1=" << xs1 );
}
}
}
if( sum_flux == 0 ) return 0;
xs1 /= sum_flux;
return xs1;
}
const int iza,
const int mt,
const AmpxLibrary& idnxs,
const std::vector<float>& mg_flux )
{
SCOPED_TIMER( "tst.check_og_xs" );
EXPECT_FALSE( field == nullptr );
EXPECT_FALSE( channeltype == nullptr );
EXPECT_FALSE( nuclide == nullptr );
Origen::SP_Channel channel = nuclide->loss( channeltype->id() );
EXPECT_FALSE( channel == nullptr );
float xs = channel->xs();
float ref_xs = explicit_calc_og_xs( iza, mt, idnxs, mg_flux );
EXPECT_FLOAT_EQ( ref_xs, xs ) << "iza=" << iza << " mt=" << mt;
return xs;
}
{
if( !success ) return;
// Import cross sections, only allowing n,gamma.
std::vector<int> allowed_mts( 1, 102 );
rr->set_allowed_mts( allowed_mts );
rr->set_max_isomeric_state( 1 );
// load reaction resource (mapping to storage_bounds)
Origen::ReactionResourceIO rr_io;
success = success && rr_io.load( *rr, Origen::jeff252_filepath, opts );
// Set reaction resource on adapter.
tsa.set_reaction_resource( rr );
// Get a constant flux.
std::vector<float> mg_flux( 252, 1.0 );
// Create a transition system.
// Update a transition system (will not actually be able to update
// anything because system is empty).
tsa.update( &ts, mg_flux );
// Now actually extend the system to include the proper nuclides
// and channels.
tsa.extend( &ts );
// get neutron field
EXPECT_FALSE( field == nullptr );
// Update and check we have some nuclides now.
tsa.update( &ts, mg_flux );
EXPECT_EQ( 1057u, ts.num_nuclides() ) << ts.to_string() << std::endl;
Origen::SP_ChannelType n_fiss = ts.channel_type( field, 18 );
EXPECT_TRUE( n_fiss == nullptr );
SCOPED_TRACE( "Testing 1g cross sections with constant flux" );
{
tsa.update( &ts, mg_flux );
float n0002003_mt102 = check_og_xs( ts, 2003, 102, *idnxs, mg_flux );
float n0008016_mt102 = check_og_xs( ts, 8016, 102, *idnxs, mg_flux );
float n0092235_mt102 = check_og_xs( ts, 92235, 102, *idnxs, mg_flux );
float n1076191_mt102 = check_og_xs( ts, 1076191, 102, *idnxs, mg_flux );
// make sure the 1g cross sections don't change when the flux
// magnitude changes
mg_flux.clear();
mg_flux.resize( 252, 2.0 );
tsa.update( &ts, mg_flux );
EXPECT_FLOAT_EQ( n0002003_mt102,
check_og_xs( ts, 2003, 102, *idnxs, mg_flux ) );
EXPECT_FLOAT_EQ( n0008016_mt102,
check_og_xs( ts, 8016, 102, *idnxs, mg_flux ) );
EXPECT_FLOAT_EQ( n0092235_mt102,
check_og_xs( ts, 92235, 102, *idnxs, mg_flux ) );
EXPECT_FLOAT_EQ( n1076191_mt102,
check_og_xs( ts, 1076191, 102, *idnxs, mg_flux ) );
}
SCOPED_TRACE( "Testing 1g cross sections with non-constant flux" );
{
mg_flux.clear();
mg_flux.resize( 100, 1.0 );
mg_flux.resize( 252, 2.0 );
tsa.update( &ts, mg_flux );
float n0002003_mt102 = check_og_xs( ts, 2003, 102, *idnxs, mg_flux );
float n0008016_mt102 = check_og_xs( ts, 8016, 102, *idnxs, mg_flux );
float n0092235_mt102 = check_og_xs( ts, 92235, 102, *idnxs, mg_flux );
float n1076191_mt102 = check_og_xs( ts, 1076191, 102, *idnxs, mg_flux );
// make sure the 1g cross sections don't change when the flux
// magnitude changes
mg_flux.clear();
mg_flux.resize( 100, 2.0 );
mg_flux.resize( 252, 4.0 );
tsa.update( &ts, mg_flux );
EXPECT_FLOAT_EQ( n0002003_mt102,
check_og_xs( ts, 2003, 102, *idnxs, mg_flux ) );
EXPECT_FLOAT_EQ( n0008016_mt102,
check_og_xs( ts, 8016, 102, *idnxs, mg_flux ) );
EXPECT_FLOAT_EQ( n0092235_mt102,
check_og_xs( ts, 92235, 102, *idnxs, mg_flux ) );
EXPECT_FLOAT_EQ( n1076191_mt102,
check_og_xs( ts, 1076191, 102, *idnxs, mg_flux ) );
}
}
TEST_F( TransitionSystemAdapterTesterJeff252, ThreeGroupStorage )
{
if( !success ) return;
// group boundaries
std::vector<float> three_group_bounds;
three_group_bounds.push_back( 20.e6 ); // 20 MeV
three_group_bounds.push_back( 2.5e4 ); // 25 keV
three_group_bounds.push_back( 4.1e0 ); // 4.1 eV
three_group_bounds.push_back( 1.e-5 ); // 1e-5 eV
rr->set_storage_bounds( three_group_bounds );
ASSERT_EQ( 4u, rr->storage_bounds().size() );
// load reaction resource (mapping to storage_bounds)
Origen::ReactionResourceIO rr_io;
success = success && rr_io.load( *rr, Origen::jeff252_filepath, opts );
// Set reaction resource on adapter.
tsa.set_reaction_resource( rr );
// Get the flux in the storage group structure (mapped from the library
// group flux so it will be consistent.)
std::vector<float> mg_flux( tsa.reaction_resource()->storage_flux() );
std::vector<float> library_flux( tsa.reaction_resource()->library_flux() );
// check sums
float ref_sum = ScaleSTL::sum( library_flux );
EXPECT_NEAR( ref_sum, ScaleSTL::sum( mg_flux ), ref_sum * 1e-5 );
// Create a transition system.
tsa.extend( &ts );
tsa.update( &ts, mg_flux );
ASSERT_EQ( 1637u, ts.num_nuclides() ) << ts.to_string() << std::endl;
SCOPED_TRACE( "Testing 1g cross sections with 3-group flux" );
{
check_og_xs( ts, 2003, 102, *idnxs, library_flux );
check_og_xs( ts, 8016, 102, *idnxs, library_flux );
check_og_xs( ts, 92235, 102, *idnxs, library_flux );
check_og_xs( ts, 1076191, 102, *idnxs, library_flux );
}
if( false )
{
std::ofstream ofs( "tsxs1.out" );
printTransitionSystemXs( ts, ofs );
}
}
TEST_F( TransitionSystemAdapterTesterJeff252, FissionYields )
{
if( !success ) return;
// load reaction resource (mapping to storage_bounds)
Origen::ReactionResourceIO rr_io;
success = success && rr_io.load( *rr, Origen::jeff252_filepath, opts );
// Create an adapter.
rr->import_fission_yields( yr );
// Set reaction resource on adapter.
tsa.set_reaction_resource( rr );
// Create a transition system.
// extend the system
tsa.extend( &ts );
// get neutron field
EXPECT_FALSE( field == nullptr );
// get all the stuff we need
Origen::SP_ChannelType fission = ts.channel_type( field, 18 );
ASSERT_TRUE( fission != nullptr );
Origen::SP_Nuclide u235 = ts.nuclide( 92235 );
ASSERT_TRUE( u235 != nullptr );
Origen::SP_Channel u235_fission = u235->loss( fission->id() );
ASSERT_TRUE( u235_fission != nullptr );
ASSERT_EQ( fission, u235_fission->channel_type() );
Origen::SP_Nuclide xe135 = ts.nuclide( 54135 );
ASSERT_TRUE( xe135 != nullptr );
Origen::SP_Transition t_xe135 = u235_fission->transition( xe135->id() );
ASSERT_TRUE( t_xe135 != nullptr );
Origen::SP_Nuclide i135 = ts.nuclide( 53135 );
ASSERT_TRUE( i135 != nullptr );
Origen::SP_Transition t_i135 = u235_fission->transition( i135->id() );
ASSERT_TRUE( t_i135 != nullptr );
Origen::SP_Nuclide nd149 = ts.nuclide( 60149 );
ASSERT_TRUE( nd149 != nullptr );
Origen::SP_Transition t_nd149 = u235_fission->transition( nd149->id() );
ASSERT_TRUE( t_nd149 != nullptr );
Origen::SP_Nuclide pm149 = ts.nuclide( 61149 );
ASSERT_TRUE( pm149 != nullptr );
Origen::SP_Transition t_pm149 = u235_fission->transition( pm149->id() );
ASSERT_TRUE( t_pm149 != nullptr );
Origen::SP_Nuclide sm149 = ts.nuclide( 62149 );
ASSERT_TRUE( sm149 != nullptr );
Origen::SP_Transition t_sm149 = u235_fission->transition( sm149->id() );
ASSERT_TRUE( t_sm149 != nullptr );
// setup a zero flux vector
std::vector<float> mg_flux( 252, 1e-20 );
tsa.update( &ts, mg_flux );
if( print )
{
std::cout << u235_fission->to_string() << std::endl;
}
// perturb each group and check results for some important
// yields (uncomment std::cout to make a table)
if( print )
{
std::cout << "g"
<< " "
<< "yld_xe135"
<< " "
<< "yld_i135"
<< " "
<< "yld_nd149"
<< " "
<< "yld_pm149"
<< " "
<< "yld_sm149" << std::scientific << std::setprecision( 4 )
<< std::endl;
}
size_t ngroups = 5;
for( size_t g = 1; g <= ngroups; ++g )
{
std::stringstream time;
time << "tst.yield_trial(" << std::setw( 3 ) << std::setfill( '0' ) << g
<< ")";
SCOPED_TIMER( time.str() );
mg_flux[g - 1] = 1e20; // put a spike in this group
tsa.update( &ts, mg_flux );
double yld_xe135 = t_xe135->yield();
double yld_i135 = t_i135->yield();
double yld_nd149 = t_nd149->yield();
double yld_pm149 = t_pm149->yield();
double yld_sm149 = t_sm149->yield();
if( print )
{
std::cout << g << " " << yld_xe135 << " " << yld_i135 << " "
<< yld_nd149 << " " << yld_pm149 << " " << yld_sm149
<< std::endl;
}
EXPECT_LT( 2.0e-2, yld_i135 ) << "I-135 yield should be >2% for group "
<< g;
EXPECT_GT( 4.0e-2, yld_i135 ) << "I-135 yield should be <4% for group "
<< g;
EXPECT_LT( 1.0e-3, yld_xe135 )
<< "Xe-135 yield should be >0.1% for group " << g;
EXPECT_GT( 6.0e-3, yld_xe135 )
<< "Xe-135 yield should be <0.6% for group " << g;
EXPECT_LT( 1.0e-5, yld_nd149 )
<< "Nd-149 yield should be >0.001% for group " << g;
EXPECT_GT( 7.0e-4, yld_nd149 )
<< "Nd-149 yield should be <0.07% for group " << g;
EXPECT_LT( 1.0e-8, yld_pm149 )
<< "Pm-149 yield should be >1e-8 for group " << g;
EXPECT_GT( 1.4e-5, yld_pm149 )
<< "Pm-149 yield should be <1.4e-5 for group " << g;
EXPECT_LT( 1.0e-13, yld_sm149 )
<< "Sm-149 yield should be >1e-13 for group " << g;
EXPECT_GT( 3.0e-8, yld_sm149 )
<< "Sm-149 yield should be <3e-8 for group " << g;
mg_flux[g - 1] = 1e-20;
}
if( false )
{
std::ofstream ofs( "tsxs2.out" );
printTransitionSystemXs( ts, ofs );
}
}
TEST_F( TransitionSystemAdapterTesterJeff252, LibraryBuilder )
{
if( !success ) return;
// Create a transition system.
// variables for Origen library builder
std::vector<int> num_decay_parents;
std::vector<int> num_parents;
std::vector<int> parent_positions;
std::vector<int> transition_ids;
std::vector<double> decay_constants;
std::vector<double> fission_xs;
std::vector<double> loss_xs;
std::vector<double> transition_coeff;
std::vector<float> n_production_xs;
std::vector<float> flux;
std::vector<std::vector<Origen::Transition::WP>> reaction_gain_sequence;
// get nuclides
Origen::Vec_Int sizzzaaa_list =
// Create a decay adapter and extend the transition system with decay data
// Build basic decay library from extended transition system
{
SCOPED_TIMER( "create_decay_library" );
bld.create_decay_transitions( &num_decay_parents,
&num_parents,
&parent_positions,
&transition_ids,
&decay_constants,
&transition_coeff,
sizzzaaa_list,
true,
ts );
}
// include EVERY reaction
bld.set_minimum_reaction_coeff( 0.0 );
// Create a reaction adapter.
Origen::ReactionResourceIO rr_io;
success = success && rr_io.load( *rr, Origen::jeff252_filepath, opts );
Origen::TransitionSystemAdapter_AmpxN tsa_reac;
rr->import_fission_yields( yr );
tsa_reac.set_reaction_resource( rr );
tsa_reac.extend( &ts );
std::vector<float> mg_flux = rr->library_flux();
tsa_reac.update( &ts, mg_flux );
// build reaction transitions
{
SCOPED_TIMER( "create_reaction_transitions" );
bld.create_reaction_transitions( &num_decay_parents,
&num_parents,
&parent_positions,
&transition_ids,
&transition_coeff,
&fission_xs,
&loss_xs,
&n_production_xs,
&flux,
&reaction_gain_sequence,
sizzzaaa_list,
true,
ts );
}
EXPECT_EQ( 64079u, reaction_gain_sequence.size() );
// Create this library.
Origen::Library lib;
ScaleUtils::IO::DB opts;
bld.create_library( &lib,
num_decay_parents,
num_parents,
parent_positions,
transition_ids,
decay_constants,
fission_xs,
loss_xs,
transition_coeff,
n_production_xs,
flux,
sizzzaaa_list );
// check some n_production_xs/loss_xs (like nubar)
std::list<std::pair<int, float>> reference_kinf = {
{1014, 0.00779756}, // 20092238
{1010, 1.92759}, // 20092235
{1068, 3.41061}, // 20097248
{1069, 3.41061}, // 21097248
{1083, 4.21324}, // 20099252
{1032, 1.83269}, // 20094239
{1034, 2.16627}, // 20094241
};
for( auto index_kinf : reference_kinf )
{
auto index = index_kinf.first;
auto kinf = index_kinf.second;
EXPECT_NEAR(
kinf, n_production_xs[index] / loss_xs[index], kinf * 1.e-4 );
}
// std::cerr<<opts<<std::endl;
// do this loop some times to simulate the real timing breakdown
const bool debug = false;
{
SCOPED_TIMER( "update_reaction_transitions" );
{
mg_flux[7] = 90.0;
tsa_reac.update( &ts, mg_flux );
size_t tind = 0; // transition index (includes decay slots)
size_t rind = 0; // reaction index
for( size_t nind = 0; nind < sizzzaaa_list.size(); ++nind )
{
tind += num_decay_parents[nind];
for( int j = num_decay_parents[nind]; j < num_parents[nind];
++j, ++rind, ++tind )
{
auto gain = reaction_gain_sequence[rind][0];
if( gain == nullptr )
{
EXPECT_FALSE( gain == nullptr )
<< "must have non-null transition for sizzzaaa="
<< sizzzaaa_list[nind]
<< " mt=" << transition_ids[tind] << std::endl;
continue;
}
else
{
transition_coeff[tind] = gain->coeff();
}
if( debug )
{
std::cout << std::setw( 9 ) << tind << std::setw( 8 )
<< std::fixed << gain->product()->id()
<< std::setw( 4 ) << std::fixed << "<--"
<< std::setw( 8 ) << std::fixed
<< gain->channel()->parent()->id()
<< std::setw( 6 ) << std::fixed
<< std::boolalpha << gain->is_byproduct()
<< std::setw( 6 ) << std::fixed
<< gain->channel()->channel_type()->id()
<< std::setw( 13 ) << std::scientific
<< gain->coeff() << std::endl;
}
}
}
}
}
}
std::string prefix = "ampx" )
{
auto sp_nxs = loadAmpx( ampxlib );
AmpxLibrary& nxs = *sp_nxs;
int num_ampx_nuclides = nxs.getNumberNuclides();
for( int j = 0; j < num_ampx_nuclides; j++ )
{
LibraryNuclide* ampx_nuclide = nxs.getNuclideAt( j );
std::cout << ampx_nuclide->getThermal2DTemp() << std::endl;
int iza = ampx_nuclide->getId();
int num_1dxs = ampx_nuclide->getNumNeutron1dData();
for( int k = 0; k < num_1dxs; k++ )
{
CrossSection1d* xs = ampx_nuclide->getNeutron1dDataAt( k );
int mt = xs->getMt();
bool partial;
int final_isomeric_state, initial_isomeric_state;
ScaleData::Utils::transform_mt(
mt, partial, final_isomeric_state, initial_isomeric_state );
std::stringstream ss;
ss << prefix << "/" << iza << "/" << mt;
if( partial )
{
ss << "/" << final_isomeric_state;
}
std::stringstream dir;
dir << "mkdir -p " << ss.str();
system( dir.str().c_str() );
ss << "/xs.txt";
std::ofstream fs( ss.str() );
auto bounds = nxs.getNeutronEnergyBounds();
for( size_t g = 0; g < bounds->getBoundsSize() - 1; ++g )
{
float eavg = 0.5 * ( bounds->getBounds()[g] +
bounds->getBounds()[g + 1] );
fs << eavg << " " << xs->getValues()[g] << std::endl;
}
fs.close();
}
}
}
TEST( Output, AmpxLib )
{
// output_ampx_library_text_files("/Users/ww5/lib/scale_dev_data/scale.rev03.xn56v7.1","ampx56v7.1");
}
#define MIRROR_VECTOR( PARENT1 , VAR1, PARENT2, ADD )\
auto new##VAR1 = PARENT1.VAR1;\
auto old##VAR1 = PARENT2.VAR1;\
for(size_t i=0; i<old##VAR1.size(); ++i )\
{\
new##VAR1.push_back( ADD + old##VAR1[i] );\
}\
PARENT1.VAR1 = new##VAR1
TEST( Output, DISABLED_DoubleMatrix )
{
// Here's what I did to make this work.
//
// First remove the DISABLED_ above so this test runs.
//
// You'll need these SCALE data files locally.
//
// cp ~/dev/scale/data/origen_library/pwr.rev03.orglib transition.def
// cp ~/dev/scale/data/origen.rev04.end7dec .
// rm -f double.f33 && \
// make -j4 && \
// ./OrigenManager_tstTransitionSystemAdapter_AmpxN.exe \
// --gtest_filter=Output.DoubleMatrix && \
// ../../../Core/obiwan view double.f33
//
// The view of the library looks good. But ORIGEN will not run
// the library because of 3 sublib hardcoding.
//
// Consider the following simple input.
// =shell
// cp ${RTNDIR}/double.f33 double
// end
// =origen
// case{
// lib{
// file="double"
// }
// time=[100] %days
// save{
// file="double.f71" steps=ALL
// }
// }
// end
//
// It bombs on:
//
// terminate called after throwing an instance of 'nemesis::assertion'
// what(): Assertion: ilite + iact + ifp == itot,
// failed in /Users/ww5/dev/scale/tip/packages/Origen/Core/dc/TransitionMatrixP.cpp:161
//
// Looking through the code, there are a lot of other places where we're
// going to hit these problems.
//
//
// I will add data in place to this library.
::Origen::Library lib;
{
ScaleUtils::IO::DB opts;
io.load(lib,"transition.def",opts);
}
// This is the decay-only data library.
// Alternatively, because the lib has decay data on it,
// I could have used that but it would have looked a little more
// convoluted.
{
ScaleUtils::IO::DB opts;
io.load(dec,"origen.rev04.end7dec",opts);
}
//
// The existing library has 3 main parts that we need
// to append to create the new library.
//
// 1. The nuclide set needs to be extended to include
// clones of sublibs 1-3 nuclides as sublib 4-6 and
// copies of all the decay data.
// 2. The transition structure needs to be extended to
// include the decay from sublib S to S+3 as well
// as all the decay physics between everything in
// sublib 4-6.
// 3. The transition coefficients need to be updated to
// include the decay constant from S --> S+3 in the
// total for a nuclide in sublib 1-3 AND in the
// transition coeff (gain) fro nuclides in
// sublib 4-6.
//
// STEP 1a: nuclide set
size_t num_nuclides;
{
num_nuclides = ids2.size();
for(size_t i=0; i<num_nuclides; ++i )
{
ids.push_back( 3*ScaleData::Utils::PIZZZAAA_P_OFFSET + ids2[i] );
}
nucset.set_ids( ids );
def.set_nuclide_set( nucset );
lib.set_definition( def );
}
// STEP 1b: copy decay data
// NOTE: if you want to add a decay from S-->S+3, you'd have to
// do it here on the vector decay_constants.
{
auto dd2 = dec.decay_data();
MIRROR_VECTOR( dd , decay_constants, dd2, 0 );
MIRROR_VECTOR( dd , recoverable_energy_values, dd2, 0 );
MIRROR_VECTOR( dd , q_fractions_from_photons, dd2, 0 );
MIRROR_VECTOR( dd , rcg_for_air, dd2, 0 );
MIRROR_VECTOR( dd , rcg_for_water, dd2, 0 );
MIRROR_VECTOR( dd , abund_of_lite_nuclides, dd2, 0 );
MIRROR_VECTOR( dd , masses, dd2, 0 );
lib.set_decay_data( dd );
}
// STEP 2 & 3: add extra transition structure and coefficients.
{
auto new_transition_ids = tr.transition_ids();
auto new_num_parents = tr.num_parents();
auto new_num_decay_parents = tr.num_decay_parents();
auto new_loss_xs = tc.loss_xs();
auto new_kappa_capture = tc.kappa_capture();
auto tr2 = dec.transition_structure();
auto old_parent_positions = tr2.parent_positions();
auto old_transition_ids = tr2.transition_ids();
auto old_num_parents = tr2.num_parents();
auto old_num_decay_parents = tr2.num_decay_parents();
auto tc2 = dec.transition_coeff_at( 0 );
auto old_matrix = tc2.matrix();
size_t t=0;
for(size_t i=0; i<num_nuclides; ++i)
{
// This is the transition (gain) from S-->S+3:
// We are *AT* the daughter and are adding a back-reference
// gain from the parent at i+1.
new_parent_positions.push_back( i+1 );
//dummy value
new_transition_ids.push_back( 0 );
//set a nonzero value to prove that we generate these isotopes
//copy over other transition data
for(size_t j=0; j<old_num_parents[i]; ++j)
{
new_parent_positions.push_back( num_nuclides+old_parent_positions[t] );
new_transition_ids.push_back( old_transition_ids[t] );
new_matrix.push_back( old_matrix[t] );
++t;
}
//record parents
new_num_parents.push_back( old_num_parents[i]+1 );
new_num_decay_parents.push_back( old_num_decay_parents[i]+1 );
//append 0 to total loss and energy generated so vectors are
//the correct size
new_loss_xs.push_back( 0.0 );
new_kappa_capture.push_back( 0.0 );
}
//overwrite old transition structure with the new
tr.set_parent_positions(new_parent_positions);
tr.set_transition_ids(new_transition_ids);
tr.set_num_parents(new_num_parents);
tr.set_num_decay_parents(new_num_decay_parents);
lib.set_transition_structure( tr );
//overwrite old transition coeff with the new
tc.set_matrix( new_matrix );
tc.set_kappa_capture( new_kappa_capture );
tc.set_loss_xs( new_loss_xs );
lib.set_transition_coeff_at( tc, 0 );
}
//finally ready to save data
{
ScaleUtils::IO::DB opts;
io.save(lib,"double.f33",opts);
}
}
