d9/df8/tst_comp_manager_8cpp-example.html

#include "Origen/Core/xf/CompManager.h"


#include "Nemesis/gtest/nemesis_gtest.hh"


#include <cmath>

#include <cstdlib>

#include <functional>

#include <memory>

#include <numeric>

#include <sstream>


#include "Nemesis/harness/DBC.hh"

#include "ScaleInput/Parser/ScaleParser.h"

#include "Standard/Interface/jdebug.h"


#include "ScaleData/Core/StdComp.h"


using namespace std;

using namespace ScaleData;

using namespace ScaleUtils;

using Origen::CompManager;


namespace ScaleData

{

class FakeStdComp : public StdComp

{

  public:

    FakeStdComp() : StdComp()

    {

        {

            // Element 5000

            Element* element = new Element();

            element->setId( 5000 );

            element->setMass( 10.811000 );

            {

                // nuclide 5010

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 5010 );

                nuclide->setMass( 10.012940 );

                nuclide->setSigmaScatter( 2.023610 );

                nuclide->setSigmaTotal( 36.998859 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 18.430941 );

            }  // end of nuclide 5010

            {

                // nuclide 5011

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 5011 );

                nuclide->setMass( 11.009310 );

                nuclide->setSigmaScatter( 5.034320 );

                nuclide->setSigmaTotal( 5.034430 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 81.569054 );

            }  // end of nuclide 5011

            add( element );

        }  // end of Element 5000


        {

            // Element 92000

            Element* element = new Element();

            element->setId( 92000 );

            element->setMass( 238.028915 );

            {

                // nuclide 92234

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 92234 );

                nuclide->setMass( 234.040955 );

                nuclide->setSigmaScatter( 12.187160 );

                nuclide->setSigmaTotal( 16.095421 );

                nuclide->setNuSigmaF( 2.669690 );

                add( nuclide );

                element->setNuclide( nuclide, 0.005310 );

            }  // end of nuclide 92234

            {

                // nuclide 92235

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 92235 );

                nuclide->setMass( 235.043930 );

                nuclide->setSigmaScatter( 11.902490 );

                nuclide->setSigmaTotal( 35.223831 );

                nuclide->setNuSigmaF( 90.231522 );

                add( nuclide );

