dc/d1a/tst_kernel_radau5_c_8cpp-example.html

/******************************************************************************

 * Unit test for the c interface to RADAU5 solver 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/radau5/src/kernel_radau5.hpp"

#include "ScaleUtils/IO/nprintf.h"


typedef std::vector<int> Vec_Int;

typedef std::vector<double> Vec_Dbl;

using ScaleUtils::IO::nprintf;


TEST( tstSolverRadau5, test )

{

    int itot, non;

    Vec_Int non0, kd, loc, ZAI;

    DoubleList DLdiag, DLa;

    double flux;

    int i, j;


    // Load a library

    Origen::Library lib;

    bool loaded = Origen::FakeFactory::Library_scale_pwr( lib );

    Insist( loaded, "load failed!" );

    Origen::TransitionMatrixP::SP trx = lib.newsp_transition_matrix_at( 0 );


    EXPECT_EQ( 0, trx->first_transition( 0 ) );

    EXPECT_EQ( 555, trx->last_transition( 0, 1.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 );

    Vec_Dbl diag( DLdiag.begin(), DLdiag.end() );

    Vec_Dbl a( DLa.begin(), DLa.end() );


    // Set initial composition and step length


    int zero_flux_step = 0, is_adjoint = 0;


    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 = 1 * 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.9935947875e-02, 3.1830272138e-06, 5.0309302491e-07, 3.0836538981e-06};


    // 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]);

      }

    */


    // the toleranses are rather high to make the running time reasonable

    // (even with just one day step)


    kernel_radau5( &*trx,

                   itot,

                   non,

                   &n0[0],

                   &n[0],

                   delta_t,

                   &diag[0],

                   &a[0],

                   &source_term[0],

                   zero_flux_step,

                   is_adjoint,

                   1.0e-10,

                   1.0e-3 );


    // 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];

        // std::cout << nprintf("%4d %.5e %.5e \n",i,n[i],n_ref[j]);

        EXPECT_NEAR( n[i], n_ref[j], n_ref[j] * 0.05 );

    }

}