tstStateSet.cpp

./Core/dc/tests/tstStateSet.cpp

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

#include "Origen/Core/dc/FakeFactory.h"
#include "Origen/Core/dc/StateSet.h"

#include "Origen/Core/Definitions.h"

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

#include "Origen/Core/fn/print.h"

#include "Origen/Core/config.h"

using namespace Origen;

TEST( StateSet, Interpolation )
{
    StateSet f;
    FakeFactory::StateSet_random1( f );

    bool print_iso = false;
    if( print_iso )
    {
        std::ofstream oldfile( "old.csv" );
        printConcentrations_csv( f, GATOMS, nullptr, "", oldfile );
    }

    size_t n = f.states_size();
    float dt = f.states_at( n - 1 ).definition().time() -
               f.states_at( 0 ).definition().time();

    int case_number = f.states_at( 0 ).definition().case_number();

    StateSet newf;
    size_t newn = 100;
    for( size_t i = 0; i < newn; ++i )
    {
        SP_State state = f.interpolate( i * dt / ( newn - 1 ), case_number );
        newf.add_states( *state );
    }

    if( print_iso )
    {
        std::ofstream newfile( "new.csv" );
        printConcentrations_csv( newf, GATOMS, nullptr, "", newfile );
    }

    Vec_Flt test;
    {
        Vec_Dbl c;
        newf.states_at( 1 ).concs().get_vals( c );
        test.assign( c.begin(), c.end() );
    }

    Vec_Dbl c1, c2;
    f.states_at( 0 ).concs().get_vals( c1 );
    f.states_at( 1 ).concs().get_vals( c2 );
    ASSERT_EQ( c1.size(), c2.size() );

    // linearly interpolate reference solution
    Vec_Flt ref( c1.size() );
    float tref = newf.states_at( 1 ).definition().time();
    float t1 = f.states_at( 0 ).definition().time();
    float t2 = f.states_at( 1 ).definition().time();
    for( size_t i = 0; i < c1.size(); ++i )
    {
        ref[i] = c1[i] + ( c2[i] - c1[i] ) * ( tref - t1 ) / ( t2 - t1 );
    }

    EXPECT_VEC_SOFTEQ( ref, test, 1e-5 );
}

TEST( StateSet, Basic )
{
    StateSet set;

    // nuclide set
    NuclideSet nuc_set;
    nuc_set.set_ids( {8016, 92235} );

    // concentrations with default units
    Concentrations concs;
    concs.set_nuclide_set( nuc_set );
    concs.set_vals( {1.5, 10.0} );

    // create first state (copies concs' internals)
    State s1;
    s1.set_concs( concs );

    // create second state with different concentrations
    // reuse first concentrations to avoid duplicating NuclideSet
    State s2;
    concs.set_vals( {16., 5.} );
    s2.set_concs( concs );

    // add states to set
    set.add_states( s1 );
    set.add_states( s2 );

    EXPECT_FLOAT_EQ( 8016,
                     set.states_at( 0 ).concs().nuclide_set().ids_at( 0 ) );
    EXPECT_FLOAT_EQ( 92235,
                     set.states_at( 0 ).concs().nuclide_set().ids_at( 1 ) );

    EXPECT_FLOAT_EQ( 1.5, set.states_at( 0 ).concs().vals_at( 0 ) );
    EXPECT_FLOAT_EQ( 10.0, set.states_at( 0 ).concs().vals_at( 1 ) );
    EXPECT_FLOAT_EQ( 16.0, set.states_at( 1 ).concs().vals_at( 0 ) );
    EXPECT_FLOAT_EQ( 5.0, set.states_at( 1 ).concs().vals_at( 1 ) );
}

TEST( StateSet, Convenience )
{
    StateSet f;

    FakeFactory::StateSet_random1( f );

    Vec_Dbl energies = f.energies();
    Vec_Dbl times = f.times();
    Vec_Dbl powers = f.powers();
    Vec_Dbl fluxes = f.fluxes();

    EXPECT_EQ( 11, f.neutron_spectra().size() );
    EXPECT_EQ( 11, f.alpha_spectra().size() );
    EXPECT_EQ( 11, f.beta_spectra().size() );
    EXPECT_EQ( 11, f.gamma_spectra().size() );
    EXPECT_EQ( 44, f.concs().size() );

    EXPECT_EQ( 44, f.energies().size() );
    EXPECT_EQ( 44, f.times().size() );
    EXPECT_EQ( 44, f.powers().size() );
    EXPECT_EQ( 44, f.fluxes().size() );

    for( size_t i = 0; i < f.states_size(); ++i )
    {
        EXPECT_EQ( energies[i], f.states_at( i ).definition().energy() );
        EXPECT_EQ( times[i], f.states_at( i ).definition().time() );
        EXPECT_EQ( powers[i], f.states_at( i ).definition().power() );
        EXPECT_EQ( fluxes[i], f.states_at( i ).definition().flux() );
    }
}

