CSMMLab / KiT-RT / #95

Committed 28 May 2025 02:06AM UTC coverage: 46.655% (-11.1%) from 57.728%

Build # #95

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

Release 2025 (#44)

* New neural models (#30)

* extend kitrt script

* added new neural models

* extend to m3 models

* added M3_2D models

* added M2 and M4 models

* configure ml optimizer for high order models

* added configuration for high order models

* add option for rotation postprcessing in neural networks. remove unneccessary python scripts

* started rotational symmetric postprocessing

* added rotational invariance postprocessing for m2 and m1 models

* fix post merge bugs

* created hohlraum mesh

* add hohlraum tes case

* add hohlraum test case

* add hohlraum cfg filees

* fixed hohlraum testcase

* add hohlraum cfg files

* changed hohlraum cfg

* changed hohlraum testcase to isotropic inflow source boundary condition

* added ghost cell bonudary for hohlraum testcase

* update readme and linesource mesh creator

* added proper scaling for linesource reference solution

* regularized newton debugging

* Data generator with reduced objective functional (#33)

* added reduced optimizer for sampling

* remove old debugging comments

* mesh acceleration (#38)

* added new ansatz (#36)

* added new ansatz

* debug new ansatz

* branch should be abandoned here

* debug new ansatz

* fix scaling error new ansatz

* fix config errors

* temporarily fixed dynamic ansatz rotation bug

* fix inheritance error for new ansatz

* mesh acceleration

* add structured hohlraum

* Mesh acc (#41)

* mesh acceleration

* added floor value for starmap moment methods

* enable accelleration

* delete minimum value starmap

* Update README.md

* Spherical harmonics nn (#40)

* added new ansatz

* debug new ansatz

* branch should be abandoned here

* debug new ansatz

* fix scaling error new ansatz

* fix config errors

* temporarily fixed dynamic ansatz rotation bug

* fix inheritance error for new ansatz

* mesh acceleration

* add structured hohlraum

* added sph... (continued)

Run Details

767 of 3019 new or added lines in 51 files covered. (25.41%)

131 existing lines in 8 files now uncovered.

4422 of 9478 relevant lines covered (46.66%)

57163.05 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

83.06

/src/velocitybasis/sphericalharmonics.cpp

#include "velocitybasis/sphericalharmonics.hpp"
#include "toolboxes/errormessages.hpp"
#include <math.h>

SphericalHarmonics::SphericalHarmonics( unsigned L_degree ) {
    _LMaxDegree = L_degree;
    _spatialDim = 3;    // Default no projection

    unsigned associatedLegendrePolySize = GlobalIdxAssLegendreP( _LMaxDegree, _LMaxDegree ) + 1;

    _aParam       = std::vector<double>( associatedLegendrePolySize, 0.0 );
    _bParam       = std::vector<double>( associatedLegendrePolySize, 0.0 );
    _assLegendreP = std::vector<double>( associatedLegendrePolySize, 0.0 );

    ComputeCoefficients();

    _YBasis = Vector( GetBasisSizeUnprojected(), 0.0 );
}

SphericalHarmonics::SphericalHarmonics( unsigned L_degree, unsigned short spatialDim ) {
    _LMaxDegree = L_degree;
    _spatialDim = (unsigned)spatialDim;

    unsigned associatedLegendrePolySize = GlobalIdxAssLegendreP( _LMaxDegree, _LMaxDegree ) + 1;

    _aParam       = std::vector<double>( associatedLegendrePolySize, 0.0 );
    _bParam       = std::vector<double>( associatedLegendrePolySize, 0.0 );
    _assLegendreP = std::vector<double>( associatedLegendrePolySize, 0.0 );

    ComputeCoefficients();

    _YBasis = Vector( GetBasisSizeUnprojected(), 0.0 );
}

unsigned SphericalHarmonics::GetBasisSize() {
    unsigned count = 0;

    switch( _spatialDim ) {
        case 1: return _LMaxDegree + 1; break;
        case 2:
            for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
                for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
                    if( ( idx_l + idx_k ) % 2 == 0 ) {
                        count++;
                    }
                }
            }
            return count;
            break;
        default:
            return GetGlobalIndexBasis( _LMaxDegree, _LMaxDegree ) + 1; /* +1, since globalIdx computes indices */
            break;
    }
}

unsigned SphericalHarmonics::GetBasisSizeUnprojected() {
    return GetGlobalIndexBasis( _LMaxDegree, _LMaxDegree ) + 1; /* +1, since globalIdx computes indices */
}

unsigned SphericalHarmonics::GetCurrDegreeSize( unsigned currDegreeL ) { return 2 * currDegreeL + 1; }

unsigned SphericalHarmonics::GetGlobalIndexBasis( int l_degree, int k_order ) {
    if( l_degree < 0 ) ErrorMessages::Error( "Negative polynomial degrees not supported.", CURRENT_FUNCTION );
    if( k_order < -l_degree || k_order > l_degree )
        ErrorMessages::Error( "Order k of spherical harmonics basis must be in bounds -l<= k <= l.", CURRENT_FUNCTION );

    return (unsigned)( k_order + l_degree /*number of previous indices in current level*/ +
                       l_degree * l_degree ) /* number of previous indices untill level l-1 */;
}

Vector SphericalHarmonics::ComputeSphericalBasis( double my, double phi, double r ) {
    // radius r is not used and always assumed to be 1

    ComputeAssLegendrePoly( my );
    ComputeYBasis( phi );

    // For 1D, just use the terms of order k=0
    if( _spatialDim == 1 ) {
        Vector YBasis1D( GetBasisSize(), 0.0 );
        for( unsigned idx_l = 0; idx_l <= _LMaxDegree; idx_l++ ) {
            YBasis1D[idx_l] = _YBasis[GetGlobalIndexBasis( idx_l, 0 )];
        }
        return r * YBasis1D;
    }
    // For 2D, just use the terms idx_l + idx_k % 2 == 0
    if( _spatialDim == 2 ) {
        Vector YBasis2D( GetBasisSize(), 0.0 );
        unsigned count = 0;
        for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
            for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
                if( ( idx_l + idx_k ) % 2 == 0 ) {
                    YBasis2D[count] = _YBasis[GetGlobalIndexBasis( idx_l, idx_k )];
                    count++;
                }
            }
        }
        return r * YBasis2D;
    }
    // Do nothing in the 3D case
    return r * _YBasis;
}

Vector SphericalHarmonics::ComputeSphericalBasisKarthesian( 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 );

    ComputeAssLegendrePoly( my );
    ComputeYBasis( phi );
    return _YBasis;
}

std::vector<double> SphericalHarmonics::GetAssLegendrePoly( const double my ) {
    ComputeAssLegendrePoly( my );
    return _assLegendreP;
}

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
    for( unsigned l_idx = 2; l_idx <= _LMaxDegree; l_idx++ ) {
        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;

            _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. ) );
        }
    }
}

