tstTransitionSystemAdapter_DecayResource.cpp

./Manager/libld/tests/tstTransitionSystemAdapter_DecayResource.cpp

#include <fstream>
#include <iostream>

#include "Nemesis/gtest/nemesis_gtest.hh"
#include "Origen/Core/Definitions.h"
#include "Origen/Core/TestPaths.h"
#include "Origen/Core/dc/FakeFactory.h"
#include "Origen/Core/dc/Library.h"
#include "Origen/Core/dc/LibraryHeader.h"
#include "Origen/Core/dc/NuclideSet.h"
#include "Origen/Core/dc/TransitionCoeff.h"
#include "Origen/Core/dc/TransitionStructure.h"
#include "Origen/Core/fn/diff.h"
#include "Origen/Core/fn/populate.h"
#include "Origen/Core/fn/print.h"
#include "Origen/Core/io/DecayResourceIO.h"
#include "Origen/Core/io/LibraryIO.h"
#include "Origen/Core/re/DecayResource.h"
#include "Origen/Core/re/defaultNuclideResource.h"
#include "Origen/Core/ts/TransitionSystem_Gen.h"
#include "Origen/Manager/libld/LibraryBuilder.h"
#include "Origen/Manager/libld/TransitionSystemAdapter_DecayResource.h"
#include "ScaleUtils/IO/DB.h"

class TSA_DecayResourceTester : public ::testing::Test
{
  public:
    void SetUp()
    {
        // Load a DecayResource from a file.
        Origen::DecayResource dr;
        Origen::DecayResourceIO dr_io;
        ScaleUtils::IO::DB opts;
        dr_io.load( dr, Origen::decay_filepath, opts );

        // Create the TransitionSystemAdapter for the decay resource and
        // populate a transition system from it.
        Origen::TransitionSystemAdapter_DecayResource tsa;
        ts = Origen::SP_TransitionSystem( new Origen::TransitionSystem_Gen() );
        tsa.extend( &*ts, true, dr );
    }

    Origen::SP_TransitionSystem ts;
};

TEST_F( TSA_DecayResourceTester, Single1003 )
{
    // output a nuclide in detail
    const bool debug = false;
    if( debug )
    {
        std::ofstream ofs( "1003.json" );
        ofs << ts->nuclide( 1003 )->to_string();
    }

    // create a list
    std::vector<int> sizzzaaa_list{10001003, 10002003};

    // creating decay library from transition system initialized in SetUp
    Origen::Library lib;
    Origen::LibraryBuilder bld;
    bld.create_decay_library( &lib, sizzzaaa_list, true, *ts );

    // check that ids made it to library
    EXPECT_EQ( 10001003, lib.definition().nuclide_set().ids_at( 0 ) );
    EXPECT_EQ( 10002003, lib.definition().nuclide_set().ids_at( 1 ) );

    // check decay constants
    EXPECT_NEAR( 1.7829e-09, lib.decay_data().decay_constants.at( 0 ), 1e-10 );
    EXPECT_EQ( 0.0, lib.decay_data().decay_constants.at( 1 ) );

    // check transition ids
    EXPECT_EQ( 1, lib.transition_structure().transition_ids_size() );
    EXPECT_EQ( -2, lib.transition_structure().transition_ids_at( 0 ) );

    // check transition coeff
    EXPECT_EQ( 1, lib.transition_coeff_size() );
    EXPECT_NEAR(
        1.7829e-09, lib.transition_coeff_at( 0 ).matrix_at( 0 ), 1e-10 );

    // output library
    if( debug )
    {
        std::ofstream ofs( "Single1003.txt" );
        ofs << lib.to_string() << std::endl;
        ScaleUtils::IO::DB opts;
        saveLibrary( lib, "Single1003.f33", opts );
    }
}

// TODO : Re-enable test after library updated!
TEST_F( TSA_DecayResourceTester, DISABLED_FullCompare )
{
    std::vector<int> sizzzaaa_list =
        Origen::FakeFactory::createNuclideIds_scale62();
    EXPECT_EQ( 2240u, sizzzaaa_list.size() );

    Origen::Library lib;
    Origen::LibraryBuilder bld;

    // creating decay library from transition system initialized in SetUp
    bld.create_decay_library( &lib, sizzzaaa_list, true, *ts );

    Origen::LibraryIO lib_io;
    ScaleUtils::IO::DB opts;

    Origen::Library lib_orig;
    bool loaded = lib_io.load( lib_orig, Origen::decaylib_filepath, opts );
    EXPECT_TRUE( loaded );

    auto diff = diffLibrary( lib_orig, lib );
    EXPECT_FALSE( diff->have_difference() )
        << "    outputting only differing lines to screen...\n"
        << "    outputting old/new libraries to lib_orig.f33/lib.f33...\n"
        << "    outputting full diff to file lib.diff.txt\n";

    if( diff->have_difference() )
    {
        diff->to_string( true ); /*diff only*/
        std::ofstream ofs( "lib.diff.txt" );
        ofs << diff->to_string();
        lib_io.save( lib_orig, "lib_orig.f33", opts );
        lib_io.save( lib, "lib.f33", opts );
    }
}

