tstTransitionMatrixP.cpp
./Core/dc/tests/tstTransitionMatrixP.cpp
#include <iostream>
#include <vector>
#include "Nemesis/gtest/nemesis_gtest.hh"
#include "Nemesis/harness/DBC.hh"
#include "Origen/Core/dc/Case.h"
#include "Origen/Core/dc/FakeFactory.h"
#include "Origen/Core/dc/Library.h"
#include "Origen/Core/dc/TransitionMatrixP.h"
#include "Standard/Interface/jdebug.h"
using namespace std;
using namespace Origen;
/*** Fold and extract a dense matrix from the TransitionMatrix ***/
double flux )
{
IntegerList* non0 = trx.get_non0_();
IntegerList* loc = trx.get_loc();
DoubleList a, d;
trx.fold_flux( flux, d, a );
std::vector<std::vector<double>> A( itot,
std::vector<double>( itot, 0.0 ) );
// Add diagonal elements
for( int i = 0; i < itot; ++i ) A[i][i] += d[i];
// Add offdiagonal elements
int j = 0;
for( int i = 0; i < itot; ++i ) // rows (daughters)
for( int p = 0; p < ( *non0 )[i]; ++p )
{
A[i][( *loc )[j]] += a[j];
j++;
}
return A;
}
{
TransitionMatrixP trx;
const int itot_loc = 11;
const int non_loc = 22;
trx.set_non( non_loc );
trx.set_itot( itot_loc );
trx.set_ilite( 5 );
trx.set_iact( 4 );
trx.set_ifp( 2 );
// cout<< "Expected value of itot is "<< itot_loc << ", value in
// TransitionMatrixP is "<< trx.get_itot() << endl;
Insist( trx.get_itot() == itot_loc,
"TransitionMatrixP.set_itot failed to set itot!" );
vector<float> dis( trx.get_itot(), 1E-2 );
vector<double> tocap( trx.get_itot(), 1E-3 );
vector<double> fiss( trx.get_iact(), 5E-4 );
vector<float> genneu( trx.get_itot(), 1E-5 );
vector<int> kd_( trx.get_itot(), 1 );
vector<int> non0_( trx.get_itot(), 2 );
vector<int> loc( trx.get_non() );
// set up a stupid structure where all parents are the
// next nuclide
int k = 0;
{
for( int j = 0; j < non0_[i]; ++j )
{
Require( k < (int)loc.size() );
loc.at( k ) = i + 1;
if( ( i + 1 ) == trx.get_itot() ) loc.at( k ) = 0;
++k;
}
}
vector<double> a( trx.get_non(), 2E-3 );
vector<int> mt( trx.get_non(), 4 );
mt[8] = -1;
mt[9] = 16;
vector<int> nucl( trx.get_itot() );
vector<int> nuc_typ( trx.get_itot() );
vector<float> abund( trx.get_ilite(), 1.0 );
nucl[0] = 10010;
nuc_typ[0] = 1;
nucl[1] = 10020;
nuc_typ[1] = 1;
nucl[2] = 561260;
nuc_typ[2] = 1;
nucl[3] = 561271;
nuc_typ[3] = 1;
nucl[4] = 651400;
nuc_typ[4] = 1;
nucl[5] = 922340;
nuc_typ[5] = 2;
nucl[6] = 922350;
nuc_typ[6] = 2;
nucl[7] = 942390;
nuc_typ[7] = 2;
nucl[8] = 952421;
nuc_typ[8] = 2;
nucl[9] = 341350;
nuc_typ[9] = 3;
nucl[10] = 371360;
nuc_typ[10] = 3;
trx.set_nucl( nucl );
trx.set_dis( dis );
trx.set_tocap( tocap );
trx.set_fiss( fiss );
trx.set_loc( loc );
trx.set_a( a );
trx.set_mt( mt );
trx.set_kd_( kd_ );
trx.set_non0_( non0_ );
trx.set_abund( abund );
trx.finish_init();
// checking get_transition_ids and filter_reaction_mts
{
int nind = 5;
Vec_Int tids;
trx.get_transition_ids( nind, tids );
for( size_t i = 0; i < tids.size(); ++i )
{
jDebugTaggedLine( "tids[" << i << "]=" << tids[i] << std::endl );
}
Insist( tids.size() == 2,
"TransitionMatrixP::get_transition_ids.size() not correct" );
Insist( tids[0] == -1, "TransitionMatrixP::tids[0]!=-1" );
Insist( tids[1] == 16, "TransitionMatrixP::tids[1]!=16" );
Vec_Int mts;
trx.filter_reaction_mts( tids.begin(), tids.end(), mts );
Insist( mts.size() == 1, "only one mts" );
Insist( mts[0] == 16, "TransitionMatrixP::mts[0]!=16" );
trx.get_possible_transition_ids( tids );
Insist( tids.size() == 3, "Only 3 possibilities" );
Insist( tids[0] == -1, "first is -1" );
Insist( tids[1] == 4, "second is 4" );
Insist( tids[2] == 16, "third is 16" );
trx.filter_reaction_mts( tids.begin(), tids.end(), mts );
Insist( mts.size() == 2, "only two mts" );
Insist( mts[0] == 4, "TransitionMatrixP::mts[0]!=4" );
Insist( mts[1] == 16, "TransitionMatrixP::mts[1]!=16" );
}
{
// cout<< "Expected value of nuc_typ[" << i << "] is "<< nuc_typ[i]
// << ", value in TransitionMatrixP is "<< nuc_typ_internal->at(i)
// << endl;
Insist( ( *nuc_typ_internal )[i] == nuc_typ[i],
"TransitionMatrixP sublibs incorrect!" );
// cout<< "Expected value of dis[" << i << "] is "<< dis[i] << ",
// value in TransitionMatrixP is "<< disRet->at(i) << endl;
Insist( dis[i] == disRet->at( i ),
"TransitionMatrixP.set_dis failed to set dis!" );
// cout<< "Expected value of tocap[" << i << "] is "<< tocap[i] <<
// ", value in TransitionMatrixP is "<< tocapRet->at(i) << endl;
Insist( tocap[i] == tocapRet->at( i ),
"TransitionMatrixP.set_tocap failed to set tocap!" );
}
{
// cout<< "Expected value of fiss[" << i << "] is "<< fiss[i] << ",
// value in TransitionMatrixP is "<< fissRet->at(i) << endl;
Insist( fiss[i] == fissRet->at( i ),
"TransitionMatrixP.set_fiss failed to set fiss!" );
}
const IntegerList* locRet = trx.get_loc();
{
// cout<< "Expected value of a[" << i << "] is "<< a[i] << ", value
// in TransitionMatrixP is "<< aRet->at(i) << endl;
Insist( a[i] == aRet->at( i ),
"TransitionMatrixP.set_a failed to set a!" );
// cout<< "Expected value of mt[" << i << "] is "<< mt[i] << ",
// value in TransitionMatrixP is "<< mtRet->at(i) << endl;
Insist( mt[i] == mtRet->at( i ),
"TransitionMatrixP.set_mt failed to set a!" );
// cout<< "Expected value of loc[" << i << "] is "<< loc[i] << ",
// value in TransitionMatrixP is "<< locRet->at(i) << endl;
Insist( loc[i] == locRet->at( i ),
"TransitionMatrixP.set_loc failed to set a!" );
}
trx.check_library();
vector<std::string> errors;
if( !pass )
{
for( size_t j = 0; j < errors.size(); ++j )
{
std::cerr << j << ". " << errors[j] << std::endl;
}
}
EXPECT_TRUE( pass );
}
TEST( TransitionMatrixP, TallyNuclides )
{
TransitionMatrixP trx;
const int itot_loc = 11;
const int non_loc = 22;
trx.set_non( non_loc );
trx.set_itot( itot_loc );
trx.set_ilite( 5 );
trx.set_iact( 4 );
trx.set_ifp( 2 );
vector<float> dis( trx.get_itot(), 1E-2 );
vector<double> tocap( trx.get_itot(), 1E-3 );
vector<double> fiss( trx.get_iact(), 5E-4 );
vector<float> genneu( trx.get_itot(), 1E-5 );
vector<int> kd_( trx.get_itot(), 1 );
vector<int> non0_( trx.get_itot(), 2 );
vector<int> loc( trx.get_non() );
// set up a stupid structure where all parents are the
// next nuclide
int k = 0;
{
for( int j = 0; j < non0_[i]; ++j )
{
Require( k < (int)loc.size() );
loc.at( k ) = i + 1;
if( ( i + 1 ) == trx.get_itot() ) loc.at( k ) = 0;
++k;
}
}
vector<double> a( trx.get_non(), 2E-3 );
vector<int> mt( trx.get_non(), 4 );
mt[8] = -1;
mt[9] = 16;
vector<int> nucl( trx.get_itot() );
vector<int> nuc_typ( trx.get_itot() );
vector<float> abund( trx.get_ilite(), 1.0 );
nucl[0] = 10010;
nucl[1] = 10020;
nucl[2] = 561260;
nucl[3] = 561271;
nucl[4] = 651400;
nucl[5] = 922340;
nucl[6] = 922350;
nucl[7] = 942390;
nucl[8] = 952421;
nucl[9] = 341350;
nucl[10] = 371360;
trx.set_nucl( nucl );
trx.set_dis( dis );
trx.set_tocap( tocap );
trx.set_fiss( fiss );
trx.set_loc( loc );
trx.set_a( a );
trx.set_mt( mt );
trx.set_kd_( kd_ );
trx.set_non0_( non0_ );
trx.set_abund( abund );
// Add kappa values to test energy tally
std::vector<double> kappa_capt = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11};
std::vector<double> kappa_fiss = {100, 110, 120, 130};
FloatList kappa_capt_fl;
for( int i = 0; i < kappa_capt.size(); i++ )
kappa_capt_fl.push_back( kappa_capt[i] );
trx.set_kappa_capture( kappa_capt_fl );
FloatList kappa_fiss_fl;
for( int i = 0; i < kappa_fiss.size(); i++ )
kappa_fiss_fl.push_back( kappa_fiss[i] );
trx.set_kappa_fission( kappa_fiss_fl );
trx.finish_init();
// Get the unmodified matrix in dense form
double flux = 1.0e14;
std::vector<std::vector<double>> A1 = GetDenseMatrix( trx, flux );
// >>> ADD AND CHECK TALLY NUCLIDE
std::vector<double> tally_dw = {1, 1e-10, 2, 20};
std::vector<int> tally_di = {0, 1, 2, 10};
std::vector<double> tally_rw = {5, 10};
std::vector<int> tally_ri = {1, 3};
trx.add_tally_nuclide( tally_dw, tally_di, tally_rw, tally_ri );
std::vector<std::vector<double>> A2 = GetDenseMatrix( trx, flux );
// Check that the physical matrix was not modified
for( int i = 0; i < A1.size(); ++i )
for( int j = 0; j < A1[0].size(); ++j ) EXPECT_EQ( A1[i][j], A2[i][j] );
// This is what we would expect to get (the flux is converted to bcm)
std::vector<double> t2 = {1,
5 * flux * 1e-24 + 1e-10,
2,
10 * flux * 1e-24,
0,
0,
0,
0,
0,
0,
20,
0};
EXPECT_EQ( 12, A2.size() );
EXPECT_EQ( 12, A2[11].size() );
EXPECT_VEC_EQ( t2, A2[11] );
// >>> ADD AND CHECK ENERGY RELEASED TALLY
trx.add_energy_tally( flux );
std::vector<std::vector<double>> A3 = GetDenseMatrix( trx, flux );
// Check that the physical matrix or the other tallies were not modified
for( int i = 0; i < A2.size(); ++i )
for( int j = 0; j < A2[0].size(); ++j ) EXPECT_EQ( A2[i][j], A3[i][j] );
// Check the energy tally
std::vector<double> t3 = {
1E-3 * flux / 1e24 * kappa_capt[0],
1E-3 * flux / 1e24 * kappa_capt[1],
1E-3 * flux / 1e24 * kappa_capt[2],
1E-3 * flux / 1e24 * kappa_capt[3],
1E-3 * flux / 1e24 * kappa_capt[4],
1E-3 * flux / 1e24 * kappa_capt[5] +
5E-4 * flux / 1e24 * ( kappa_fiss[0] - kappa_capt[5] ),
1E-3 * flux / 1e24 * kappa_capt[6] +
5E-4 * flux / 1e24 * ( kappa_fiss[1] - kappa_capt[6] ),
1E-3 * flux / 1e24 * kappa_capt[7] +
5E-4 * flux / 1e24 * ( kappa_fiss[2] - kappa_capt[7] ),
1E-3 * flux / 1e24 * kappa_capt[8] +
5E-4 * flux / 1e24 * ( kappa_fiss[3] - kappa_capt[8] ),
1E-3 * flux / 1e24 * kappa_capt[9],
1E-3 * flux / 1e24 * kappa_capt[10],
0,
0};
EXPECT_EQ( 13, A3.size() );
EXPECT_EQ( 13, A3[12].size() );
for( int i = 0; i < t3.size(); ++i )
if( std::abs( t3[i] - A3[12][i] ) > t3[i] * 1e-6 )
{
printf( "Mismatch in element i=%d of power tally\n", i );
EXPECT_FLOAT_EQ( t3[i], A3[12][i] );
}
// >>> ADD AND CHECK NUMBER OF NEUTRON INDUCED FISSIONS TALLY
trx.add_num_fissions_tally();
std::vector<std::vector<double>> A4 = GetDenseMatrix( trx, flux );
// Check that the physical matrix or the other tallies were not modified
for( int i = 0; i < A3.size(); ++i )
for( int j = 0; j < A3[0].size(); ++j ) EXPECT_EQ( A3[i][j], A4[i][j] );
// Check the fission tally
std::vector<double> t4 = {0,
0,
0,
0,
0,
5E-4 * flux / 1e24,
5E-4 * flux / 1e24,
5E-4 * flux / 1e24,
5E-4 * flux / 1e24,
0,
0,
0,
0,
0};
EXPECT_EQ( 14, A4.size() );
EXPECT_EQ( 14, A4[13].size() );
for( int i = 0; i < t4.size(); ++i )
if( std::abs( t4[i] - A4[13][i] ) > t4[i] * 1e-6 )
{
printf( "Mismatch in element i=%d of power tally\n", i );
EXPECT_FLOAT_EQ( t4[i], A4[12][i] );
}
// >>> ADD AND CHECK NUMBER OF NEUTRON CAPTURES TALLY
trx.add_num_captures_tally();
std::vector<std::vector<double>> A5 = GetDenseMatrix( trx, flux );
// Check that the physical matrix or the other tallies were not modified
for( int i = 0; i < A4.size(); ++i )
for( int j = 0; j < A4[0].size(); ++j ) EXPECT_EQ( A4[i][j], A5[i][j] );
// Check the fission tally
std::vector<double> t5( 11, 1E-3 * flux / 1e24 );
t5.push_back( 0 );
t5.push_back( 0 );
t5.push_back( 0 );
t5.push_back( 0 );
EXPECT_EQ( 15, A5.size() );
EXPECT_EQ( 15, A5[14].size() );
for( int i = 0; i < t5.size(); ++i )
if( std::abs( t5[i] - A5[14][i] ) > t5[i] * 1e-6 )
{
printf( "Mismatch in element i=%d of power tally\n", i );
EXPECT_FLOAT_EQ( t5[i], A5[12][i] );
}
// >>> ADD CONCENTRATION INTEGRAL TALLIES
std::vector<int> tally_zaids;
trx.add_concentration_tallies( tally_zaids );
std::vector<std::vector<double>> A6 = GetDenseMatrix( trx, flux );
// Check that the physical matrix or the other tallies were not modified
for( int i = 0; i < A5.size(); ++i )
for( int j = 0; j < A5[0].size(); ++j ) EXPECT_EQ( A5[i][j], A6[i][j] );
EXPECT_VEC_EQ( nucl, tally_zaids );
EXPECT_EQ( 15 + itot_loc, A6.size() );
for( int i = 0; i < itot_loc; i++ )
{
// Check the i:th concentration tally
std::vector<double> t6( 15 + itot_loc, 0.0 );
t6[i] = 1;
EXPECT_EQ( 15 + itot_loc, A6[15].size() );
for( int j = 0; j < t6.size(); ++j )
if( A6[15 + i][j] != t6[j] )
{
printf( "Mismatch in element j=%d of the i=%d:th conc. tally\n",
j,
i );
EXPECT_EQ( t6[j], A6[15 + i][j] );
}
}
}
// use print_mts_iterator to get new reference values
{
for( size_t k = 0; k < mts.size(); k++ )
{
cout << mts[k];
if( k < ( mts.size() - 1 ) )
{
cout << ",";
}
}
return;
}
TEST( TransitionMatrixP, XsChanger )
{
Origen::Library lib;
TransitionMatrixP::SP trx = lib.newsp_transition_matrix_at( 0 );
// Test these MTs
IntegerList mts;
trx->find_all_mts( mts );
// print_mts_iterator(mts);
const int mts_reference[23] = {16, 17, 18, 22, 23, 24, 25, 28,
29, 32, 33, 34, 37, 51, 102, 103,
104, 105, 106, 107, 108, 111, 112};
for( size_t k = 0; k < mts.size(); k++ )
{
EXPECT_EQ( mts_reference[k], mts[k] ) << "mt values on library changed";
}
// fget the ninds of u5 and u8
IntegerList* nucl = trx->get_nucl();
int u235 = trx->find_nuclide_guess( 922350 );
int u238 = trx->find_nuclide_guess( 922380 );
EXPECT_EQ( 1008, u235 );
EXPECT_EQ( 1012, u238 );
// Check parents.
const int pmts_reference[9] = {-1, -1, -1, 102, 16, 17, 37, 103, 107};
const int nucl_reference[9] = {
922351, 932350, 942390, 922340, 922360, 922370, 922380, 932350, 942380};
IntegerList parents, pmts, tindsp;
trx->find_parents( u235, parents, pmts, tindsp );
EXPECT_EQ( 9, pmts.size() );
for( int k = 0; k < 9; k++ )
{
EXPECT_EQ( pmts_reference[k], pmts[k] ) << " at k=" << k;
EXPECT_EQ( nucl_reference[k], nucl->at( parents[k] ) ) << " at k=" << k;
}
{
IntegerList dds, ddm, tdd;
trx->find_daughters( 960, dds, ddm, tdd );
}
// Check daughters.
const int dmts_reference[9] = {-1, 107, -1, 107, 103, 37, 17, 16, 51};
const int nucl_reference2[9] = {
20040, 20040, 902310, 902320, 912350, 922320, 922330, 922340, 922351};
IntegerList daughters, dmts, tindsd;
trx->find_daughters( u235, daughters, dmts, tindsd );
for( int k = 0; k < 9; k++ )
{
EXPECT_EQ( dmts_reference[k], dmts[k] ) << " at k=" << k;
EXPECT_EQ( nucl_reference2[k], nucl->at( daughters[k] ) ) << " at k="
<< k;
}
for( size_t k = 10; k < dmts.size(); k++ )
{
EXPECT_EQ( 18, dmts[k] );
}
trx->find_daughters( ipb203, daughters, dmts, tindsd );
// find n,gamma in Pb203
size_t j;
bool found = false;
for( j = 0; j < dmts.size(); ++j )
{
if( dmts[j] == 102 )
{
found = true;
break;
}
}
EXPECT_TRUE( found );
if( found )
{
EXPECT_EQ( 102, dmts[j] );
}
double signg = trx->get_xs( ipb203, 102 );
EXPECT_NEAR( 1.0, signg / 0.78850060701370239, 1e-5 );
// find Hg-199m daughter
found = false;
for( j = 0; j < daughters.size(); ++j )
{
if( nucl->at( daughters[j] ) == 801991 )
{
found = true;
break;
}
}
EXPECT_TRUE( found );
EXPECT_EQ( 22, dmts[j] );
// Check fast transition access with just u-235 and u-238
{
IntegerList nind_list, transition_ids;
Integer2dList tinds_2d;
DoubleList new_xs;
nind_list.push_back( u235 );
transition_ids.push_back( 18 );
new_xs.push_back( 34.3434 );
nind_list.push_back( u238 );
transition_ids.push_back( 102 );
new_xs.push_back( 6.0102 );
trx->setup_fast_transition_access(
nind_list, transition_ids, tinds_2d );
EXPECT_LT( 0, tinds_2d.size() );
trx->fast_set_xs( nind_list, transition_ids, tinds_2d, new_xs );
tinds_2d.deepDelete(); // has pointers inside
}
// Check fast transition access with all actinides
// increasing them by 10%
{
IntegerList nind_list, transition_ids;
Integer2dList tinds_2d;
DoubleList new_xs;
int lb = trx->get_ilite();
int ub = lb + trx->get_iact();
for( int n = lb; n < ub; n++ )
{
// fission xs
nind_list.push_back( n );
transition_ids.push_back( 18 );
double old_fiss = trx->get_xs( n, 18 );
new_xs.push_back( old_fiss );
// capture xs
nind_list.push_back( n );
transition_ids.push_back( 102 );
double old_ngam = trx->get_xs( n, 102 );
new_xs.push_back( old_ngam );
}
// Increase old xs by 10% to test change
DoubleList old_xs;
for( size_t n = 0; n < new_xs.size(); n++ )
{
old_xs.push_back( new_xs[n] );
new_xs[n] *= 1.1;
}
// SLOW ACCESS
for( size_t k = 0; k < nind_list.size(); k++ )
{
trx->set_xs( nind_list[k], transition_ids[k], new_xs[k] );
}
// check that the change took place
double target1 = 1.1;
for( size_t k = 0; k < nind_list.size(); k++ )
{
double myold = old_xs[k];
if( myold > 0 )
{
double mynew = trx->get_xs( nind_list[k], transition_ids[k] );
EXPECT_NEAR( 0.0, std::abs( mynew / myold - target1 ), 1e-5 );
}
}
// FAST SETUP
trx->setup_fast_transition_access(
nind_list, transition_ids, tinds_2d );
// Increase by 20% to test change
for( size_t n = 0; n < new_xs.size(); n++ )
{
old_xs[n] = new_xs[n];
new_xs[n] *= 1.2;
}
// FAST ACCESS
trx->fast_set_xs( nind_list, transition_ids, tinds_2d, new_xs );
// check that the change took place
double target2 = 1.2;
for( size_t k = 0; k < nind_list.size(); k++ )
{
double myold = old_xs[k];
if( myold > 0 )
{
double mynew = trx->get_xs( nind_list[k], transition_ids[k] );
EXPECT_NEAR( 0.0, std::abs( mynew / myold - target2 ), 1e-5 );
}
}
tinds_2d.deepDelete();
}
}