                element->setNuclide( nuclide, 0.711366 );

            }  // end of nuclide 92235

            {

                // nuclide 92238

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 92238 );

                nuclide->setMass( 238.050797 );

                nuclide->setSigmaScatter( 12.326360 );

                nuclide->setSigmaTotal( 14.627080 );

                nuclide->setNuSigmaF( 0.659700 );

                add( nuclide );

                element->setNuclide( nuclide, 99.283325 );

            }  // end of nuclide 92238

            add( element );

        }  // end of Element 92000


        {

            // Element 8000

            Element* element = new Element();

            element->setId( 8000 );

            element->setMass( 15.999400 );

            // Will be running the following  'add(nuclide);' for each nuclide

            {

                // nuclide 8016

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 8016 );

                nuclide->setMass( 15.994910 );

                nuclide->setSigmaScatter( 3.886960 );

                nuclide->setSigmaTotal( 3.886960 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 99.729004 );

            }  // end of nuclide 8016

            {

                // nuclide 8017

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 8017 );

                nuclide->setMass( 16.999130 );

                nuclide->setSigmaScatter( 3.740000 );

                nuclide->setSigmaTotal( 3.745010 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 0.040374 );

            }  // end of nuclide 8017

            {

                // nuclide 8018

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 8018 );

                nuclide->setMass( 17.999161 );

                nuclide->setSigmaScatter( 3.790000 );

                nuclide->setSigmaTotal( 3.790000 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 0.230623 );

            }  // end of nuclide 8018

            add( element );

        }  // end of Element 8000


        {

            // Element 1000

            Element* element = new Element();

            element->setId( 1000 );

            element->setMass( 1.007940 );

            // Will be running the following  'add(nuclide);' for each nuclide

            {

                // nuclide 1001

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 1001 );

                nuclide->setMass( 1.007830 );

                nuclide->setSigmaScatter( 20.380871 );

                nuclide->setSigmaTotal( 20.387819 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 99.977020 );

            }  // end of nuclide 1001

            {

                // nuclide 1002

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 1002 );

                nuclide->setMass( 2.014100 );

                nuclide->setSigmaScatter( 3.394860 );

                nuclide->setSigmaTotal( 3.394870 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 0.022980 );

            }  // end of nuclide 1002

            add( element );

        }  // end of Element 1000


        {

            // Element 47000

            Element* element = new Element();

            element->setId( 47000 );

            element->setMass( 107.868202 );

            // Will be running the following  'add(nuclide);' for each nuclide

            {

                // nuclide 47107

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 47107 );

                nuclide->setMass( 106.905098 );

                nuclide->setSigmaScatter( 6.530190 );

                nuclide->setSigmaTotal( 8.737310 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 51.376183 );

            }  // end of nuclide 47107

            {

                // nuclide 47109

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 47109 );

                nuclide->setMass( 108.904747 );

                nuclide->setSigmaScatter( 7.038600 );

                nuclide->setSigmaTotal( 10.780460 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 48.623821 );

            }  // end of nuclide 47109

            add( element );

        }  // end of Element 47000


        {

            // Element 49000

            Element* element = new Element();

            element->setId( 49000 );

            element->setMass( 114.818001 );

            // Will be running the following  'add(nuclide);' for each nuclide

            {

                // nuclide 49113

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 49113 );

                nuclide->setMass( 112.904060 );

                nuclide->setSigmaScatter( 4.732010 );

                nuclide->setSigmaTotal( 8.292040 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 4.218485 );

            }  // end of nuclide 49113

            {

                // nuclide 49115

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 49115 );

                nuclide->setMass( 114.903877 );

                nuclide->setSigmaScatter( 4.958740 );

                nuclide->setSigmaTotal( 7.681160 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 95.781517 );

            }  // end of nuclide 49115

            add( element );

        }  // end of Element 49000


        {

            // Element 48000

            Element* element = new Element();

            element->setId( 48000 );

            element->setMass( 112.411003 );

            // Will be running the following  'add(nuclide);' for each nuclide

            {

                // nuclide 48110

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 48110 );

                nuclide->setMass( 109.903000 );

                nuclide->setSigmaScatter( 3.782650 );

                nuclide->setSigmaTotal( 4.442520 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 12.211275 );

            }  // end of nuclide 48110

            {

                // nuclide 48111

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 48111 );

                nuclide->setMass( 110.904182 );

                nuclide->setSigmaScatter( 4.629640 );

                nuclide->setSigmaTotal( 6.828030 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 12.628360 );

            }  // end of nuclide 48111

            {

                // nuclide 48112

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 48112 );

                nuclide->setMass( 111.902763 );

                nuclide->setSigmaScatter( 4.412760 );

                nuclide->setSigmaTotal( 4.882320 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 24.020782 );

            }  // end of nuclide 48112

            {

                // nuclide 48113

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 48113 );

                nuclide->setMass( 112.904404 );

                nuclide->setSigmaScatter( 5.824120 );

                nuclide->setSigmaTotal( 7.657330 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 12.273576 );

            }  // end of nuclide 48113

            {

                // nuclide 48114

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 48114 );

                nuclide->setMass( 113.903358 );

                nuclide->setSigmaScatter( 4.068730 );

                nuclide->setSigmaTotal( 4.524480 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 29.111269 );

            }  // end of nuclide 48114

            {

                // nuclide 48116

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 48116 );

                nuclide->setMass( 115.904762 );

                nuclide->setSigmaScatter( 4.475170 );

                nuclide->setSigmaTotal( 4.621490 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 7.722753 );

            }  // end of nuclide 48116

            {

                // nuclide 48106

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 48106 );

                nuclide->setMass( 105.906464 );

                nuclide->setSigmaScatter( 6.988970 );

                nuclide->setSigmaTotal( 7.762620 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 1.177664 );

            }  // end of nuclide 48106

            {

                // nuclide 48108

                Nuclide* nuclide = new Nuclide();

                nuclide->setId( 48108 );

                nuclide->setMass( 107.904182 );

                nuclide->setSigmaScatter( 4.777420 );

                nuclide->setSigmaTotal( 5.024270 );

                nuclide->setNuSigmaF( 0.000000 );

                add( nuclide );

                element->setNuclide( nuclide, 0.854314 );

            }  // end of nuclide 48108

            add( element );

        }  // end of Element 48000


        {

            // Element 1095242

            Nuclide* nuclide = new Nuclide();

            nuclide->setId( 1095242 );

            nuclide->setMass( 242.059509 );

            nuclide->setSigmaScatter( 11.695570 );

            nuclide->setSigmaTotal( 67.769852 );

            nuclide->setNuSigmaF( 1494.928467 );

            add( nuclide );

        }  // end of Element 1095242


        {

            // Element 95242

            Nuclide* nuclide = new Nuclide();

            nuclide->setId( 95242 );

            nuclide->setMass( 242.059555 );

            nuclide->setSigmaScatter( 10.692500 );

            nuclide->setSigmaTotal( 33.912220 );

            nuclide->setNuSigmaF( 457.578094 );

            add( nuclide );

        }  // end of Element 95242

    }

};

}  // end namespace ScaleData


// Use a google macro

#define EXPECT_EMPTY_STR( action_and_result_str ) \

    EXPECT_EQ( "", action_and_result_str )


TEST( CompManager, test_water )

{

    // Main testing procedure.

    FakeStdComp fakeStdComp;

    Origen::CompManager mCompMan = CompManager( fakeStdComp );


    // evaluate elemental H2O

    EXPECT_EMPTY_STR(

        mCompMan.initComp( "elemental_h2o", "FORM", "ABS", true ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "elemental_h2o", "H", 2.0 ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "elemental_h2o", "O", 1.0 ) );

    EXPECT_EMPTY_STR( mCompMan.finalizeComp( "elemental_h2o" ) );

    std::cout << mCompMan.toString( "elemental_h2o" ) << std::endl;


    // double H2O

    EXPECT_EMPTY_STR( mCompMan.initComp( "double_h2o", "FORM", "ABS", true ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "double_h2o", "elemental_h2o", 1.5 ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "double_h2o", "elemental_h2o", 1.5 ) );

    EXPECT_EMPTY_STR( mCompMan.finalizeComp( "double_h2o" ) );

    jDebugLine( mCompMan.toString( "double_h2o" ) );


    {

        std::map<int, double> map_out;


        EXPECT_EMPTY_STR( mCompMan.getNumMap( "elemental_h2o", map_out ) );


        double dens = 1.0;

        EXPECT_EMPTY_STR(

            mCompMan.getConcMap( "elemental_h2o", dens, map_out ) );


        EXPECT_EMPTY_STR( mCompMan.getNumMap( "double_h2o", map_out ) );


        EXPECT_EMPTY_STR( mCompMan.getConcMap( "double_h2o", dens, map_out ) );

    }


    // evaluate simple H2O

    EXPECT_EMPTY_STR( mCompMan.initComp( "simple_h2o", "FORM", "ABS", true ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "simple_h2o", "1001", 2.0 ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "simple_h2o", "8016", 1.0 ) );

    EXPECT_EMPTY_STR( mCompMan.finalizeComp( "simple_h2o" ) );

    jDebugLine( mCompMan.toString( "simple_h2o" ) );


    // borated H2O

    double ppmbor = 1300;

    EXPECT_EMPTY_STR( mCompMan.initComp( "borated_h2o", "WT", "PPM", true ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "borated_h2o", "elemental_h2o", 1e6 ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "borated_h2o", "B", ppmbor ) );

    EXPECT_EMPTY_STR( mCompMan.finalizeComp( "borated_h2o" ) );

    jDebugLine( mCompMan.toString( "borated_h2o" ) );


    // test ppm boron

    {

        std::map<int, double> map_out;

        mCompMan.getWtMap( "borated_h2o", map_out );

        double ppmbor_check = 1e6 * ( map_out[5010] + map_out[5011] ) /

                              ( map_out[1001] + map_out[1002] + map_out[8016] +

                                map_out[8017] + map_out[8018] );

        EXPECT_SOFTEQ( ppmbor, ppmbor_check, 5.e-7 );

    }


    // natural UO2

    EXPECT_EMPTY_STR( mCompMan.initComp( "UO2[nat]", "NUM", "ABS", true ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "UO2[nat]", "U", 1.0 ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "UO2[nat]", "O", 2.0 ) );

    EXPECT_EMPTY_STR( mCompMan.finalizeComp( "UO2[nat]" ) );

    jDebugLine( mCompMan.toString( "UO2[nat]" ) );


    // AIC control rods

    EXPECT_EMPTY_STR( mCompMan.initComp( "aic", "WT", "PCT", false ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "aic", "Ag", 80.0 ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "aic", "In", 15.0 ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "aic", "Cd", 5.0 ) );

    EXPECT_EMPTY_STR( mCompMan.finalizeComp( "aic" ) );

    jDebugLine( mCompMan.toString( "aic" ) );


    // enriched UO2

    EXPECT_EMPTY_STR( mCompMan.initComp( "U[2.5]", "WT", "PCT", false ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "U[2.5]", "U-235", 2.5 ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "U[2.5]", "U-238", -100.0 ) );

    EXPECT_EMPTY_STR( mCompMan.initComp( "UO2[2.5]", "NUM", "ABS", true ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "UO2[2.5]", "U[2.5]", 1.0 ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "UO2[2.5]", "O", 2.0 ) );

    EXPECT_EMPTY_STR( mCompMan.finalizeComp( "UO2[2.5]" ) );

    jDebugLine( mCompMan.toString( "UO2[2.5]" ) );


    // using metastable

    EXPECT_EMPTY_STR( mCompMan.initComp( "myAm", "WT", "PCT", false ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "myAm", "Am-242m", 2.5 ) );

    EXPECT_EMPTY_STR( mCompMan.addPart( "myAm", "Am-242", -100.0 ) );

    EXPECT_EMPTY_STR( mCompMan.finalizeComp( "myAm" ) );

    jDebugLine( mCompMan.toString( "myAm" ) );


    std::map<int, double> map_out;


    EXPECT_EMPTY_STR( mCompMan.getNumMap( "borated_h2o", map_out ) );


    double dens = 1.0;

    EXPECT_EMPTY_STR( mCompMan.getConcMap( "borated_h2o", dens, map_out ) );


    dens = 10.16;

    EXPECT_EMPTY_STR( mCompMan.getConcMap( "aic", dens, map_out ) );


    // test getDataContainer

    DataContainer dcCompositionList;

    EXPECT_EMPTY_STR( mCompMan.getDC( &dcCompositionList, false ) );


    EXPECT_EQ( 16, mCompMan.getNumComps() );


    // test getName

    EXPECT_EQ( "elemental_h2o", mCompMan.getCompName( 1 ) );


    // test ppm boron by number (instead of by weight)

    {

        std::map<int, double> map_out;

        mCompMan.getNumMap( "borated_h2o", map_out );

        double ppmbor_check = 1e6 * ( map_out[5010] + map_out[5011] ) /

                              ( map_out[1001] + map_out[1002] + map_out[8016] +

                                map_out[8017] + map_out[8018] );

        EXPECT_SOFTEQ( 722.098, ppmbor_check, 5.e-7 );

    }


    dcCompositionList.deepDelete();

}


TEST( CompManager, FailNotValidNuclide )

{

    // Main testing procedure.

    FakeStdComp fakeStdComp;

    Origen::CompManager mCompMan = CompManager( fakeStdComp );

    EXPECT_EMPTY_STR( mCompMan.initComp( "x", "WT", "PCT", false ) );

    EXPECT_EQ(

        "part='u2' --> 92002 is not a valid NUCLIDE in this StdCompLibrary",

        mCompMan.addPart( "x", "u2", 100 ) );

    EXPECT_EQ(

        "cannot finalize composition name=x unless number fractions >=0 and "

        "sum to 1.0 (sum=0.00000e+00)",

        mCompMan.finalizeComp( "x" ) );

}


TEST( CompManager, Fail8001 )

{

    // Main testing procedure.

    FakeStdComp fakeStdComp;

    Origen::CompManager mCompMan = CompManager( fakeStdComp );

    EXPECT_EMPTY_STR( mCompMan.initComp( "x", "WT", "PCT", false ) );

    // 8008 is first allowed "oxygen" by ZA rules that A>=Z because A=Z+N

    EXPECT_EQ( "part='8001' is not a valid NUCLIDE/ELEMENT id",

               mCompMan.addPart( "x", "8001", 100 ) );

    EXPECT_EQ(

        "cannot finalize composition name=x unless number fractions >=0 and "

        "sum to 1.0 (sum=0.00000e+00)",

        mCompMan.finalizeComp( "x" ) );

}