tstGridView.cpp

./Core/dc/tests/tstGridView.cpp

#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <iostream>

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

#include "Origen/Core/dc/FakeFactory.h"
#include "Origen/Core/dc/GridData_TransitionCoeff.h"
#include "Origen/Core/dc/GridView_State.h"
#include "Origen/Core/dc/GridView_StateInfo.h"
#include "Origen/Core/dc/GridView_TransitionCoeff.h"
#include "Origen/Core/dc/Library.h"

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

#include "Origen/Core/config.h"

using namespace Origen;

TEST( GridView_StateInfo, Basic )
{
    StateSet f;

    FakeFactory::StateSet_random1( f );

    // get two state info
    SCP_StateInfo info1 = f.states_at( 0 ).scp_definition();
    SCP_StateInfo info2 = f.states_at( 1 ).scp_definition();

    float t1 = info1->time();
    float t2 = info2->time();
    float bu1 = info1->energy();
    float bu2 = info2->energy();

    // make views
    GridView_StateInfo view1( info1 );
    GridView_StateInfo view2( info2 );

    // matches itself
    EXPECT_TRUE( view1.match( *info1 ) );

    // matches second one
    EXPECT_TRUE( view1.match( *info2 ) );

    // check view of object 1
    {
        GridView_StateInfo::GD_t::I_t nr = view1.nr();
        ASSERT_EQ( 5, nr );

        const GridView_StateInfo::GD_t::Vec_I& rextents = view1.rextents();
        ASSERT_EQ( 1, rextents[0] );
        ASSERT_EQ( 1, rextents[1] );

        const GridView_StateInfo::GD_t::Vec_Rptr& r = view1.r();
        ASSERT_EQ( t1, r[0][0] );
        ASSERT_EQ( bu1, r[1][0] );

        // const StateInfo& obj = view1.obj();
        // EXPECT_TRUE(obj==(*info1));
    }

    // check view of object 2
    {
        GridView_StateInfo::GD_t::I_t nr = view2.nr();
        ASSERT_EQ( 5, nr );

        const GridView_StateInfo::GD_t::Vec_I& rextents = view2.rextents();
        ASSERT_EQ( 1, rextents[0] );
        ASSERT_EQ( 1, rextents[1] );

        const GridView_StateInfo::GD_t::Vec_Rptr& r = view2.r();
        ASSERT_EQ( t2, r[0][0] );
        ASSERT_EQ( bu2, r[1][0] );

        // const StateInfo& obj = view2.obj();
        // EXPECT_TRUE(obj==(*info1));
    }

    // build a phony StateSet just from headers
    SP_State state1( new State() );
    state1->set_definition( *info1 );

    SP_State state2( new State() );
    state2->set_definition( *info2 );

    StateSet state_set;

    state_set.add_states( *state1 );
    state_set.add_states( *state2 );

    ASSERT_EQ( info1->case_number(), info2->case_number() );

    SCP_StateInfo info =
        state_set.interpolate_info( 0.5 * ( t1 + t2 ), info1->case_number() );

    ASSERT_FALSE( info == nullptr );

    EXPECT_FLOAT_EQ( 0.5 * ( bu1 + bu2 ), info->energy() );
    EXPECT_FLOAT_EQ( 0.5 * ( t1 + t2 ), info->time() );
}