void SphericalHarmonics::ComputeAssLegendrePoly( const double my ) {
    const double sintheta = std::sqrt( 1. - my * my );
    double temp           = std::sqrt( .5 / M_PI );

    _assLegendreP[GlobalIdxAssLegendreP( 0, 0 )] = temp;

    if( _LMaxDegree > 0 ) {
        const double SQRT3                           = std::sqrt( 3 );       // 1.7320508075688772935
        _assLegendreP[GlobalIdxAssLegendreP( 1, 0 )] = my * SQRT3 * temp;    // 1.224744871391589
        const double SQRT3DIV2                       = -std::sqrt( 3. / 2. );

        temp                                         = SQRT3DIV2 * sintheta * temp;
        _assLegendreP[GlobalIdxAssLegendreP( 1, 1 )] = temp;

        for( unsigned l_idx = 2; l_idx <= _LMaxDegree; l_idx++ ) {
            for( unsigned k_idx = 0; k_idx < l_idx - 1; k_idx++ ) {
                _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] =
                    _aParam[GlobalIdxAssLegendreP( l_idx, k_idx )] *
                    ( my * _assLegendreP[GlobalIdxAssLegendreP( l_idx - 1, k_idx )] +
                      _bParam[GlobalIdxAssLegendreP( l_idx, k_idx )] * _assLegendreP[GlobalIdxAssLegendreP( l_idx - 2, k_idx )] );
            }
            _assLegendreP[GlobalIdxAssLegendreP( l_idx, l_idx - 1 )] = my * std::sqrt( 2 * ( l_idx - 1 ) + 3 ) * temp;

            temp = -std::sqrt( 1.0 + 0.5 / l_idx ) * sintheta * temp;

            _assLegendreP[GlobalIdxAssLegendreP( l_idx, l_idx )] = temp;
        }
    }
}

