tstYieldResource.cpp

./Core/re/tests/tstYieldResource.cpp

#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <istream>
#include <sstream>
#include <string>
#include <vector>

#include "Origen/Core/Definitions.h"
#include "Origen/Core/dc/FakeFactory.h"
#include "Origen/Core/fn/print.h"
#include "Origen/Core/re/YieldResource.h"

#include "ScaleSTL/Functions.h"

#include "Nemesis/gtest/nemesis_gtest.hh"

TEST( YieldResource, BasicAccessors )
{
    Origen::YieldResource yr;

    // Add a nuclide.
    yr.add_fissionable( 92235 );
    {
        Origen::YieldParent* yp = yr.fissionable( 92235 );
        EXPECT_TRUE( yp != nullptr );
    }

    // Reference from adding and getting the same.
    {
        Origen::YieldParent& yp = yr.add_fissionable( 1092235 );
        Origen::YieldParent* yp2 = yr.fissionable( 1092235 );
        ASSERT_TRUE( yp2 != nullptr );
        EXPECT_TRUE( &yp == yp2 );
    }

    // Can add a data set to fissionable nuclide.
    {
        Origen::YieldParent* yp1 = yr.fissionable( 92235 );
        ASSERT_TRUE( yp1 != nullptr );
        Origen::YieldParent& yp = *yp1;

        {
            Origen::Vec_Int id = {54135, 1054135, 49102};
            Origen::Vec_Dbl val = {0.1, 0.2, 0.3};
            yp.add_data( 0.0623, id, val );
        }

        int ienergy = yp.energy_index( 0.0623 );
        ASSERT_FALSE( ienergy == -1 );

        Origen::YieldData& ys = yp.data_at( ienergy );

        EXPECT_FLOAT_EQ( 0.0623, ys.energy() );

        EXPECT_EQ( 3, ys.num_yields() );

        int k1 = yp.fp_index( 54135 );
        EXPECT_NE( -1, k1 );

        int k2 = yp.fp_index( 1054135 );
        EXPECT_NE( -1, k2 );

        int k3 = yp.fp_index( 49102 );
        EXPECT_NE( -1, k3 );

        int k4 = yp.fp_index( 50102 );
        EXPECT_EQ( -1, k4 );
    }

    // Can add another data set with different yields.
    {
        Origen::YieldParent* yp1 = yr.fissionable( 92235 );
        ASSERT_TRUE( yp1 != nullptr );
        Origen::YieldParent& yp = *yp1;

        {
            Origen::Vec_Int id = {50102};
            Origen::Vec_Dbl val = {0.4};
            yp.add_data( 2e4, id, val );
        }

        int ienergy = yp.energy_index( 2e4 );
        ASSERT_FALSE( ienergy == -1 );

        Origen::YieldData& ys = yp.data_at( ienergy );

        EXPECT_FLOAT_EQ( 2e4, ys.energy() );

        EXPECT_EQ( 4, ys.num_yields() );

        int k1 = yp.fp_index( 54135 );
        EXPECT_NE( -1, k1 );

        int k2 = yp.fp_index( 1054135 );
        EXPECT_NE( -1, k2 );

        int k3 = yp.fp_index( 49102 );
        EXPECT_NE( -1, k3 );

        int k4 = yp.fp_index( 50102 );
        EXPECT_NE( -1, k4 );

        // should now be able to get all 4 yield values for both energies
        EXPECT_FLOAT_EQ( 0.0, ys.yield_at( k1 ) );
        EXPECT_FLOAT_EQ( 0.0, ys.yield_at( k2 ) );
        EXPECT_FLOAT_EQ( 0.0, ys.yield_at( k3 ) );
        EXPECT_FLOAT_EQ( 0.4, ys.yield_at( k4 ) );

        // for the other energy we should also be able to get all 4
        int ienergy2 = yp.energy_index( 0.0623 );
        ASSERT_FALSE( ienergy2 == -1 );

        Origen::YieldData& ys2 = yp.data_at( ienergy2 );

        EXPECT_FLOAT_EQ( 0.1, ys2.yield_at( k1 ) );
        EXPECT_FLOAT_EQ( 0.2, ys2.yield_at( k2 ) );
        EXPECT_FLOAT_EQ( 0.3, ys2.yield_at( k3 ) );
        EXPECT_FLOAT_EQ( 0.0, ys2.yield_at( k4 ) );
    }

    // Can add a new fission product and interpolate
    {
        Origen::YieldParent* yp1 = yr.fissionable( 92235 );
        ASSERT_TRUE( yp1 != nullptr );

        int k = yp1->add_fp( 51103 );
        ASSERT_NE( 5, k );
        ASSERT_EQ( 2, yp1->num_energies() );

        EXPECT_FLOAT_EQ( 0.0, yp1->data_at( 0 ).yield_at( k ) );
        EXPECT_FLOAT_EQ( 0.0, yp1->data_at( 1 ).yield_at( k ) );

        ASSERT_TRUE( yp1->fp_index( 49102 ) != -1 );

        float v0 = yp1->interp_yield( 49102, 0.0623 );
        EXPECT_FLOAT_EQ( 0.3, v0 );

        float v1 = yp1->interp_yield( 51103, 2e4 );
        EXPECT_FLOAT_EQ( 0.0, v1 );

        float v2 = yp1->interp_yield( 49102, 1e4 );
        float w1 = ( 2e4 - 1e4 ) / ( 2e4 - 0.0623 );
        float w2 = 1.0 - w1;
        float v2_ref = w1 * 0.0 + w2 * 0.3;
        EXPECT_NEAR( v2_ref, v2, 0.0001 );
    }
}