TEST( GridView_StateInfo, WithTagManager )
{
    StateSet f;

    FakeFactory::StateSet_random1( f );

    // add tag managers to state set
    {
        TagManager tagman1;
        TagManager tagman2;

        tagman1.setInterpTag( "boron", 600 );
        tagman2.setInterpTag( "boron", 1200 );
        State state1 = f.states_at( 0 );
        State state2 = f.states_at( 1 );
        StateInfo info1 = state1.definition();
        StateInfo info2 = state2.definition();

        info1.set_tag_manager( tagman1 );
        info2.set_tag_manager( tagman2 );
        state1.set_definition( info1 );
        state2.set_definition( info2 );
        f.set_states_at( state1, 0 );
        f.set_states_at( state2, 1 );
    }

    SCP_StateInfo info1 = f.states_at( 0 ).scp_definition();
    SCP_StateInfo info2 = f.states_at( 1 ).scp_definition();

    // make views
    GridView_StateInfo view1( info1 );
    GridView_StateInfo view2( info2 );

    // matches itself
    EXPECT_TRUE( view1.match( *info1 ) );

    // matches second one
    EXPECT_TRUE( view1.match( *info2 ) ) << info1->toString() << "\n"
                                         << info2->toString();

    // 6 responses with addition of tag manager
    ASSERT_EQ( 6, view1.nr() );
    ASSERT_EQ( 6, view2.nr() );

    // build a phony StateSet just from headers
    State state1;
    state1.set_definition( *info1 );

    State state2;
    state2.set_definition( *info2 );

    StateSet state_set;

    state_set.add_states( state1 );
    state_set.add_states( state2 );

    ASSERT_EQ( info1->case_number(), info2->case_number() );

    float t1 = info1->time();
    float t2 = info2->time();
    SP_StateInfo info =
        state_set.interpolate_info( 0.5f * ( t1 + t2 ), info1->case_number() );

    ASSERT_FALSE( info == nullptr );

    EXPECT_FLOAT_EQ( 900, info->tag_manager().getInterpTag( "boron" ) );
}

TEST( GridView_State, Basic )
{
    StateSet f;

    FakeFactory::StateSet_random1( f );

    // get two state
    SCP_State state1 = f.scp_states_at( 0 );
    SCP_State state2 = f.scp_states_at( 1 );

    // get two state info
    SCP_StateInfo info1 = state1->scp_definition();
    SCP_StateInfo info2 = state2->scp_definition();

    // make views
    GridView_State view1( state1 );
    GridView_State view2( state2 );

    // matches itself
    EXPECT_TRUE( view1.match( *state1 ) );

    // matches second one
    EXPECT_TRUE( view1.match( *state2 ) );

    // 6 responses with addition of tag manager
    ASSERT_EQ( 1, view1.nr() );
    ASSERT_EQ( 1, view2.nr() );

    // build a new state set
    StateSet state_set;

    state_set.add_states( *state1 );
    state_set.add_states( *state2 );

    ASSERT_EQ( info1->case_number(), info2->case_number() );

    float t1 = info1->time();
    float t2 = info2->time();
    SP_State state = state_set.interpolate(
        0.5 * ( t1 + t2 ), state1->definition().case_number() );

    ASSERT_FALSE( state == nullptr );

    // check interpolation
    double c1 = state1->concs().vals_at( 0 );
    double c2 = state2->concs().vals_at( 0 );
    double c = state->concs().vals_at( 0 );

    EXPECT_FLOAT_EQ( 0.5 * ( c1 + c2 ), c );
}

