tstCRAMcommons.cpp
./Solver/cram/src/tests/tstCRAMcommons.cpp
/******************************************************************************
* Common functions for the tstCRAM* tests. *
* *
*****************************************************************************/
#include "tstCRAMcommons.hpp"
#include <cmath>
using std::printf;
using std::fprintf;
/******************************************************************************
* rand() is not quaranteed to be platform independent. While the tests should*
* work with any random sequence, one mixes random numbers with an externally *
* calculated reference solution and others involve some tolerances that are *
* essentially arbitrary and might not quite work for some specific randoms. *
* *
* parameters are from numerical recipes (actually from wikipedia) *
*****************************************************************************/
{
long long int const a = 1664525LL, c = 1013904223LL, m = 1013904223LL;
long long int static val = 1337762451LL;
if( seed > 0 ) val = seed;
val = ( a * val + c ) % m;
return (double)val / (double)( m - 1 );
}
/******************************************************************************
* Print debug/dev output about the accuracy of the results. *
* *
* This does not (yet?) perform any actual tests on the results, so this part *
* is always passed in the gtest sense. *
* *
* Note that the sign of the difference is ignored. The signs of the results *
* are also ignored for the purpose of binning. *
* *
* Input: *
* itot number of nuclides *
* n[itot] pointer to an array with the resultsx) *
* n_ref[itot] pointer to array with reference results *
* *
*****************************************************************************/
{
double diff, ref_adens;
int afrac_bin;
int i, j;
// bins for relative error and atomic fraction (2D)
const int n_error_bins = 19;
double error_limits[n_error_bins + 1] = {
0.0, 1e-16, 1e-15, 1e-14, 1e-13, 1e-12, 1e-10, 1e-9, 1e-8, 1e-7,
1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e-0, 2e0, 1e2, 1e99};
const int n_afrac_bins = 9;
double afrac_limits[n_afrac_bins + 1] = {
0.0, 1e-50, 1e-40, 1e-30, 1e-25, 1e-20, 1e-15, 1e-10, 1e-5, 1};
int larger[n_error_bins][n_afrac_bins] = {{0}};
int afrac_bin_totals[n_afrac_bins] = {0};
// calculate total atomic density (in ref)
ref_adens = 0.0;
for( i = 0; i < itot; i++ ) ref_adens += std::fabs( n_ref[i] );
for( i = 0; i < itot; i++ )
{
// find atomic fraction bin for nuclide i
for( afrac_bin = 0; afrac_bin < n_afrac_bins; afrac_bin++ )
if( std::fabs( n_ref[i] / ref_adens ) <
afrac_limits[afrac_bin + 1] )
break;
afrac_bin_totals[afrac_bin]++;
// find error bin for nuclide i
double dmin_double = std::numeric_limits<double>::denorm_min();
if( n_ref[i] != 0.0 )
diff = std::fabs( n[i] / n_ref[i] - 1 );
else if( n[i] == 0 )
diff = 0.0;
else if( n[i] == dmin_double )
diff =
1; // sign that the result from cram was culled, mark as 1 even
// if the real rel dif is infinite
else
diff = 1e99;
for( j = 0; j < n_error_bins - 1; j++ )
if( diff < error_limits[j + 1] ) break;
larger[j][afrac_bin]++;
}
// Write an overview type output about the accuracy
printf( "\n log10(afrac) | " );
for( i = 0; i < n_afrac_bins; i++ )
printf( " % 4d", (int)std::round( std::log10( afrac_limits[i + 1] ) ) );
printf( "\n" );
for( j = n_error_bins - 1; j >= 0; j-- )
{
printf( "%.1e-%.1e|", error_limits[j], error_limits[j + 1] );
for( i = 0; i < n_afrac_bins; i++ ) printf( " %4d", larger[j][i] );
printf( "\n" );
}
printf( "\n" );
}
/******************************************************************************
* Write a file with all the input needed for kernel_cram and the results it *
* generated. *
* *
* Arguments are the same as for kernel_cram() except that n is now input, *
* reporter is not passed and there are also two additional arguments: *
* comment is written to the file. DO NOT ADD NEWLINE *
* (Hint: make it the library name and folding flux) *
* filepath is path to the file where to write (will be overwritten) *
* *
*****************************************************************************/
double* n0,
double** source_term,
int source_order,
double delta_t,
int itot,
int non,
int* non0,
int* kd,
double* diag,
double* a,
int* loc,
int zero_flux_step,
int is_adjoint,
int cram_order,
int internal_steps,
int* ZAI,
const char* filepath,
const char* comment )
{
int i, j;
FILE* fp;
fp = std::fopen( filepath, "w" );
// get the time
std::time_t time_now = std::time( nullptr );
// write the header
fprintf( fp,
"// Written by tstCRAM_WriteBlessedResults() at %s",
std::asctime( std::localtime( &time_now ) ) );
fprintf( fp,
"// Presumably there was high certainty that the solver worked "
"back then.\n" );
fprintf( fp, "// Comment: %s\n", comment );
fprintf( fp, "//\n" );
// write scalars
fprintf( fp,
"// itot non zero_flux_step is_adjoint cram_order "
"internal_steps source_order delta_t\n" );
fprintf( fp,
" %5d %6d %1d %1d %3d %3d "
" %7d %.17e\n",
itot,
non,
zero_flux_step,
is_adjoint,
cram_order,
internal_steps,
source_order,
delta_t );
// write [itot] arrays
fprintf( fp,
"// ZAI[i] n[i] n_0[i] "
"non0[i] kd[i] diag[i] for i=0... itot-1\n" );
for( i = 0; i < itot; i++ )
fprintf( fp,
" %7d % .17e % .17e %7i %7i % .17e\n",
ZAI[i],
n[i],
n0[i],
non0[i],
kd[i],
diag[i] );
// write source term
if( source_order >= 0 )
{
fprintf(
fp,
"// source_term[i][j] for i=0... itot-1, j=0... source_order\n" );
for( i = 0; i < itot; i++ )
{
fprintf( fp, " " );
for( j = 0; j < source_order + 1; j++ )
fprintf( fp, " % .17e", source_term[i][j] );
fprintf( fp, "\n" );
}
}
// write [non] arrays
fprintf( fp,
"// loc[i] a[i] for i=0... non-1 (a=offdiag in "
"kernel_cram)\n" );
for( i = 0; i < non; i++ ) fprintf( fp, " %6i %.17e\n", loc[i], a[i] );
fprintf( fp, "\n" );
std::fclose( fp );
}
/******************************************************************************
* Read the file generated by WriteBlessedResults() *
* *
* Arguments are the same as for WriteBlessedResults() except they are now *
* pointers to scalars or vectors (because: Look mama, with c++ features!). *
* The comment is also not recovered. *
* *
*****************************************************************************/
Vec_Dbl* n0,
double*** source_term,
int* source_order,
double* delta_t,
int* itot,
int* non,
Vec_Int* non0,
Vec_Int* kd,
Vec_Dbl* diag,
Vec_Dbl* a,
Vec_Int* loc,
int* zero_flux_step,
int* is_adjoint,
int* cram_order,
int* internal_steps,
Vec_Int* ZAI,
const char* filepath )
{
int i, j;
std::string line;
// open the file
std::ifstream infile( filepath );
if( infile.fail() )
{
std::cerr << "Cannot open " << filepath << "for reading.\n";
std::exit( EXIT_FAILURE );
}
// skip header
std::getline( infile, line );
std::getline( infile, line );
std::getline( infile, line );
std::getline( infile, line );
// read: itot non zero_flux_step is_adjoint cram_order source_order delta_t
std::getline( infile, line );
std::getline( infile, line );
std::istringstream iss( line );
if( !( iss >> *itot >> *non >> *zero_flux_step >> *is_adjoint >>
*cram_order >> *internal_steps >> *source_order >> *delta_t ) )
{
std::cerr << "Unexpected line in " << filepath << ":\n";
std::cerr << " " << line << "\n";
std::exit( EXIT_FAILURE );
}
// resize the vectors so we can stream directly to the elements
( *n ).resize( *itot );
( *n0 ).resize( *itot );
( *diag ).resize( *itot );
( *a ).resize( *non );
( *non0 ).resize( *itot );
( *kd ).resize( *itot );
( *ZAI ).resize( *itot );
( *loc ).resize( *non );
// read: ZAI[i], n[i], n0[i], source_term[i], non0[i], kd[i], diag[i]
std::getline( infile, line );
for( i = 0; i < *itot; i++ )
{
std::getline( infile, line );
std::istringstream iss( line );
if( !( iss >> ( *ZAI )[i] >> ( *n )[i] >> ( *n0 )[i] >> ( *non0 )[i] >>
( *kd )[i] >> ( *diag )[i] ) )
{
std::cerr << "Unexpected line in " << filepath << ":\n";
std::cerr << " " << line << "\n";
std::exit( EXIT_FAILURE );
}
}
// read source term
std::getline( infile, line );
if( source_term >= 0 )
{
*source_term = (double**)std::malloc( *itot * sizeof( double* ) );
for( i = 0; i < *itot; i++ )
( *source_term )[i] =
(double*)std::calloc( *source_order + 1, sizeof( double ) );
for( i = 0; i < *itot; i++ )
{
std::getline( infile, line );
std::istringstream iss( line );
for( j = 0; j < ( *source_order ) + 1; j++ )
if( !( iss >> ( *source_term )[i][j] ) )
{
std::cerr << "Unexpected line in " << filepath << ":\n";
std::cerr << " " << line << "\n";
std::exit( EXIT_FAILURE );
}
}
}
else
source_term = nullptr;
// read: loc[i], a[i]
std::getline( infile, line );
for( i = 0; i < *non; i++ )
{
std::getline( infile, line );
std::istringstream iss( line );
if( !( iss >> ( *loc )[i] >> ( *a )[i] ) )
{
std::cerr << "Unexpected line in " << filepath << ":\n";
std::cerr << " " << line << "\n";
std::exit( EXIT_FAILURE );
}
}
}