// This function uses GTEST macros to test various sizes on the library
// for consistency.
void gtestLibrarySizes( const Origen::Library& lib )
{
    // get various components
    const Origen::LibraryHeader& def = lib.definition();
    const Origen::NuclideSet& nucset = def.nuclide_set();
    const Origen::DecayData& dd = lib.decay_data();
    const Origen::TransitionStructure& tstruct = lib.transition_structure();

    // get reference numbers of stuff
    size_t total = lib.total_nuclides();
    size_t num_lt = nucset.num_in_sublib( Origen::SUBLIB_1LT );
    size_t num_ac = nucset.num_in_sublib( Origen::SUBLIB_2AC );
    size_t num_fp = nucset.num_in_sublib( Origen::SUBLIB_3FP );
    ASSERT_LT( 0, num_lt );
    ASSERT_LT( 0, num_ac );
    ASSERT_LT( 0, num_fp );

    // test number of nuclides
    EXPECT_EQ( total, nucset.ids_size() );
    EXPECT_EQ( total, nucset.total_nuclides() );
    EXPECT_EQ( total, def.total_nuclides() );

    // test numbers in decay data
    EXPECT_EQ( total, dd.decay_constants.size() );
    EXPECT_EQ( total, dd.recoverable_energy_values.size() );
    EXPECT_EQ( total, dd.q_fractions_from_photons.size() );
    EXPECT_EQ( total, dd.rcg_for_air.size() );
    EXPECT_EQ( total, dd.rcg_for_water.size() );
    EXPECT_EQ( num_lt, dd.abund_of_lite_nuclides.size() );
    EXPECT_LE( 0, dd.fissionable_nuclide_ids.size() );
    EXPECT_EQ( total, dd.masses.size() );

    // test number of transitions
    size_t non = tstruct.total_transitions();
    EXPECT_EQ( non, tstruct.transition_ids_size() );
    EXPECT_EQ( non, tstruct.parent_positions_size() );
    EXPECT_EQ( total, tstruct.num_parents_size() );
    EXPECT_EQ( total, tstruct.num_decay_parents_size() );

    size_t nfluxgrp = 0;
    if( lib.transition_coeff_size() > 0 )
    {
        nfluxgrp = lib.transition_coeff_at( 0 ).flux_size();
    }
    for( size_t i = 0; i < lib.transition_coeff_size(); ++i )
    {
        const Origen::TransitionCoeff& tcoeff = lib.transition_coeff_at( i );
        EXPECT_EQ( non, tcoeff.matrix_size() );
        EXPECT_EQ( total, tcoeff.loss_xs_size() );
        EXPECT_EQ( num_ac, tcoeff.fission_xs_size() );
        EXPECT_EQ( num_ac + num_lt, tcoeff.neutron_yields_size() );
        EXPECT_EQ( nfluxgrp, tcoeff.flux_size() );
    }
}

// This function uses GTEST macros to test that masses are something
// that's close but not equal to the mass number.
void gtestLibraryMasses( const Origen::Library& lib )
{
    // get various components
    const Origen::LibraryHeader& def = lib.definition();
    const Origen::NuclideSet& nucset = def.nuclide_set();
    const Origen::DecayData& dd = lib.decay_data();

    EXPECT_EQ( nucset.total_nuclides(), dd.masses.size() );

    for( size_t i = 0; i < dd.masses.size(); ++i )
    {
        int izzzaaa = nucset.izzzaaa_at( i );
        int a = nucset.a_at( i );
        double da = std::abs( dd.masses[i] - a );

        // check that it's not zero UNLESS it's C-12
        if( izzzaaa != 6012 && izzzaaa != 1) EXPECT_LT( 1e-5, da ) << nucset.izzzaaa_at( i );
        if( izzzaaa == 6012 && izzzaaa != 1) EXPECT_GT( 1e-6, da ) << nucset.izzzaaa_at( i );

        // check that it's not more than 1
        if (izzzaaa != 1) EXPECT_GT( 1.0, da );
    }
}

