tstLibrary.cpp

./Core/dc/tests/tstLibrary.cpp

/*
 * File:   tstLibrary.cpp
 * Author: Brenden T. Mervin, mervinbt@ornl.gov
 * Created on Thu Apr  4 15:32:59 EDT 2013
 */

#include <stdio.h>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <string>
#include <vector>

#include "Origen/Core/config.h"
#include "Origen/Core/dc/FakeFactory.h"
#include "Origen/Core/dc/GridView_TransitionCoeff.h"
#include "Origen/Core/dc/Library.h"
#include "Origen/Core/dc/SerialRegistrar.h"
#include "Origen/Core/dc/TransitionCoeff.h"
#include "Origen/Core/io/LibraryIO.h"
#include "Standard/Interface/SerialFactory.h"
#include "Standard/Interface/SerialStream.h"

#include "ScaleUtils/Math/Interp1D_CubicSpline_SF.h"
#include "ScaleUtils/Math/interpwts.h"

#include "Nemesis/gtest/nemesis_gtest.hh"
#include "Nemesis/harness/DBC.hh"

using namespace Origen;
using namespace Scale;

class LibraryTest : public ::testing::Test
{
  protected:
    virtual void SetUp()
    {
        // nuclide set
        NuclideSet nucSet;
        nucSet.set_ids_legacy( 2, 2, 1, {1, 2, 235, 128, 135} );

        // header
        LibraryHeader header;
        header.set_nuclide_set( nucSet );
        header.set_month_made( 3 );
        header.set_day_made( 29 );
        header.set_year_made( 2013 );
        header.set_title( "TestHeader" );
        lib.set_definition( header );

        // tag manager
        TagManager tagMan;
        tagMan.setIdTag( "Lattice type", "Westinghouse 17x17" );
        tagMan.setIdTag( "Pin location", "Corner" );
        tagMan.setIdTag( "Fuel type", "UO2" );
        tagMan.setInterpTag( "Enrichment", 4.5 );
        tagMan.setInterpTag( "Number of Gadolinium Rods", 8 );
        lib.set_tag_manager( tagMan );

        // transition structure
        TransitionStructure tStruct;
        FakeFactory::TransitionStructure_random1( tStruct, header );
        lib.set_transition_structure( tStruct );

        // decay data
        DecayData decayDat;
        FakeFactory::DecayData_random1( decayDat, header );
        lib.set_decay_data( decayDat );

        // transition coeff
        TransitionCoeff transCoeff;
        FakeFactory::TransitionCoeff_random1( transCoeff, header );
        lib.set_transition_coeff_at( transCoeff, 0 );
    }

    Library lib;
};

/*
 * Test serialization of the Library.
 */
TEST_F( LibraryTest, Initialize )
{
    ASSERT_TRUE( lib.scp_definition() != nullptr );

    //\todo Add some checks on the library definition values here

    // Verify tag manager set correctly
    ASSERT_TRUE( lib.scp_tag_manager() != nullptr );

    // Verify transition_coeff set correctly
    ASSERT_TRUE( lib.scp_transition_coeff_at( 0 ) != nullptr );

    // Verify decay data is set correctly
    ASSERT_TRUE( lib.scp_decay_data() != nullptr );

    // Verify the transitionStructure is set correctly
    ASSERT_TRUE( lib.scp_transition_structure() != nullptr );

    // set up the serialization factory
    SerialStream ss;
    Standard::SerialFactory sf;
    SerialRegistrar( &sf );
    ss.setFactory( &sf );

    // serialize
    int error = lib.serialize( &ss );
    ASSERT_EQ( error, 0 ) << "Library serialization failed" << std::endl;

    // make a new lib to test the serialization
    Library new_lib;

    // deserialize
    error = new_lib.deserialize( &ss );
    ASSERT_EQ( error, 0 );

    // Check that the deserialized Library matches the original
    EXPECT_TRUE( lib == new_lib )
        << "Deserialized library does not match original" << std::endl;

    // Test to/from JSON
    Json::Value oldJSON = lib.to_json();
    Json::Value newJSON = new_lib.to_json();
    EXPECT_TRUE( newJSON == oldJSON );
}

