Commit 865bc46e authored by steffen.schotthoefer's avatar steffen.schotthoefer
Browse files

harmonics unit test added; PN solver validation test added; deleted safety...

harmonics unit test added; PN solver validation test added; deleted safety outputs for PN Solver; logger does not output stuff on screen (need to fix); Flux Jacobians for MN implemented
parent 08963ca1
Pipeline #94240 failed with stages
in 27 minutes and 54 seconds
......@@ -4,8 +4,6 @@
#include "settings/globalconstants.h"
#include "settings/typedef.h"
#include "toolboxes/errormessages.h"
#include <iostream>
#include <string>
class QuadratureBase
{
......@@ -27,6 +25,12 @@ class QuadratureBase
* @returns double result: result of the quadrature rule */
double Integrate( double( f )( double x0, double x1, double x2 ) );
/*! @brief Integrates vector valued f(x,y,z) with the quadrature. Each dimension is integrated by itself.
* @param : double(f)( double x0, double x1, double x2 ) : density function that depends on a three spatial dimensions.
* @param : len : lenght of vector
* @returns double result: result of the quadrature rule (vector valued) */
std::vector<double> Integrate( std::vector<double>( f )( double x0, double x1, double x2 ), unsigned len );
// Quadrature Hub
/*! @brief Creates a quadrature rule with a given name and a given order.
* @param: std::string name: Name of the quadrature rule
......
#ifndef MNSOLVER_H
#define MNSOLVER_H
#include <pnsolver.h>
#include <sphericalharmonics.h>
#include "solverbase.h"
#include "sphericalharmonics.h"
class MNSolver : Solver
{
......@@ -13,37 +13,24 @@ class MNSolver : Solver
*/
MNSolver( Config* settings );
/*! @brief MNSolver destructor */
~MNSolver();
/**
* @brief Solve functions runs main time loop
*/
void Solve() override;
private:
/*! @brief: Total number of equations in the system */
unsigned _nTotalEntries;
/*! @brief: Max Order of Moments */
unsigned short _nMaxMomentsOrder;
/*! @brief: Absorbtion coefficient for all energies */
VectorVector _sigmaA;
/*! @brief: Class to compute and store current spherical harmonics basis */
SphericalHarmonics _basis;
/*! @brief: Diagonal of the scattering matrix (its a diagonal matrix by construction) */
Vector _scatterMatDiag;
/*! @brief: Diagonals of the system matrices (each sysmatric is a diagonal matrix by construction)
layout: Nx2 (in 2d), N = size of system */
VectorVector _A;
///*! @brief: Diagonal of the system matrix in x direction (its a diagonal matrix by construction) */
// Vector _Ax;
//
///*! @brief: Diagonal of the system matrix in y direction (its a diagonal matrix by construction) */
// Vector _Ay;
//
///*! @brief: Diagonal of the system matrix in z direction (its a diagonal matrix by construction) */
// Vector _Az;
unsigned _nTotalEntries; /*! @brief: Total number of equations in the system */
unsigned short _nMaxMomentsOrder; /*! @brief: Max Order of Moments */
VectorVector _sigmaA; /*! @brief: Absorbtion coefficient for all energies */
SphericalHarmonics _basis; /*! @brief: Class to compute and store current spherical harmonics basis */
Vector _scatterMatDiag; /*! @brief: Diagonal of the scattering matrix (its a diagonal matrix by construction) */
VectorVector _A; /*! @brief: Diagonals of the system matrices (each sysmatric is a diagonal matrix by construction)
layout: Nx2 (in 2d), N = size of system */
QuadratureBase* _quadrature; /*! @brief: Quadrature rule to compute flux Jacobian */
/*! @brief : computes the global index of the moment corresponding to basis function (l,k)
* @param : degree l, it must hold: 0 <= l <=_nq
......
......@@ -30,7 +30,7 @@ class Solver
VectorVector _sigmaT; // total cross section for all energies
std::vector<VectorVector> _Q; // external source term
Matrix _scatteringKernel; // scattering kernel for the quadrature
VectorVector _quadPoints; // quadrature points, dim(_quadPoints) = (_nTimeSteps,spatialDim)
VectorVector _quadPoints; // quadrature points, dim(_quadPoints) = (_nSystem,spatialDim)
Vector _weights; // quadrature weights, dim(_weights) = (_NCells)
std::vector<BOUNDARY_TYPE> _boundaryCells; // boundary type for all cells, dim(_boundary) = (_NCells)
std::vector<double> _solverOutput; // PROTOTYPE: Outputfield for solver
......@@ -55,6 +55,8 @@ class Solver
*/
Solver( Config* settings );
~Solver();
/**
* @brief Create constructor
* @param settings stores all needed information
......
......@@ -30,6 +30,12 @@ class SphericalHarmonics
*/
std::vector<double> ComputeSphericalBasis( double my, double phi );
/*! @brief : Computes all N = L² +2L basis functions at point (x, y, z) on the unit sphere
* @param : x,y,z = coordinates on unit sphere
* @return : vector of basis functions at point (x,y,z) with size N = L² +2L
*/
std::vector<double> ComputeSphericalBasis( double x, double y, double z );
private:
/*! @brief: maximal degree of the spherical harmonics basis (this is "L" in the comments)*/
unsigned _LMaxDegree;
......
......@@ -13,62 +13,20 @@
int main( int argc, char** argv ) {
MPI_Init( &argc, &argv );
// std::string filename = ParseArguments( argc, argv );
//
// // Load Settings from File
// Config* config = new Config( filename );
//
// // Print input file and run info to file
// PrintLogHeader( filename );
//
// // Build solver
// Solver* solver = Solver::Create( config );
//
// // Run solver and export
// solver->Solve();
// solver->Save();
//
std::string filename = ParseArguments( argc, argv );
SphericalHarmonics testBase( 2 );
double d_my = 2. / 50.;
double d_phi = 2. * M_PI / 50.;
std::vector<double> erg( 4, 0.0 );
// Load Settings from File
Config* config = new Config( filename );
std::ofstream outputFile;
std::string filename = "harmonicTest.csv";
outputFile.open( filename );
// Print input file and run info to file
PrintLogHeader( filename );
outputFile << "my"
<< ","
<< "phi"
<< ","
<< "v0"
<< ","
<< "v1"
<< ","
<< "v2"
<< ","
<< "v3"
<< ","
<< "v4"
<< ","
<< "v5"
<< ","
<< "v6"
<< ","
<< "v7"
<< ","
<< "v8" << std::endl;
// Build solver
Solver* solver = Solver::Create( config );
for( int my_idx = -25; my_idx <= 25; my_idx++ ) {
for( int phi_idx = 0; phi_idx <= 50; phi_idx++ ) {
erg = testBase.ComputeSphericalBasis( my_idx * d_my, phi_idx * d_phi );
outputFile << my_idx * d_my << "," << phi_idx * d_phi << "," << erg[0] << "," << erg[1] << "," << erg[2] << "," << erg[3] << "," << erg[4]
<< "," << erg[5] << "," << erg[6] << "," << erg[7] << "," << erg[8] << std::endl;
}
}
outputFile.close();
// Run solver and export
solver->Solve();
solver->Save();
MPI_Finalize();
return EXIT_SUCCESS;
......
......@@ -20,7 +20,7 @@ QuadratureBase* QuadratureBase::CreateQuadrature( QUAD_NAME name, unsigned order
}
double QuadratureBase::Integrate( double( f )( double x0, double x1, double x2 ) ) {
double result = 0;
double result = 0.0;
for( unsigned i = 0; i < _nq; i++ ) {
double x = _points[i][0];
double y = _points[i][1];
......@@ -31,6 +31,23 @@ double QuadratureBase::Integrate( double( f )( double x0, double x1, double x2 )
return result;
}
std::vector<double> QuadratureBase::Integrate( std::vector<double>( f )( double x0, double x1, double x2 ), unsigned len ) {
std::vector<double> result( len, 0.0 );
std::vector<double> funcEval( len, 0.0 );
for( unsigned i = 0; i < _nq; i++ ) {
double x = _points[i][0];
double y = _points[i][1];
double z = _points[i][2];
double w = _weights[i];
funcEval = f( x, y, z );
for( unsigned idx_len = 0; idx_len < len; idx_len++ ) {
result[idx_len] += w * funcEval[idx_len];
}
}
return result;
}
double QuadratureBase::SumUpWeights() { return sum( _weights ); }
void QuadratureBase::PrintWeights() {
......
#include <mnsolver.h>
#include "solvers/mnsolver.h"
MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder( settings->GetMaxMomentOrder() ), _basis( _nMaxMomentsOrder ) {
// Is this good (fast) code using a constructor list?
_nTotalEntries = GlobalIndex( _nMaxMomentsOrder, int( _nMaxMomentsOrder ) ) + 1;
_quadrature = QuadratureBase::CreateQuadrature( settings->GetQuadName(), settings->GetQuadOrder() );
// transform sigmaT and sigmaS in sigmaA.
_sigmaA = VectorVector( _nEnergies, Vector( _nCells, 0 ) ); // Get rid of this extra vektor!
......@@ -26,6 +26,8 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder(
ComputeSystemMatrices();
}
MNSolver::~MNSolver() { delete _quadrature; }
int MNSolver::GlobalIndex( int l, int k ) const {
int numIndicesPrevLevel = l * l; // number of previous indices untill level l-1
int prevIndicesThisLevel = k + l; // number of previous indices in current level
......@@ -33,13 +35,28 @@ int MNSolver::GlobalIndex( int l, int k ) const {
}
void MNSolver::ComputeSystemMatrices() {
for( int l_idx = 0; l_idx <= _nMaxMomentsOrder; l_idx++ ) {
for( unsigned k_idx = -l_idx; k_idx <= l_idx; k_idx++ ) {
_A[idx_sys] = Vector( 2, 0 );
// compute entries. =<Omega_i*m^(l,k)>
// Initialize entries of _A
for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
Vector entry( 2, 0.0 );
_A[idx_sys] = entry;
}
// Compute Integrals _A[idx_system][i]=<Omega_i*m^(l,k)>
std::vector<double> funcEval( _nTotalEntries, 0.0 );
for( unsigned i = 0; i < _nq; i++ ) {
double x = _quadrature->GetPoints()[i][0];
double y = _quadrature->GetPoints()[i][1];
double z = _quadrature->GetPoints()[i][2];
double w = _quadrature->GetWeights()[i];
funcEval = _basis.ComputeSphericalBasis( x, y, z ); // = m_l^k
for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
_A[idx_system][0] += w * x * funcEval[idx_system];
_A[idx_system][1] += w * y * funcEval[idx_system];
}
}
}
void MNSolver::Solve() {
int rank;
......@@ -85,12 +102,10 @@ void MNSolver::Solve() {
// store flux contribution on psiNew_sigmaS to save memory
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][l] == _nCells )
psiNew[idx_cell][idx_system] +=
_g->Flux( _quadPoints[idx_system], _psi[idx_cell][idx_system], _psi[idx_cell][idx_system], _normals[idx_cell][l] );
_g->Flux( _A[idx_system], _psi[idx_cell][idx_system], _psi[idx_cell][idx_system], _normals[idx_cell][l] );
else
psiNew[idx_cell][idx_system] += _g->Flux( _quadPoints[idx_system],
_psi[idx_cell][idx_system],
_psi[_neighbors[idx_cell][l]][idx_system],
_normals[idx_cell][l] );
psiNew[idx_cell][idx_system] += _g->Flux(
_A[idx_system], _psi[idx_cell][idx_system], _psi[_neighbors[idx_cell][l]][idx_system], _normals[idx_cell][l] );
}
// time update angular flux with numerical flux and total scattering cross section
psiNew[idx_cell][idx_system] = _psi[idx_cell][idx_system] -
......@@ -103,27 +118,29 @@ void MNSolver::Solve() {
// TODO: figure out a more elegant way
// add external source contribution
if( _Q.size() == 1u ) { // constant source for all energies
if( _Q[0][idx_cell].size() == 1u ) // isotropic source
psiNew[idx_cell] += _dE * _Q[0][idx_cell][0];
else
psiNew[idx_cell] += _dE * _Q[0][idx_cell];
}
else {
if( _Q[0][idx_cell].size() == 1u ) // isotropic source
psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell][0];
else
psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell];
}
// if( _Q.size() == 1u ) { // constant source for all energies
// if( _Q[0][idx_cell].size() == 1u ) // isotropic source
// psiNew[idx_cell] += _dE * _Q[0][idx_cell][0];
// else
// psiNew[idx_cell] += _dE * _Q[0][idx_cell];
//}
// else {
// if( _Q[0][idx_cell].size() == 1u ) // isotropic source
// psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell][0];
// else
// psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell];
//}
}
_psi = psiNew;
for( unsigned i = 0; i < _nCells; ++i ) {
fluxNew[i] = dot( _psi[i], _weights );
_solverOutput[i] = fluxNew[i];
_psi = psiNew;
double mass = 0.0;
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
fluxNew[idx_cell] = _psi[idx_cell][0]; // zeroth moment is raditation densitiy we are interested in
_solverOutput[idx_cell] = _psi[idx_cell][0];
mass += _psi[idx_cell][0] * _areas[idx_cell];
}
// Save( n );
dFlux = blaze::l2Norm( fluxNew - fluxOld );
fluxOld = fluxNew;
if( rank == 0 ) log->info( "{:03.8f} {:01.5e}", _energies[n], dFlux );
if( rank == 0 ) log->info( "{:03.8f} {:01.5e}", _energies[idx_energy], dFlux );
}
}
......@@ -38,33 +38,33 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
// Fill System Matrices
ComputeSystemMatrices();
std::cout << "System Matrix Set UP!" << std::endl;
// std::cout << "System Matrix Set UP!" << std::endl;
// Compute Decomposition in positive and negative (eigenvalue) parts of flux jacobians
ComputeFluxComponents();
// Compute diagonal of the scatter matrix (it's a diagonal matrix)
ComputeScatterMatrix();
std::cout << "scatterMatrix : " << _scatterMatDiag << "\n";
// std::cout << "scatterMatrix : " << _scatterMatDiag << "\n";
// AdaptTimeStep();
if( settings->GetCleanFluxMat() ) CleanFluxMatrices();
std::cout << "--------\n";
std::cout << "_Ax :\n" << _Ax << "\n "; // _AxP \n" << _AxPlus << "\n _AxM \n" << _AxMinus << "\n";
std::cout << "_Ay :\n" << _Ay << "\n "; //_AyP \n" << _AyPlus << "\n _AyM \n" << _AyMinus << "\n";
std::cout << "_Az :\n" << _Az << "\n "; //_AzP \n" << _AzPlus << "\n _AzM \n" << _AzMinus << "\n";
std::cout << "_AxPlus :\n" << _AxPlus << "\n "; // _AxP \n" << _AxPlus << "\n _AxM \n" << _AxMinus << "\n";
std::cout << "_AyPlus :\n" << _AyPlus << "\n "; //_AyP \n" << _AyPlus << "\n _AyM \n" << _AyMinus << "\n";
std::cout << "_AzPlus :\n" << _AzPlus << "\n ";
std::cout << "_AxMinus :\n" << _AxMinus << "\n "; // _AxP \n" << _AxPlus << "\n _AxM \n" << _AxMinus << "\n";
std::cout << "_AyMinus :\n" << _AyMinus << "\n "; //_AyP \n" << _AyPlus << "\n _AyM \n" << _AyMinus << "\n";
std::cout << "_AzMinus :\n" << _AzMinus << "\n ";
std::cout << "_nCells: " << _nCells << "\n";
// std::cout << "--------\n";
// std::cout << "_Ax :\n" << _Ax << "\n "; // _AxP \n" << _AxPlus << "\n _AxM \n" << _AxMinus << "\n";
// std::cout << "_Ay :\n" << _Ay << "\n "; //_AyP \n" << _AyPlus << "\n _AyM \n" << _AyMinus << "\n";
// std::cout << "_Az :\n" << _Az << "\n "; //_AzP \n" << _AzPlus << "\n _AzM \n" << _AzMinus << "\n";
//
// std::cout << "_AxPlus :\n" << _AxPlus << "\n "; // _AxP \n" << _AxPlus << "\n _AxM \n" << _AxMinus << "\n";
// std::cout << "_AyPlus :\n" << _AyPlus << "\n "; //_AyP \n" << _AyPlus << "\n _AyM \n" << _AyMinus << "\n";
// std::cout << "_AzPlus :\n" << _AzPlus << "\n ";
//
// std::cout << "_AxMinus :\n" << _AxMinus << "\n "; // _AxP \n" << _AxPlus << "\n _AxM \n" << _AxMinus << "\n";
// std::cout << "_AyMinus :\n" << _AyMinus << "\n "; //_AyP \n" << _AyPlus << "\n _AyM \n" << _AyMinus << "\n";
// std::cout << "_AzMinus :\n" << _AzMinus << "\n ";
//
// std::cout << "_nCells: " << _nCells << "\n";
}
void PNSolver::Solve() {
......@@ -151,13 +151,15 @@ void PNSolver::Solve() {
for( unsigned i = 0; i < _nCells; ++i ) {
fluxNew[i] = _psi[i][0]; // zeroth moment is raditation densitiy we are interested in
_solverOutput[i] = _psi[i][0];
mass += _psi[i][0];
mass += _psi[i][0] * _areas[i];
}
dFlux = blaze::l2Norm( fluxNew - fluxOld );
fluxOld = fluxNew;
if( rank == 0 ) log->info( "{:03.8f} {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
if( rank == 0 ) {
// std::cout << "PseudoTime: " << _energies[idx_energy] << " | Flux Change: " << dFlux << " | Mass: " << mass << "\n";
log->info( "{:03.8f} {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
}
Save( idx_energy );
}
}
......@@ -322,16 +324,16 @@ void PNSolver::ComputeFluxComponents() {
}
// Compute Spectral Radius
std::cout << "Eigenvalues x direction " << eigenValuesX << "\n";
std::cout << "Eigenvalues y direction " << eigenValuesY << "\n";
std::cout << "Eigenvalues z direction " << eigenValues << "\n";
std::cout << "Spectral Radius X " << blaze::max( blaze::abs( eigenValuesX ) ) << "\n";
std::cout << "Spectral Radius Y " << blaze::max( blaze::abs( eigenValuesY ) ) << "\n";
std::cout << "Spectral Radius Z " << blaze::max( blaze::abs( eigenValues ) ) << "\n";
// std::cout << "Eigenvalues x direction " << eigenValuesX << "\n";
// std::cout << "Eigenvalues y direction " << eigenValuesY << "\n";
// std::cout << "Eigenvalues z direction " << eigenValues << "\n";
//
// std::cout << "Spectral Radius X " << blaze::max( blaze::abs( eigenValuesX ) ) << "\n";
// std::cout << "Spectral Radius Y " << blaze::max( blaze::abs( eigenValuesY ) ) << "\n";
// std::cout << "Spectral Radius Z " << blaze::max( blaze::abs( eigenValues ) ) << "\n";
_combinedSpectralRadius = blaze::max( blaze::abs( eigenValues + eigenValuesX + eigenValuesY ) );
std::cout << "Spectral Radius combined " << _combinedSpectralRadius << "\n";
// std::cout << "Spectral Radius combined " << _combinedSpectralRadius << "\n";
}
void PNSolver::ComputeScatterMatrix() {
......@@ -372,7 +374,7 @@ void PNSolver::Save() const {
}
std::vector<std::vector<double>> scalarField( 1, flux );
std::vector<std::vector<std::vector<double>>> results{ scalarField };
ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( _nEnergies ), results, fieldNames, _mesh );
ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
}
void PNSolver::Save( int currEnergy ) const {
......
......@@ -25,6 +25,8 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
_settings->SetNQuadPoints( _nq );
ScatteringKernel* k = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), quad );
_scatteringKernel = k->GetScatteringKernel();
// delete QuadratureBase
delete quad;
// set time step
_dE = ComputeTimeStep( settings->GetCFL() );
......@@ -49,6 +51,11 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
_solverOutput.resize( _nCells ); // Currently only Flux
}
Solver::~Solver() {
delete _mesh;
delete _problem;
}
double Solver::ComputeTimeStep( double cfl ) const {
double maxEdge = 1000.;
for( unsigned j = 0; j < _nCells; ++j ) {
......
......@@ -22,9 +22,23 @@ std::vector<double> SphericalHarmonics::ComputeSphericalBasis( double my, double
return _YBasis;
}
std::vector<double> SphericalHarmonics::ComputeSphericalBasis( double x, double y, double z ) {
// transform (x,y,z) into (my,phi)
double my = z;
double phi = 0.0;
if( y >= 0 )
phi = acos( x );
else
phi = acos( -x ) + M_PI;
ComputeAssLegendrePoly( my );
ComputeYBasis( phi );
return _YBasis;
}
void SphericalHarmonics::ComputeCoefficients() {
// m in paper is here denoted by k
double ls = 0.0; // l^2
double lm1s = 0.0; // (l-1)^2
double ks = 0.0; // k^2
......@@ -32,8 +46,7 @@ void SphericalHarmonics::ComputeCoefficients() {
ls = l_idx * l_idx;
lm1s = ( l_idx - 1 ) * ( l_idx - 1 );
for( unsigned k_idx = 0; k_idx < l_idx - 1; k_idx++ ) {
ks = k_idx * k_idx;
ks = k_idx * k_idx;
_aParam[GlobalIdxAssLegendreP( l_idx, k_idx )] = std::sqrt( ( 4 * ls - 1. ) / ( ls - ks ) );
_bParam[GlobalIdxAssLegendreP( l_idx, k_idx )] = -std::sqrt( ( lm1s - ks ) / ( 4 * lm1s - 1. ) );
}
......
This source diff could not be displayed because it is too large. You can view the blob instead.
% ---- File specifications ----
OUTPUT_DIR = ../../result
OUTPUT_FILE = rtsn_test_linesource_PN
LOG_DIR = ../../result/logs
MESH_FILE = linesource.su2
% ---- Solver specifications ----
CFL_NUMBER = 0.5
TIME_FINAL = 0.5
PROBLEM = LINESOURCE
SOLVER = PN_SOLVER
% ---- Boundary Conditions ----
BC_NEUMANN = ( void )
QUAD_ORDER = 2
This diff is collapsed.
......@@ -63,20 +63,20 @@ TEST_CASE( "linesource_SN", "[testcases]" ) {
}
}
// TEST_CASE( "linesource_PN", "[testcases]" ) {
// char config_file_name[MAX_STRING_SIZE] = "../tests/input/linesource_PN.cfg";
//
// Config* config = new Config( config_file_name );
// Solver* solver = Solver::Create( config );
// solver->Solve();
// solver->Save();
//
// auto test = readVTKFile( "../result/rtsn_test_linesource_PN.vtk" );
// auto reference = readVTKFile( "../tests/input/linesource_reference_PN.vtk" );
//
// double eps = 1e-3;
// REQUIRE( test.size() == reference.size() );
// for( unsigned i = 0; i < test.size(); ++i ) {
// REQUIRE( std::fabs( test[i] - reference[i] ) < eps );
// }
//}
TEST_CASE( "linesource_PN", "[testcases]" ) {
char config_file_name[MAX_STRING_SIZE] = "../tests/input/linesource_PN.cfg";
Config* config = new Config( config_file_name );
Solver* solver = Solver::Create( config );
solver->Solve();
solver->Save();
auto test = readVTKFile( "../result/rtsn_test_linesource_PN.vtk" );
auto reference = readVTKFile( "../tests/input/rtsn_test_linesource_PN_reference.vtk" );
double eps = 1e-3;
REQUIRE( test.size() == reference.size() );
for( unsigned i = 0; i < test.size(); ++i ) {
REQUIRE( std::fabs( test[i] - reference[i] ) < eps );
}
}
#include "catch.hpp"
#include "solvers/sphericalharmonics.h"
#include <fstream>
#include <sstream>
TEST_CASE( "test the spherical harmonics basis computation", "[spherical harmonics]" ) {
std::string text_line;
std::ifstream case_file;
SphericalHarmonics testBase( 2 );
double my = 0.0;
double phi = 0.0;
std::vector<double> values( 9, 0.0 );
std::vector<double> result( 4, 0.0 );
case_file.open( "harmonicBasis_reference.csv", std::ios::in );
// Remove title line
getline( case_file, text_line );
while( getline( case_file, text_line ) ) {
// give line to stringstream
std::stringstream ss( text_line );
// Read values
ss >> my >> phi >> values[0] >> values[1] >> values[2] >> values[3] >> values[4] >> values[5] >> values[6] >> values[7] >> values[8];
result = testBase.ComputeSphericalBasis( my, phi );
for( unsigned idx = 0; idx < 9; idx++ ) {
REQUIRE( std::fabs( result[idx] - values[idx] ) < std::numeric_limits<double>::epsilon() );
}
}
case_file.close();
}