class LibraryBuilderTester : public ::testing::Test
{
  public:
    void SetUp()
    {
        //
        // load decay resource
        //
        ScaleUtils::IO::DB opts;
        std::string decay_resource( Origen::decay_filepath );
        dr = Origen::SP_DecayResource(
            Origen::loadDecayResource( decay_resource, opts ) );
        Insist( opts.numError() == 0, "decay resource had errors" );

        //
        // load nuclide resource
        //
        nr = std::make_shared<Origen::NuclideResource>();
        Origen::defaultNuclideResource( &*nr );

        //
        // set on library builder
        //
        bld.set_decay_resource( dr );
        bld.set_nuclide_resource( nr );

        // create a transition system adapter to create a general transition
        // system from this data
        // see TransitionSystem and friends to_string methods
        Origen::TransitionSystemAdapter_DecayResource tsa;
        ts = std::make_shared<Origen::TransitionSystem_Gen>();
        tsa.extend( &*ts, true, *dr );
    }

    Origen::SP_NuclideResource nr;
    Origen::SP_DecayResource dr;
    Origen::SP_TransitionSystem ts;
    Origen::LibraryBuilder bld;
};

TEST_F( LibraryBuilderTester, CreateNewEnd7Dec )
{
    // list of ids
    std::vector<int> sizzzaaa_list =
        Origen::FakeFactory::createNuclideIds_scale62();

    // build an ORIGEN library with the standard set of SCALE 6.2 ids
    Origen::Library lib;
    bld.create_decay_library( &lib, sizzzaaa_list, true, *ts );

    // check sizes and then write to disk
    ScaleUtils::IO::DB opts;
    opts.set( "fileFormat", "s62b" );
    bool saved = saveLibrary( lib, "end7dec", opts );
    EXPECT_TRUE( saved );
}