// Test the Library cloning capability
TEST( Library, Clone )
{
    SP_Library fakeLib = SP_Library( new Library() );
    FakeFactory::Library_random1( *fakeLib );
    TagManager tagMan;
    tagMan.setIdTag( "Cheese", "Stinky" );
    tagMan.setInterpTag( "StinkFactor", 4.0 );
    SP_Library cloneLib = fakeLib->clone();

    // Check that everything copied correctly
    ASSERT_TRUE( fakeLib->definition() == cloneLib->definition() );
    ASSERT_TRUE( fakeLib->tag_manager() == cloneLib->tag_manager() );
    ASSERT_TRUE( fakeLib->decay_data() == cloneLib->decay_data() );
    ASSERT_TRUE( fakeLib->transition_structure() ==
                 cloneLib->transition_structure() );

    // Check each individual transition coefficient
    ASSERT_TRUE( fakeLib->transition_coeff_size() ==
                 cloneLib->transition_coeff_size() );
    for( size_t i = 0; i < fakeLib->transition_coeff_size(); ++i )
    {
        TransitionCoeff fakeTC = fakeLib->transition_coeff_at( i );
        SP_TransitionCoeff cloneTC = SP_TransitionCoeff( new TransitionCoeff );
        fakeLib->get_transition_coeff_at( *cloneTC, i );
        EXPECT_EQ( fakeTC.burnup(), cloneTC->burnup() );
        EXPECT_TRUE( fakeTC.approx_eq( *cloneTC ) );
    }

    // Check if clone actually makes a deep copy
    SCP_TransitionCoeff refTC = fakeLib->scp_transition_coeff_at( 0 );
    SP_TransitionCoeff tmpTC = refTC->clone();
    tmpTC->scale( 2.0 );
    SP_TransitionCoeff tmpTC2 = tmpTC->clone();
    EXPECT_FALSE( tmpTC->approx_eq( *refTC ) );
    EXPECT_TRUE( tmpTC2->approx_eq( *tmpTC ) );

    // Set a new TC on the library
    fakeLib->set_transition_coeff_at( *tmpTC, 0 );
    EXPECT_FALSE( *fakeLib == *cloneLib );
    fakeLib->set_transition_coeff_at( *refTC, 0 );
    EXPECT_TRUE( *fakeLib == *cloneLib );

    // Check that the library header has been deep copied
    LibraryHeader refLH = fakeLib->definition();
    LibraryHeader tmpLH = cloneLib->definition();
    EXPECT_NE( fakeLib->scp_definition(), cloneLib->scp_definition() );
    EXPECT_TRUE( *fakeLib == *cloneLib );

    // Check that the TransitionStructure has been deep copied
    TransitionStructure refTS = fakeLib->transition_structure();
    TransitionStructure tmpTS = cloneLib->transition_structure();
    EXPECT_TRUE( refTS == tmpTS );

    size_t itot = fakeLib->definition().nuclide_set().total_nuclides();
    std::vector<int> tmpNumParents( itot, 1 );
    tmpTS.set_num_parents( tmpNumParents );

    EXPECT_FALSE( refTS == tmpTS );
    fakeLib->set_transition_structure( tmpTS );
    EXPECT_FALSE( *fakeLib == *cloneLib );

    fakeLib->set_transition_structure( refTS );
    EXPECT_TRUE( *fakeLib == *cloneLib );

    // Check that the DecayData structure has been deep copied
    DecayData refDD = fakeLib->decay_data();
    DecayData tmpDD = cloneLib->decay_data();
    EXPECT_TRUE( refDD == tmpDD );

    tmpDD.q_fractions_from_photons = std::vector<float>( itot, 0.5 );

    EXPECT_FALSE( refDD == tmpDD );
    fakeLib->set_decay_data( tmpDD );
    EXPECT_FALSE( *fakeLib == *cloneLib );

    fakeLib->set_decay_data( refDD );
    EXPECT_TRUE( *fakeLib == *cloneLib );

    // Check that the TagManager has been deep copied
    TagManager refTM = fakeLib->tag_manager();
    TagManager tmpTM = cloneLib->tag_manager();
    EXPECT_TRUE( refTM == tmpTM );

    tmpTM.setIdTag( "Cheese", "Smelly" );

    EXPECT_FALSE( refTM == tmpTM );
    fakeLib->set_tag_manager( tmpTM );
    EXPECT_FALSE( *fakeLib == *cloneLib );

    fakeLib->set_tag_manager( refTM );
    EXPECT_TRUE( *fakeLib == *cloneLib );
}

