Commit 2d1c9684 authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

merged with jannick

parents faf089e7 f6af86e2
......@@ -12,7 +12,7 @@ set( KITRT_DEBUG_OPTIONS -Wall -Wextra -Wpedantic )
### LIBRARIES ###################################
set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH})
find_package( OpenMP REQUIRED )
......@@ -82,6 +82,8 @@ target_compile_options( ${CMAKE_PROJECT_NAME} PUBLIC "$<$<CONFIG:RELEASE>:${KITR
### BUILD UNIT TESTS ############################
include( CTest )
set( CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/ext/Catch2/contrib ${CMAKE_MODULE_PATH} )
include( Catch )
set( CATCH_INCLUDE_DIR ${CMAKE_SOURCE_DIR}/ext/Catch2/single_include/catch2 )
add_compile_definitions( BUILD_TESTING )
add_compile_definitions( TESTS_PATH="${CMAKE_SOURCE_DIR}/tests/" )
......@@ -99,6 +101,6 @@ if( CMAKE_COMPILER_IS_GNUCXX )
endif()
target_compile_options( unit_tests PUBLIC "$<$<CONFIG:RELWITHDEBINFO>:${KITRT_RELWITHDEBINFO_OPTIONS}>" )
target_compile_options( unit_tests PUBLIC "$<$<CONFIG:RELEASE>:${KITRT_RELEASE_OPTIONS}>" )
add_test( NAME unit_tests COMMAND unit_tests )
catch_discover_tests( unit_tests )
enable_testing()
#################################################
......@@ -257,7 +257,6 @@ void ICRU::ELINIT( std::vector<double> IZ, std::vector<double> STF ) {
void ICRU::DCSEL0( double E ) {
if( E < 49.9999e0 || E > 1.0001e8 ) {
std::cout << "Energy " << E << std::endl;
ErrorMessages::Error( "Outside of supported energy range", CURRENT_FUNCTION );
}
double EL = std::log( E );
......
......@@ -97,7 +97,7 @@ void MLOptimizer::initialize_python() {
PyRun_SimpleString( ( "import sys\nsys.path.append('" + pyPath + "')" ).c_str() );
}
std::cout << "Python working directory is: " << pyPath << " \n";
// std::cout << "Python working directory is: " << pyPath << " \n";
init_numpy();
}
......
......@@ -21,12 +21,11 @@ VectorVector MuscleBoneLung::SetupIC() {
return psi;
}
std::vector<double> MuscleBoneLung::GetDensity( const VectorVector& cellMidPoints ) {
std::cout<<"Length of mesh "<< cellMidPoints.size() << std::endl;
std::vector<double> densities ( 167, 1.04 ); //muscle layer
std::vector<double> bone( 167, 1.85 );
std::vector<double> lung ( 665, 0.3 ); //maybe this is just the lung tissue and after that there should be air?
densities.insert(densities.end(),bone.begin(),bone.end());
densities.insert(densities.end(),lung.begin(),lung.end());
std::cout<<"Length of densities "<< densities.size() << std::endl;
return densities; }
std::vector<double> MuscleBoneLung::GetDensity( const VectorVector& cellMidPoints ) {
std::vector<double> densities( 167, 1.04 ); // muscle layer
std::vector<double> bone( 167, 1.85 );
std::vector<double> lung( 665, 0.3 ); // maybe this is just the lung tissue and after that there should be air?
densities.insert( densities.end(), bone.begin(), bone.end() );
densities.insert( densities.end(), lung.begin(), lung.end() );
return densities;
}
......@@ -2,6 +2,7 @@
#include "catch.hpp"
#include <Python.h>
#include <filesystem>
#include <mpi.h>
int main( int argc, char** argv ) {
......@@ -12,6 +13,8 @@ int main( int argc, char** argv ) {
const int result = Catch::Session().run( argc, argv );
if( Py_IsInitialized() ) Py_Finalize();
std::filesystem::remove_all( std::string( TESTS_PATH ) + "result" );
MPI_Finalize();
return result;
}
......@@ -46,9 +46,11 @@ TEST_CASE( "SN_SOLVER", "[validation_tests]" ) {
double eps = 1e-3;
REQUIRE( test.size() == reference.size() );
bool errorWithinBounds = true;
for( unsigned i = 0; i < test.size(); ++i ) {
REQUIRE( std::fabs( test[i] - reference[i] ) < eps );
if( std::fabs( test[i] - reference[i] ) > eps ) errorWithinBounds = false;
}
REQUIRE( errorWithinBounds );
}
SECTION( "linesource" ) {
......@@ -64,9 +66,11 @@ TEST_CASE( "SN_SOLVER", "[validation_tests]" ) {
double eps = 1e-3;
REQUIRE( test.size() == reference.size() );
bool errorWithinBounds = true;
for( unsigned i = 0; i < test.size(); ++i ) {
REQUIRE( std::fabs( test[i] - reference[i] ) < eps );
if( std::fabs( test[i] - reference[i] ) > eps ) errorWithinBounds = false;
}
REQUIRE( errorWithinBounds );
}
}
......@@ -83,11 +87,13 @@ TEST_CASE( "PN_SOLVER", "[validation_tests]" ) {
auto test = readVTKFile( std::string( TESTS_PATH ) + "result/rtsn_test_checkerboard_PN.vtk" );
auto reference = readVTKFile( std::string( TESTS_PATH ) + pn_fileDir + "checkerboard_PN_reference.vtk" );
double eps = 1e-3;
double eps = 1e-3;
bool errorWithinBounds = true;
REQUIRE( test.size() == reference.size() );
for( unsigned i = 0; i < test.size(); ++i ) {
REQUIRE( std::fabs( test[i] - reference[i] ) < eps );
if( std::fabs( test[i] - reference[i] ) > eps ) errorWithinBounds = false;
}
REQUIRE( errorWithinBounds );
}
SECTION( "linesource" ) {
......@@ -101,12 +107,14 @@ TEST_CASE( "PN_SOLVER", "[validation_tests]" ) {
auto test = readVTKFile( std::string( TESTS_PATH ) + "result/rtsn_test_linesource_PN.vtk" );
auto reference = readVTKFile( std::string( TESTS_PATH ) + pn_fileDir + "linesource_PN_reference.vtk" );
double eps = 1e-3;
double eps = 1e-3;
bool errorWithinBounds = true;
REQUIRE( test.size() == reference.size() );
for( unsigned i = 0; i < test.size(); ++i ) {
REQUIRE( std::fabs( test[i] - reference[i] ) < eps );
if( std::fabs( test[i] - reference[i] ) > eps ) errorWithinBounds = false;
}
REQUIRE( errorWithinBounds );
}
}
......@@ -124,11 +132,14 @@ TEST_CASE( "MN_SOLVER", "[validation_tests]" ) {
auto test = readVTKFile( std::string( TESTS_PATH ) + "result/rtsn_test_checkerboard_MN.vtk" );
auto reference = readVTKFile( std::string( TESTS_PATH ) + mn_fileDir + "checkerboard_MN_reference.vtk" );
double eps = 1e-3;
double eps = 1e-3;
bool errorWithinBounds = true;
REQUIRE( test.size() == reference.size() );
for( unsigned i = 0; i < test.size(); ++i ) {
REQUIRE( std::fabs( test[i] - reference[i] ) < eps );
if( std::fabs( test[i] - reference[i] ) > eps ) errorWithinBounds = false;
}
REQUIRE( errorWithinBounds );
}
SECTION( "linesource" ) {
......@@ -145,12 +156,14 @@ TEST_CASE( "MN_SOLVER", "[validation_tests]" ) {
auto test = readVTKFile( std::string( TESTS_PATH ) + "result/rtsn_test_linesource_MN_Quad.vtk" );
auto reference = readVTKFile( std::string( TESTS_PATH ) + mn_fileDir + "linesource_MN_Quad_reference.vtk" );
double eps = 1e-3;
double eps = 1e-3;
bool errorWithinBounds = true;
REQUIRE( test.size() == reference.size() );
for( unsigned i = 0; i < test.size(); ++i ) {
REQUIRE( std::fabs( test[i] - reference[i] ) < eps );
if( std::fabs( test[i] - reference[i] ) > eps ) errorWithinBounds = false;
}
REQUIRE( errorWithinBounds );
}
{ // --- Maxwell Boltzmann Entropy ---
......@@ -164,11 +177,14 @@ TEST_CASE( "MN_SOLVER", "[validation_tests]" ) {
auto test = readVTKFile( std::string( TESTS_PATH ) + "result/rtsn_test_linesource_MN_MB.vtk" );
auto reference = readVTKFile( std::string( TESTS_PATH ) + mn_fileDir + "linesource_MN_MB_reference.vtk" );
double eps = 1e-3;
double eps = 1e-3;
bool errorWithinBounds = true;
REQUIRE( test.size() == reference.size() );
for( unsigned i = 0; i < test.size(); ++i ) {
REQUIRE( std::fabs( test[i] - reference[i] ) < eps );
if( std::fabs( test[i] - reference[i] ) > eps ) errorWithinBounds = false;
}
REQUIRE( errorWithinBounds );
}
}
}
......
......@@ -42,31 +42,35 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmoni
SECTION( "Test against analytical solution" ) {
std::vector<double> legendre;
Vector moment;
std::vector<bool> validLegendrePoly( 6, true );
std::vector<bool> validMoment( 9, true );
for( double my = -1.0; my < 1.0; my += 0.1 ) {
legendre = testBase.GetAssLegendrePoly( my );
REQUIRE( std::fabs( legendre[0] - P0_0( my ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( legendre[1] - P1_0( my ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( legendre[2] - P1_1( my ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( legendre[3] - P2_0( my ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( legendre[4] - P2_1( my ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( legendre[5] - P2_2( my ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
if( std::fabs( legendre[0] - P0_0( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[0] = false;
if( std::fabs( legendre[1] - P1_0( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[1] = false;
if( std::fabs( legendre[2] - P1_1( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[2] = false;
if( std::fabs( legendre[3] - P2_0( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[3] = false;
if( std::fabs( legendre[4] - P2_1( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[4] = false;
if( std::fabs( legendre[5] - P2_2( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[5] = false;
for( double phi = 0.0; phi < 2 * M_PI; phi += 0.1 ) {
moment = testBase.ComputeSphericalBasis( my, phi );
REQUIRE( std::fabs( moment[0] - Y0_0( my, phi ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( moment[1] - Y1_m1( my, phi ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( moment[2] - Y1_0( my, phi ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( moment[3] - Y1_1( my, phi ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( moment[4] - Y2_m2( my, phi ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( moment[5] - Y2_m1( my, phi ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( moment[6] - Y2_0( my, phi ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( moment[7] - Y2_1( my, phi ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
REQUIRE( std::fabs( moment[8] - Y2_2( my, phi ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
if( std::fabs( moment[0] - Y0_0( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[0] = false;
if( std::fabs( moment[1] - Y1_m1( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[1] = false;
if( std::fabs( moment[2] - Y1_0( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[2] = false;
if( std::fabs( moment[3] - Y1_1( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[3] = false;
if( std::fabs( moment[4] - Y2_m2( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[4] = false;
if( std::fabs( moment[5] - Y2_m1( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[5] = false;
if( std::fabs( moment[6] - Y2_0( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[6] = false;
if( std::fabs( moment[7] - Y2_1( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[7] = false;
if( std::fabs( moment[8] - Y2_2( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[8] = false;
}
}
REQUIRE( std::all_of( validLegendrePoly.begin(), validLegendrePoly.end(), []( bool v ) { return v; } ) );
REQUIRE( std::all_of( validMoment.begin(), validMoment.end(), []( bool v ) { return v; } ) );
}
// Remove title line
......@@ -86,6 +90,7 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmoni
getline( case_file, text_line );
bool errorWithinBounds = true;
while( getline( case_file, text_line ) ) {
// give line to stringstream
......@@ -97,9 +102,10 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmoni
result = testBase.ComputeSphericalBasis( my, phi );
for( unsigned idx = 0; idx < 9; idx++ ) {
REQUIRE( std::fabs( result[idx] - values[idx] ) < 1e2 * std::numeric_limits<double>::epsilon() );
if( std::fabs( result[idx] - values[idx] ) > 1e2 * std::numeric_limits<double>::epsilon() ) errorWithinBounds = false;
}
}
REQUIRE( errorWithinBounds );
case_file.close();
}
......@@ -128,24 +134,22 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmoni
}
}
bool errorWithinBounds = true;
bool orthogonality = true;
for( unsigned idx_row = 0; idx_row < 9; idx_row++ ) {
for( unsigned idx_col = 0; idx_col < 9; idx_col++ ) {
if( idx_row == idx_col ) {
if( std::fabs( results( idx_row, idx_col ) - 1.0 ) > 1e2 * std::numeric_limits<double>::epsilon() ) {
std::cout << "Error in entry (" << idx_row << ", " << idx_col << ") with result." << results( idx_row, idx_col ) << std::endl;
}
// Orthogonality
REQUIRE( std::fabs( results( idx_row, idx_col ) - 1.0 ) < 1e2 * std::numeric_limits<double>::epsilon() );
if( std::fabs( results( idx_row, idx_col ) - 1.0 ) > 1e2 * std::numeric_limits<double>::epsilon() ) orthogonality = false;
}
else {
if( std::fabs( results( idx_row, idx_col ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) {
std::cout << "Error in entry (" << idx_row << ", " << idx_col << ") with result." << results( idx_row, idx_col ) << std::endl;
}
// Normality
REQUIRE( std::fabs( results( idx_row, idx_col ) ) < 1e2 * std::numeric_limits<double>::epsilon() );
if( std::fabs( results( idx_row, idx_col ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) errorWithinBounds = false;
}
}
}
REQUIRE( errorWithinBounds );
REQUIRE( orthogonality );
}
SECTION( "test parity - carthesian coordinates" ) {
......@@ -158,6 +162,7 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmoni
double x, y, z;
bool errorWithinBounds = true;
for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
x = quad.GetPoints()[idx_quad][0];
y = quad.GetPoints()[idx_quad][1];
......@@ -177,10 +182,11 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmoni
else
result = moment2[idx_sys] + moment1[idx_sys];
REQUIRE( std::fabs( result ) < 1e2 * std::numeric_limits<double>::epsilon() );
if( std::fabs( result ) > 1e2 * std::numeric_limits<double>::epsilon() ) errorWithinBounds = false;
}
}
}
REQUIRE( errorWithinBounds );
}
SECTION( "test parity - polar coordinates" ) {
......@@ -192,6 +198,7 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmoni
double result = 0.;
// // test in polar coordinates
bool errorWithinBounds = true;
for( double my = -1.0; my < 1.0; my += 0.1 ) {
for( double phi = 0.0; phi < 2 * M_PI; phi += 0.1 ) {
moment2 = testBase.ComputeSphericalBasis( my, phi );
......@@ -206,10 +213,11 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmoni
else
result = moment2[idx_sys] + moment1[idx_sys];
REQUIRE( std::fabs( result ) < 1e2 * std::numeric_limits<double>::epsilon() );
if( std::fabs( result ) > 1e2 * std::numeric_limits<double>::epsilon() ) errorWithinBounds = false;
}
}
}
}
REQUIRE( errorWithinBounds );
}
}
......@@ -47,11 +47,13 @@ TEST_CASE( "convert image data to grayscale matrix", "[image I/O]" ) {
return x.size() == refMatrix[0].size();
} ) ); // consistency check if all columns of the read-in file have equal length
bool matricesEqual = true;
for( unsigned i = 0; i < gsImage.rows(); ++i ) {
for( unsigned j = 0; j < gsImage.columns(); ++j ) {
REQUIRE( refMatrix[i][j] == gsImage( i, j ) ); // all values match
if( refMatrix[i][j] != gsImage( i, j ) ) matricesEqual = false; // all values match
}
}
REQUIRE( matricesEqual );
}
SECTION( "interpolation of grayscale matrix onto the generated mesh" ) {
......
......@@ -34,9 +34,11 @@ TEST_CASE( "interpolation tests", "[interpolation]" ) {
Interpolation interp( x, y, Interpolation::linear );
Vector res = interp( xq );
REQUIRE( res.size() == ref.size() );
bool errorWithinBounds = true;
for( unsigned i = 0; i < res.size(); ++i ) {
REQUIRE( std::fabs( res[i] - ref[i] ) < 1e-6 );
if( std::fabs( res[i] - ref[i] ) > 1e-6 ) errorWithinBounds = false;
}
REQUIRE( errorWithinBounds );
}
SECTION( "loglinear" ) {
......@@ -47,9 +49,11 @@ TEST_CASE( "interpolation tests", "[interpolation]" ) {
Interpolation interp( x, y, Interpolation::loglinear );
Vector res = interp( xq );
REQUIRE( res.size() == ref.size() );
bool errorWithinBounds = true;
for( unsigned i = 0; i < res.size(); ++i ) {
REQUIRE( std::fabs( res[i] - ref[i] ) < 1e-6 );
if( std::fabs( res[i] - ref[i] ) > 1e-6 ) errorWithinBounds = false;
}
REQUIRE( errorWithinBounds );
}
SECTION( "cubic" ) {
......@@ -60,8 +64,10 @@ TEST_CASE( "interpolation tests", "[interpolation]" ) {
Interpolation interp( x, y, Interpolation::cubic );
Vector res = interp( xq );
REQUIRE( res.size() == ref.size() );
bool errorWithinBounds = true;
for( unsigned i = 0; i < res.size(); ++i ) {
REQUIRE( std::fabs( res[i] - ref[i] ) < 1e-6 );
if( std::fabs( res[i] - ref[i] ) > 1e-6 ) errorWithinBounds = false;
}
REQUIRE( errorWithinBounds );
}
}
......@@ -18,10 +18,13 @@ TEST_CASE( "test all scattering kernels", "[kernel]" ) {
auto weights = quad->GetWeights();
Isotropic kernel( quad );
Matrix scatteringMatrix = kernel.GetScatteringKernel();
bool errorWithinBounds = true;
for( unsigned i = 0; i < scatteringMatrix.rows(); ++i ) {
for( unsigned j = 0; j < scatteringMatrix.columns(); ++j ) {
REQUIRE( std::fabs( scatteringMatrix( i, j ) - ( weights[j] / ( 4 * M_PI ) ) ) < std::numeric_limits<double>::epsilon() );
if( std::fabs( scatteringMatrix( i, j ) - ( weights[j] / ( 4 * M_PI ) ) ) > std::numeric_limits<double>::epsilon() )
errorWithinBounds = false;
}
}
REQUIRE( errorWithinBounds );
}
}
......@@ -21,9 +21,10 @@ TEST_CASE( "unit mesh tests", "[mesh]" ) {
}
SECTION( "neighbor and faces are sorted equally" ) {
auto n = mesh->GetNormals();
auto neighbors = mesh->GetNeighbours();
double eps = 1e-7;
auto n = mesh->GetNormals();
auto neighbors = mesh->GetNeighbours();
double eps = 1e-7;
bool errorWithinBounds = true;
for( unsigned i = 0; i < mesh->GetNumCells(); ++i ) {
for( unsigned j = 0; j < mesh->GetNumNodesPerCell(); ++j ) {
unsigned pos;
......@@ -32,29 +33,35 @@ TEST_CASE( "unit mesh tests", "[mesh]" ) {
for( unsigned k = 0; k < neighbors[nID].size(); ++k ) {
if( neighbors[nID][k] == i ) pos = k;
}
REQUIRE( blaze::l2Norm( n[i][j] + n[nID][pos] ) < eps );
if( blaze::l2Norm( n[i][j] + n[nID][pos] ) > eps ) errorWithinBounds = false;
}
}
REQUIRE( errorWithinBounds );
}
SECTION( "sum over all normals yields zero" ) {
auto n = mesh->GetNormals();
double eps = 1e-7;
auto n = mesh->GetNormals();
double eps = 1e-7;
bool errorWithinBounds = true;
for( unsigned i = 0; i < mesh->GetNumCells(); ++i ) {
Vector sum( 2, 0.0 );
for( unsigned j = 0; j < mesh->GetNumNodesPerCell(); ++j ) {
sum += n[i][j];
}
REQUIRE( blaze::l2Norm( sum ) < eps );
if( blaze::l2Norm( sum ) > eps ) errorWithinBounds = false;
}
REQUIRE( errorWithinBounds );
}
SECTION( "mesh does not have any unassigned faces" ) {
auto neighbors = mesh->GetNeighbours();
auto boundaryType = mesh->GetBoundaryTypes();
auto neighbors = mesh->GetNeighbours();
auto boundaryType = mesh->GetBoundaryTypes();
bool noUnassignedFaces = true;
for( unsigned i = 0; i < mesh->GetNumCells(); ++i ) {
REQUIRE( ( neighbors[i].size() == mesh->GetNumNodesPerCell() ||
( neighbors[i].size() < mesh->GetNumNodesPerCell() && boundaryType[i] != BOUNDARY_TYPE::NONE ) ) );
if( !( neighbors[i].size() == mesh->GetNumNodesPerCell() ||
( neighbors[i].size() < mesh->GetNumNodesPerCell() && boundaryType[i] != BOUNDARY_TYPE::NONE ) ) )
noUnassignedFaces = false;
}
REQUIRE( noUnassignedFaces );
}
}
......@@ -30,6 +30,9 @@ TEST_CASE( "unit numericalflux tests", "[numericalflux]" ) {
Vector psiLVect{ 1.0, 2.0, 3.0, 4.0 };
Vector psiRVect{ 3.0, 4.0, 1.0, 2.0 };
bool fluxCancelingScalar = true;
bool fluxCancelingMatrix = true;
for( unsigned l = 0; l < 10; ++l ) {
Vector n( 2, 0.0 );
std::default_random_engine generator;
......@@ -56,16 +59,18 @@ TEST_CASE( "unit numericalflux tests", "[numericalflux]" ) {
// ---- Test all fluxes ----
// a) Upwinding methods - scalar
REQUIRE( std::fabs( g.Flux( omega, psiL, psiR, n ) - ( -g.Flux( omega, psiR, psiL, -n ) ) ) <
1e2 * std::numeric_limits<double>::epsilon() );
if( std::fabs( g.Flux( omega, psiL, psiR, n ) - ( -g.Flux( omega, psiR, psiL, -n ) ) ) > 1e2 * std::numeric_limits<double>::epsilon() )
fluxCancelingScalar = false;
// b) Upwinding methods - Systems MatrixFlux, 4x4 Matrices, i.e. P1 case;
resultFluxPlus = g.Flux( AxP, AxM, AyP, AyM, AyP, AyM, psiLVect, psiRVect, n );
resultFluxMinus = g.Flux( AxP, AxM, AyP, AyM, AyP, AyM, psiRVect, psiLVect, -n );
REQUIRE( blaze::norm( resultFluxPlus + resultFluxMinus ) < 1e2 * std::numeric_limits<double>::epsilon() );
if( blaze::norm( resultFluxPlus + resultFluxMinus ) > 1e2 * std::numeric_limits<double>::epsilon() ) fluxCancelingMatrix = false;
}
REQUIRE( fluxCancelingScalar );
REQUIRE( fluxCancelingMatrix );
}
SECTION( "test consistency" ) {
......@@ -75,6 +80,9 @@ TEST_CASE( "unit numericalflux tests", "[numericalflux]" ) {
Vector resultFluxOrig( 4, 0 );
Vector psiLVect{ 1.0, 2.0, 3.0, 4.0 };
bool consitencyScalar = true;
bool consitencyMatrix = true;
for( unsigned l = 0; l < 10; ++l ) {
std::default_random_engine generator;
std::normal_distribution<double> distribution( -10.0, 10.0 );
......@@ -102,14 +110,16 @@ TEST_CASE( "unit numericalflux tests", "[numericalflux]" ) {
// a) Upwinding methods - scalar
REQUIRE( std::fabs( g.Flux( omega, psi, psi, n ) - inner * psi ) < 1e2 * std::numeric_limits<double>::epsilon() );
if( std::fabs( g.Flux( omega, psi, psi, n ) - inner * psi ) > 1e2 * std::numeric_limits<double>::epsilon() ) consitencyScalar = false;
// b) Upwinding methods - MatrixFlux, 4x4 Matrices, i.e. P1 case;
resultFluxOrig = ( n[0] * ( AxP + AxM ) * psiLVect + +n[1] * ( AyP + AyM ) * psiLVect );
resultFlux = g.Flux( AxP, AxM, AyP, AyM, AyP, AyM, psiLVect, psiLVect, n );
REQUIRE( blaze::norm( resultFlux - resultFluxOrig ) < 1e2 * std::numeric_limits<double>::epsilon() );
if( blaze::norm( resultFlux - resultFluxOrig ) > 1e2 * std::numeric_limits<double>::epsilon() ) consitencyMatrix = false;
}
REQUIRE( consitencyScalar );
REQUIRE( consitencyMatrix );
}
}
......@@ -50,7 +50,6 @@ TEST_CASE( "checking angular integral of scattering cross sections to be unity",
for( unsigned a = 1; a < angles.size(); ++a ) {
integral_sXS[e] += 0.5 * ( angles[a] - angles[a - 1] ) * ( sXS[e][a] + sXS[e][a - 1] );
}
std::cout << integral_sXS[e] << " " << tXS[e] << std::endl;
REQUIRE( std::fabs( integral_sXS[e] - tXS[e] ) / std::fabs( tXS[e] ) < 0.05 );
}
*/
......
......@@ -114,6 +114,8 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
// Set quadOrder
config->SetQuadOrder( quadratureorder );
bool errorWithinBounds = true;
QuadratureBase* Q = QuadratureBase::CreateQuadrature( config );
VectorVector points = Q->GetPoints();
for( unsigned i = 0; i < Q->GetNq(); i++ ) {
......@@ -126,7 +128,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
lowAccuracyTesting );
printf( "Computed result %.15f", norm( points[i] ) );
}
REQUIRE( approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) );
if( !approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) ) errorWithinBounds = false;
}
// Special case for Gauss Legendre with half weights
......@@ -147,11 +149,12 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
lowAccuracyTesting );
printf( "Computed result %.15f", norm( points[i] ) );
}
REQUIRE( approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) );
if( !approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) ) errorWithinBounds = false;
}
config->SetSNAllGaussPts( true );
}
REQUIRE( errorWithinBounds );
}
}
}
......
Markdown is supported
0% or