exNEAMS.cpp
./Manager/mzd/tests/exNEAMS.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/TransitionMatrixP.h"
#include "Origen/Core/io/LibraryIO.h"
#include "Origen/Core/xf/MultiZoneDepleter.h"
#include "Origen/Core/xf/Solver_Fake.h"
#include "Origen/Solver/matrex/Solver_matrex.h"
#include "ScaleUtils/IO/DB.h"
#include "Standard/Interface/AbstractWriter.h"
#include "Standard/Interface/BasicIOWriter.h"
#include "Standard/Interface/Communicator.h"
#include "Standard/Interface/jdebug.h"
using ScaleUtils::IO::nprintf;
using namespace Origen;
// GLOBAL COMMUNICATOR
Standard::Communicator world;
TEST( Example, 1CreationOfAMaterial )
{
// Get a general 2237-nuclide ORIGEN library.
FakeFactory::Library_scale_pwr( *lib );
// Extract a transition matrix from it. We will use this PWR transition
// matrix for all materials in this test.
SP_TransitionMatrix trx( lib->newsp_transition_matrix_at( 0 ) );
// Create a single 3.6% enriched UO2 material.
SP_Material mat;
{
std::vector<int> ids; // nuclide ids in IZZZAAA format
std::vector<double> numden; // atoms/barn-cm
double volume = 5.6; // cm^3
int id = 1234; // material id
std::string name = "material_1234"; // material name
// Use the FakeFactory to get a reasonable initial composition.
FakeFactory::vera_uox_e360( ids, numden );
// Initialize the material.
// Set the beginning of step number densities.
mat->set_numden_bos( numden, ids );
}
}
TEST( Example, 2SimpleDepletionByFlux )
{
// Get a general 2237-nuclide ORIGEN library.
FakeFactory::Library_scale_pwr( *lib );
// Extract a transition matrix from it. We will use this PWR transition
// matrix for all materials in this test.
SP_TransitionMatrix trx( lib->newsp_transition_matrix_at( 0 ) );
// Create a single 3.6% enriched UO2 material.
SP_Material mat;
{
std::vector<int> ids; // nuclide ids in IZZZAAA format
std::vector<double> numden; // atoms/barn-cm
double volume = 5.6; // cm^3
int id = 1234; // material id
std::string name = "material_1234"; // material name
// Use the FakeFactory to get a reasonable initial composition.
FakeFactory::vera_uox_e360( ids, numden );
// Initialize the material.
// Set the beginning of step number densities.
mat->set_numden_bos( numden, ids );
}
// Above this line this is the same as 1CreationOfMaterial
//--------------------------------------------------------------------------
// Create depletion times/fluxes.
// 0 --> 3 days at 1e14 n/cm^2s
// 3 --> 250 days at 1e14 n/cm^2s
// 250 --> 500 days at 1e14 n/cm^2s
// 500 --> 515 days at 0 (decay)
// 515 --> 530 days at 0 (decay)
std::vector<double> time{0, 3, 250, 500, 515, 530};
std::vector<double> flux{1, 1, 1, 0, 0};
for( size_t j = 0; j < time.size(); ++j )
time[j] *= 86400.0; // scale to seconds
for( size_t j = 0; j < flux.size(); ++j )
flux[j] *= 1e14; // scale to n/cm^2s
size_t nsteps = flux.size();
// Create a Bell depletion/decay solver.
Solver_matrex slv;
// Enter step loop.
for( size_t j = 0; j < nsteps; ++j )
{
// Add a new time step to this material.
double dt = time[j + 1] - time[j];
mat->add_step( dt );
mat->set_flux( flux[j] );
mat->set_transition_matrix( trx );
// Solve the step using only data on the material, solve takes
// the bos vector, time, flux, and pointer to eos.
SP_Vec_Dbl n0 = mat->amount_bos();
SP_Vec_Dbl n1 = mat->amount_eos();
slv.set_transition_matrix( &*mat->transition_matrix() );
slv.solve( *n0, mat->flux(), mat->dt(), &*n1 );
slv.clear();
}
// Inspect results stored in material.
if( false ) std::cout << mat->to_string() << std::endl;
}
TEST( Example, 3MultiZoneDepletionByFlux )
{
// Get a general 2237-nuclide ORIGEN library.
FakeFactory::Library_scale_pwr( *lib );
// Extract a transition matrix from it. We will use this PWR transition
// matrix for all materials in this test.
SP_TransitionMatrix trx( lib->newsp_transition_matrix_at( 0 ) );
// Create materials.
size_t num_mats = 4;
{
std::vector<int> ids;
std::vector<double> numden;
Origen::FakeFactory::vera_uox_e360( ids, numden );
int id = 0;
for( size_t m = 0; m < num_mats; ++m )
{
++id;
std::string name = nprintf( "mat%d", id );
double volume = 5.6;
Origen::SP_Material mat(
mat->set_numden_bos( numden, ids );
matmap->add( mat );
}
}
// Create depletion times/fluxes.
// cycle1-------------->
// deplete------>decay->
std::vector<double> time{0, 3, 250, 500, 515, 530};
std::vector<double> flux{1, 1, 1, 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] *= 1e14;
size_t nsteps = flux.size();
// Create multi-zone depleter.
// Set solver to be used.
mzd.set_solver( slv );
// Set material map.
mzd.set_materials( matmap );
// Initialize number of substeps and whether to deplete by power.
// Used in mzd.update(dc) below.
mzd.set_substep_dtrel( {1.0} );
mzd.unset_constraint();
// Enter step loop.
for( size_t j = 0; j < nsteps; ++j )
{
// Add new step data for each material.
mzd.add_step( time[j + 1] - time[j] );
for( auto it = matmap->begin(); it != matmap->end(); ++it )
{
// int id = it->first;
SP_Material mat = it->second;
double rflux = 1.0; // rel. flux in this material
double mat_flux = flux[j] * rflux; // mult. by global flux
mat->set_flux( mat_flux );
mat->set_transition_matrix( trx );
}
// Execute solve for all materials.
mzd.execute();
}
// Have a look at the materials in the map.
if( false ) std::cout << matmap->to_string() << std::endl;
}
TEST( Example, 4MultiZoneDepletionByPower )
{
// Get a general 2237-nuclide ORIGEN library.
FakeFactory::Library_scale_pwr( *lib );
// Extract a transition matrix from it. We will use this PWR transition
// matrix for all materials in this test.
SP_TransitionMatrix trx( lib->newsp_transition_matrix_at( 0 ) );
// Create materials.
size_t num_mats = 4;
{
std::vector<int> ids;
std::vector<double> numden;
Origen::FakeFactory::vera_uox_e360( ids, numden );
int id = 0;
for( size_t m = 0; m < num_mats; ++m )
{
++id;
std::string name = nprintf( "mat%d", id );
double volume = 5.6;
Origen::SP_Material mat(
mat->set_numden_bos( numden, ids );
matmap->add( mat );
}
}
// Create depletion times/fluxes.
// cycle1-------------->
// deplete------>decay->
std::vector<double> time{0, 3, 250, 500, 515, 530};
std::vector<double> flux{1, 1, 1, 0, 0};
for( size_t j = 0; j < time.size(); ++j ) time[j] *= 86400.0;
size_t nsteps = flux.size();
// Create multi-zone depleter.
// Set solver to be used.
mzd.set_solver( slv );
// Set material map.
mzd.set_materials( matmap );
// Initialize number of substeps and whether to deplete by power.
mzd.set_constraint_its_max( 1 );
mzd.set_power_constraint( 7.5e3 );
mzd.set_substep_dtrel( {1.0} ); // single substep
// Enter step loop.
for( size_t j = 0; j < nsteps; ++j )
{
jDebugLine( "# initiating step=" << j );
// Add new step data for each material and update time.
mzd.add_step( time[j + 1] - time[j] );
for( auto it = matmap->begin(); it != matmap->end(); ++it )
{
// int id = it->first;
SP_Material mat = it->second;
double rflux = 1.0; // rel. flux in this material
double mat_flux = flux[j] * rflux; // mult. by global factor
mat->set_flux( mat_flux );
mat->set_transition_matrix( trx );
}
// Execute step solve for all materials.
mzd.execute();
}
// Check the reference solve.
SP_Material mat = matmap->begin()->second;
size_t i = lib->definition().nuclide_set().lookup_index( 30070172 );
EXPECT_FLOAT_EQ( 0.0000000e+00, mat->amount_at( 0 )->at( i ) );
EXPECT_FLOAT_EQ( 1.2144864e-14, mat->amount_at( 1 )->at( i ) );
EXPECT_FLOAT_EQ( 1.1990722e-11, mat->amount_at( 2 )->at( i ) );
EXPECT_FLOAT_EQ( 2.5051308e-11, mat->amount_at( 3 )->at( i ) );
EXPECT_FLOAT_EQ( 2.5380880e-11, mat->amount_at( 4 )->at( i ) );
EXPECT_FLOAT_EQ( 2.5393699e-11, mat->amount_at( 5 )->at( i ) );
// Have a look at the materials in the map.
if( false ) std::cout << matmap->to_string() << std::endl;
}
TEST( Example, DISABLED_5MultiZoneWithTransportCoupling )
{
std::string library_path;
std::string jeff_path;
std::string yields_path;
// Initialize transition matrix updater by paths.
tmu.load( library_path, jeff_path, yields_path );
// Get a general 2237-nuclide ORIGEN library.
// Create 16 materials.
size_t num_mats = 4;
{
std::vector<int> ids;
std::vector<double> numden;
Origen::FakeFactory::vera_uox_e360( ids, numden );
int id = 0;
for( size_t m = 0; m < num_mats; ++m )
{
++id;
std::string name = nprintf( "mat%d", id );
double volume = 5.6;
Origen::SP_Material mat(
mat->set_numden_bos( numden, ids );
matmap->add( mat );
}
}
// Create depletion times.
std::vector<double> time{0, 3, 250, 500};
for( size_t j = 0; j < time.size(); ++j ) time[j] *= 86400.0;
size_t nsteps = time.size() - 1;
// Create multi-zone depleter.
// Set solver to be used.
mzd.set_solver( slv );
// Set material map.
mzd.set_materials( matmap );
// Initialize.
mzd.set_constraint_its_max( 1 );
mzd.set_power_constraint( 7.5e3 );
mzd.set_substep_dtrel( {1.0} ); // single substep
// Enter step loop.
for( size_t j = 0; j < nsteps; ++j )
{
// Add new step data for each material.
mzd.add_step( time[j + 1] - time[j] );
for( auto it = matmap->begin(); it != matmap->end(); ++it )
{
// int id = it->first;
SP_Material mat = it->second;
// In a real simulation, this flux come from the transport solution.
Vec_Flt mg_flux( 252, 1.0 );
// Now we can update specific one group cross sections to
// supplement/override data set from just the flux.
Vec_Int nuclide_list = {92235, 94239};
Vec_Int num_mts = {1, 2};
Vec_Int mt_list = {102, 102, 18};
Vec_Flt xs = {17.0, 9., 100.};
// Get a new transition matrix.
tmu.get_transition_matrix( &*trx );
mat->set_flux( rflux );
mat->set_transition_matrix( trx );
}
// Execute solve for all materials.
mzd.execute();
}
// Have a look at the materials in the map.
if( false ) std::cout << matmap->to_string() << std::endl;
}
TEST( Example, DISABLED_6MultiZoneWithAbsorptionTreatment )
{
std::string library_path;
std::string jeff_path;
std::string yields_path;
// Initialize transition matrix updater by paths.
tmu.load( library_path, jeff_path, yields_path );
// Get a general 2237-nuclide ORIGEN library.
SP_Library lib = tmu.sp_library();
// Create 16 materials.
size_t num_mats = 4;
{
std::vector<int> ids;
std::vector<double> numden;
Origen::FakeFactory::vera_uox_e360( ids, numden );
int id = 0;
for( size_t m = 0; m < num_mats; ++m )
{
++id;
std::string name = nprintf( "mat%d", id );
double volume = 5.6;
Origen::SP_Material mat(
mat->set_numden_bos( numden, ids );
matmap->add( mat );
}
}
// Create depletion times.
std::vector<double> time{0, 3, 250, 500};
for( size_t j = 0; j < time.size(); ++j ) time[j] *= 86400.0;
size_t nsteps = time.size() - 1;
// Create multi-zone depleter.
// Set solver to be used.
mzd.set_solver( slv );
// Set material map.
mzd.set_materials( matmap );
mzd.set_constraint_its_max( 1 );
mzd.set_power_constraint( 7.5e3 );
mzd.set_substep_dtrel( {1.0} ); // single substep
// Enter step loop.
for( size_t j = 0; j < nsteps; ++j )
{
// Add new step data for each material.
mzd.add_step( time[j + 1] - time[j] );
for( auto it = matmap->begin(); it != matmap->end(); ++it )
{
// int id = it->first;
SP_Material mat = it->second;
// In a real simulation, this flux come from the transport solution.
Vec_Flt mg_flux( 252, 1.0 );
// This method performs a special fission/absorption xs only
// update. However, most of the important self-shielding effects
// can be taken into account with this approach.
Vec_Int nuclide_list = {92235, 94239};
Vec_Flt fis_xs = {40.0, 100.0};
Vec_Flt abs_xs = {57.0, 109.0};
"IZZZAAA", nuclide_list, fis_xs, abs_xs );
// Get a new transition matrix.
tmu.get_transition_matrix( &*trx );
mat->set_flux( rflux );
mat->set_transition_matrix( trx );
}
// Execute solve for all materials.
mzd.execute();
}
// Have a look at the materials in the map.
if( false ) std::cout << matmap->to_string() << std::endl;
}
TEST( Example, 7MultiZoneWithMPI )
{
// Get a general 2237-nuclide ORIGEN library.
FakeFactory::Library_scale_pwr( *lib );
// Extract a transition matrix from it. We will use this PWR transition
// matrix for all materials in this test.
SP_TransitionMatrix trx( lib->newsp_transition_matrix_at( 0 ) );
// Create materials.
size_t num_mats = 4;
{
std::vector<int> ids;
std::vector<double> numden;
Origen::FakeFactory::vera_uox_e360( ids, numden );
int id = 0;
for( size_t m = 0; m < num_mats; ++m )
{
++id;
std::string name = nprintf( "mat%d", id );
double volume = 5.6;
Origen::SP_Material mat(
mat->set_numden_bos( numden, ids );
matmap->add( mat );
}
}
// Create depletion times/fluxes.
// cycle1-------------->
// deplete------>decay->
std::vector<double> time{0, 3, 250, 500, 515, 530};
std::vector<double> flux{1, 1, 1, 0, 0};
for( size_t j = 0; j < time.size(); ++j ) time[j] *= 86400.0;
size_t nsteps = flux.size();
// Create multi-zone depleter.
// Set solver to be used.
mzd.set_solver( slv );
// Set material map.
mzd.set_materials( matmap );
// Set communicator.
mzd.set_constraint_its_max( 1 );
mzd.set_power_constraint( 7.5e3 );
mzd.set_substep_dtrel( {1.0} ); // single substep
// Enter step loop.
for( size_t j = 0; j < nsteps; ++j )
{
// Add new step data for each material.
mzd.add_step( time[j + 1] - time[j] );
for( auto it = matmap->begin(); it != matmap->end(); ++it )
{
// int id = it->first;
SP_Material mat = it->second;
double rflux = 1.0; // rel. flux in this material
double mat_flux = flux[j] * rflux; // mult. by global factor
mat->set_flux( mat_flux );
mat->set_transition_matrix( trx );
}
if( flux[j] == 0.0 ) mzd.set_power_constraint( 0.0 );
// Execute solve for all materials.
mzd.execute();
// The status object is useful for getting parallel info.
}
}
TEST( Example, 8ExternallyDriven ) {}