// Test out 1-D interpolation over burnups on a Library
TEST( Library, Basic_Interp )
{
    SP_Library fakeLib = SP_Library( new Library );
    FakeFactory::Library_random1( *fakeLib );

    // replace burnups to be realistic and easy to interp
    std::vector<float> burnups = {1000., 2000., 3000.};
    for( size_t i = 0; i < burnups.size(); ++i )
    {
        TransitionCoeff tc = fakeLib->transition_coeff_at( i );
        tc.set_burnup( burnups[i] );
        fakeLib->set_transition_coeff_at( tc, i );
        EXPECT_FLOAT_EQ( burnups[i],
                         fakeLib->transition_coeff_at( i ).burnup() );
    }

    const int numR = 8;  // 1 burnup + 7 TC vectors = 8 interpolable responses

    std::vector<SP_TransitionCoeff> origTC( burnups.size() );
    for( size_t i = 0; i < burnups.size(); ++i )
    {
        origTC[i] = SP_TransitionCoeff( new TransitionCoeff );
        fakeLib->get_transition_coeff_at( *origTC[i], i );
        GridView_TransitionCoeff view( origTC[i] );

        // Check that the view matches itself
        EXPECT_TRUE( view.match( *origTC[i] ) );

        // Check that we see 8 responses (1 burnup + 7 TC vectors)
        ASSERT_EQ( numR, view.nr() );
    }
    std::vector<float> burnupNew = {1000., 1500., 2000., 2500., 3000.};

    // Set a function pointer to the interpolator
    Library::InterpFunction intF =
        &( ScaleUtils::Math::interpwts_lin<std::vector<float>, float> );
    fakeLib->set_interpolator( intF );
    SP_Library libNew = fakeLib->interpolate( burnupNew );

    // New library has the correct number of burnups
    ASSERT_EQ( libNew->transition_coeff_size(), burnupNew.size() );

    // New library burnups are equal to burnups in burnupNew

    for( size_t i = 0; i < burnupNew.size(); i++ )
    {
        SP_TransitionCoeff tcoeff = SP_TransitionCoeff( new TransitionCoeff );
        libNew->get_transition_coeff_at( *tcoeff, i );

        // New Library has correct burnups
        EXPECT_EQ( tcoeff->burnup(), burnupNew[i] );

        std::vector<float> weights( 2 );
        SP_TransitionCoeff lesser, greater;

        // Vectors in Transition Coefficients interpolated correctly
        if( burnupNew[i] == burnups[int( i / 2 )] )
        {
            // interpolated burnups that match an original burnup (1000, 2000,
            // 3000)
            lesser = greater = origTC[i / 2];
            weights = {1.0, 0.0};
        }
        else
        {
            // so can Assert that positions 2,4 ==
            // (0.5*tc1+0.5*tc2),(0.5*tc2+0.5*tc3)
            // which is equivalent to 2,4  ==
            // (0.5*pos1+0.5*pos3),(0.5*pos3+0.5*pos5)
            lesser = origTC[( i - 1 ) / 2];
            greater = origTC[( i + 1 ) / 2];
            weights = {0.5, 0.5};
        }
        std::vector<SP_TransitionCoeff> interpCoeffs = {lesser, greater};

        // Create a fresh TC for comparison
        SP_TransitionCoeff tcComp = tcoeff->clone();
        tcComp->scale( 0.0 );

        // Add the scaled contributions from each TC on the interpolation grid
        for( size_t i = 0; i != weights.size(); ++i )
        {
            SP_TransitionCoeff tcClone = interpCoeffs[i]->clone();
            tcClone->scale( weights[i] );
            *tcComp += *tcClone;
            tcClone.reset();
        }

        // Check that the interpolated TC at this burnup matches what we expect
        EXPECT_FLOAT_EQ( tcComp->burnup(), tcoeff->burnup() );
        EXPECT_TRUE( tcComp->approx_eq( *tcoeff ) );
    }  // end loop over burnups
}