void timing_smart_pointer( int repeat, Origen::StateSet& state_set )
{
    SCP_NuclideSet nuclide_set;
    SCP_Concentrations concs;
    SCP_State state;
    SCP_Vec_Dbl vals;
    SCP_Vec_Int ids;
    int units;
    int isum;
    double dsum;
    for( int r = 0; r < repeat; ++r )
    {
        isum = 0;
        dsum = 0.0;
        for( size_t i = 0; i < state_set.states_size(); ++i )
        {
            state = state_set.scp_states_at( i );
            concs = state->scp_concs();
            units = concs->units();
            vals = concs->scp_vals();
            nuclide_set = concs->scp_nuclide_set();
            ids = nuclide_set->scp_ids();
            for( size_t j = 0; j < vals->size(); ++j )
            {
                isum += ids->at( j );
                dsum += vals->at( j );
            }
        }
        // std::cout << "isum="<<isum<<"\n";
        // std::cout << "dsum="<<dsum<<"\n";
    }
    (void)units;
    (void)isum;
    (void)dsum;
}

void timing_ref_at( int repeat, Origen::StateSet& state_set )
{
    int units;
    int isum;
    double dsum;
    for( int r = 0; r < repeat; ++r )
    {
        isum = 0;
        dsum = 0.0;
        for( size_t i = 0; i < state_set.states_size(); ++i )
        {
            const State& state = state_set.states_at( i );
            const Concentrations& concs = state.concs();
            units = concs.units();
            const Vec_Dbl& vals = concs.vals();
            const NuclideSet& nuclide_set = concs.nuclide_set();
            const Vec_Int& ids = nuclide_set.ids();
            for( size_t j = 0; j < vals.size(); ++j )
            {
                isum += ids.at( j );
                dsum += vals.at( j );
            }
        }
        // std::cout << "isum="<<isum<<"\n";
        // std::cout << "dsum="<<dsum<<"\n";
    }
    (void)units;
    (void)isum;
    (void)dsum;
}

void timing_ref( int repeat, Origen::StateSet& state_set )
{
    int units;
    int isum;
    double dsum;
    for( int r = 0; r < repeat; ++r )
    {
        isum = 0;
        dsum = 0.0;
        for( size_t i = 0; i < state_set.states_size(); ++i )
        {
            const State& state = state_set.states_at( i );
            const Concentrations& concs = state.concs();
            units = concs.units();
            const Vec_Dbl& vals = concs.vals();
            const NuclideSet& nuclide_set = concs.nuclide_set();
            const Vec_Int& ids = nuclide_set.ids();
            for( size_t j = 0; j < vals.size(); ++j )
            {
                isum += ids[j];
                dsum += vals[j];
            }
        }
        // std::cout << "isum="<<isum<<"\n";
        // std::cout << "dsum="<<dsum<<"\n";
    }
    (void)units;
    (void)isum;
    (void)dsum;
}

#include <ctime>
float cpu_time()
{
    clock_t t;
    t = clock();
    return ( ( (float)t ) / CLOCKS_PER_SEC );
}

TEST( StateSet, Timing )
{
    StateSet state_set;
    FakeFactory::StateSet_random1( state_set );
    int repeat = 1000;

    {
        float t1 = cpu_time();
        timing_smart_pointer( repeat, state_set );
        float t2 = cpu_time();
        std::cout << ScaleUtils::IO::nprintf(
            "%20s %8.3f\n", "SmartPointer", ( t2 - t1 ) * 1000 / repeat );
    }

    {
        float t1 = cpu_time();
        timing_ref_at( repeat, state_set );
        float t2 = cpu_time();
        std::cout << ScaleUtils::IO::nprintf(
            "%20s %8.3f\n", "Reference_at", ( t2 - t1 ) * 1000 / repeat );
    }

    {
        float t1 = cpu_time();
        timing_ref( repeat, state_set );
        float t2 = cpu_time();
        std::cout << ScaleUtils::IO::nprintf(
            "%20s %8.3f\n", "Reference[]", ( t2 - t1 ) * 1000 / repeat );
    }
}