// Test out basic 2-D interpolation using TagManager with
// GridView_TransitionCoeff
TEST( GridView_TransitionCoeff, DISABLED_ND_Grid_Interp )
{
    using ScaleUtils::Math::interpwts_lin;

    // Create a 2-D grid of fake libraries associated with fake enrichment /
    // void fractions
    std::vector<SP_Library> fakeLibs;
    GridData_TransitionCoeff InterpGrid;
    std::vector<size_t> libOrder;
    std::vector<TagManager> fakeTagMans;
    const int numR = 8;  // 1 burnup + 7 TC vectors = 8 interpolable responses

    std::vector<float> voidFracs = {0.0, 40.0, 80.0};
    std::vector<float> enrich = {0.5, 1.0, 1.5, 3.5, 5.0};
    std::vector<std::string> dimNames = {"void", "enrich"};

    for( size_t i = 0; i != voidFracs.size(); ++i )
    {
        for( size_t j = 0; j != enrich.size(); ++j )
        {
            // Create a set of arbitrary CX's & tag our fake library
            SP_Library fakeLib = SP_Library( new Library() );
            FakeFactory::Library_random1( *fakeLib );
            SP_TagManager fakeTags = SP_TagManager( new TagManager() );
            fakeTags->setInterpTag( dimNames[0], voidFracs[i] );
            fakeTags->setInterpTag( dimNames[1], enrich[j] );
            fakeLib->set_tag_manager( *fakeTags );

            // Give the library somewhat more realistic burnups
            for( size_t k = 0; k != fakeLib->transition_coeff_size(); ++k )
            {
                ASSERT_TRUE( fakeLib->has_transition_coeff_at( k ) );
                TransitionCoeff tc = fakeLib->transition_coeff_at( k );
                tc.set_burnup( 5000.0 * k );
                fakeLib->set_transition_coeff_at( tc, k );
            }
            fakeLibs.push_back( fakeLib );
        }
    }
    // Now extract all of our tagManager objects from the Libraries,
    // similar to how we would for a real N-D interpolation
    for( auto& lib : fakeLibs )
    {
        fakeTagMans.push_back( lib->tag_manager() );
    }

    // Each TransitionCoeff has 8 interpolable responses => nr=8
    InterpGrid.setx( fakeTagMans, dimNames, numR );
    libOrder = InterpGrid.get_tm_indexMap();

    // check that xgrid is the same as what we expect it to be
    std::vector<float> voidGrid;
    std::vector<float> enrichGrid;
    voidGrid.assign( InterpGrid.x( 0 ),
                     InterpGrid.x( 0 ) + InterpGrid.nx( 0 ) );
    enrichGrid.assign( InterpGrid.x( 1 ),
                       InterpGrid.x( 1 ) + InterpGrid.nx( 1 ) );

    ASSERT_EQ( InterpGrid.nx( 0 ), 3 );  // number of void xgrid points
    ASSERT_EQ( InterpGrid.nx( 1 ), 5 );  // number of enrichment xgrid points

    EXPECT_TRUE( voidGrid == voidFracs );  // void xgrid = void points
    EXPECT_TRUE( enrichGrid ==
                 enrich );  // enrichment xgrid = enrichment points

    ASSERT_EQ(
        InterpGrid.nr(),
        numR );  // should see 1 burnup + 7 transition coefficient responses

    // Now re-order the Libraries based on libOrder
    std::vector<SP_Library> gridOrderedLibs;
    for( auto& idx : libOrder )
    {
        Check( idx < fakeLibs.size() );
        gridOrderedLibs.push_back( fakeLibs[idx] );
    }

    // Check that the library order is correct
    EXPECT_FLOAT_EQ( gridOrderedLibs[0]->tag_manager().getInterpTag( "enrich" ),
                     0.5 );
    EXPECT_FLOAT_EQ( gridOrderedLibs[0]->tag_manager().getInterpTag( "void" ),
                     0.0 );
    EXPECT_FLOAT_EQ( gridOrderedLibs[3]->tag_manager().getInterpTag( "enrich" ),
                     3.5 );
    EXPECT_FLOAT_EQ( gridOrderedLibs[3]->tag_manager().getInterpTag( "void" ),
                     0.0 );
    EXPECT_FLOAT_EQ( gridOrderedLibs[7]->tag_manager().getInterpTag( "enrich" ),
                     1.5 );
    EXPECT_FLOAT_EQ( gridOrderedLibs[7]->tag_manager().getInterpTag( "void" ),
                     40.0 );
    EXPECT_FLOAT_EQ(
        gridOrderedLibs[11]->tag_manager().getInterpTag( "enrich" ), 1.0 );
    EXPECT_FLOAT_EQ( gridOrderedLibs[11]->tag_manager().getInterpTag( "void" ),
                     80.0 );
    EXPECT_FLOAT_EQ(
        gridOrderedLibs[14]->tag_manager().getInterpTag( "enrich" ), 5.0 );
    EXPECT_FLOAT_EQ( gridOrderedLibs[14]->tag_manager().getInterpTag( "void" ),
                     80.0 );

    // Now try a 2-D interpolation
    // Note that we don't even need responses here, yet!
    const float enrInterp = 2.5;
    const float voidInterp = 20.0;
    GridData_TransitionCoeff::Vec_X wt_enr;
    GridData_TransitionCoeff::Vec_X wt_void;

    // Calculate the grid interpolation weights
    interpwts_lin( enrich, enrInterp, &wt_enr );
    interpwts_lin( voidFracs, voidInterp, &wt_void );
    GridData_TransitionCoeff::Vec_Xptr wts = {&wt_void[0], &wt_enr[0]};

    ASSERT_EQ( wt_enr.size(), enrich.size() );
    ASSERT_EQ( wt_void.size(), voidFracs.size() );

    // Now add the responses to the grid via GridView objects
    // Do this for the first burnup value on the Library
    for( size_t bu = 0; bu < gridOrderedLibs[0]->transition_coeff_size(); ++bu )
    {
        std::vector<GridView_TransitionCoeff> gridResponses;
        for( auto& lib : gridOrderedLibs )
        {
            GridView_TransitionCoeff gvTC( lib->scp_transition_coeff_at( bu ) );
            gridResponses.push_back( gvTC );
        }
        // Set the response vector for this Grid from the GridView
        InterpGrid.setr( gridResponses );

        // Finally, interpolate on the calculated weights and return a new
        // TransitionCoeff via
        // GridData_TransitionCoeff::interpolate_transition_coeff
        GridData_TransitionCoeff::Vec_Z z = InterpGrid.rgridsum( wts );
        SP_TransitionCoeff interpTC =
            InterpGrid.interpolate_transition_coeff( wts );

        // Check that the burnup of the interpolated TC is the same as the
        // original
        // (i.e., as we have not interpolated on burnup)
        EXPECT_FLOAT_EQ( gridOrderedLibs[0]->transition_coeff_at( bu ).burnup(),
                         interpTC->burnup() );

        // Grab up the TCs we're interpolating on
        std::vector<SP_TransitionCoeff> gridTCs;
        gridTCs.push_back( gridOrderedLibs[2]
                               ->transition_coeff_at( bu )
                               .clone() );  // ENR 1.5 / VOID 0
        gridTCs.push_back( gridOrderedLibs[3]
                               ->transition_coeff_at( bu )
                               .clone() );  // ENR 3.5 / VOID 0
        gridTCs.push_back( gridOrderedLibs[7]
                               ->transition_coeff_at( bu )
                               .clone() );  // ENR 1.5 / VOID 40
        gridTCs.push_back( gridOrderedLibs[8]
                               ->transition_coeff_at( bu )
                               .clone() );  // ENR 3.5 / VOID 40

        const std::vector<float> tcWeights = {0.25, 0.25, 0.25, 0.25};

        // Create a fresh TC for comparison
        SP_TransitionCoeff tcComp = gridTCs[0]->clone();
        tcComp->scale( 0.0 );

        // Add the scaled contributions from each TC on the interpolation grid
        ASSERT_EQ( gridTCs.size(), tcWeights.size() );
        for( size_t i = 0; i != gridTCs.size(); ++i )
        {
            SP_TransitionCoeff tcClone = gridTCs[i]->clone();
            tcClone->scale( tcWeights[i] );
            *tcComp += *tcClone;
            tcClone.reset();
        }
        EXPECT_FLOAT_EQ( tcComp->burnup(), interpTC->burnup() );
        EXPECT_TRUE( tcComp->approx_eq( *interpTC ) );

    }  // end loop over burnups
}