void SphericalHarmonics::ComputeYBasis( const double phi ) {
    for( unsigned l_idx = 0; l_idx <= _LMaxDegree; l_idx++ ) {
        _YBasis[GetGlobalIndexBasis( l_idx, 0 )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, 0 )] * 0.5 * M_SQRT2;    // M_SQRT2 = sqrt(2)
    }

    // helper constants
    double c1 = 1.0;
    double c2 = std::cos( phi );
    double s1 = 0.0;
    double s2 = -std::sin( phi );
    double tc = 2.0 * c2;

    double s = 0.0;
    double c = 0.0;

    for( unsigned k_idx = 1; k_idx <= _LMaxDegree; k_idx++ ) {
        s  = tc * s1 - s2;    // addition theorem
        c  = tc * c1 - c2;    // addition theorem
        s2 = s1;
        s1 = s;
        c2 = c1;
        c1 = c;
        for( unsigned l_idx = k_idx; l_idx <= _LMaxDegree; l_idx++ ) {
            _YBasis[GetGlobalIndexBasis( l_idx, -k_idx )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] * s;
            _YBasis[GetGlobalIndexBasis( l_idx, k_idx )]  = _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] * c;
        }
    }
}

1	#include "velocitybasis/sphericalharmonics.hpp"
2	#include "toolboxes/errormessages.hpp"
3	#include <math.h>
4
5	SphericalHarmonics::SphericalHarmonics( unsigned L_degree ) {	14✔
6	_LMaxDegree = L_degree;	14✔
7	_spatialDim = 3; // Default no projection	14✔
8
9	unsigned associatedLegendrePolySize = GlobalIdxAssLegendreP( _LMaxDegree, _LMaxDegree ) + 1;	14✔
10
11	_aParam = std::vector<double>( associatedLegendrePolySize, 0.0 );	14✔
12	_bParam = std::vector<double>( associatedLegendrePolySize, 0.0 );	14✔
13	_assLegendreP = std::vector<double>( associatedLegendrePolySize, 0.0 );	14✔
14
15	ComputeCoefficients();	14✔
16
17	_YBasis = Vector( GetBasisSizeUnprojected(), 0.0 );	14✔
18	}	14✔
19
20	SphericalHarmonics::SphericalHarmonics( unsigned L_degree, unsigned short spatialDim ) {	36✔
21	_LMaxDegree = L_degree;	36✔
22	_spatialDim = (unsigned)spatialDim;	36✔
23
24	unsigned associatedLegendrePolySize = GlobalIdxAssLegendreP( _LMaxDegree, _LMaxDegree ) + 1;	36✔
25
26	_aParam = std::vector<double>( associatedLegendrePolySize, 0.0 );	36✔
27	_bParam = std::vector<double>( associatedLegendrePolySize, 0.0 );	36✔
28	_assLegendreP = std::vector<double>( associatedLegendrePolySize, 0.0 );	36✔
29
30	ComputeCoefficients();	36✔
31
32	_YBasis = Vector( GetBasisSizeUnprojected(), 0.0 );	36✔
33	}	36✔
34
35	unsigned SphericalHarmonics::GetBasisSize() {	36✔
36	unsigned count = 0;	36✔
37
38	switch( _spatialDim ) {	36✔
39	case 1: return _LMaxDegree + 1; break;	×
NEW 40	case 2:	×
NEW 41	for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {	×
NEW 42	for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {	×
NEW 43	if( ( idx_l + idx_k ) % 2 == 0 ) {	×
NEW 44	count++;	×
45	}
46	}
47	}
NEW 48	return count;	×
49	break;
50	default:	36✔
51	return GetGlobalIndexBasis( _LMaxDegree, _LMaxDegree ) + 1; /* +1, since globalIdx computes indices */	36✔
52	break;
53	}
54	}
55
56	unsigned SphericalHarmonics::GetBasisSizeUnprojected() {	50✔
57	return GetGlobalIndexBasis( _LMaxDegree, _LMaxDegree ) + 1; /* +1, since globalIdx computes indices */	50✔
58	}
59
60	unsigned SphericalHarmonics::GetCurrDegreeSize( unsigned currDegreeL ) { return 2 * currDegreeL + 1; }	×
61
62	unsigned SphericalHarmonics::GetGlobalIndexBasis( int l_degree, int k_order ) {	125,660✔
63	if( l_degree < 0 ) ErrorMessages::Error( "Negative polynomial degrees not supported.", CURRENT_FUNCTION );	125,660✔
64	if( k_order < -l_degree \|\| k_order > l_degree )	125,660✔
65	ErrorMessages::Error( "Order k of spherical harmonics basis must be in bounds -l<= k <= l.", CURRENT_FUNCTION );	×
66
67	return (unsigned)( k_order + l_degree /number of previous indices in current level/ +	125,660✔
68	l_degree * l_degree ) /* number of previous indices untill level l-1 */;	125,660✔
69	}
70
71	Vector SphericalHarmonics::ComputeSphericalBasis( double my, double phi, double r ) {	9,610✔
72	// radius r is not used and always assumed to be 1
73
74	ComputeAssLegendrePoly( my );	9,610✔
75	ComputeYBasis( phi );	9,610✔
76
77	// For 1D, just use the terms of order k=0
78	if( _spatialDim == 1 ) {	9,610✔
79	Vector YBasis1D( GetBasisSize(), 0.0 );	×
80	for( unsigned idx_l = 0; idx_l <= _LMaxDegree; idx_l++ ) {	×
81	YBasis1D[idx_l] = _YBasis[GetGlobalIndexBasis( idx_l, 0 )];	×
82	}
83	return r * YBasis1D;	×
84	}
85	// For 2D, just use the terms idx_l + idx_k % 2 == 0
86	if( _spatialDim == 2 ) {	9,610✔
87	Vector YBasis2D( GetBasisSize(), 0.0 );	×
88	unsigned count = 0;	×
89	for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {	×
90	for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {	×
NEW 91	if( ( idx_l + idx_k ) % 2 == 0 ) {	×
92	YBasis2D[count] = _YBasis[GetGlobalIndexBasis( idx_l, idx_k )];	×
93	count++;	×
94	}
95	}
96	}
97	return r * YBasis2D;	×
98	}
99	// Do nothing in the 3D case
100	return r * _YBasis;	19,220✔
101	}
102
103	Vector SphericalHarmonics::ComputeSphericalBasisKarthesian( double x, double y, double z ) {	1,154✔
104
105	// transform (x,y,z) into (my,phi)
106	double my = z;	1,154✔
107	double phi = 0.0;	1,154✔
108
109	if( y >= 0 )	1,154✔
110	phi = acos( x );	578✔
111	else
112	phi = 2 * M_PI - acos( x );	576✔
113
114	ComputeAssLegendrePoly( my );	1,154✔
115	ComputeYBasis( phi );	1,154✔
116	return _YBasis;	1,154✔
117	}
118
119	std::vector<double> SphericalHarmonics::GetAssLegendrePoly( const double my ) {	42✔
120	ComputeAssLegendrePoly( my );	42✔
121	return _assLegendreP;	42✔
122	}
123
124	void SphericalHarmonics::ComputeCoefficients() {	50✔
125	// m in paper is here denoted by k
126	double ls = 0.0; // l^2	50✔
127	double lm1s = 0.0; // (l-1)^2	50✔
128	double ks = 0.0; // k^2	50✔
129	for( unsigned l_idx = 2; l_idx <= _LMaxDegree; l_idx++ ) {	118✔
130	ls = l_idx * l_idx;	68✔
131	lm1s = ( l_idx - 1 ) * ( l_idx - 1 );	68✔
132	for( unsigned k_idx = 0; k_idx < l_idx - 1; k_idx++ ) {	196✔
133	ks = k_idx * k_idx;	128✔
134
135	_aParam[GlobalIdxAssLegendreP( l_idx, k_idx )] = std::sqrt( ( 4 * ls - 1. ) / ( ls - ks ) );	128✔
136	_bParam[GlobalIdxAssLegendreP( l_idx, k_idx )] = -std::sqrt( ( lm1s - ks ) / ( 4 * lm1s - 1. ) );	128✔
137	}
138	}
139	}	50✔
140
141	void SphericalHarmonics::ComputeAssLegendrePoly( const double my ) {	10,806✔
142	const double sintheta = std::sqrt( 1. - my * my );	10,806✔
143	double temp = std::sqrt( .5 / M_PI );	10,806✔
144
145	_assLegendreP[GlobalIdxAssLegendreP( 0, 0 )] = temp;	10,806✔
146
147	if( _LMaxDegree > 0 ) {	10,806✔
148	const double SQRT3 = std::sqrt( 3 ); // 1.7320508075688772935	10,806✔
149	_assLegendreP[GlobalIdxAssLegendreP( 1, 0 )] = my * SQRT3 * temp; // 1.224744871391589	10,806✔
150	const double SQRT3DIV2 = -std::sqrt( 3. / 2. );	10,806✔
151
152	temp = SQRT3DIV2 * sintheta * temp;	10,806✔
153	_assLegendreP[GlobalIdxAssLegendreP( 1, 1 )] = temp;	10,806✔
154
155	for( unsigned l_idx = 2; l_idx <= _LMaxDegree; l_idx++ ) {	21,548✔
156	for( unsigned k_idx = 0; k_idx < l_idx - 1; k_idx++ ) {	21,484✔
157	_assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] =	10,742✔
158	_aParam[GlobalIdxAssLegendreP( l_idx, k_idx )] *	10,742✔
159	( my * _assLegendreP[GlobalIdxAssLegendreP( l_idx - 1, k_idx )] +	10,742✔
160	_bParam[GlobalIdxAssLegendreP( l_idx, k_idx )] * _assLegendreP[GlobalIdxAssLegendreP( l_idx - 2, k_idx )] );	10,742✔
161	}
162	_assLegendreP[GlobalIdxAssLegendreP( l_idx, l_idx - 1 )] = my * std::sqrt( 2 * ( l_idx - 1 ) + 3 ) * temp;	10,742✔
163
164	temp = -std::sqrt( 1.0 + 0.5 / l_idx ) * sintheta * temp;	10,742✔
165
166	_assLegendreP[GlobalIdxAssLegendreP( l_idx, l_idx )] = temp;	10,742✔
167	}
168	}
169	}	10,806✔
170
171	void SphericalHarmonics::ComputeYBasis( const double phi ) {	10,764✔
172	for( unsigned l_idx = 0; l_idx <= _LMaxDegree; l_idx++ ) {	42,992✔
173	_YBasis[GetGlobalIndexBasis( l_idx, 0 )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, 0 )] * 0.5 * M_SQRT2; // M_SQRT2 = sqrt(2)	32,228✔
174	}
175
176	// helper constants
177	double c1 = 1.0;	10,764✔
178	double c2 = std::cos( phi );	10,764✔
179	double s1 = 0.0;	10,764✔
180	double s2 = -std::sin( phi );	10,764✔
181	double tc = 2.0 * c2;	10,764✔
182
183	double s = 0.0;	10,764✔
184	double c = 0.0;	10,764✔
185
186	for( unsigned k_idx = 1; k_idx <= _LMaxDegree; k_idx++ ) {	32,228✔
187	s = tc * s1 - s2; // addition theorem	21,464✔
188	c = tc * c1 - c2; // addition theorem	21,464✔
189	s2 = s1;	21,464✔
190	s1 = s;	21,464✔
191	c2 = c1;	21,464✔
192	c1 = c;	21,464✔
193	for( unsigned l_idx = k_idx; l_idx <= _LMaxDegree; l_idx++ ) {	53,628✔
194	_YBasis[GetGlobalIndexBasis( l_idx, -k_idx )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] * s;	32,164✔
195	_YBasis[GetGlobalIndexBasis( l_idx, k_idx )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] * c;	32,164✔
196	}
197	}
198	}	10,764✔

CSMMLab / KiT-RT / #95

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous