tstKernelMatrexC.cpp

./Solver/matrex/src/tests/tstKernelMatrexC.cpp

/******************************************************************************
 * Unit test for the c interface to old MATREX kernel
 * Author: Aarno Isotalo
 *
 *****************************************************************************/

#include <iostream>
#include "Nemesis/gtest/nemesis_gtest.hh"
#include "Origen/Core/dc/FakeFactory.h"
#include "Origen/Core/dc/TransitionMatrixP.h"
#include "Origen/Core/io/LibraryIO.h"
#include "Origen/Solver/matrex/src/kernel_matrex.hpp"
typedef std::vector<int> Vec_Int;

TEST( tstKernelMatrexC, test )
{
    int itot, non;
    Vec_Int non0, kd, loc, ZAI;
    DoubleList DLdiag, DLa;
    double flux;
    int i, j;

    // Load a library.
    Origen::Library lib;
    ASSERT_TRUE( Origen::FakeFactory::Library_scale_pwr( lib ) );
    Origen::TransitionMatrixP::SP trx = lib.newsp_transition_matrix_at( 0 );

    itot = trx->get_itot();
    non = trx->get_non();

    // Convert the Lists to vectors

    non0.resize( itot );
    kd.resize( itot );
    loc.resize( non );
    ZAI.resize( itot );

    copy( ( trx->get_non0_() )->begin(),
          ( trx->get_non0_() )->end(),
          non0.begin() );
    copy( ( trx->get_kd_() )->begin(), ( trx->get_kd_() )->end(), kd.begin() );
    copy( ( trx->get_loc() )->begin(), ( trx->get_loc() )->end(), loc.begin() );
    copy(
        ( trx->get_nucl() )->begin(), ( trx->get_nucl() )->end(), ZAI.begin() );

    // Fold library

    flux = 3.5e14;
    trx->fold_flux( flux, DLdiag, DLa );
    Origen::Vec_Dbl diag( DLdiag.begin(), DLdiag.end() );
    Origen::Vec_Dbl a( DLa.begin(), DLa.end() );

    // Set initial composition and step length

    int zero_flux_step = 0, is_adjoint = 0;

    Origen::Vec_Dbl n( itot ), n0( itot ), source_term( itot );
    n0[trx->find_nuclide_guess( 922350 )] = 0.05;
    n0[trx->find_nuclide_guess( 922380 )] = 0.95;
    double delta_t = 100 * 60 * 60 * 24;

    // solve with CRAM for reference results (The results are hardcoded because
    // CRAM might be turned off)

    int const num_ref = 4;
    int i_ref[num_ref] = {trx->find_nuclide_guess( 922350 ),
                          trx->find_nuclide_guess( 942390 ),
                          trx->find_nuclide_guess( 541350 ),
                          trx->find_nuclide_guess( 551370 )};
    double n_ref[num_ref] = {
        4.3984488451e-02, 1.9475224164e-03, 5.8150312205e-07, 3.0786326668e-04};

    // This allows the ref. results to be recalculated
    /*
    std::vector<double> n_cram(itot);
    kernel_cram (&n_cram[0], &n0[0], nullptr, -1, delta_t, itot,
                 non, &non0[0], &kd[0], &diag[0], &a[0], &loc[0],
    zero_flux_step,
                 is_adjoint, 16, &ZAI[0]);
    for (j=0; j<num_ref; j++)
      {
        i = i_ref[j];
        printf("%d %.10e\n",ZAI[i], n_cram[i]);
      }
    */

    kernel_matrex( &*trx,
                   itot,
                   non,
                   &n0[0],
                   &n[0],
                   delta_t,
                   &diag[0],
                   &a[0],
                   &source_term[0],
                   zero_flux_step,
                   is_adjoint,
                   21,
                   100,
                   0.0 );

    // check a few nuclides. This only checks that the solver actually ran so
    // there is no need for rigorous analysis

    for( j = 0; j < num_ref; j++ )
    {
        i = i_ref[j];
        EXPECT_NEAR( n[i], n_ref[j], n_ref[j] * 0.05 );
    }
}