TEST( YieldResource, AllEnergy )
{
    Origen::YieldResource yr;

    // first nuclide
    Origen::YieldParent& yp1 = yr.add_fissionable( 92235 );
    {
        Origen::Vec_Int id = {54135, 1054135, 49102};
        Origen::Vec_Dbl val = {0.1, 0.2, 0.3};
        yp1.add_data( 0.0623, id, val );
    }
    {
        Origen::Vec_Int id = {54135, 1054135, 49102};
        Origen::Vec_Dbl val = {0.15, 0.25, 0.35};
        yp1.add_data( 1e6, id, val );
    }
    {
        Origen::Vec_Int id = {54135, 1054135, 49102};
        Origen::Vec_Dbl val = {0.2, 0.3, 0.4};
        yp1.add_data( 2e7, id, val );
    }

    // second nuclide
    Origen::YieldParent& yp2 = yr.add_fissionable( 94239 );
    {
        Origen::Vec_Int id = {50102};
        Origen::Vec_Dbl val = {0.33};
        yp2.add_data( 1e7, id, val );
    }

    Origen::Vec_Dbl energies;
    yr.get_all_energies( &energies );

    ASSERT_EQ( 4u, energies.size() );
    EXPECT_FLOAT_EQ( 2e7, energies[0] );
    EXPECT_FLOAT_EQ( 1e7, energies[1] );
    EXPECT_FLOAT_EQ( 1e6, energies[2] );
    EXPECT_FLOAT_EQ( 0.0623, energies[3] );
}

class YieldResourceTester : public ::testing::Test, public Origen::YieldResource
{
  protected:
    virtual void SetUp()
    {
        bool status = Origen::FakeFactory::YieldResource_scale( *this );
        ASSERT_TRUE( status );
    }
};

TEST_F( YieldResourceTester, IOFunctions )
{
    Origen::printInterpolatedYields( *this, 92235, 1.0e4, 0, std::cout );
    Origen::printYields( *this, std::cout );
    Origen::printYieldEnergies( *this, std::cout );
}

