Commit e3bdadf1 authored by steffen.schotthoefer's avatar steffen.schotthoefer
Browse files

fix #57 and #58. GLquad is now correct and tested for SphericalHarmonics....

fix #57 and #58. GLquad is now correct and tested for SphericalHarmonics. SphericalHarmonics normality unit test implemented. MN_Solver works (unvalidated) (solves #8). Changed tolarance for quadrature rules, due to weird behaviour of MC quadrature. Its now e-12 instead of e-15
parent 37925060
Pipeline #96611 failed with stages
in 6 minutes and 44 seconds
......@@ -5,6 +5,7 @@
class QGaussLegendreTensorized : public QuadratureBase
{
// Implementation is done accordingly to Kendall Atkinson 1981, Australian Matematical Society.
private:
double Pythag( const double a, const double b );
std::pair<Vector, Matrix> ComputeEigenValTriDiagMatrix( const Matrix& mat );
......
......@@ -22,53 +22,53 @@ double testFunc2( double x, double y, double z ) { return x + y + z; }
int main( int argc, char** argv ) {
MPI_Init( &argc, &argv );
QMonteCarlo quad( 1000 );
int maxMomentDegree = 2;
// test quadrature;
std::cout << " Integration test:" << quad.Integrate( testFunc2 ) << "\n";
SphericalHarmonics testBase( maxMomentDegree );
double my, phi, w, resTest = 0.0;
Vector moment = testBase.ComputeSphericalBasis( 0, 1, 0 );
// 9 basis moments if degree = 2
// unsigned nTotalMoments = moment.size();
// QGaussLegendreTensorized quad( 8 );
// int maxMomentDegree = 8;
//
Vector results( moment.size(), 0.0 );
int idx_basis = 0, idx_basisM = 0;
for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
my = quad.GetPointsSphere()[idx_quad][0];
phi = quad.GetPointsSphere()[idx_quad][1];
w = quad.GetWeights()[idx_quad];
moment = testBase.ComputeSphericalBasis( my, phi );
double funcEval = testFunc( my, phi );
// std::cout << "f(my= " << my << ", phi= " << phi / M_PI << "*pi) = " << funcEval << " w " << w << "\n";
resTest += w * funcEval;
std::cout << "my, phi " << my << " , " << phi << " ::\t mom5 " << moment[5] << " mom7 \t" << moment[7] << "\n";
// results[5] += w * moment[5] * moment[5];
// results[7] += w * moment[7] * moment[7];
//
for( int idx_l = 0; idx_l <= maxMomentDegree; idx_l++ ) {
idx_basis = testBase.GlobalIdxBasis( idx_l, 0 );
results[idx_basis] += w * moment[idx_basis] * moment[idx_basis];
for( int idx_k = 1; idx_k <= idx_l; idx_k++ ) {
idx_basis = testBase.GlobalIdxBasis( idx_l, idx_k );
idx_basisM = testBase.GlobalIdxBasis( idx_l, -idx_k );
results[idx_basis] += w * moment[idx_basis] * moment[idx_basis];
results[idx_basisM] += w * moment[idx_basisM] * moment[idx_basisM];
}
}
}
std::cout << " testInt " << resTest << "\n";
std::cout << " Squared Integration:\n " << results << "\n";
//// test quadrature;
// std::cout << " Integration test:" << quad.Integrate( testFunc2 ) << "\n";
//
// SphericalHarmonics testBase( maxMomentDegree );
//
// double my, phi, w, resTest = 0.0;
// Vector moment = testBase.ComputeSphericalBasis( 0, 1, 0 );
//// 9 basis moments if degree = 2
//// unsigned nTotalMoments = moment.size();
////
// Vector results( moment.size(), 0.0 );
//
// int idx_basis = 0, idx_basisM = 0;
//
// for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
// my = quad.GetPointsSphere()[idx_quad][0];
// phi = quad.GetPointsSphere()[idx_quad][1];
// w = quad.GetWeights()[idx_quad];
// moment = testBase.ComputeSphericalBasis( my, phi );
//
// double funcEval = testFunc( my, phi );
// std::cout << "f(my= " << my << ", phi= " << phi / M_PI << "*pi) = " << funcEval << " w " << w << "\n";
// resTest += w * funcEval;
//
// // std::cout << "my, phi " << my << " , " << phi << " ::\t mom5 " << moment[5] << " mom7 \t" << moment[7] << "\n";
// // results[5] += w * moment[5] * moment[5];
// // results[7] += w * moment[7] * moment[7];
// //
// for( int idx_l = 0; idx_l <= maxMomentDegree; idx_l++ ) {
// idx_basis = testBase.GlobalIdxBasis( idx_l, 0 );
//
// results[idx_basis] += w * moment[idx_basis] * moment[idx_basis];
//
// for( int idx_k = 1; idx_k <= idx_l; idx_k++ ) {
// idx_basis = testBase.GlobalIdxBasis( idx_l, idx_k );
// idx_basisM = testBase.GlobalIdxBasis( idx_l, -idx_k );
//
// results[idx_basis] += w * moment[idx_basis] * moment[idx_basis];
// results[idx_basisM] += w * moment[idx_basisM] * moment[idx_basisM];
// }
// }
//}
// std::cout << " testInt " << resTest << "\n";
// std::cout << " Squared Integration:\n " << results << "\n";
//
// Normalized integration
// Vector resultsNormal( moment.size(), 0.0 );
......@@ -122,20 +122,20 @@ int main( int argc, char** argv ) {
// myfile.close();
//
// std::string filename = ParseArguments( argc, argv );
//
// // CD 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 );
// CD 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();
MPI_Finalize();
return EXIT_SUCCESS;
......
......@@ -79,6 +79,7 @@ VectorVector LineSource_PN::SetupIC() {
psi[j][0] = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) );
}
// Debugging jump test case
// for( unsigned j = 0; j < cellMids.size(); ++j ) {
// if( cellMids[j][0] > -0.2 && cellMids[j][1] > -0.2 && cellMids[j][1] < 0.2 && cellMids[j][0] < 0.2 )
// psi[j][0] = 1.0;
......
......@@ -28,9 +28,6 @@ void QGaussLegendreTensorized::SetPointsAndWeights() {
nodes1D[i] = evSys.first[i];
weights1D[i] = 2 * std::pow( evSys.second( 0, i ), 2 );
}
for( unsigned i = 0; i < _order; ++i ) {
nodes1D[i] = ( nodes1D[i] + 1.0 ) * 0.5;
}
// sort nodes increasingly and also reorder weigths for consistency
std::vector<unsigned> sortOrder( nodes1D.size() );
......@@ -48,8 +45,8 @@ void QGaussLegendreTensorized::SetPointsAndWeights() {
phi[i] = ( i + 0.5 ) * M_PI / _order;
}
// only use the spheres upper half
unsigned range = std::floor( _order / 2.0 );
unsigned range = std::floor( _order / 2.0 ); // comment (steffen): why do we only need half of the points: => In 2D we would count everything
// twice. (not wrong with scaling, but expensive)
// resize points and weights
_points.resize( _nq );
......@@ -71,7 +68,7 @@ void QGaussLegendreTensorized::SetPointsAndWeights() {
_points[j * ( 2 * _order ) + i][2] = nodes1D[j];
_pointsSphere[j * ( 2 * _order ) + i][0] = nodes1D[j];
_pointsSphere[j * ( 2 * _order ) + i][0] = phi[i];
_pointsSphere[j * ( 2 * _order ) + i][1] = phi[i];
_weights[j * ( 2 * _order ) + i] = 2.0 * M_PI / _order * weights1D[j];
}
......
......@@ -76,4 +76,6 @@ void QLookupQuadrature::SetPointsAndWeights() {
_weights[idx] /= sumWeights;
_weights[idx] *= 4.0 * PI_NUMBER;
}
// Transform _points to _pointsSphere
}
......@@ -53,30 +53,11 @@ void MNSolver::ComputeMoments() {
double my, phi, w;
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
phi = _quadrature->GetPointsSphere()[idx_quad][0];
my = _quadrature->GetPointsSphere()[idx_quad][1];
my = _quadrature->GetPointsSphere()[idx_quad][0];
phi = _quadrature->GetPointsSphere()[idx_quad][1];
_moments[idx_quad] = _basis.ComputeSphericalBasis( my, phi );
}
Vector normalizationFactor( _moments[0].size(), 0.0 );
// Normalize Moments
// for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
// w = _quadrature->GetWeights()[idx_quad];
//
// for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
// normalizationFactor[idx_sys] += w * _moments[idx_quad][idx_sys] * _moments[idx_quad][idx_sys];
// }
//}
// std::cout << "norm" << normalizationFactor << "\n";
//
// for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
// for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
// _moments[idx_quad][idx_sys] /= sqrt( normalizationFactor[idx_sys] );
// }
//}
}
Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
......@@ -87,18 +68,20 @@ Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
// std::cout << "--------\n";
// std::cout << " alphas curr cell" << _alpha[idx_cell] << " alphas neigh cell:\n";
bool good = _alpha[idx_cell] == 0.0;
// bool good = _alpha[_neighbors[idx_cell][idx_neigh]] == 0.0;
int idx_good = 0;
for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
// std::cout << _alpha[idx_neigh] << "\n";
if( _alpha[idx_neigh] == 0.0 ) idx_good++;
}
if( idx_good == 2 && good ) {
std::cout << "Candidate\n";
}
// int idx_good = 0;
// for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
// std::cout << _alpha[_neighbors[idx_cell][idx_neigh]] << "\n";
// if( _alpha[_neighbors[idx_cell][idx_neigh]] == 0.0 ) idx_good++;
//}
// if( idx_good == 2 && good ) {
// std::cout << "Candidate\n";
//}
Vector flux( _nTotalEntries, 0.0 );
Vector omega( 3, 0.0 );
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
w = _quadrature->GetWeights()[idx_quad];
......@@ -107,24 +90,29 @@ Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
entropyL = dot( _alpha[idx_cell], _moments[idx_quad] );
// _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _moments[idx_quad] ) );
omega = _quadrature->GetPoints()[idx_quad];
// std::cout << "omega " << omega << " ||omega|| " << omega[0] * omega[0] + omega[1] * omega[1] + omega[2] * omega[2] << " entropyL " <<
// entropyL
// << "\n";
for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
// Store fluxes in psiNew, to save memory
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neigh] == _nCells )
entropyR = entropyL;
else {
entropyR = dot( _alpha[idx_neigh], _moments[idx_quad] );
entropyR = dot( _alpha[_neighbors[idx_cell][idx_neigh]], _moments[idx_quad] );
//_entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_neigh], _moments[idx_quad] ) );
}
entropyFlux += _g->Flux( _quadrature->GetPoints()[idx_quad], entropyL, entropyR, _normals[idx_cell][idx_neigh] );
// std::cout << " entropyR " << entropyR << " flux " << _g->Flux( omega, entropyL, entropyR, _normals[idx_cell][idx_neigh] ) << " \n";
entropyFlux += _g->Flux( omega, entropyL, entropyR, _normals[idx_cell][idx_neigh] );
}
// std::cout << "\n entropy Flux at cell [" << idx_cell << "] and quad [" << idx_quad << "]: " << entropyFlux << "\n";
// std::cout << "\n entropy Flux at cell [" << idx_cell << "] and quad [" << idx_quad << "]: " << entropyFlux << " with moment\n"
// << _moments[idx_quad] << "\n";
// integrate
flux += _moments[idx_quad] * ( w * entropyFlux );
}
// std::cout << "\n integrate entropy Flux at cell [" << idx_cell << "] and : " << flux << "\n";
// std::cout << std::endl;
// std::cout << " integrate entropy Flux at cell [" << idx_cell << "] : " << flux << "\n";
return flux;
}
......
......@@ -122,31 +122,31 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical harmoni
}
}
// SECTION( "test normality" ) {
//
// QGaussLegendreTensorized quad( 6 );
//
// double x, y, z, w;
// Vector moment = testBase.ComputeSphericalBasis( 0, 1, 0 );
// // 9 basis moments if degree = 2
//
// Vector results( moment.size(), 0.0 );
//
// for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
// x = quad.GetPoints()[idx_quad][0];
// y = quad.GetPoints()[idx_quad][1];
// z = quad.GetPoints()[idx_quad][2];
// w = quad.GetWeights()[idx_quad];
// moment = testBase.ComputeSphericalBasis( x, y, z );
//
// for( unsigned idx_sys = 0; idx_sys < 9; idx_sys++ ) {
// results[idx_sys] += w * moment[idx_sys] * moment[idx_sys];
// }
// }
// for( unsigned idx_sys = 0; idx_sys < 9; idx_sys++ ) {
// REQUIRE( std::fabs( results[idx_sys] - 1 ) < std::numeric_limits<double>::epsilon() );
// }
//}
SECTION( "test normality" ) {
QGaussLegendreTensorized quad( 6 );
double my, phi, w;
Vector moment = testBase.ComputeSphericalBasis( 0, 1, 0 );
// 9 basis moments if degree = 2
Vector results( moment.size(), 0.0 );
for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
my = quad.GetPointsSphere()[idx_quad][0];
phi = quad.GetPointsSphere()[idx_quad][1];
w = quad.GetWeights()[idx_quad];
moment = testBase.ComputeSphericalBasis( my, phi );
for( unsigned idx_sys = 0; idx_sys < maxMomentDegree; idx_sys++ ) {
results[idx_sys] += w * moment[idx_sys] * moment[idx_sys];
}
}
for( unsigned idx_sys = 0; idx_sys < maxMomentDegree; idx_sys++ ) {
REQUIRE( std::fabs( results[idx_sys] - 1 ) < 1e2 * std::numeric_limits<double>::epsilon() );
}
}
SECTION( "test parity - carthesian coordinates" ) {
......
......@@ -15,7 +15,7 @@ std::vector<std::vector<int>> quadratureorders = {
};
bool approxequal( double a, double b, bool lowAccuracy = false ) {
double tol = 1e-15; // For computed quadrature weights
double tol = 1e-12; // For computed quadrature weights
if( lowAccuracy ) tol = 5e-5; // Mainly for Lookup Quadratures
return std::abs( a - b ) < tol;
}
......@@ -36,7 +36,7 @@ TEST_CASE( "Quadrature weights sum to 4*pi.", "[quadrature]" ) {
quadratureorder,
std::abs( Q->SumUpWeights() - 4 * M_PI ),
lowAccuracyTesting );
printf( "Computed result %.15f", Q->SumUpWeights() );
printf( "Computed result %.15f \n", Q->SumUpWeights() );
}
REQUIRE( approxequal( Q->SumUpWeights(), 4 * M_PI, lowAccuracyTesting ) );
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment