tstSolver_cram.cpp
./Solver/cram/tests/tstSolver_cram.cpp
/**************************************************************************/ /**
* \file tstSolver_cram.cpp
* \brief Unit tests for the class Solver_cram
* \author Aarno Isotalo & William Wieselquist
*
*****************************************************************************/
#include "Nemesis/gtest/nemesis_gtest.hh"
#include "Origen/Core/config.h"
#include "Origen/Core/dc/TransitionMatrixP.h"
#include "Origen/Core/io/LibraryIO.h"
#include "Origen/Solver/cram/Solver_cram.h"
#include "ScaleUtils/IO/DB.h"
// For Will's MultiZoneDepleter test
#include "Origen/Core/fn/print.h"
#include "Origen/Core/xf/MultiZoneDepleter.h"
#include "Origen/Core/dc/FakeFactory.h"
#include "Standard/Interface/Communicator.h"
/**************************************************************************/ /**
* \brief Test the check_opts() function.
*****************************************************************************/
TEST( Solver_cram, check_opts )
{
Origen::Solver_cram solver;
ScaleUtils::IO::DB opts;
// Get the default options
solver.get_defaults( opts );
EXPECT_TRUE( opts.has<std::string>( "solver" ) );
EXPECT_TRUE( opts.has<int>( "cramOrder" ) );
EXPECT_TRUE( opts.has<int>( "cramInternalSubsteps" ) );
EXPECT_TRUE( opts.has<bool>( "isAdjoint" ) );
EXPECT_TRUE( opts.has<int>( "cramRemoveNegatives" ) );
EXPECT_TRUE( opts.has<double>( "cramResultsCutoff" ) );
EXPECT_TRUE( opts.has<bool>( "tallyEnergyReleased" ) );
EXPECT_TRUE( opts.has<bool>( "tallyNumFissions" ) );
EXPECT_TRUE( opts.has<bool>( "tallyNumCaptures" ) );
EXPECT_TRUE( opts.has<bool>( "tallyAveConcentrations" ) );
// Check the default options
// *** Test allowing missing arguments
solver.get_defaults( opts );
opts.remove<int>( "cramOrder" );
EXPECT_FALSE( solver.check_opts( opts, false ) );
opts.clearErrorStack();
EXPECT_TRUE( solver.check_opts( opts, true ) );
solver.get_defaults( opts );
opts.remove<int>( "cramInternalSubsteps" );
EXPECT_FALSE( solver.check_opts( opts, false ) );
opts.clearErrorStack();
EXPECT_TRUE( solver.check_opts( opts, true ) );
solver.get_defaults( opts );
opts.remove<bool>( "isAdjoint" );
EXPECT_FALSE( solver.check_opts( opts, false ) );
opts.clearErrorStack();
EXPECT_TRUE( solver.check_opts( opts, true ) );
solver.get_defaults( opts );
opts.remove<int>( "cramRemoveNegatives" );
EXPECT_FALSE( solver.check_opts( opts, false ) );
opts.clearErrorStack();
EXPECT_TRUE( solver.check_opts( opts, true ) );
solver.get_defaults( opts );
opts.remove<double>( "cramResultsCutoff" );
EXPECT_FALSE( solver.check_opts( opts, false ) );
opts.clearErrorStack();
EXPECT_TRUE( solver.check_opts( opts, true ) );
solver.get_defaults( opts );
opts.remove<bool>( "tallyEnergyReleased" );
EXPECT_FALSE( solver.check_opts( opts, false ) );
opts.clearErrorStack();
EXPECT_TRUE( solver.check_opts( opts, true ) );
solver.get_defaults( opts );
opts.remove<bool>( "tallyNumFissions" );
EXPECT_FALSE( solver.check_opts( opts, false ) );
opts.clearErrorStack();
EXPECT_TRUE( solver.check_opts( opts, true ) );
solver.get_defaults( opts );
opts.remove<bool>( "tallyNumCaptures" );
EXPECT_FALSE( solver.check_opts( opts, false ) );
opts.clearErrorStack();
EXPECT_TRUE( solver.check_opts( opts, true ) );
solver.get_defaults( opts );
opts.remove<bool>( "tallyAveConcentrations" );
EXPECT_FALSE( solver.check_opts( opts, false ) );
opts.clearErrorStack();
EXPECT_TRUE( solver.check_opts( opts, true ) );
solver.get_defaults( opts );
opts.remove<std::string>( "solver" );
EXPECT_FALSE( solver.check_opts( opts, false ) );
opts.clearErrorStack();
EXPECT_TRUE( solver.check_opts( opts, true ) );
// *** Test incorrect input
solver.get_defaults( opts, true );
opts.set<int>( "cramInternalSubsteps", -1 );
EXPECT_FALSE( solver.check_opts( opts, true ) );
opts.clearErrorStack();
solver.get_defaults( opts, true );
opts.set<int>( "cramOrder", 13 );
EXPECT_FALSE( solver.check_opts( opts, true ) );
opts.clearErrorStack();
solver.get_defaults( opts, true );
opts.set<std::string>( "solver", "matrex" );
EXPECT_FALSE( solver.check_opts( opts, true ) );
opts.clearErrorStack();
// Print all errors for manual checking of their formatting (kills the test)
if( false )
{
ScaleUtils::IO::DB dummy;
solver.check_opts( dummy, false );
dummy.set<int>( "cramOrder", 13 );
dummy.set<int>( "cramInternalSubsteps", -1 );
dummy.set<std::string>( "solver", "matrex" );
dummy.set<bool>( "isAdjoint", false );
dummy.set<int>( "cramRemoveNegatives", 5 );
dummy.set<double>( "cramResultsCutoff", -1.0 );
solver.check_opts( dummy, false );
}
}
/**************************************************************************/ /**
* \brief Test the set_opts() function.
*****************************************************************************/
TEST( Solver_cram, set_opts )
{
Origen::Solver_cram solver;
ScaleUtils::IO::DB opts;
// Set some options.
opts.set<int>( "cramOrder", 14 );
opts.set<int>( "cramInternalSubsteps", 5 );
opts.set<bool>( "isAdjoint", true );
opts.set<int>( "cramRemoveNegatives", 0 );
opts.set<double>( "cramResultsCutoff", 1e-50 );
opts.set<bool>( "tallyEnergyReleased", true );
opts.set<bool>( "tallyNumFissions", true );
opts.set<bool>( "tallyNumCaptures", true );
opts.set<bool>( "tallyAveConcentrations", true );
solver.set_opts( opts );
// Check that the options were correctly set.
ScaleUtils::IO::DB status = solver.status();
EXPECT_EQ( opts.get<int>( "cramOrder" ), status.get<int>( "cramOrder" ) );
EXPECT_EQ( opts.get<int>( "cramInternalSubsteps" ),
status.get<int>( "cramInternalSubsteps" ) );
EXPECT_EQ( opts.get<bool>( "isAdjoint" ), status.get<bool>( "isAdjoint" ) );
EXPECT_EQ( opts.get<int>( "cramRemoveNegatives" ),
status.get<int>( "cramRemoveNegatives" ) );
EXPECT_EQ( opts.get<double>( "cramResultsCutoff" ),
status.get<double>( "cramResultsCutoff" ) );
EXPECT_EQ( opts.get<bool>( "tallyEnergyReleased" ),
status.get<bool>( "tallyEnergyReleased" ) );
EXPECT_EQ( opts.get<bool>( "tallyNumFissions" ),
status.get<bool>( "tallyNumFissions" ) );
EXPECT_EQ( opts.get<bool>( "tallyNumCaptures" ),
status.get<bool>( "tallyNumCaptures" ) );
EXPECT_EQ( opts.get<bool>( "tallyAveConcentrations" ),
status.get<bool>( "tallyAveConcentrations" ) );
EXPECT_EQ( "cram", status.get<std::string>( "solver" ) );
}
/**************************************************************************/ /**
* \brief Test solving.
*****************************************************************************/
TEST( Solver_cram, solve )
{
Origen::Solver_cram solver;
// Load a library.
Origen::Library lib;
ASSERT_TRUE( Origen::FakeFactory::Library_scale_pwr( lib ) );
// Number of nuclides.
int itot = trx->get_itot();
// Step length.
double dt = 100 * 60 * 60 * 24;
// Set the transition matrix.
solver.set_transition_matrix( &*trx );
// Set source term.
Origen::Vec_Dbl src( itot, 0.0 );
src[trx->find_nuclide_guess( 922350 )] = 0.05e-7;
src[trx->find_nuclide_guess( 922380 )] = 0.20e-7;
src[trx->find_nuclide_guess( 932390 )] = 0.05e-7;
src[trx->find_nuclide( 1, 350790 )] = 1e-7;
src[trx->find_nuclide( 1, 350810 )] = 1e-20;
src[trx->find_nuclide( 1, 350780 )] = 1e-7;
src[trx->find_nuclide( 1, 350800 )] = 1e-20;
solver.set_source( &src );
// Initial composition.
Origen::Vec_Dbl n0( itot, 0.0 );
n0[trx->find_nuclide_guess( 922350 )] = 0.05;
n0[trx->find_nuclide_guess( 922380 )] = 0.95;
// Final composition.
Origen::Vec_Dbl n( itot, 0.0 );
// Solve.
EXPECT_TRUE( solver.solve( n0, 3.5e14, dt, &n ) );
// Check the results. Reference values were calculated with CRAM in
// tstCRAM.FreshFuelDepletionWithSource
// (Origen/Manager/Wrapper/tests/tstCRAM.cpp)
// and extracted by adding the below lines. This only needs to test that
// CRAM was called properly.
/*
printf("%.14e\n",n[trx->find_nuclide_guess(922350)]);
printf("%.14e\n",n[trx->find_nuclide_guess(942390)]);
printf("%.14e\n",n[trx->find_nuclide(3,541350)]);
printf("%.14e\n",n[trx->find_nuclide(3,430990)]);
*/
EXPECT_FLOAT_EQ( n[trx->find_nuclide_guess( 922350 )], 0.085176535 );
EXPECT_FLOAT_EQ( n[trx->find_nuclide_guess( 942390 )], 0.037553333 );
EXPECT_FLOAT_EQ( n[trx->find_nuclide( 3, 541350 )], 2.0385921e-06 );
EXPECT_FLOAT_EQ( n[trx->find_nuclide( 3, 430990 )], 0.00062933267 );
// Check that removal gets set.
double u236 = n[trx->find_nuclide_guess( 922360 )];
// Set removal.
Origen::Vec_Dbl lambda( itot, 0.0 );
lambda[trx->find_nuclide_guess( 922360 )] = 1; // U-236 disappears fast
solver.set_continuous_loss( &lambda );
// Solve.
EXPECT_TRUE( solver.solve( n0, 3.5e14, dt, &n ) );
// Check that U236 amount decreased.
EXPECT_TRUE( n[trx->find_nuclide_guess( 922360 )] < 1e-5 * u236 );
}
/**************************************************************************/ /**
* \brief Test the depletion tallies (mostly energy released)
*****************************************************************************/
TEST( Solver_cram, tallies )
{
Origen::Solver_cram solver;
// Load a library.
Origen::Library lib;
ASSERT_TRUE( Origen::FakeFactory::Library_scale_pwr( lib ) );
Origen::TransitionMatrixP::SP trx = lib.newsp_transition_matrix_at( 0 );
// Set the transition matrix.
solver.set_transition_matrix( &*trx );
// Turn the tallies on.
ScaleUtils::IO::DB opts;
opts.set<bool>( "tallyEnergyReleased", true );
opts.set<bool>( "tallyNumFissions", true );
opts.set<bool>( "tallyNumCaptures", true );
opts.set<bool>( "tallyAveConcentrations", true );
solver.set_opts( opts );
// Step length & neutron flux.
double dt = 1 * 60 * 60 * 24;
double flux = 3.5e14;
// Initial composition
Origen::Vec_Dbl n0( trx->get_itot(), 0.0 );
n0[trx->find_nuclide_guess( 922350 )] = 0.05;
n0[trx->find_nuclide_guess( 922380 )] = 0.95;
// Solve
Origen::Vec_Dbl n( trx->get_itot(), 0.0 );
EXPECT_TRUE( solver.solve( n0, 3.5e14, dt, &n ) );
// Get the tally results
const ScaleUtils::IO::DB& status = solver.status();
double energy_tally = status.get<double>( "energyReleased" );
double fission_tally = status.get<double>( "numFissions" );
double capture_tally = status.get<double>( "numCaptures" );
// Get the total power at BOS and EOS
double BOS_power = trx->get_reaction_power( n0, 1, flux ) / 1e24;
double EOS_power = trx->get_reaction_power( n, 1, flux ) / 1e24;
// Check that the energy tally is close. It should for such short step
double energy_ref = dt * ( BOS_power + EOS_power ) / 2;
EXPECT_NEAR( energy_ref, energy_tally, energy_ref * 1e-3 );
// Check that the other tallies returned something
EXPECT_GT( fission_tally, 0.0 );
EXPECT_LT( fission_tally, 1.0 );
EXPECT_GT( capture_tally, 0.0 );
EXPECT_LT( capture_tally, 1.0 );
// Check that the U238 Average concentration is reasnable
// This is a one day depletion
double u238 = status.get<double>( "ave_922380" );
EXPECT_LT( u238, 0.94999 );
EXPECT_GT( u238, 0.949 );
}
/******************************************************************************
* The below part actually tests MultiZoneDepleter as much as Solver_matrex
*
*****************************************************************************/
using ScaleUtils::IO::nprintf;
// GLOBAL COMMUNICATOR
Standard::Communicator world;
{
size_t ntimes = mat->ntimes();
std::vector<int> ids, user_to_master, master_to_user;
unsigned int w = 13;
unsigned int pr = 3;
mat->mapping_sets( ids, &user_to_master, &master_to_user );
std::ofstream summary( file + "_summary.txt" );
summary << std::setw( 5 ) << "index" << std::setw( 11 ) << "time(d)"
<< std::setw( w ) << "power(W)" << std::setw( w ) << "mass(g)"
<< std::setw( w ) << "act(Bq)" << std::setw( w ) << "xs(1/cm)"
<< std::endl;
for( size_t p = 0; p < ntimes; ++p )
{
std::vector<double> numden( ids.size() );
Origen::SP_DoubleList amount = mat->amount_at( p );
double volume = mat->cold_volume();
mat->get_numden_at( &numden, p, ids );
summary << std::fixed << std::setw( 5 ) << ( p + 1 )
<< std::setprecision( 3 ) << std::scientific << std::setw( 11 )
<< mat->time_at( p ) / 86400. << std::setw( w )
<< ( p == 0 ? 0.0 : mat->power_over( p - 1 ) );
// output number densities
{
std::stringstream ss;
ss << file << "_nd_" << std::setw( 2 ) << std::setfill( '0' ) << p
<< ".txt";
std::ofstream ofs( ss.str() );
for( size_t i = 0; i < ids.size(); ++i )
{
ofs << std::fixed << std::setw( 3 ) << ( i + 1 )
<< std::setw( 9 ) << ids[i] << std::scientific
<< std::setprecision( 5 ) << std::setw( 15 ) << numden[i]
<< std::endl;
}
ofs.close();
}
// output number densities
{
std::stringstream ss;
std::vector<double> numden_all( mat->total_nuclides() );
mat->get_numden_at( &numden_all, p );
ss << file << "_nda_" << std::setw( 2 ) << std::setfill( '0' ) << p
<< ".txt";
std::ofstream ofs( ss.str() );
for( size_t i = 0; i < numden_all.size(); ++i )
{
ofs << std::scientific << std::setprecision( 5 )
<< std::setw( 15 ) << numden_all[i] << std::endl;
}
ofs.close();
}
// output masses in grams
{
std::stringstream ss;
ss << file << "_m_" << std::setw( 2 ) << std::setfill( '0' ) << p
<< ".txt";
std::ofstream ofs( ss.str() );
double mass_sum = 0.0;
for( size_t i = 0; i < ids.size(); ++i )
{
int m = user_to_master[i];
double mass =
m < 0 ? 0.0 : mat->library()->decay_data().masses[m];
double val = numden[i] * volume * mass / 0.60221413;
ofs << std::fixed << std::setw( 3 ) << ( i + 1 )
<< std::setw( 9 ) << ids[i] << std::scientific
<< std::setprecision( 5 ) << std::setw( 15 ) << val
<< std::endl;
mass_sum += val;
}
summary << std::setprecision( pr ) << std::scientific
<< std::setw( w ) << mass_sum;
ofs.close();
}
// output activities
{
std::stringstream ss;
ss << file << "_bq_" << std::setw( 2 ) << std::setfill( '0' ) << p
<< ".txt";
std::ofstream ofs( ss.str() );
double bq_sum = 0.0;
for( size_t i = 0; i < ids.size(); ++i )
{
int m = user_to_master[i];
double lambda =
m < 0 ? 0.0
: mat->library()->decay_data().decay_constants[m];
double val = 0.60221413e+23 * numden[i] * volume * lambda;
ofs << std::fixed << std::setw( 3 ) << ( i + 1 )
<< std::setw( 9 ) << ids[i] << std::scientific
<< std::setprecision( 5 ) << std::setw( 15 ) << val
<< std::endl;
bq_sum += val;
}
summary << std::setprecision( pr ) << std::scientific
<< std::setw( w ) << bq_sum;
ofs.close();
}
// output macros
{
std::stringstream ss;
ss << file << "_xs_" << std::setw( 2 ) << std::setfill( '0' ) << p
<< ".txt";
std::ofstream ofs( ss.str() );
for( size_t i = 0; i < ids.size(); ++i )
{
int m = user_to_master[i];
double sigma =
m < 0
? 0.0
: mat->library()->transition_coeff_at( 0 ).loss_xs().at(
m );
double val = numden[i] * sigma;
ofs << std::fixed << std::setw( 3 ) << ( i + 1 )
<< std::setw( 9 ) << ids[i] << std::scientific
<< std::setprecision( 5 ) << std::setw( 15 ) << val
<< std::endl;
xs_sum += val;
}
summary << std::setprecision( pr ) << std::scientific
<< std::setw( w ) << xs_sum;
ofs.close();
}
summary << std::endl;
}
}
/**************************************************************************/ /**
* \brief Test through the MultiZoneDepleter class.
*****************************************************************************/
TEST( CycleDepletion, SingleMaterial )
{
// Get a transition matrix and library.
Origen::SP_Library lib = std::make_shared<Origen::Library>();
{
ASSERT_TRUE( Origen::FakeFactory::Library_scale_pwr( *lib ) );
trx = lib->newsp_transition_matrix_at( 0 );
}
jDebugTaggedLine( "finished SetUp" );
// Setup materials.
int id1 = 1;
// double assembly_mass = 0.4; //MTU
double total_mass = 0.4 / 0.8815; // MT
double volume = total_mass * 1e6 / 10.2; // 10.2 g/cc
// Create first begining-of-step number density.
Origen::SP_Material mat1 = std::make_shared<Origen::Material>(
"GENERAL", lib, "mat1", id1, volume );
std::vector<int> ids;
std::vector<double> numden;
Origen::FakeFactory::vera_uox_e360( ids, numden );
mat1->set_numden_bos( numden, ids );
// Associate materials with map.
Origen::SP_MaterialMap matmap = std::make_shared<Origen::MaterialMap>();
matmap->add( mat1 );
// Create depletion times/fluxes.
// cycle1--------------> cycle2--------------->
// deplete------>decay-> deplete----->decay--->
std::vector<double> time{
0, 3, 250, 500, 515, 530, 533, 780, 1030, 1045, 1060};
std::vector<double> flux{4, 4, 4, 0, 0, 3, 3, 3, 0, 0};
for( size_t j = 0; j < time.size(); ++j ) time[j] *= 86400.0;
for( size_t j = 0; j < flux.size(); ++j ) flux[j] *= 1.0e14;
size_t nsteps = flux.size();
// Create step depleter.
DataContainer dc;
// set solver
mzd.set_solver( slv );
// set materials
mzd.set_materials( matmap );
// Enter step loop.
for( size_t j = 0; j < nsteps; ++j )
{
std::stringstream ss;
ss << "loop(" << std::setw( 3 ) << std::setfill( '0' ) << ( j + 1 )
<< ")";
SCOPED_TIMER( ss.str() );
// Add new step data.
double stepTime = time[j + 1] - time[j];
mzd.add_step( stepTime );
mat1->set_flux( flux[j] );
mat1->set_transition_matrix( trx );
// Update data container.
dc.addInt( "numSubsteps", 1 );
dc.addDouble( "stepTime", stepTime );
dc.addBool( "depleteByPower", false );
mzd.update( dc );
// Step in time.
mzd.execute();
EXPECT_FLOAT_EQ( stepTime, mat1->dt() );
EXPECT_FLOAT_EQ( flux[j], mat1->flux() );
// Step output.
/*
std::cout << std::fixed <<
" step="<< std::setw(3) << (j+1)<<
std::scientific << std::setprecision(2) <<
" dt="<<mat1->dt()<<
" flux="<< mat1->flux()<<
" power="<< mat1->power()<<
std::endl;
*/
}
// Print nuclides of interest.
print_summary_1( mat1, "base" );
std::ofstream ofs( "matmap.txt" );
ofs << matmap->to_string();
ofs.close();
// timing table
Origen::printTimingReport( std::cout );
// wrapup
dc.deepDelete();
}
