tstMaterial.cpp
./Core/dc/tests/tstMaterial.cpp
#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
#include "Nemesis/gtest/nemesis_gtest.hh"
#include "Nemesis/harness/DBC.hh"
#include "Origen/Core/dc/FakeFactory.h"
#include "Origen/Core/dc/Material.h"
#include "Origen/Core/dc/MaterialMap.h"
#include "Origen/Core/xf/Solver_Fake.h"
#include "ScaleSTL/Functions.h"
#include "ScaleUtils/IO/Utils.h"
using std::string;
namespace Origen
{
namespace test
{
TEST( Material, Setup )
{
ASSERT_TRUE( created );
// initialize material
double volume = 0.55 * 0.55 * 3.14;
string name = "FUEL.1";
int id = 1001;
string library_type = "random1";
// check initial conditions
ASSERT_EQ( 1, mat->ntimes() );
ASSERT_EQ( 0, mat->nsteps() );
ASSERT_EQ( 0, mat->point_bos() );
ASSERT_EQ( -1, mat->point_eos() );
ASSERT_EQ( 0.0, mat->time_bos() );
ASSERT_EQ( -1.0, mat->time_eos() );
ASSERT_EQ( 0.0, mat->time_at( 0 ) );
// ASSERT_TRUE(mat->time_at(1)==0.0); //should throw
ASSERT_TRUE( mat->amount_bos() != nullptr );
ASSERT_TRUE( mat->amount_eos() == nullptr );
ASSERT_TRUE( mat->transition_matrix() == nullptr );
// check variables
ASSERT_EQ( library_type, mat->library_type() );
ASSERT_EQ( library, mat->library() );
ASSERT_EQ( name, mat->name() );
ASSERT_EQ( id, mat->id() );
ASSERT_EQ( volume, mat->cold_volume() );
// check number of nuclides
size_t total_nuclides = library->definition().total_nuclides();
EXPECT_EQ( total_nuclides, mat->amount_bos()->size() );
// check equal to zero
for( size_t i = 0; i < total_nuclides; ++i )
{
EXPECT_EQ( 0.0, mat->amount_bos()->at( i ) );
}
std::vector<double> numden( mat->total_nuclides() );
mat->get_numden_bos( &numden );
ASSERT_EQ( total_nuclides, numden.size() );
for( size_t i = 0; i < total_nuclides; ++i )
{
EXPECT_EQ( 0.0, numden[i] );
}
}
TEST( Material, Accessors )
{
// reference bos data
Vec_Int ref_ids;
Vec_Dbl ref_vals;
FakeFactory::vera_uox_e360( ref_ids, ref_vals );
// create library
string library_type = "pwr.*.orglib";
ASSERT_TRUE( created );
// initialize material
double volume = 0.55 * 0.55 * 3.14;
string name = "FUEL.1";
int id = 1001;
mat->set_numden_bos( ref_vals, ref_ids );
// check masses
EXPECT_FLOAT_EQ( 8.58616, mat->initial_hm_mass() );
EXPECT_FLOAT_EQ( 9.74114, mat->initial_mass() );
// get number densities
{
Vec_Dbl vals( ref_ids.size() );
mat->get_numden_bos( &vals, ref_ids );
for( size_t i = 0; i < ref_ids.size(); ++i )
{
EXPECT_FLOAT_EQ( ref_vals[i], vals[i] );
}
}
// get a full transition matrix
SP_TransitionMatrixP trx_m1 =
mat->library()->newsp_transition_matrix_at( 0 );
double flux_m1 = 1.2;
// add a step
mat->allocate_step();
mat->set_transition_matrix( trx_m1 );
mat->set_flux( flux_m1 );
mat->set_dt( 65.0 );
// check ability to change time step length
ASSERT_TRUE( mat->time_eos() == 65.0 );
mat->set_dt( 123.456 );
// check conditions after a step is added
ASSERT_EQ( 2, mat->ntimes() );
ASSERT_EQ( 1, mat->nsteps() );
ASSERT_EQ( 0, mat->point_bos() );
ASSERT_EQ( 1, mat->point_eos() );
ASSERT_EQ( 0.0, mat->time_bos() );
ASSERT_EQ( 123.456, mat->time_eos() );
ASSERT_EQ( 0.0, mat->time_at( 0 ) );
ASSERT_EQ( 123.456, mat->time_at( 1 ) );
// check variables
ASSERT_EQ( trx_m1, mat->transition_matrix() );
ASSERT_TRUE( mat->amount_bos() != nullptr );
ASSERT_TRUE( mat->amount_eos() != nullptr );
// set eos amount
{
Vec_Int ids1( 10 );
Vec_Dbl vals1( 10 );
// 3 - fission product
// 2 - actinide
// 1 - activation product
// SIZZZAAA //final SIZZZAAA
ids1[0] = 92235;
vals1[0] = 0.90e-3; // 20092235
ids1[1] = 92238;
vals1[1] = 0.99e-2; // 20092238
ids1[2] = 64156;
vals1[2] = 0.10e-2; // 30064156
ids1[3] = 48114;
vals1[3] = 0.20e-2; // 30048114
ids1[4] = 1050121;
vals1[4] = 0.30e-2; // 31050121
ids1[5] = 20082210;
vals1[5] = 0.40e-2; // 20082210
ids1[6] = 20082208;
vals1[6] = 0.50e-2; // 20082208
ids1[7] = 30046105;
vals1[7] = 0.60e-2; // 30046105
ids1[8] = 10048114;
vals1[8] = 0.70e-2; // 10048114
ids1[9] = 30064158;
vals1[9] = 0.80e-2; // 30064158
mat->set_numden_eos( vals1, ids1 );
// get the numden using the full SIZZZAAA
Vec_Int ids2( 10 );
Vec_Dbl vals2( 10 );
ids2[0] = 20092235;
ids2[1] = 20092238;
ids2[2] = 30064156;
ids2[3] = 30048114;
ids2[4] = 31050121;
ids2[5] = 20082210;
ids2[6] = 20082208;
ids2[7] = 30046105;
ids2[8] = 10048114;
ids2[9] = 30064158;
mat->get_numden_eos( &vals2, ids2 );
for( size_t i = 0; i < ids2.size(); ++i )
{
EXPECT_FLOAT_EQ( vals1[i], vals2[i] );
}
// test being able to extract two sublibs worth
Vec_Int ids3( 1 );
Vec_Dbl vals3( 1 );
ids3[0] = 48114; // will get both 10048114+30048114=0.90e-2
mat->get_numden_eos( &vals3, ids3 );
EXPECT_FLOAT_EQ( 0.90e-2, vals3[0] );
}
}
{
// reference bos data
Vec_Int ref_ids;
Vec_Dbl ref_vals;
FakeFactory::vera_uox_e360( ref_ids, ref_vals );
// create library
string library_type = "pwr.*.orglib";
bool created = FakeFactory::Library_scale_pwr( *library );
ASSERT_TRUE( created );
MaterialMap matmap;
// initialize materials
double total_volume = 0.0;
for( size_t j = 1; j < 10; ++j )
{
double volume = j * 0.55 * 0.55 * 3.14;
total_volume += volume;
string name = ScaleUtils::IO::nprintf( "FUEL.%d", j );
int id = 1000 + j;
SP_Material mat(
new Material( library_type, library, name, id, volume ) );
mat->set_numden_bos( ref_vals, ref_ids );
matmap.add( mat );
}
EXPECT_FLOAT_EQ( total_volume, matmap.cold_volume() );
EXPECT_FLOAT_EQ( 10.255454 * total_volume, matmap.initial_mass() );
EXPECT_FLOAT_EQ( 10.255454 * total_volume * 0.88143277,
matmap.initial_hm_mass() );
}
{
SP_Material mat = std::make_shared<Material>(
Vec_Int ref_ids;
Vec_Dbl ref_vals;
FakeFactory::vera_uox_e360( ref_ids, ref_vals );
mat->set_numden_bos( ref_vals, ref_ids );
mat->set_solver( std::make_shared<Origen::Solver_Fake>() );
return mat;
}
TEST( Material, Flush )
{
auto mat = newspMaterial_local();
Vec_Dbl dtrel( 4, 0.25 ); // relative distribution of time over the step
Vec_Dbl flux, power;
size_t ref_nsteps = 10;
size_t ref_npoints = ref_nsteps + 1;
Vec_Dbl ref_fluence( 1, 0.0 );
Vec_Dbl ref_flux;
Vec_Dbl ref_dt;
for( size_t i = 0; i < ref_nsteps; ++i )
{
double dt = 86400. * 30 * ( 1 + i / 2 );
ref_dt.push_back( dt );
ref_flux.push_back( f );
ref_fluence.push_back( ref_fluence.back() + dt * f );
mat->add_step( dt );
mat->set_flux( f );
mat->set_transition_matrix(
mat->library()->newsp_transition_matrix_at( 0 ) );
mat->solve( dtrel, &flux, &power );
EXPECT_FLOAT_EQ( dt, mat->dt() );
EXPECT_FLOAT_EQ( f, mat->flux() );
EXPECT_FLOAT_EQ( ref_fluence[i], mat->fluence_bos() );
EXPECT_FLOAT_EQ( ref_fluence[i + 1], mat->fluence_eos() );
}
// check stuff before flush
EXPECT_EQ( ref_nsteps, mat->nsteps() );
for( size_t i = 0; i < ref_nsteps; ++i )
{
EXPECT_TRUE( mat->transition_matrix_over( i ) != nullptr );
EXPECT_FLOAT_EQ( ref_dt[i], mat->dt_over( i ) );
EXPECT_FLOAT_EQ( ref_flux[i], mat->flux_over( i ) );
EXPECT_FLOAT_EQ( ref_fluence[i], mat->fluence_at( i ) );
EXPECT_FLOAT_EQ( ref_fluence[i + 1], mat->fluence_at( i + 1 ) );
}
// check amounts
EXPECT_EQ( ref_npoints, mat->ntimes() );
for( size_t p = 0; p < ref_npoints; ++p )
{
EXPECT_TRUE( mat->amount_at( p ) != nullptr );
}
//
// FLUSHING
//
// now flush and check that big data is cleared but scalar data remains
mat->flush_oldest();
// check point data (bos and eos for last will be non-null)
EXPECT_TRUE( mat->amount_bos() != nullptr );
EXPECT_TRUE( mat->amount_eos() != nullptr );
for( size_t p = 0; p < ref_npoints; ++p )
{
if( p >= ref_npoints - 2 )
{
EXPECT_TRUE( mat->amount_at( p ) != nullptr );
}
else
{
EXPECT_TRUE( mat->amount_at( p ) == nullptr );
}
}
// check interval data
for( size_t i = 0; i < ref_nsteps; ++i )
{
if( i >= ref_nsteps - 1 )
{
EXPECT_TRUE( mat->transition_matrix_over( i ) != nullptr );
}
else
{
EXPECT_TRUE( mat->transition_matrix_over( i ) == nullptr );
}
// flux/fluence should be same as above check
EXPECT_FLOAT_EQ( ref_dt[i], mat->dt_over( i ) );
EXPECT_FLOAT_EQ( ref_flux[i], mat->flux_over( i ) );
EXPECT_FLOAT_EQ( ref_fluence[i], mat->fluence_at( i ) );
EXPECT_FLOAT_EQ( ref_fluence[i + 1], mat->fluence_at( i + 1 ) );
}
}
TEST( Material, PowerBurnup )
{
auto mat = newspMaterial_local();
Vec_Dbl dtrel( 2, 0.50 ); // relative distribution of time over the step
Vec_Dbl flux, power;
double initial_hm = mat->initial_hm_mass();
double hm0 = 45.197456;
EXPECT_FLOAT_EQ( hm0, initial_hm ); // grams
mat->add_step( 86400. ); // 1 day
mat->set_power( 300.0 ); // 300 Watts
mat->set_transition_matrix(
mat->library()->newsp_transition_matrix_at( 0 ) );
mat->solve( dtrel, &flux, &power );
double bu1 = 300.0 / hm0;
EXPECT_FLOAT_EQ( 300.0, mat->power() );
EXPECT_FLOAT_EQ( 0.0, mat->burnup_bos() );
EXPECT_FLOAT_EQ( bu1, mat->burnup_eos() );
mat->add_step( 86400. * 10000 ); // 10000 days
mat->set_power( 200.0 ); // 200 Watts
mat->set_transition_matrix(
mat->library()->newsp_transition_matrix_at( 0 ) );
mat->solve( dtrel, &flux, &power );
double bu2 = bu1 + 200.0 * 10000 / hm0;
EXPECT_FLOAT_EQ( 200.0, mat->power() );
EXPECT_FLOAT_EQ( bu1, mat->burnup_bos() );
EXPECT_FLOAT_EQ( bu2, mat->burnup_eos() );
// verify we're storing these
EXPECT_FLOAT_EQ( 300.0, mat->power_over( 0 ) );
EXPECT_FLOAT_EQ( bu1, mat->burnup_at( 1 ) );
}
TEST( Material, NuclideRemoval )
{
// mat --> removal
// necessary output
Vec_Dbl flux( 4, 0.0 );
Vec_Dbl power( 4, 0.0 );
// input
Vec_Dbl dtrel( 4, 0.25 );
// mix1
auto mat = newspMaterial_local();
// do first step
mat->add_step( 10 );
mat->set_power( 300 ); // 300 W
mat->set_transition_matrix(
mat->library()->newsp_transition_matrix_at( 0 ) );
mat->solve( dtrel, &flux, &power );
// do 2nd step (decay) and check concs
mat->add_step( 0.01 );
mat->set_power( 0 ); // decay
mat->set_transition_matrix(
mat->library()->newsp_transition_matrix_at( 0 ) );
Vec_Int removal_ids{54135}; // why doesn't this work with 54000 only 54135?
Vec_Dbl removal_vals{10.};
mat->set_removal( r );
mat->solve( dtrel, &flux, &power );
// concentrations we're interested in
Concentrations concs_w( NuclideSet( Vec_Int{54135} ) );
mat->get_concs_at( &concs_w, 2 );
// re-do 2nd step (decay) without removal
mat->flush_newest();
mat->add_step( 0.01 );
mat->set_power( 0 ); // decay
mat->set_transition_matrix(
mat->library()->newsp_transition_matrix_at( 0 ) );
mat->set_removal( nullptr );
mat->solve( dtrel, &flux, &power );
// concentrations we're interested in
Concentrations concs_wo( NuclideSet( Vec_Int{54135} ) );
mat->get_concs_at( &concs_wo, 2 );
EXPECT_LT( concs_w.vals_at( 0 ), concs_wo.vals_at( 0 ) );
}
#if 0
TEST( Material, RemovalToFeed )
{
// mat1 --> removal --> mat2
//mix1
auto mat1= newspMaterial_local();
mat1->add_step( 86400. ); // 1 day
mat1->set_power( 300.0 ); // 300 Watts
mat1->set_transition_matrix(
mat1->library()->newsp_transition_matrix_at( 0 ) );
Vec_Int removal_ids{54000};
Vec_Dbl removal_vals{1e-4};
mat1->set_removal_over( r, 0 );
//mix2
SP_Material mat2 = std::make_shared<Material>("scale_pwr",FakeFactory::newspLibrary_scale_pwr(),"mix2",2,5.0);
mat2->add_step( 86400. ); // 1 day
mat2->set_power( 0.0 );
mat2->set_transition_matrix(
mat2->library()->newsp_transition_matrix_at( 0 ) );
//necessary output
Vec_Dbl flux, power;
//test loop with varying substeps
for( size_t i=0; i<4; ++i)
{
//test power of 2 # of substeps
size_t n=2**i;
Vec_Dbl dtrel( n, 1.0/n );
mat1->solve( dtrel, &flux, &power );
//transform removal over step to feed rate
auto f2 = mat1->average_removal_rate_over( 0 );
mat2->set_feed_rate_over( f2, 0 );
mat2->solve( dtrel, &flux, &power );
}
}
#endif
} // namespace test
} // namespace Origen