TEST_F( YieldResourceTester, Removal )
{
    // check existence of nuclides we will use in testing
    ASSERT_TRUE( this->fissionable( 92235 ) != nullptr );
    ASSERT_TRUE( this->fissionable( 94239 ) != nullptr );

    // uncomment to print initial yields
    // std::ofstream orig("orig.txt");
    // Origen::printYieldEnergies(*this,orig);

    // get Cs135 yields before
    double p9_cs135_yield = this->interp_yield( 94239, 53135, 0.5e6 );
    double p0_cs135_yield = this->interp_yield( 94240, 53135, 0.5e6 );

    // do some simple accessing to verify we can remove something
    int iu235 = this->fissionable_index( 92235 );
    EXPECT_TRUE( iu235 != -1 );

    int num_fissionables = this->num_fissionables();

    bool found = this->remove_fissionable( 92235 );
    EXPECT_TRUE( found );

    iu235 = this->fissionable_index( 92235 );
    EXPECT_TRUE( iu235 == -1 );

    EXPECT_EQ( num_fissionables - 1, this->num_fissionables() );

    // get Cs135 yields after
    EXPECT_FLOAT_EQ( 0.0, this->interp_yield( 92235, 53135, 0.5e6 ) );
    EXPECT_FLOAT_EQ( p9_cs135_yield,
                     this->interp_yield( 94239, 53135, 0.5e6 ) );
    EXPECT_FLOAT_EQ( p0_cs135_yield,
                     this->interp_yield( 94240, 53135, 0.5e6 ) );

    // create a system with only Pu239, Pu240, Am242m in it
    Origen::YieldParent* pu239 = this->fissionable( 94239 );
    ASSERT_TRUE( pu239 != nullptr );

    Origen::YieldParent* pu240 = this->fissionable( 94240 );
    ASSERT_TRUE( pu240 != nullptr );

    Origen::YieldParent* am242m = this->fissionable( 1095242 );
    ASSERT_TRUE( am242m != nullptr );

    Origen::Vec_Int ids2 = {94239, 94240, 1095242};
    this->set_fissionable_ids( ids2 );
    EXPECT_EQ( 3, this->num_fissionables() );

    // need to get new pointers after modifying system
    pu239 = this->fissionable( 94239 );
    ASSERT_TRUE( pu239 != nullptr );
    pu240 = this->fissionable( 94240 );
    ASSERT_TRUE( pu240 != nullptr );
    am242m = this->fissionable( 1095242 );
    ASSERT_TRUE( am242m != nullptr );

    ASSERT_TRUE( pu239 != nullptr );
    ASSERT_EQ( 3, pu239->num_energies() );

    Origen::Vec_Dbl energies;
    pu239->get_energies( &energies );
    EXPECT_EQ( 3u, energies.size() );

    // remove middle energy and see if we get close to the same (5% tolerance)
    pu239->remove_data_at( 1 );  // remove middle energy
    EXPECT_NEAR( p9_cs135_yield,
                 pu239->interp_yield( 53135, 0.5e6 ),
                 0.05 * p9_cs135_yield );
    // new reference
    p9_cs135_yield = pu239->interp_yield( 53135, 0.5e6 );

    Origen::Vec_Dbl new_energies;
    pu239->get_energies( &new_energies );
    ASSERT_EQ( 2u, new_energies.size() );

    EXPECT_FLOAT_EQ( energies[0], new_energies[0] );
    EXPECT_FLOAT_EQ( energies[2], new_energies[1] );

    found = pu239->remove_fp( 56135 );
    EXPECT_TRUE( found );
    EXPECT_FALSE( pu239->has_fp( 56135 ) );

    // keep only the I/Xe mass chain
    Origen::Vec_Int ids3 = {52135, 53135, 54135, 55135, 56135};
    pu239->set_fp_ids( ids3 );
    EXPECT_EQ( 5, pu239->num_fps() );

    // this was added back in the set_fp_ids above but should now be zero yield
    EXPECT_TRUE( pu239->has_fp( 56135 ) );
    EXPECT_EQ( 0.0, pu239->interp_yield( 0.5e6, 56135 ) );

    // still should get this yield the same
    EXPECT_FLOAT_EQ( p9_cs135_yield,
                     this->interp_yield( 94239, 53135, 0.5e6 ) );

    // Now pu239 has only 5 fission products but everything else has many more.
    // We are going to restrict all fission products.
    Origen::Vec_Int ids4 = {
        1003, 2003, 2004, 52135, 53135, 54135, 55135, 56135};
    this->set_fp_ids( ids4 );

    // don't need new pointers
    EXPECT_EQ( 8, pu239->num_fps() );
    EXPECT_EQ( 8, pu240->num_fps() );
    EXPECT_EQ( 8, am242m->num_fps() );

    EXPECT_FLOAT_EQ( p9_cs135_yield, pu239->interp_yield( 53135, 0.5e6 ) );
    EXPECT_FLOAT_EQ( p0_cs135_yield,
                     this->interp_yield( 94240, 53135, 0.5e6 ) );

    // uncomment to print final yields
    // std::ofstream final("final.txt");
    // Origen::printYields(*this,final);

    // final clear
    this->clear();
    EXPECT_EQ( 0, this->num_fissionables() );
}