Commit 45a834c7 authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

added monomial basis on unit sphere, added unit tests for this basis


Former-commit-id: 444eeaa2
parent 596994ca
......@@ -83,6 +83,9 @@ class Config
// Scattering Kernel
KERNEL_NAME _kernelName; /*!< @brief Scattering Kernel Name*/
// Spherical Basis
SPHERICAL_BASIS_NAME _sphericalBasisName; /*!< @brief: Name of the basis on the unit sphere */
// Optimizer
OPTIMIZER_NAME _entropyOptimizerName; /*!< @brief Choice of optimizer */
double _optimizerEpsilon; /*!< @brief termination criterion epsilon for Newton Optmizer */
......@@ -273,6 +276,8 @@ class Config
// Scattering Kernel
KERNEL_NAME inline GetKernelName() const { return _kernelName; }
// Basis name
SPHERICAL_BASIS_NAME inline GetSphericalBasisName() const { return _sphericalBasisName; }
// Output Structure
std::vector<VOLUME_OUTPUT> inline GetVolumeOutput() { return _volumeOutput; }
unsigned short inline GetNVolumeOutput() { return _nVolumeOutput; }
......
......@@ -141,4 +141,10 @@ enum SCALAR_OUTPUT { ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT };
inline std::map<std::string, SCALAR_OUTPUT> ScalarOutput_Map{
{ "ITER", ITER }, { "MASS", MASS }, { "RMS_FLUX", RMS_FLUX }, { "VTK_OUTPUT", VTK_OUTPUT }, { "CSV_OUTPUT", CSV_OUTPUT } };
// Spherical Basis Name
enum SPHERICAL_BASIS_NAME { SPHERICAL_HARMONICS, SPHERICAL_MONIMIALS };
inline std::map<std::string, SPHERICAL_BASIS_NAME> SphericalBasis_Map{ { "SPHERICAL_HARMONICS", SPHERICAL_HARMONICS },
{ "SPHERICAL_MONIMIALS", SPHERICAL_MONIMIALS } };
#endif // GLOBAL_CONSTANTS_H
......@@ -69,7 +69,7 @@ class SphericalHarmonics : public SphericalBasisBase
*/
std::vector<double> _assLegendreP;
/*! @brief sperical harmonic basis functions of
/*! @brief: spherical harmonic basis functions of
* degree 0 <= l <= L and order -l <= k <= l
* length : N = L² +2L
*/
......
......@@ -283,6 +283,10 @@ void Config::SetConfigOptions() {
/*! @brief Scattering kernel \n DESCRIPTION: Describes used scattering kernel \n DEFAULT KERNEL_Isotropic \ingroup Config */
AddEnumOption( "KERNEL", _kernelName, Kernel_Map, KERNEL_Isotropic );
/*! @brief Spherical Basis \n DESCRIPTION: Describes the chosen set of basis functions for on the unit sphere (e.g. for Moment methods) \n DEFAULT
* SPHERICAL_HARMONICS \ingroup Config */
AddEnumOption( "SPHERICAL_BASIS", _sphericalBasisName, SphericalBasis_Map, SPHERICAL_HARMONICS );
// Output related options
/*! @brief Volume output \n DESCRIPTION: Describes output groups to write to vtk \ingroup Config */
AddEnumListOption( "VOLUME_OUTPUT", _nVolumeOutput, _volumeOutput, VolOutput_Map );
......
......@@ -29,6 +29,7 @@ QuadratureBase* QuadratureBase::Create( Config* settings ) {
case QUAD_Product: return new QProduct( settings );
default: ErrorMessages::Error( "Creator for the chose quadrature does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
}
return nullptr;
}
QuadratureBase* QuadratureBase::Create( QUAD_NAME name, unsigned quadOrder ) {
......@@ -45,6 +46,7 @@ QuadratureBase* QuadratureBase::Create( QUAD_NAME name, unsigned quadOrder ) {
case QUAD_Product: return new QProduct( quadOrder );
default: ErrorMessages::Error( "Creator for the chose quadrature does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
}
return nullptr;
}
double QuadratureBase::Integrate( double( f )( double x0, double x1, double x2 ) ) {
......
/*!
* @file sphericalbasisbase.cpp
* @brief Base Class to handle basis classes on the unit sphere
* @author S. Schotthöfer
*
*/
#include "toolboxes/sphericalbasisbase.h"
#include "common/config.h"
#include "toolboxes/errormessages.h"
#include "toolboxes/sphericalharmonics.h"
#include "toolboxes/sphericalmonomials.h"
SphericalBasisBase* SphericalBasisBase::Create( Config* settings ) {
SPHERICAL_BASIS_NAME name = settings->GetSphericalBasisName();
unsigned maxMomentDegree = settings->GetMaxMomentDegree();
switch( name ) {
case SPHERICAL_HARMONICS: return new SphericalHarmonics( maxMomentDegree );
case SPHERICAL_MONIMIALS: return new SphericalMonomials( maxMomentDegree );
default: ErrorMessages::Error( "Creator for the chosen basis does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
}
return nullptr;
}
/*!
* @file sphericalmonomials.cpp
* @brief Class for efficient computation of a spherical monomial basis
* @author S. Schotthöfer
*/
#include "toolboxes/sphericalmonomials.h"
SphericalMonomials::SphericalMonomials( unsigned L_degree ) {
_LMaxDegree = L_degree;
_spatialDim = 3; // Spatial dimension 2 is hardcoded
unsigned basisSize = ComputeBasisSize( L_degree );
_YBasis = Vector( basisSize );
}
Vector SphericalMonomials::ComputeSphericalBasis( double my, double phi ) {
unsigned idx_vector = 0;
// go over all degrees of polynomials
for( unsigned idx_degree = 0; idx_degree <= _LMaxDegree; idx_degree++ ) {
// elem = Omega_x^a+ Omega_y^b +Omega_z^c : a+b+c = idx_degree
double omega_X = Omega_x( my, phi );
double omega_Y = Omega_y( my, phi );
double omega_Z = Omega_z( my );
double tmp = .0;
for( unsigned a = 0; a <= idx_degree; a++ ) {
for( unsigned b = 0; b <= idx_degree - a; b++ ) {
unsigned c = idx_degree - a - b; // c uniquely defined
tmp = Power( omega_X, a ) * Power( omega_Y, b ) * Power( omega_Z, c );
_YBasis[idx_vector] = Power( omega_X, a ) * Power( omega_Y, b ) * Power( omega_Z, c );
idx_vector++;
}
}
}
return _YBasis;
}
Vector SphericalMonomials::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 = 2 * M_PI - acos( x );
return ComputeSphericalBasis( my, phi );
}
/*! @brief: Computes the length of the basis vector for a given max degree and
* spatial dimension dim. len of a single oder: (degree + _spatialDim -1) over (degree)
* @return: lenght of whole basis */
unsigned SphericalMonomials::ComputeDimensionSize( unsigned degree ) {
return Factorial( degree + _spatialDim - 1 ) / ( Factorial( degree ) * Factorial( _spatialDim - 1 ) );
}
unsigned SphericalMonomials::ComputeBasisSize( unsigned degree ) {
unsigned basisLen = 0;
for( unsigned idx_degree = 0; idx_degree <= degree; idx_degree++ ) {
basisLen += ComputeDimensionSize( idx_degree );
}
return basisLen;
}
unsigned SphericalMonomials::Factorial( unsigned n ) { return ( n == 1 || n == 0 ) ? 1 : Factorial( n - 1 ) * n; }
double SphericalMonomials::Omega_x( double my, double phi ) { return sqrt( 1 - my * my ) * sin( phi ); }
double SphericalMonomials::Omega_y( double my, double phi ) { return sqrt( 1 - my * my ) * cos( phi ); }
double SphericalMonomials::Omega_z( double my ) { return my; }
double SphericalMonomials::Power( double basis, unsigned exponent ) {
if( exponent == 0 ) return 1.0;
double result = basis;
for( unsigned i = 1; i < exponent; i++ ) {
result = result * basis;
}
return result;
}
......@@ -3,6 +3,7 @@
#include "common/config.h"
#include "quadratures/qgausslegendretensorized.h"
#include "toolboxes/sphericalharmonics.h"
#include "toolboxes/sphericalmonomials.h"
#include <fstream>
#include <iostream>
......@@ -28,7 +29,7 @@ double P2_0( double my ) { return sqrt( 5 / ( 8 * M_PI ) ) * ( 3 * my * my - 1 )
double P2_1( double my ) { return -1 * sqrt( 15 / ( 4 * M_PI ) ) * my * sqrt( 1 - my * my ); }
double P2_2( double my ) { return sqrt( 15 / ( 16 * M_PI ) ) * ( 1 - my * my ); }
TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmonics]" ) {
TEST_CASE( "test spherical harmonics basis ", "[spherical_harmonics]" ) {
std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/solvers/unit_harmonics.cfg";
......@@ -221,3 +222,57 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmoni
REQUIRE( errorWithinBounds );
}
}
double Omega_x( double my, double phi ) { return sqrt( 1 - my * my ) * sin( phi ); }
double Omega_y( double my, double phi ) { return sqrt( 1 - my * my ) * cos( phi ); }
double Omega_z( double my ) { return my; }
double SphericalMonomial_0( double /* my */, double /* phi */ ) { return 1; }
double SphericalMonomial_1( double my, double /*phi*/ ) { return Omega_z( my ); } // omega_z
double SphericalMonomial_2( double my, double phi ) { return Omega_y( my, phi ); } // omega_y
double SphericalMonomial_3( double my, double phi ) { return Omega_x( my, phi ); } // omega_x
double SphericalMonomial_4( double my, double /*phi*/ ) { return Omega_z( my ) * Omega_z( my ); } // omega_z^2
double SphericalMonomial_5( double my, double phi ) { return Omega_y( my, phi ) * Omega_z( my ); } // omega_y*omega_z
double SphericalMonomial_6( double my, double phi ) { return Omega_y( my, phi ) * Omega_y( my, phi ); } // omega_y^2
double SphericalMonomial_7( double my, double phi ) { return Omega_x( my, phi ) * Omega_z( my ); } // omega_x*omega_z
double SphericalMonomial_8( double my, double phi ) { return Omega_x( my, phi ) * Omega_y( my, phi ); } // omega_x*omega_y
double SphericalMonomial_9( double my, double phi ) { return Omega_x( my, phi ) * Omega_x( my, phi ); } // omega_x^2
TEST_CASE( "test spherical monomial basis", "[spherical_monomials]" ) {
unsigned maxMomentDegree = 2; //==> 6+3+1 basis functions
SphericalMonomials testBase( maxMomentDegree );
SECTION( "Test against analytical solution" ) {
Vector moment;
std::vector<bool> validMoment( 10, true );
for( double my = -1.0; my < 1.0; my += 0.1 ) {
for( double phi = 0.0; phi < 2 * M_PI; phi += 0.1 ) {
moment = testBase.ComputeSphericalBasis( my, phi );
std::cout << "(algo|real): (" << moment[0] << " | " << SphericalMonomial_0( my, phi ) << ")" << std::endl;
std::cout << "(algo|real): (" << moment[1] << " | " << SphericalMonomial_1( my, phi ) << ")" << std::endl;
std::cout << "(algo|real): (" << moment[2] << " | " << SphericalMonomial_2( my, phi ) << ")" << std::endl;
std::cout << "(algo|real): (" << moment[3] << " | " << SphericalMonomial_3( my, phi ) << ")" << std::endl;
std::cout << "(algo|real): (" << moment[4] << " | " << SphericalMonomial_4( my, phi ) << ")" << std::endl;
std::cout << "(algo|real): (" << moment[5] << " | " << SphericalMonomial_5( my, phi ) << ")" << std::endl;
std::cout << "(algo|real): (" << moment[6] << " | " << SphericalMonomial_6( my, phi ) << ")" << std::endl;
std::cout << "(algo|real): (" << moment[7] << " | " << SphericalMonomial_7( my, phi ) << ")" << std::endl;
std::cout << "(algo|real): (" << moment[8] << " | " << SphericalMonomial_8( my, phi ) << ")" << std::endl;
std::cout << "(algo|real): (" << moment[9] << " | " << SphericalMonomial_9( my, phi ) << ")" << std::endl;
if( std::fabs( moment[0] - SphericalMonomial_0( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[0] = false;
if( std::fabs( moment[1] - SphericalMonomial_1( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[1] = false;
if( std::fabs( moment[2] - SphericalMonomial_2( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[2] = false;
if( std::fabs( moment[3] - SphericalMonomial_3( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[3] = false;
if( std::fabs( moment[4] - SphericalMonomial_4( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[4] = false;
if( std::fabs( moment[5] - SphericalMonomial_5( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[5] = false;
if( std::fabs( moment[6] - SphericalMonomial_6( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[6] = false;
if( std::fabs( moment[7] - SphericalMonomial_7( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[7] = false;
if( std::fabs( moment[8] - SphericalMonomial_8( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[8] = false;
if( std::fabs( moment[9] - SphericalMonomial_9( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[8] = false;
}
}
REQUIRE( std::all_of( validMoment.begin(), validMoment.end(), []( bool v ) { return v; } ) );
}
}
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