// Test out 1-D interpolation over burnups on a Library using SplineFit
TEST( Library, Basic_Interp1D )
{
    using ScaleUtils::Math::Numeric::Interp1D_CubicSpline_SF;
    SP_Library fakeLib = SP_Library( new Library );
    FakeFactory::Library_random1( *fakeLib );

    // replace burnups to be realistic and easy to interp
    std::vector<float> burnups = {1000., 2000., 3000., 4000., 5000., 6000.};
    for( size_t i = 0; i < burnups.size(); ++i )
    {
        TransitionCoeff tc;
        FakeFactory::TransitionCoeff_random1( tc, fakeLib->definition() );
        tc.set_burnup( burnups[i] );
        fakeLib->set_transition_coeff_at( tc, i );
        EXPECT_FLOAT_EQ( burnups[i],
                         fakeLib->transition_coeff_at( i ).burnup() );
    }

    std::vector<SP_TransitionCoeff> origTC( burnups.size() );
    for( size_t i = 0; i < burnups.size(); ++i )
    {
        origTC[i] = SP_TransitionCoeff( new TransitionCoeff );
        fakeLib->get_transition_coeff_at( *origTC[i], i );
    }

    std::vector<double> burnupNew = {1000., 1500., 2000., 2500., 3000.};
    SP_Library libNew = fakeLib->interpolate_Interp1D( burnupNew );

    // New library has the correct number of burnups
    ASSERT_EQ( libNew->transition_coeff_size(), burnupNew.size() );

    // New library burnups are equal to burnups in burnupNew

    for( size_t i = 0; i < burnupNew.size(); i++ )
    {
        SP_TransitionCoeff tcoeff = SP_TransitionCoeff( new TransitionCoeff );
        libNew->get_transition_coeff_at( *tcoeff, i );

        // New Library has correct burnups
        EXPECT_EQ( tcoeff->burnup(), burnupNew[i] );

        std::vector<float> weights( 2 );
        SP_TransitionCoeff lesser, greater;

        // Vectors in Transition Coefficients interpolated correctly
        if( burnupNew[i] == burnups[int( i / 2 )] )
        {
            // interpolated burnups that match an original burnup (1000, 2000,
            // 3000)
            lesser = greater = origTC[i / 2];
            weights = {1.0, 0.0};
        }
        else
        {
            // so can Assert that positions 2,4 ==
            // (0.5*tc1+0.5*tc2),(0.5*tc2+0.5*tc3)
            // which is equivalent to 2,4  ==
            // (0.5*pos1+0.5*pos3),(0.5*pos3+0.5*pos5)
            lesser = origTC[( i - 1 ) / 2];
            greater = origTC[( i + 1 ) / 2];
            weights = {0.5, 0.5};
        }
        std::vector<SP_TransitionCoeff> interpCoeffs = {lesser, greater};

        // Create a fresh TC for comparison
        SP_TransitionCoeff tcComp = tcoeff->clone();
        tcComp->scale( 0.0 );

        // Add the scaled contributions from each TC on the interpolation grid
        for( size_t i = 0; i != weights.size(); ++i )
        {
            SP_TransitionCoeff tcClone = interpCoeffs[i]->clone();
            tcClone->scale( weights[i] );
            *tcComp += *tcClone;
            tcClone.reset();
        }

        // Check that the interpolated TC at this burnup matches what we expect
        EXPECT_FLOAT_EQ( tcComp->burnup(), tcoeff->burnup() );
        EXPECT_TRUE( tcComp->approx_eq( *tcoeff ) );
    }  // end loop over burnups
}

TEST( Library, Interp1D_Thinner )
{
    // "Re-grid" a Library to a desired set of burnups using interpolation

    Library baseLib;
    Library refLib;
    Library scaledLib;
    TransitionCoeff baseTC;

    Vec_Flt base_bu_list = {0.0, 1.0, 10.0, 20.0};
    Vec_Flt ref_bu_list = {0.0, 5.0, 15.0, 20.0};

    // Setup the base & reference libraries from the same base TransitionCoeff
    FakeFactory::TransitionCoeff_random2( baseTC );
    baseTC.set_burnup( base_bu_list[0] );
    baseLib.add_transition_coeff( baseTC );
    refLib.add_transition_coeff( baseTC );

    // Create a linearly-scaled series of TCs on the Library, based on the
    // burnup grid
    FakeFactory::Library_linear1D( baseLib, base_bu_list, 1.5 );
    FakeFactory::Library_linear1D( refLib, ref_bu_list, 1.5 );

    // Set a function pointer to the interpolator
    Library::InterpFunction intF =
        &( ScaleUtils::Math::interpwts_lin<std::vector<float>, float> );
    baseLib.set_interpolator( intF );
    baseLib.regrid( ref_bu_list );
    ASSERT_EQ( refLib.transition_coeff_size(),
               baseLib.transition_coeff_size() );

    for( size_t i = 0; i != refLib.transition_coeff_size(); ++i )
    {
        SCP_TransitionCoeff tcNew = baseLib.scp_transition_coeff_at( i );
        SCP_TransitionCoeff tcRef = refLib.scp_transition_coeff_at( i );

        EXPECT_FLOAT_EQ( tcRef->burnup(), tcNew->burnup() )
            << "Burnups mismatched for Library " << i << " / "
            << refLib.transition_coeff_size() << std::endl;
        EXPECT_TRUE( tcRef->approx_eq( *tcNew ) )
            << "Library " << i << " / " << refLib.transition_coeff_size()
            << " is does not match expected Library" << std::endl;
    }
}

TEST( Library, FortranTestConsistency )
{
    Library lib;
    FakeFactory::Library_random1( lib, 333 );
    EXPECT_FLOAT_EQ( 1.50642300E-03, lib.decay_data().decay_constants.at( 0 ) );
    EXPECT_FLOAT_EQ( 8.69777426E-03,
                     lib.transition_coeff_at( 0 ).matrix().at( 0 ) );
}

// Testing the fakeFactory function to a create a 1D linearly scaled library of
// TCs
TEST( FakeFactory, linear1D )
{
    SP_Library baseLib = SP_Library( new Library );
    SP_TransitionCoeff baseTC = SP_TransitionCoeff( new TransitionCoeff );
    FakeFactory ff;
    ff.TransitionCoeff_random2( *baseTC );
    baseTC->set_burnup( 1000 );
    baseLib->add_transition_coeff( *baseTC );

    std::vector<float> bu_list = {1000, 2000, 3000};
    float scaleBU = 3;
    ff.Library_linear1D( *baseLib, bu_list, scaleBU );

    ASSERT_TRUE( baseLib->transition_coeff_size() == 3 );
    EXPECT_FLOAT_EQ( baseLib->transition_coeff_at( 0 ).burnup(), 1000 );
    EXPECT_FLOAT_EQ( baseLib->transition_coeff_at( 1 ).burnup(), 2000 );
    EXPECT_FLOAT_EQ( baseLib->transition_coeff_at( 2 ).burnup(), 3000 );

    std::vector<float> weights = {1.0, 2.0, 3.0};

    SP_TransitionCoeff tmpTC = SP_TransitionCoeff( new TransitionCoeff );
    baseLib->get_transition_coeff_at( *tmpTC, 0 );

    SP_TransitionCoeff scaledTC = tmpTC->clone();
    scaledTC->scale( weights[0] );
    EXPECT_FLOAT_EQ( scaledTC->burnup(), bu_list[0] );
    EXPECT_TRUE( baseLib->transition_coeff_at( 0 ).approx_eq( *scaledTC ) );

    scaledTC = tmpTC->clone();
    scaledTC->scale( weights[1] );
    EXPECT_FLOAT_EQ( scaledTC->burnup(), bu_list[1] );
    EXPECT_TRUE( baseLib->transition_coeff_at( 1 ).approx_eq( *scaledTC ) );

    scaledTC = tmpTC->clone();
    scaledTC->scale( weights[2] );
    EXPECT_FLOAT_EQ( scaledTC->burnup(), bu_list[2] );
    EXPECT_TRUE( baseLib->transition_coeff_at( 2 ).approx_eq( *scaledTC ) );
}

