tstTransitionMatrixUpdater.cpp

./Manager/libld/tests/tstTransitionMatrixUpdater.cpp

#include <memory.h>
#include <fstream>
#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/ReactionResourceIO.h"
#include "Origen/Core/io/YieldResourceIO.h"
#include "Origen/Core/ts/Channel.h"
#include "Origen/Manager/libld/TransitionMatrixUpdater.h"
#include "Origen/Manager/libld/TransitionSystemAdapter_AmpxN.h"

namespace Origen
{
namespace test
{

TEST( TransitionMatrixUpdater, Load )
{
    Origen::TransitionMatrixUpdater tmu;
    tmu.load( Origen::pwrlib_filepath,
              Origen::jeff252_filepath,
              Origen::yields_filepath );

    if( !tmu.initialized() )
    {
        tmu.clear();  // clears errors
        return;
    }
    // tmu.mg_database()->print(std::cout);
}

TEST( TransitionMatrixUpdater, LoadUpdate )
{
    Origen::TransitionMatrixUpdater tmu;
    tmu.load( Origen::pwrlib_filepath,
              Origen::jeff252_filepath,
              Origen::yields_filepath );

    if( !tmu.initialized() )
    {
        tmu.clear();  // clears errors
        return;
    }

    Origen::Vec_Flt mg_flux( 252, 1.0 );
    tmu.update_mg_flux( mg_flux );
}

TEST( TransitionMatrixUpdater, LoadUpdateGet )
{
    Origen::TransitionMatrixUpdater tmu;
    tmu.load( Origen::pwrlib_filepath,
              Origen::jeff252_filepath,
              Origen::yields_filepath );

    if( !tmu.initialized() )
    {
        tmu.clear();  // clears errors
        return;
    }

    Origen::Vec_Flt mg_flux( 252, 1.0 );
    tmu.update_mg_flux( mg_flux );

    Origen::TransitionMatrixP trx;
    tmu.get_transition_matrix( &trx );
}

TEST( TransitionMatrixUpdater, LoadUpdateGetClear )
{
    Origen::TransitionMatrixUpdater tmu;
    tmu.load( Origen::pwrlib_filepath,
              Origen::jeff252_filepath,
              Origen::yields_filepath );

    if( !tmu.initialized() )
    {
        tmu.clear();  // clears errors
        return;
    }

    Origen::Vec_Flt mg_flux( 252, 1.0 );
    tmu.update_mg_flux( mg_flux );

    Origen::TransitionMatrixP trx;
    tmu.get_transition_matrix( &trx );

    tmu.clear();
}

// for a uniform flux=1
std::map<int, std::map<int, double>> reference_sb47g_constants()
{
    std::map<int, std::map<int, double>> r;

    r[92235][16] = 0.25605130195617676;
    r[92235][17] = 0.014363746158778667;
    r[92235][18] = 1.4366118907928467;
    r[92235][102] = 0.023061279207468033;
    r[92238][16] = 0.38620460033416748;
    r[92238][17] = 0.062143206596374512;
    r[92238][18] = 0.74442368745803833;
    r[92238][102] = 0.01585003174841404;
    r[94240][16] = 0.10783249884843826;
    r[94240][17] = 0.027231113985180855;
    r[94240][18] = 1.8240015506744385;
    r[94240][102] = 0.022608127444982529;

    return r;
}

// for a uniform flux=1
std::map<int, std::map<int, double>> reference_ub47g_constants()
{
    std::map<int, std::map<int, double>> r;

    r[92235][16] = 0.3402692973613739;
    r[92235][17] = 0.038585104048252106;
    r[92235][18] = 1.6015549898147583;
    r[92235][102] = 0.022057455033063889;
    r[92238][16] = 0.57465583086013794;
    r[92238][17] = 0.16693426668643951;
    r[92238][18] = 0.84844231605529785;
    r[92238][102] = 0.015442271716892719;
    r[94240][16] = 0.18764506280422211;
    r[94240][17] = 0.073150500655174255;
    r[94240][18] = 1.9431796073913574;
    r[94240][102] = 0.021005298942327499;

    return r;
}

TEST( TransitionMatrixUpdater, StorageBounds47g )
{
    // get 47 group bounds
    Origen::Vec_Dbl mpact_47g_dbl;
    Origen::FakeFactory::energy_bounds_mpact_47g( mpact_47g_dbl );
    Origen::Vec_Flt mpact_47g;
    mpact_47g.assign( mpact_47g_dbl.begin(), mpact_47g_dbl.end() );

    // uniform flux spectrum
    Origen::Vec_Flt mg_flux( 47 );
    Origen::FakeFactory::getUniformFlux( &mg_flux, mpact_47g );

    // set storage bounds and load the data
    Origen::TransitionMatrixUpdater tmu( mpact_47g );
    tmu.load( Origen::pwrlib_filepath,
              Origen::jeff252_filepath,
              Origen::yields_filepath );

    if( !tmu.initialized() )
    {
        tmu.clear();  // clears errors
        return;
    }

    jDebugBlock( {
        std::ofstream ofs( "mg.txt" );
        tmu.mg_database()->print( ofs );
    } );

    // update flux
    tmu.update_mg_flux( mg_flux );

    Origen::TransitionMatrixP trx;
    tmu.get_transition_matrix( &trx );
    // spot check values
    for( auto x1 : reference_sb47g_constants() )
    {
        int id = x1.first;
        for( auto x2 : x1.second )
        {
            int mt = x2.first;
            double ref = x2.second;
            EXPECT_NEAR(
                ref,
                tmu.transition_system()->nuclide( id )->loss( mt )->xs(),
                1e-5 * ref );
        }
    }

    jDebugBlock( {
        std::ofstream ofs( "StorageBounds47g.trx.txt" );
        ofs << trx.toQString( 3, "F", 0, 0 ).toStdString() << std::endl;
    } );

    tmu.clear();
}

TEST( TransitionMatrixUpdater, UserBounds47g )
{
    // use default storage bounds and load the data
    Origen::TransitionMatrixUpdater tmu;
    tmu.load( Origen::pwrlib_filepath,
              Origen::jeff252_filepath,
              Origen::yields_filepath );

    if( !tmu.initialized() )
    {
        tmu.clear();  // clears errors
        return;
    }

    jDebugBlock( {
        std::ofstream ofs( "mg.txt" );
        tmu.mg_database()->print( ofs );
    } );

    // get 47g user bounds
    Origen::Vec_Dbl mpact_47g_dbl;
    Origen::FakeFactory::energy_bounds_mpact_47g( mpact_47g_dbl );
    Origen::Vec_Flt mpact_47g;
    mpact_47g.assign( mpact_47g_dbl.begin(), mpact_47g_dbl.end() );
    tmu.set_user_bounds( mpact_47g );

    // uniform flux spectrum
    Origen::Vec_Flt mg_flux( 47 );
    Origen::FakeFactory::getUniformFlux( &mg_flux, mpact_47g );

    // update
    tmu.update_mg_flux( mg_flux );

    Origen::TransitionMatrixP trx;
    tmu.get_transition_matrix( &trx );

    // spot check values
    for( auto x1 : reference_ub47g_constants() )
    {
        int id = x1.first;
        for( auto x2 : x1.second )
        {
            int mt = x2.first;
            double ref = x2.second;
            EXPECT_NEAR(
                ref,
                tmu.transition_system()->nuclide( id )->loss( mt )->xs(),
                1e-5 * ref );
        }
    }

    jDebugBlock( {
        std::ofstream ofs( "UserBounds47g.trx.txt" );
        ofs << trx.toQString( 3, "F", 0, 0 ).toStdString() << std::endl;
    } );

    tmu.clear();
}

TEST( TransitionMatrixUpdater, FetchNuclides )
{
    Origen::TransitionMatrixUpdater tmu;
    tmu.load( Origen::pwrlib_filepath,
              Origen::jeff252_filepath,
              Origen::yields_filepath );

    if( !tmu.initialized() )
    {
        tmu.clear();  // clears errors
        return;
    }

    // some checks on providing access to skeleton
    Origen::SP_TransitionMatrixP base_trx = tmu.base_transition_matrix();
    ASSERT_TRUE( base_trx != nullptr );
    EXPECT_EQ( 2237, base_trx->get_itot() );
    std::vector<Origen::SP_Nuclide> ts_nuclides;

    // Test ability to fetch when FP have +5000 AAAI
    {
        std::string id_form = "ZZZAAAI+FP5000";
        std::vector<int> n_list;
        n_list.push_back( 922350 );
        n_list.push_back( 636510 );
        n_list.push_back( 10020 );
        n_list.push_back( 922351 );

        tmu.fetch_ts_nuclides( &ts_nuclides, id_form, n_list );
        ASSERT_EQ( 4u, ts_nuclides.size() );
        for( size_t i = 0; i < ts_nuclides.size(); ++i )
        {
            ASSERT_TRUE( ts_nuclides[i] != nullptr );
        }
        EXPECT_EQ( 92235, ts_nuclides[0]->id() );
        EXPECT_EQ( 63151, ts_nuclides[1]->id() );
        EXPECT_EQ( 1002, ts_nuclides[2]->id() );
        EXPECT_EQ( 1092235, ts_nuclides[3]->id() );
    }

    // Test ability to fetch ZZZAAAI
    {
        std::string id_form = "ZZZAAAI";
        std::vector<int> n_list;
        n_list.push_back( 922350 );
        n_list.push_back( 631510 );
        n_list.push_back( 10020 );
        n_list.push_back( 922351 );

        tmu.fetch_ts_nuclides( &ts_nuclides, id_form, n_list );
        ASSERT_EQ( 4u, ts_nuclides.size() );
        for( size_t i = 0; i < ts_nuclides.size(); ++i )
        {
            ASSERT_TRUE( ts_nuclides[i] != nullptr );
        }
        EXPECT_EQ( 92235, ts_nuclides[0]->id() );
        EXPECT_EQ( 63151, ts_nuclides[1]->id() );
        EXPECT_EQ( 1002, ts_nuclides[2]->id() );
        EXPECT_EQ( 1092235, ts_nuclides[3]->id() );
    }

    // Test ability to fetch IZZZAAA
    {
        std::string id_form = "IZZZAAA";
        std::vector<int> n_list;
        n_list.push_back( 92235 );
        n_list.push_back( 63151 );
        n_list.push_back( 1002 );
        n_list.push_back( 1092235 );

        tmu.fetch_ts_nuclides( &ts_nuclides, id_form, n_list );
        ASSERT_EQ( 4u, ts_nuclides.size() );
        for( size_t i = 0; i < ts_nuclides.size(); ++i )
        {
            ASSERT_TRUE( ts_nuclides[i] != nullptr );
        }
        EXPECT_EQ( 92235, ts_nuclides[0]->id() );
        EXPECT_EQ( 63151, ts_nuclides[1]->id() );
        EXPECT_EQ( 1002, ts_nuclides[2]->id() );
        EXPECT_EQ( 1092235, ts_nuclides[3]->id() );
    }
}

TEST( TransitionMatrixUpdater, LoadWithUserReactions )
{
    /*
     Try to lad the libraries adding in the following user-defined
     reaction channels:
        H-1 (n,g)     // A reaction that already exists
        H-1 (n,2n)    // A reaction that does not exist (yes its silly)
        H-2  ---      // Nuclide with zero MTs is allowed
        B-11 (n,2n)   // A reaction with nothing special in it.
        B-11 (n,na)   // A reaction which (at the time of writing this) is
                      // in .orglib but not in the MG data.
        U-235 (n,f)   // This should not break fission.
    */

    std::vector<int> user_nuclide_list = {1001, 1002, 5011, 92235};
    std::vector<int> user_num_mts = {2, 0, 2, 1};
    std::vector<int> user_mt_list = {102, 16, 16, 22, 18};

    // ZAI corresponding to each element of the user_mt_list.
    std::vector<int> expanded_zai = {10010, 10010, 50110, 50110, 922350};

    Origen::TransitionMatrixUpdater tmu;
    tmu.load( Origen::pwrlib_filepath,
              Origen::jeff252_filepath,
              Origen::yields_filepath,
              user_nuclide_list,
              user_num_mts,
              user_mt_list );

    if( !tmu.initialized() )
    {
        tmu.clear();  // clears errors
        return;
    }

    // Updating MG flux once is required even if you intend to set one-group.
    Origen::Vec_Flt mg_flux( 252, 1.0 );
    tmu.update_mg_flux( mg_flux );

    // Get the cross-sections.
    Origen::Vec_Flt xs( user_mt_list.size(), 0.0 );

    tmu.get_xs( "IZZZAAA", user_nuclide_list, user_num_mts, user_mt_list, &xs );

    // New reactions are zero, the pre-existing ones with mg data non-zero.
    EXPECT_TRUE( xs[0] != 0.0 );
    EXPECT_TRUE( xs[1] == 0.0 );
    EXPECT_TRUE( xs[2] != 0.0 );
    EXPECT_TRUE( xs[3] == 0.0 );
    EXPECT_TRUE( xs[2] != 0.0 );

    // Set the cross-sections.
    for( int i = 0; i < xs.size(); ++i ) xs[i] = ( i + 1 ) * 1000;

    tmu.update_user_xs(
        "IZZZAAA", user_nuclide_list, user_num_mts, user_mt_list, xs );

    // Data was copied so this won't affect the values in tmu.
    for( int i = 0; i < xs.size(); ++i ) xs[i] = 0;

    // Get the cross-sections.
    tmu.get_xs( "IZZZAAA", user_nuclide_list, user_num_mts, user_mt_list, &xs );

    for( int i = 0; i < xs.size(); ++i )
        EXPECT_DOUBLE_EQ( ( i + 1 ) * 1000, xs[i] );

    // >>>  Check that the altered values made it to the TransitionMatrix
    Origen::TransitionMatrixP trx;
    tmu.get_transition_matrix( &trx );

    DoubleList diag, offdiag;
    trx.fold_flux( 1.0e24,
                   diag,
                   offdiag );  // The 1e24 flux causes the elements
    // to be in xs units.

    /************************************************************************/
    // The get_xs function in transition matrix is broken and intended to be
    // removed so it should not be used. Instead get the actual transition
    // matrix and check values in it.

    // Convert the Lists to vectors
    int itot = trx.get_itot();
    int non = trx.get_non();
    std::vector<int> non0( itot ), kd( itot ), loc( non ), ZAI( itot );

    copy( trx.get_non0_()->begin(), trx.get_non0_()->end(), non0.begin() );
    copy( trx.get_kd_()->begin(), trx.get_kd_()->end(), kd.begin() );
    copy( trx.get_loc()->begin(), trx.get_loc()->end(), loc.begin() );
    copy( trx.get_nucl()->begin(), trx.get_nucl()->end(), ZAI.begin() );

    // >>> Convert Origen form to dense matrix

    double** B = (double**)calloc( itot, sizeof( double* ) );
    for( int i = 0; i < itot; i++ )
        B[i] = (double*)calloc( itot, sizeof( double ) );

    // Add diagonal elements
    for( int i = 0; i < itot; i++ ) B[i][i] += diag[i];

    // Add offdiagonal elements
    {
        int j = 0;
        for( int i = 0; i < itot; i++ )  // rows (parents)
            for( int p = 0; p < non0[i]; p++ )
            {
                B[i][loc[j]] += offdiag[j];
                j++;
            }
    }

    // >>> Now B is the transition matrix in a dense format.

    int idx_H1 = trx.find_nuclide( 1, 10010 );
    int idx_B11 = trx.find_nuclide( 1, 50110 );
    int idx_U235 = trx.find_nuclide( 2, 922350 );

    // The massive flux we set means that the xs should dominate diagonal

    EXPECT_NEAR( B[idx_H1][idx_H1], -3000, 0.001 * 3000 );
    EXPECT_NEAR( B[idx_B11][idx_B11], -7000, 0.001 * 7000 );
    EXPECT_NEAR( B[idx_U235][idx_U235], -5000, 0.01 * 5000 );

    // Check the daughters

    EXPECT_FLOAT_EQ( B[trx.find_nuclide( 1, 10020 )][idx_H1], 1000 );
    EXPECT_FLOAT_EQ( B[trx.find_nuclide( 1, 50100 )][idx_B11], 3000 );
    EXPECT_FLOAT_EQ( B[trx.find_nuclide( 1, 30070 )][idx_B11], 4000 );

    // For fission, check that the total production from U-235 is close

    double tot_yield = 0.0;
    for( int i = 0; i < itot; ++i )
        if( i != idx_U235 ) tot_yield += B[i][idx_U235];
    EXPECT_NEAR( tot_yield, 10000, 0.01 * 10000 );

    /************************************************************************/
    // Finally, check that we did not break the multigroup functionality

    for( int i = 25; i < 252; ++i ) mg_flux[i] = 0;  // softer spectrum

    tmu.update_mg_flux( mg_flux );
    tmu.get_transition_matrix( &trx );
    trx.fold_flux( 1.0e24, diag, offdiag );

    copy( trx.get_non0_()->begin(), trx.get_non0_()->end(), non0.begin() );
    copy( trx.get_kd_()->begin(), trx.get_kd_()->end(), kd.begin() );
    copy( trx.get_loc()->begin(), trx.get_loc()->end(), loc.begin() );
    copy( trx.get_nucl()->begin(), trx.get_nucl()->end(), ZAI.begin() );

    /*** Convert Origen form to dense matrix ***/

    double** C = (double**)calloc( itot, sizeof( double* ) );
    for( int i = 0; i < itot; i++ )
        C[i] = (double*)calloc( itot, sizeof( double ) );

    // Add diagonal elements
    for( int i = 0; i < itot; i++ ) C[i][i] += diag[i];

    // Add offdiagonal elements
    {
        int j = 0;
        for( int i = 0; i < itot; i++ )  // rows (parents)
            for( int p = 0; p < non0[i]; p++ )
            {
                C[i][loc[j]] += offdiag[j];
                j++;
            }
    }

    // >>> Now C is the new (softer spectrum) transition Matrix in dence form.
    // Check that things changed...

    EXPECT_TRUE( B[idx_H1][idx_H1] != C[idx_H1][idx_H1] );
    EXPECT_TRUE( B[idx_B11][idx_B11] != C[idx_B11][idx_B11] );
    EXPECT_TRUE( B[idx_U235][idx_U235] != C[idx_U235][idx_U235] );

    // B-11 n,2n should change too as it has MG data
    EXPECT_TRUE( B[trx.find_nuclide( 1, 50100 )][idx_B11] !=
                 C[trx.find_nuclide( 1, 50100 )][idx_B11] );

    // B-11 n,na should not change as it has no MG data
    EXPECT_FLOAT_EQ( C[trx.find_nuclide( 1, 30070 )][idx_B11], 4000 );
    EXPECT_TRUE( -C[idx_B11][idx_B11] > 4000 );
    EXPECT_TRUE( -C[idx_B11][idx_B11] < 4050 );

    // Check that the fission yields still update.
    tot_yield = 0.0;
    for( int i = 0; i < itot; ++i )
        if( i != idx_U235 ) tot_yield += C[i][idx_U235];
    EXPECT_NEAR( tot_yield, -2 * C[idx_U235][idx_U235], 0.1 * tot_yield );

    // Free memory.

    for( int i = 0; i < itot; i++ ) free( B[i] );
    free( B );

    for( int i = 0; i < itot; i++ ) free( C[i] );
    free( C );

    tmu.clear();
}

std::string int2str( int i )
{
    std::stringstream ss;
    ss << i;
    return ss.str();
}

TEST( TransitionMatrixUpdater, Example )
{
    //     nemesis::Timer t = nemesis::Timer();
    //     std::vector<std::string> errors;

    Origen::TransitionSystemAdapter_AmpxN tsa;
    {
        ScaleUtils::IO::DB opts;
        //
        // Load AMPX infinitely dilute library into memory.
        //
        // TIC();
        auto rr = std::make_shared<Origen::ReactionResource>();
        {
            Origen::ReactionResourceIO rr_io;
            rr_io.load( *rr, Origen::jeff252_filepath, opts );
        }
        ASSERT_TRUE( rr != nullptr ) << opts.writeErrorStack();
        // TOC("loaded AMPX Infinitely Dilute library");
        // mempause std::cin.get(); //0.06s 16 MB

        //
        // Load the yield resource.
        //
        // TIC();
        auto yr = std::make_shared<Origen::YieldResource>();
        {
            Origen::YieldResourceIO yr_io;
            yr_io.load( *yr, Origen::yields_filepath, opts );
        }
        ASSERT_TRUE( yr != nullptr ) << opts.writeErrorStack();
        // TOC("loaded fission yield resource");
        // mempause std::cin.get(); //0.04s 18 MB

        // TIC();
        rr->import_fission_yields( *yr );
        // TOC("loaded internal YR");
        // mempause std::cin.get(); //0.03s 35 MB

        //
        // Import .
        //
        // TIC();
        tsa.set_reaction_resource( rr );
        // TOC("loaded internal AMPX");
        // mempause std::cin.get(); //0.02s 31 MB
    }

    // TIC();
    // TransitionSystemMole ts; //use the Mole to print info while updating
    Origen::TransitionSystem_Gen ts;
    // TOC("created empty transition system");
    // mempause std::cin.get(); //0.00s 25 MB

    // TIC();
    tsa.extend( &ts );
    // TOC("extended transition system");
    // mempause std::cin.get(); //0.04s 34 MB

    // TIC();
    tsa.extend( &ts );
    // TOC("re-extended transition system");
    // mempause std::cin.get(); //0.05s 34 MB

    // tsa.print(); //print tsa internals

    // TIC();
    // std::vector<double> mg_flux0 = TestHelper::get_test_flux(1.0,252);
    std::vector<float> mg_flux( 252, 1.0 );
    // mg_flux.assign(mg_flux0.begin(),mg_flux0.end());
    // TOC("setup test data");
    // mempause std::cin.get(); //0.00s 35 MB

    {
        // TIC();
        // ts.header(); //uncomment when using the Mole
        float sumflux = tsa.update( &ts, mg_flux );
        EXPECT_FLOAT_EQ( 252, sumflux );

        // ts.footer(); //uncomment when using the Mole
        jDebugLine( "ID sumflux=" << sumflux );
        // TOC("updated transition system with internal AMPX data");
        // mempause std::cin.get(); //0.04s 35 MB
    }

    if( false )
    {
        //
        // Load AMPX self-shielded library into memory.
        //
        // TIC();
        AmpxWorkingLibrary ssnxs;
        ssnxs.setFileName( "ft04f001" );
        {
            int rtn = ssnxs.open();
            ASSERT_TRUE( rtn != 0 );
        }
        // TOC("loaded AMPX Self-Shielded library");
        // mempause std::cin.get(); //0.01s 35 MB

        // TIC();
        std::ofstream os;
        os.open( "ssnxs.txt" );
        tsa.print( ssnxs, os );
        os.close();
        // TOC("wrote file={ft04f001} to ssnxs.txt");

        int mix_id = 25;

        // TIC();
        tsa.extend_by_mix( &ts, ssnxs, mix_id );
        // TOC("extended transition system for mix_id="<<int2str(mix_id));
        // mempause std::cin.get(); //0.03s 36 MB

        // TIC();
        // ts.header(); //uncomment when using the Mole
        float sumflux = tsa.update_by_mix( &ts, ssnxs, mix_id );
        EXPECT_FLOAT_EQ( 252, sumflux );

        // ts.footer(); //uncomment when using the Mole
        jDebugLine( "SS sumflux=" << sumflux );
        // TOC("updated transition system with internal AMPX data and additional
        // library");
        // mempause std::cin.get(); //0.0s 36 MB
    }
}

}}