TEST_F( LibraryBuilderTester, AddDecayModes )
{
    // show contents
    // look for cs141 modes
    // 1009   551410    cs141 2.790e-02 3.198e+00    0.534   beta-(0) 9.997e-01
    // beta-,neutron(0) 3.500e-04
    // Origen::printDecayTable(*dr,std::cout);

    // get decay parent of interest
    // note ids in table are ZZZAAAI!!!
    // internally ORIGEN now uses IZZZAAA
    Origen::DecayParent& dp = dr->get_parent_map()[55141];
    Insist( dp.num_decay_modes() == 2, "did not find right nuclide" );

    // get decay channel of interest
    // const int Utils::DECAY_NONE                = 0;
    // const int Utils::DECAY_GAMMA               = 1;
    // const int Utils::DECAY_BETA_MINUS          = 2;
    // const int Utils::DECAY_BETA_PLUS           = 3;
    // const int Utils::DECAY_ISOMERIC_TRANSITION = 4;
    // const int Utils::DECAY_ALPHA               = 5;
    // const int Utils::DECAY_NEUTRON             = 6;
    // const int Utils::DECAY_SPONTANEOUS_FISSION = 7;
    // const int Utils::DECAY_PROTON              = 8;
    // const int Utils::DECAY_UNKNOWN             = 9;

    Origen::DecayChannel& dc = dp.get_channel_map()[26];
    // std::cout << "existing 55141[26] branching
    // ratio="<<dc.get_branch_ratio()<<std::endl;
    double dconst = dp.get_decay_constant();

    // change branching ratio for beta-,n
    double new_br = 5e-4;  // would typically read from file
    double delta_br = dc.get_branch_ratio() - new_br;
    dc.set_branch_ratio( new_br );

    // add new decay modes and set branching ratio for beta-,2n
    Origen::DecayChannel& dc_2n = dp.get_channel_map()[266];  // creating new
    double new_br_2n = 6e-4;
    delta_br -= new_br_2n;  // increment delta
    dc_2n.set_branch_ratio( new_br_2n );
    dc_2n.set_yield_type( Origen::DecayResource::ISOMERIC_YIELDS );
    dc_2n.get_yield_map()[0] = 0.95;  // 95% to ground
    dc_2n.get_yield_map()[1] = 0.05;  // 5% to first metastable

    // know what I'm dealing with so I will hard-code the renormalization of
    // other branch ratio
    Origen::DecayChannel& dc_other = dp.get_channel_map()[2];
    double other_br = dc_other.get_branch_ratio() + delta_br;
    dc_other.set_branch_ratio( other_br );

    // Origen::printDecayTable(*dr,std::cout);

    // final result should be:
    // 1009   551410    cs141 2.790e-02 3.198e+00    0.534   beta-(0) 9.989e-01
    // beta-,neutron(0) 5.000e-04  beta-,neutron,neutron(0) 5.700e-04
    // beta-,neutron,neutron(1) 3.000e-05

    // calculate reference values
    double coeff_2 = dconst * dp.get_channel_map()[2].get_branch_ratio();
    double coeff_26 = dconst * dp.get_channel_map()[26].get_branch_ratio();
    double coeff_266 =
        dconst * dp.get_channel_map()[266].get_branch_ratio() * 0.95;

    // create a transition system adapter to create a general transition system
    // from this data
    // see TransitionSystem and friends to_string methods
    Origen::TransitionSystemAdapter_DecayResource tsa;
    ts = std::make_shared<Origen::TransitionSystem_Gen>();
    tsa.extend( &*ts, true, *dr );

    std::vector<int> sizzzaaa_list =
        Origen::FakeFactory::createNuclideIds_scale62();

    // build an ORIGEN library with the standard set of SCALE 6.2 ids
    Origen::Library lib;
    bld.create_decay_library( &lib, sizzzaaa_list, true, *ts );

    // check sizes and then write to disk
    {
        SCOPED_TRACE( "in-memory" );
        gtestLibrarySizes( lib );
        gtestLibraryMasses( lib );
        ScaleUtils::IO::DB opts;
        bool saved = saveLibrary(
            lib,
            "decay.f33",
            opts );  // library can be used in a couple calculation now
        EXPECT_TRUE( saved );
    }

    // read from disk and then check sizes
    Origen::Library lib2;
    {
        ScaleUtils::IO::DB opts;
        Origen::LibraryIO io;
        bool loaded = io.load( lib2, "decay.f33", opts );
        EXPECT_TRUE( loaded );
        SCOPED_TRACE( "on-disk" );
        gtestLibrarySizes( lib2 );
        gtestLibraryMasses( lib2 );
    }

    // now look through library and see that we have some values that we
    // modified
    {
        const auto& nucset = lib2.definition().nuclide_set();
        size_t i_cs141 = nucset.lookup_index( 30055141 );
        ASSERT_GT( nucset.total_nuclides(), i_cs141 );
        ASSERT_EQ( 55141, nucset.izzzaaa_at( i_cs141 ) );
        const auto& tids = lib2.transition_structure().transition_ids();
        const auto& numdp = lib2.transition_structure().num_decay_parents();
        const auto& ppos = lib2.transition_structure().parent_positions();
        const auto& tcoeff = lib2.transition_coeff_at( 0 ).matrix();
        size_t j = 0;
        for( size_t i = 0; i < numdp.size(); ++i )
        {
            for( size_t k = 0; k < numdp[i]; ++k )
            {
                int tid = tids[j];
                int nind = ppos[j] - 1;
                if( nind == (int)i_cs141 )
                {
                    // daughter of cs141: 30056139 by mode=-266 with
                    // coeff=1.59056e-05
                    // daughter of cs141: 30056140 by mode=-26 with
                    // coeff=1.39523e-05
                    // daughter of cs141: 30056141 by mode=-2 with
                    // coeff=0.0278739
                    if( nucset.izzzaaa_at( i ) == 56139 )
                    {
                        EXPECT_EQ( -266, tid );
                        EXPECT_NEAR( coeff_266, tcoeff[j], 1e-4 * coeff_266 );
                    }
                    else if( nucset.izzzaaa_at( i ) == 56140 )
                    {
                        EXPECT_EQ( -26, tid );
                        EXPECT_NEAR( coeff_26, tcoeff[j], 1e-4 * coeff_26 );
                    }
                    else if( nucset.izzzaaa_at( i ) == 56141 )
                    {
                        EXPECT_EQ( -2, tid );
                        EXPECT_NEAR( coeff_2, tcoeff[j], 1e-4 * coeff_2 );
                    }
                    else if(nucset.izzzaaa_at(i) == 1)
                    {
                        // -266 and -26 also produce neutrons now, need to check?
                    }
                    else
                    {
                        EXPECT_TRUE( false )
                            << "daughter of cs141: " << nucset.sizzzaaa_at( i )
                            << " by mode=" << tid << " with coeff=" << tcoeff[j]
                            << " not expected";
                    }
                }
                ++j;
            }
        }
    }
}