// Testing the fakefactory function to create linear libraries
TEST( FakeFactory, linear2D )
{
    SP_Library baseLib = SP_Library( new Library );
    SP_TransitionCoeff baseTC = SP_TransitionCoeff( new TransitionCoeff );
    FakeFactory ff;
    ff.TransitionCoeff_random2( *baseTC );
    baseTC->set_burnup( 1000 );
    baseLib->add_transition_coeff( *baseTC );

    // make params for fake library to be made at
    std::vector<float> bu_list = {1000, 2000, 3000};
    std::vector<float> enr_list = {0.10, 0.20, 0.30};
    float scaleBU = 3;
    float scaleEnr = 3;
    std::string name = "Enrichment";

    std::vector<SP_Library> libVec = ff.Library_linear2D(
        *baseLib, bu_list, enr_list, scaleBU, scaleEnr, name );
    // check to see if libraries in vector are scaled correctly
    ASSERT_TRUE( libVec.size() == 3 );  // correct number of libraries

    ASSERT_TRUE( libVec[0]->transition_coeff_size() ==
                 3 );  // correct number of transCoeff
    ASSERT_TRUE( libVec[1]->transition_coeff_size() == 3 );
    ASSERT_TRUE( libVec[2]->transition_coeff_size() == 3 );

    for( size_t i = 0; i < 3; ++i )
    {
        for( size_t j = 0; j < 3; ++j )
        {
            // correct burnups
            EXPECT_FLOAT_EQ( libVec[i]->transition_coeff_at( j ).burnup(),
                             bu_list[j] );
        }
        // correct interptags
        EXPECT_TRUE( libVec[i]->tag_manager().hasTag( name ) );
        ASSERT_DOUBLE_EQ( double( enr_list[i] ),
                          libVec[i]->tag_manager().getInterpTag( name ) );
    }

    // correct scaling on TransitionCoeffs
    std::vector<float> wgts = {1.0, 2.0, 3.0, 9.0};

    SP_TransitionCoeff tempTC;

    tempTC = baseTC->clone();
    tempTC->scale( wgts[0] );
    EXPECT_TRUE( libVec[0]->transition_coeff_at( 0 ).approx_eq( *tempTC ) );

    tempTC = baseTC->clone();
    tempTC->scale( wgts[1] );
    tempTC->set_burnup( bu_list[1] );
    EXPECT_FLOAT_EQ( tempTC->burnup(),
                     libVec[0]->transition_coeff_at( 1 ).burnup() );
    EXPECT_TRUE( libVec[0]->transition_coeff_at( 1 ).approx_eq( *tempTC ) );

    tempTC = baseTC->clone();
    tempTC->scale( wgts[2] );
    tempTC->set_burnup( bu_list[2] );
    EXPECT_FLOAT_EQ( tempTC->burnup(),
                     libVec[0]->transition_coeff_at( 2 ).burnup() );
    EXPECT_TRUE( libVec[0]->transition_coeff_at( 2 ).approx_eq( *tempTC ) );

    tempTC = baseTC->clone();
    tempTC->scale( wgts[1] );
    tempTC->set_burnup( bu_list[0] );
    EXPECT_FLOAT_EQ( tempTC->burnup(),
                     libVec[1]->transition_coeff_at( 0 ).burnup() );
    EXPECT_TRUE( libVec[1]->transition_coeff_at( 0 ).approx_eq( *tempTC ) );

    tempTC = baseTC->clone();
    tempTC->scale( wgts[2] );
    tempTC->set_burnup( bu_list[0] );
    EXPECT_FLOAT_EQ( tempTC->burnup(),
                     libVec[2]->transition_coeff_at( 0 ).burnup() );
    EXPECT_TRUE( libVec[2]->transition_coeff_at( 0 ).approx_eq( *tempTC ) );

    tempTC = baseTC->clone();
    tempTC->scale( wgts[3] );
    tempTC->set_burnup( bu_list[2] );
    EXPECT_FLOAT_EQ( tempTC->burnup(),
                     libVec[2]->transition_coeff_at( 2 ).burnup() );
    EXPECT_TRUE( libVec[2]->transition_coeff_at( 2 ).approx_eq( *tempTC ) );
}