Skip to content

GitLab

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

pn solver test repaired, on the hunt for the mn flux bug, spherical basis under construction

parent 3b14b274
Expand all Hide whitespace changes
Inline Side-by-side
code/include/entropies/entropybase.h
View file @ a1bc511e
......@@ -36,6 +36,11 @@ class EntropyBase
* @returns: value of dual of the derivative at z */
virtual double EntropyPrimeDual( double y ) = 0;
/*! @brief: computes the hessian of the dual entropy functional
* @param: y = point where the hessian should be evaluated;
* @returns: value of the hessian at y */
// virtual double EntropyHessianDual( double y ) = 0;
/*! @brief: checks, if value is in domain of entropy */
virtual bool CheckDomain( double z ) = 0;
};
......
code/include/entropies/quadraticentropy.h
View file @ a1bc511e
......@@ -19,6 +19,8 @@ class QuadraticEntropy : public EntropyBase
inline double EntropyPrimeDual( double y ) override { return y; }
// inline double EntropyHessianDual( double y ) override { return y; }
inline bool CheckDomain( double z ) override { return std::isnan( z ); }
};
......
code/include/solvers/mnsolver.h
View file @ a1bc511e
......@@ -29,13 +29,14 @@ class MNSolver : public Solver
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 */
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 */
VectorVector _moments; /*! @brief: Moment Vector pre-computed at each quadrature point: dim= _nq x _nTotalEntries */
Vector _scatterMatDiag; /*! @brief: Diagonal of the scattering matrix (its a diagonal matrix by construction) */
QuadratureBase* _quadrature; /*! @brief: Quadrature rule to compute flux Jacobian */
QuadratureBase* _quadrature; /*! @brief: Quadrature rule to compute flux Jacobian */ // Shift this to SolverBase!
EntropyBase* _entropy; /*! @brief: Class to handle entropy functionals */
VectorVector _alpha; /*! @brief: Lagrange Multipliers for Minimal Entropy problem for each gridCell
......@@ -49,14 +50,13 @@ class MNSolver : public Solver
*/
int GlobalIndex( int l, int k ) const;
// /*! @brief : function for computing and setting up the System Matrices.
// The System Matrices are diagonal, so there is no need to diagonalize*/
// void ComputeSystemMatrices();
/*! @brief : Construct flux by computing the Moment of the FVM discretization at the interface of cell
* @param : idx_cell = current cell id
* @param : idx_neighbor = neighbor cell id
* @returns : flux for all moments at interface of idx_cell,idx_neighbor */
Vector ConstructFlux( unsigned idx_cell, unsigned idx_neighbor );
/*! @brief : Pre-Compute Moments at all quadrature points. */
void ComputeMoments();
};
#endif // MNSOLVER_H
code/python/harmonicTest.csv
View file @ a1bc511e
This diff is collapsed.
code/python/harmonicTest2.csv deleted 100644 → 0
View file @ 3b14b274
This diff is collapsed.
code/python/plotBasis.py
View file @ a1bc511e
......@@ -19,6 +19,18 @@ X = np.zeros((ndimGrid, ndimGrid))
Y = np.zeros((ndimGrid, ndimGrid))
Z = np.zeros((ndimGrid, ndimGrid))
X0 = np.zeros((ndimGrid, ndimGrid))
Y0 = np.zeros((ndimGrid, ndimGrid))
Z0 = np.zeros((ndimGrid, ndimGrid))
X1 = np.zeros((ndimGrid, ndimGrid))
Y1 = np.zeros((ndimGrid, ndimGrid))
Z1 = np.zeros((ndimGrid, ndimGrid))
X2 = np.zeros((ndimGrid, ndimGrid))
Y2 = np.zeros((ndimGrid, ndimGrid))
Z2 = np.zeros((ndimGrid, ndimGrid))
Harm0 = np.zeros((ndimGrid, ndimGrid))
Harm1 = np.zeros((ndimGrid, ndimGrid))
Harm2 = np.zeros((ndimGrid, ndimGrid))
......@@ -65,6 +77,19 @@ with open('harmonicTest.csv') as csv_file:
Harm8[idx_row, idx_col] = row[10]
line_count += 1
# radius version
X0[idx_row, idx_col] = X[idx_row, idx_col] * Harm0[idx_row, idx_col]
Y0[idx_row, idx_col] = Y[idx_row, idx_col] * Harm0[idx_row, idx_col]
Z0[idx_row, idx_col] = Z[idx_row, idx_col] * Harm0[idx_row, idx_col]
X1[idx_row, idx_col] = X[idx_row, idx_col] * Harm1[idx_row, idx_col]
Y1[idx_row, idx_col] = Y[idx_row, idx_col] * Harm1[idx_row, idx_col]
Z1[idx_row, idx_col] = Z[idx_row, idx_col] * Harm1[idx_row, idx_col]
X2[idx_row, idx_col] = X[idx_row, idx_col] * Harm2[idx_row, idx_col]
Y2[idx_row, idx_col] = Y[idx_row, idx_col] * Harm2[idx_row, idx_col]
Z2[idx_row, idx_col] = Z[idx_row, idx_col] * Harm2[idx_row, idx_col]
print(f'Processed {line_count} lines.')
# setup colors
......
code/python/plotsphericalbasis.py
View file @ a1bc511e
import matplotlib.pyplot as plt
from matplotlib import cm, colors
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import sph_harm
phi = np.linspace(0, np.pi, 100)
theta = np.linspace(0, 2*np.pi, 100)
phi, theta = np.meshgrid(phi, theta)
# nur fuer den Seiteneffekt: plt.gca(projection = '3d') funktioniert sonst nicht
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
theta_1d = np.linspace(0, np.pi, 91) # 2 GRAD Schritte
phi_1d = np.linspace(0, 2 * np.pi, 181) # 2 GRAD Schritte
theta_2d, phi_2d = np.meshgrid(theta_1d, phi_1d)
xyz_2d = np.array([np.sin(theta_2d) * np.sin(phi_2d),
np.sin(theta_2d) * np.cos(phi_2d),
np.cos(theta_2d)])
colormap = cm.ScalarMappable(cmap=plt.get_cmap("cool"))
colormap.set_clim(-.45, .45)
limit = .5
def show_Y_lm(l, m):
print("Y_%i_%i" % (l, m)) # zeigen, dass was passiert
plt.figure()
ax = plt.gca(projection="3d")
plt.title("$Y^{%i}_{%i}$" % (m, l))
Y_lm = sph_harm(m, l, phi_2d, theta_2d)
r = np.abs(Y_lm.real) * xyz_2d
ax.plot_surface(r[0], r[1], r[2],
facecolors=colormap.to_rgba(Y_lm.real),
rstride=2, cstride=2)
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_zlim(-limit, limit)
#ax.set_aspect("equal")
# ax.set_axis_off()
# The Cartesian coordinates of the unit sphereslack
x = np.sin(phi) * np.cos(theta)
y = np.sin(phi) * np.sin(theta)
z = np.cos(phi)
m, l = 2, 3
# Vorsicht: diese Schleifen erzeugen 16 plots (in 16 Fenstern)!
for l in range(0, 4):
for m in range(-l, l + 1):
show_Y_lm(l, m)
# Calculate the spherical harmonic Y(l,m) and normalize to [0,1]
fcolors = sph_harm(m, l, theta, phi).real
fmax, fmin = fcolors.max(), fcolors.min()
fcolors = (fcolors - fmin)/(fmax - fmin)
show_Y_lm(l=5, m=0)
show_Y_lm(l=5, m=4)
show_Y_lm(l=6, m=6)
# Set the aspect ratio to 1 so our sphere looks spherical
fig = plt.figure(figsize=plt.figaspect(1.))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, rstride=1, cstride=1, facecolors=cm.seismic(fcolors))
# Turn off the axis planes
ax.set_axis_off()
plt.show()
\ No newline at end of file
code/src/fluxes/upwindflux.cpp
View file @ a1bc511e
......@@ -26,6 +26,7 @@ double UpwindFlux::Flux( const Vector& Omega, double psiL, double psiR, const Ve
* @param n : Normal vector at the edge between left and right control volume
* @return resultFlux: Vector with resulting flux.
*/
Vector UpwindFlux::Flux( const Matrix AxPlus,
const Matrix AxMinus,
const Matrix AyPlus,
......
code/src/main.cpp
View file @ a1bc511e
......@@ -10,23 +10,69 @@
#include <iostream>
#include <string>
// ----
#include "quadratures/qgausslegendretensorized.h"
#include "quadratures/qmontecarlo.h"
#include "solvers/sphericalharmonics.h"
int main( int argc, char** argv ) {
MPI_Init( &argc, &argv );
std::string filename = ParseArguments( argc, argv );
QGaussLegendreTensorized quad( 4 );
SphericalHarmonics testBase( 2 );
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];
// std::cout << idx_quad << ": " << results[0] << "\n";
}
}
std::cout << "moment integration:\n " << results << "\n";
// Load Settings from File
Config* config = new Config( filename );
// Normalized integration
Vector resultsNormal( moment.size(), 0.0 );
// Print input file and run info to file
PrintLogHeader( filename );
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 );
// Build solver
Solver* solver = Solver::Create( config );
for( unsigned idx_sys = 0; idx_sys < 9; idx_sys++ ) {
moment[idx_sys] /= results[idx_sys];
resultsNormal[idx_sys] += w * moment[idx_sys] * moment[idx_sys];
// std::cout << idx_quad << ": " << results[0] << "\n";
}
}
std::cout << "moment integration:\n " << resultsNormal << "\n";
// Run solver and export
solver->Solve();
solver->Save();
// 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();
MPI_Finalize();
return EXIT_SUCCESS;
......
code/src/problems/linesource.cpp
View file @ a1bc511e
......@@ -72,18 +72,18 @@ VectorVector LineSource_PN::SetupIC() {
VectorVector cellMids = _mesh->GetCellMidPoints();
// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments are zero.
double t = 3.2e-4; // pseudo time for gaussian smoothing (Approx to dirac impulse)
for( unsigned j = 0; j < cellMids.size(); ++j ) {
double x = cellMids[j][0];
double y = cellMids[j][1];
psi[j][0] = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) );
}
// double t = 3.2e-4; // pseudo time for gaussian smoothing (Approx to dirac impulse)
// for( unsigned j = 0; j < cellMids.size(); ++j ) {
// if( cellMids[j][0] < 0 && cellMids[j][1] > 0 )
// psi[j][0] = 1.0; // / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) ) / ( 4 * M_PI );
// else
// psi[j][0] = 0.0;
// double x = cellMids[j][0];
// double y = cellMids[j][1];
// psi[j][0] = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) );
//}
for( unsigned j = 0; j < cellMids.size(); ++j ) {
if( cellMids[j][0] < 0 && cellMids[j][1] > 0 )
psi[j][0] = 1.0;
else
psi[j][0] = 0.0;
}
return psi;
}
code/src/problems/problembase.cpp
View file @ a1bc511e
......@@ -3,7 +3,10 @@
#include "problems/electronrt.h"
#include "problems/linesource.h"
ProblemBase::ProblemBase( Config* settings, Mesh* mesh ) : _settings( settings ), _mesh( mesh ) {}
ProblemBase::ProblemBase( Config* settings, Mesh* mesh ) {
_settings = settings;
_mesh = mesh;
}
ProblemBase::~ProblemBase() {}
......
code/src/solvers/mnsolver.cpp
View file @ a1bc511e
#include "solvers/mnsolver.h"
#include "toolboxes/textprocessingtoolbox.h"
//#include <chrono>
#include <mpi.h>
MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder( settings->GetMaxMomentDegree() ), _basis( _nMaxMomentsOrder ) {
......@@ -7,6 +8,7 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder(
_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!
......@@ -32,6 +34,10 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder(
for( unsigned idx_cells = 0; idx_cells < _nTotalEntries; idx_cells++ ) {
_alpha[idx_cells] = Vector( _nTotalEntries, 0.0 );
}
// Initialize and Pre-Compute Moments at quadrature points
_moments = VectorVector( _nq );
ComputeMoments();
}
MNSolver::~MNSolver() {
......@@ -49,32 +55,64 @@ int MNSolver::GlobalIndex( int l, int k ) const {
Vector MNSolver::ConstructFlux( unsigned idx_cell, unsigned idx_neigh ) {
// ---- Integration of Moment of flux ----
double x, y, z, w, entropy, velocity;
Vector moment( _nTotalEntries, 0.0 );
double x, y, w, entropy, velocity;
// double z;
Vector flux( _nTotalEntries, 0.0 );
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
x = _quadrature->GetPoints()[idx_quad][0];
y = _quadrature->GetPoints()[idx_quad][1];
z = _quadrature->GetPoints()[idx_quad][2];
w = _quadrature->GetWeights()[idx_quad];
moment = _basis.ComputeSphericalBasis( x, y, z );
x = _quadrature->GetPoints()[idx_quad][0];
y = _quadrature->GetPoints()[idx_quad][1];
// z = _quadrature->GetPoints()[idx_quad][2];
w = _quadrature->GetWeights()[idx_quad];
// --- get upwinding direction (2D)
// DOUBLE CHECK IF WE NEED DIMENSION SPLITTIN!!!
// DOUBLE CHECK IF WE NEED DIMENSION SPLITTING!!!
velocity = _normals[idx_cell][idx_neigh][0] * x + _normals[idx_cell][idx_neigh][1] * y;
// Perform upwinding
( velocity > 0 ) ? ( entropy = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], moment ) ) )
: ( entropy = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_neigh], moment ) ) );
( velocity > 0 ) ? ( entropy = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _moments[idx_quad] ) ) )
: ( entropy = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_neigh], _moments[idx_quad] ) ) );
// integrate
flux += moment * ( w * velocity * entropy );
flux += _moments[idx_quad] * ( w * velocity * entropy );
}
return flux;
}
void MNSolver::ComputeMoments() {
double x, y, z, w;
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
x = _quadrature->GetPoints()[idx_quad][0];
y = _quadrature->GetPoints()[idx_quad][1];
z = _quadrature->GetPoints()[idx_quad][2];
_moments[idx_quad] = _basis.ComputeSphericalBasis( x, y, z );
}
Vector results( _moments[0].size(), 0.0 );
// Normalize Moments
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
x = _quadrature->GetPoints()[idx_quad][0];
y = _quadrature->GetPoints()[idx_quad][1];
z = _quadrature->GetPoints()[idx_quad][2];
w = _quadrature->GetWeights()[idx_quad];
for( unsigned idx_sys = 0; idx_sys < 9; idx_sys++ ) {
results[idx_sys] += w * _moments[idx_quad][idx_sys] * _moments[idx_quad][idx_sys];
}
}
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
for( unsigned idx_sys = 0; idx_sys < 9; idx_sys++ ) {
_moments[idx_quad][idx_sys] /= results[idx_sys];
}
}
}
void MNSolver::Solve() {
int rank;
......@@ -97,6 +135,10 @@ void MNSolver::Solve() {
if( rank == 0 ) log->info( "{:10} {:10}", "t", "dFlux" );
if( rank == 0 ) log->info( "{:03.8f} {:01.5e} {:01.5e}", -1.0, dFlux, mass1 );
// Time measurement
// auto start = chrono::steady_clock::now();
// auto end = chrono::steady_clock::now();
// Loop over energies (pseudo-time of continuous slowing down approach)
for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
......@@ -119,18 +161,22 @@ void MNSolver::Solve() {
for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
// Store fluxes in psiNew, to save memory
psiNew[idx_cell] += ConstructFlux( idx_cell, idx_neigh );
// std::cout << " psiNew[" << idx_cell << "] " << psiNew[idx_cell] << "\n";
}
// std::cout << " TOTAL psiNew[" << idx_cell << "] " << psiNew[idx_cell] << "\n";
// ------ FVM Step ------
// NEED TO VECTORIZE
for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
psiNew[idx_cell][idx_system] = _psi[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _psi[idx_cell][idx_system] *
( _sigmaA[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
psiNew[idx_cell][idx_system] =
_psi[idx_cell][idx_system] - ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system]; /* cell averaged flux */
//- _dE * _psi[idx_cell][idx_system] *
// ( _sigmaA[idx_energy][idx_cell] /* absorbtion influence */
// + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
}
_psi = psiNew;
......
code/src/solvers/pnsolver.cpp
View file @ a1bc511e
......@@ -105,7 +105,7 @@ void PNSolver::Solve() {
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue; // Dirichlet cells stay at IC, farfield assumption
// Reset temporary variable psiNew
for( int idx_lDegree = 0; idx_lDegree <= int( _nq ); idx_lDegree++ ) {
for( int idx_lDegree = 0; idx_lDegree <= int( _LMaxDegree ); idx_lDegree++ ) {
for( int idx_kOrder = -idx_lDegree; idx_kOrder <= idx_lDegree; idx_kOrder++ ) {
idx_system = unsigned( GlobalIndex( idx_lDegree, idx_kOrder ) );
psiNew[idx_cell][idx_system] = 0.0;
......@@ -132,7 +132,7 @@ void PNSolver::Solve() {
}
// time update angular flux with numerical flux and total scattering cross section
for( int idx_lOrder = 0; idx_lOrder <= int( _nq ); idx_lOrder++ ) {
for( int idx_lOrder = 0; idx_lOrder <= int( _LMaxDegree ); idx_lOrder++ ) {
for( int idx_kOrder = -idx_lOrder; idx_kOrder <= idx_lOrder; idx_kOrder++ ) {
idx_system = unsigned( GlobalIndex( idx_lOrder, idx_kOrder ) );
......@@ -168,7 +168,7 @@ void PNSolver::ComputeSystemMatrices() {
unsigned idx_row = 0;
// loop over columns of A
for( int idx_lOrder = 0; idx_lOrder <= int( _nq ); idx_lOrder++ ) { // index of legendre polynom
for( int idx_lOrder = 0; idx_lOrder <= int( _LMaxDegree ); idx_lOrder++ ) { // index of legendre polynom
for( int idx_kOrder = -idx_lOrder; idx_kOrder <= idx_lOrder; idx_kOrder++ ) { // second index of legendre function
idx_row = unsigned( GlobalIndex( idx_lOrder, idx_kOrder ) );
......@@ -344,7 +344,7 @@ void PNSolver::ComputeScatterMatrix() {
}
}
double PNSolver::LegendrePoly( double x, int l ) {
double PNSolver::LegendrePoly( double x, int l ) { // Legacy. TO BE DELETED
// Pre computed low order polynomials for faster computation
switch( l ) {
case 0: return 1;
......
code/src/solvers/sphericalharmonics.cpp
View file @ a1bc511e
......@@ -46,7 +46,8 @@ 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. ) );
}
......@@ -55,7 +56,7 @@ void SphericalHarmonics::ComputeCoefficients() {
void SphericalHarmonics::ComputeAssLegendrePoly( const double my ) {
const double sintheta = std::sqrt( 1. - my * my );
double temp = std::sqrt( .5 / M_PI );
double temp = 1.0; // std::sqrt( .5 / M_PI );
_assLegendreP[GlobalIdxAssLegendreP( 0, 0 )] = temp;
......
code/tests/test_harmonics.cpp
View file @ a1bc511e
#include "catch.hpp"
#include "quadratures/qgausslegendretensorized.h"
#include "solvers/sphericalharmonics.h"
#include <fstream>
......@@ -22,19 +23,75 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical harmoni
// Remove title line
getline( case_file, text_line );
while( getline( case_file, text_line ) ) {
SECTION( "test to reference solution" ) {
// give line to stringstream
std::stringstream ss( text_line );
while( getline( case_file, 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];
// give line to stringstream
std::stringstream ss( text_line );
result = testBase.ComputeSphericalBasis( my, phi );
// Read values
ss >> my >> phi >> values[0] >> values[1] >> values[2] >> values[3] >> values[4] >> values[5] >> values[6] >> values[7] >> values[8];
for( unsigned idx = 0; idx < 9; idx++ ) {
REQUIRE( std::fabs( result[idx] - values[idx] ) < std::numeric_limits<double>::epsilon() );
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();
SECTION( "test orthogonality" ) {
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 = 1; idx_sys < 9; idx_sys++ ) {
results[idx_sys] += w * moment[idx_sys - 1] * moment[idx_sys];
}
}
for( unsigned idx_sys = 0; idx_sys < 9; idx_sys++ ) {
REQUIRE( std::fabs( results[idx_sys] ) < std::numeric_limits<double>::epsilon() );
}
}
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() );