test_harmonics.cpp 10 KB
Newer Older
1
2

#include "catch.hpp"
3
#include "common/config.h"
4
#include "quadratures/qgausslegendretensorized.h"
5
#include "toolboxes/sphericalharmonics.h"
6
7

#include <fstream>
8
#include <iostream>
9
#include <sstream>
10

11
double Y0_0( double, double ) { return sqrt( 1 / ( 4 * M_PI ) ); }
12
13

double Y1_m1( double my, double phi ) { return -sqrt( 3 / ( 4 * M_PI ) ) * sqrt( 1 - my * my ) * sin( phi ); }
14
double Y1_0( double my, double /*phi*/ ) { return sqrt( 3 / ( 4 * M_PI ) ) * my; }
15
16
17
18
double Y1_1( double my, double phi ) { return -sqrt( 3 / ( 4 * M_PI ) ) * sqrt( 1 - my * my ) * cos( phi ); }

double Y2_m2( double my, double phi ) { return sqrt( 15 / ( 16 * M_PI ) ) * ( 1 - my * my ) * sin( 2 * phi ); }
double Y2_m1( double my, double phi ) { return -1 * sqrt( 15 / ( 4 * M_PI ) ) * my * sqrt( 1 - my * my ) * sin( phi ); }
19
double Y2_0( double my, double /*phi*/ ) { return sqrt( 5 / ( 16 * M_PI ) ) * ( 3 * my * my - 1 ); }
20
21
double Y2_1( double my, double phi ) { return -1 * sqrt( 15 / ( 4 * M_PI ) ) * my * sqrt( 1 - my * my ) * cos( phi ); }
double Y2_2( double my, double phi ) { return sqrt( 15 / ( 16 * M_PI ) ) * ( 1 - my * my ) * cos( 2 * phi ); }
22

23
double P0_0( double /*my*/ ) { return sqrt( 1 / ( 2 * M_PI ) ); }
24
25
26
27
28
29
30
double P1_0( double my ) { return sqrt( 3 / ( 2 * M_PI ) ) * my; }
double P1_1( double my ) { return -sqrt( 3 / ( 4 * M_PI ) ) * sqrt( 1 - my * my ); }

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

31
TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmonics]" ) {
32

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
33
    std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/solvers/unit_harmonics.cfg";
34
35
36
37

    // Load Settings from File
    Config* config = new Config( filename );

38
    unsigned maxMomentDegree = 2;
39

40
    SphericalHarmonics testBase( maxMomentDegree );
41

42
43
44
    SECTION( "Test against analytical solution" ) {
        std::vector<double> legendre;
        Vector moment;
45
46
        std::vector<bool> validLegendrePoly( 6, true );
        std::vector<bool> validMoment( 9, true );
47
48
49
50
        for( double my = -1.0; my < 1.0; my += 0.1 ) {

            legendre = testBase.GetAssLegendrePoly( my );

51
52
53
54
55
56
            if( std::fabs( legendre[0] - P0_0( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[0] = false;
            if( std::fabs( legendre[1] - P1_0( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[1] = false;
            if( std::fabs( legendre[2] - P1_1( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[2] = false;
            if( std::fabs( legendre[3] - P2_0( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[3] = false;
            if( std::fabs( legendre[4] - P2_1( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[4] = false;
            if( std::fabs( legendre[5] - P2_2( my ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validLegendrePoly[5] = false;
57
58
59
60

            for( double phi = 0.0; phi < 2 * M_PI; phi += 0.1 ) {
                moment = testBase.ComputeSphericalBasis( my, phi );

61
62
63
64
65
66
67
68
69
                if( std::fabs( moment[0] - Y0_0( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[0] = false;
                if( std::fabs( moment[1] - Y1_m1( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[1] = false;
                if( std::fabs( moment[2] - Y1_0( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[2] = false;
                if( std::fabs( moment[3] - Y1_1( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[3] = false;
                if( std::fabs( moment[4] - Y2_m2( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[4] = false;
                if( std::fabs( moment[5] - Y2_m1( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[5] = false;
                if( std::fabs( moment[6] - Y2_0( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[6] = false;
                if( std::fabs( moment[7] - Y2_1( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[7] = false;
                if( std::fabs( moment[8] - Y2_2( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[8] = false;
70
71
            }
        }
72
73
        REQUIRE( std::all_of( validLegendrePoly.begin(), validLegendrePoly.end(), []( bool v ) { return v; } ) );
        REQUIRE( std::all_of( validMoment.begin(), validMoment.end(), []( bool v ) { return v; } ) );
74
75
    }

76
77
    // Remove title line

78
    SECTION( "test to reference solution" ) {
79

80
81
82
        std::string text_line;
        std::ifstream case_file;

83
84
85
86
87
88
        double my  = 0.0;
        double phi = 0.0;

        Vector values( 9, 0.0 );
        Vector result( 4, 0.0 );

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
89
        case_file.open( "unit_test/solvers/harmonicBasis_reference.csv", std::ios::in );
90
91
92

        getline( case_file, text_line );

93
        bool errorWithinBounds = true;
94
        while( getline( case_file, text_line ) ) {
95

96
97
            // give line to stringstream
            std::stringstream ss( text_line );
98

99
100
            // Read values
            ss >> my >> phi >> values[0] >> values[1] >> values[2] >> values[3] >> values[4] >> values[5] >> values[6] >> values[7] >> values[8];
101

102
103
104
            result = testBase.ComputeSphericalBasis( my, phi );

            for( unsigned idx = 0; idx < 9; idx++ ) {
105
                if( std::fabs( result[idx] - values[idx] ) > 1e2 * std::numeric_limits<double>::epsilon() ) errorWithinBounds = false;
106
            }
107
        }
108
        REQUIRE( errorWithinBounds );
109
        case_file.close();
110
    }
111

112
113
    SECTION( "test orthonormality - spherical coordinates" ) {
        // Caution: Integration only works with spherical coordinates!
114

115
        QGaussLegendreTensorized quad( config );
116
117
118
119
120

        double my, phi, w;
        Vector moment = testBase.ComputeSphericalBasis( 0, 1, 0 );
        // 9 basis moments if degree = 2

121
        Matrix results( moment.size(), moment.size(), 0.0 );
122
123

        for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
124
125
126
            my  = quad.GetPointsSphere()[idx_quad][0];
            phi = quad.GetPointsSphere()[idx_quad][1];
            // z      = quad.GetPoints()[idx_quad][2];
127
128
129
            w      = quad.GetWeights()[idx_quad];
            moment = testBase.ComputeSphericalBasis( my, phi );

130
131
132
133
            for( unsigned idx_row = 0; idx_row < 9; idx_row++ ) {
                for( unsigned idx_col = 0; idx_col < 9; idx_col++ ) {
                    results( idx_row, idx_col ) += w * moment[idx_row] * moment[idx_col];
                }
134
135
136
            }
        }

137
138
        bool errorWithinBounds = true;
        bool orthogonality     = true;
139
140
141
142
        for( unsigned idx_row = 0; idx_row < 9; idx_row++ ) {
            for( unsigned idx_col = 0; idx_col < 9; idx_col++ ) {
                if( idx_row == idx_col ) {
                    // Orthogonality
143
                    if( std::fabs( results( idx_row, idx_col ) - 1.0 ) > 1e2 * std::numeric_limits<double>::epsilon() ) orthogonality = false;
144
145
146
                }
                else {
                    // Normality
147
                    if( std::fabs( results( idx_row, idx_col ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) errorWithinBounds = false;
148
149
                }
            }
150
        }
151
152
        REQUIRE( errorWithinBounds );
        REQUIRE( orthogonality );
153
    }
154
155
156
157
158
159
160

    SECTION( "test parity - carthesian coordinates" ) {

        Vector moment1 = testBase.ComputeSphericalBasis( 0, 0 );
        Vector moment2 = testBase.ComputeSphericalBasis( 0, 0 );

        // Parity in carthesian coordinates
161
        QGaussLegendreTensorized quad( config );
162
163
164

        double x, y, z;

165
        bool errorWithinBounds = true;
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
        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];
            moment1 = testBase.ComputeSphericalBasis( x, y, z );
            moment2 = testBase.ComputeSphericalBasis( -x, -y, -z );

            int idx_sys;
            double result = 0.;

            for( int l_idx = 0; l_idx <= int( maxMomentDegree ); l_idx++ ) {
                for( int k_idx = -l_idx; k_idx <= l_idx; k_idx++ ) {
                    idx_sys = testBase.GlobalIdxBasis( l_idx, k_idx );

                    if( l_idx % 2 == 0 )
                        result = moment2[idx_sys] - moment1[idx_sys];
                    else
                        result = moment2[idx_sys] + moment1[idx_sys];

185
                    if( std::fabs( result ) > 1e2 * std::numeric_limits<double>::epsilon() ) errorWithinBounds = false;
186
187
188
                }
            }
        }
189
        REQUIRE( errorWithinBounds );
190
191
192
193
194
195
196
197
198
199
200
    }

    SECTION( "test parity - polar coordinates" ) {

        Vector moment1 = testBase.ComputeSphericalBasis( 0, 0 );
        Vector moment2 = testBase.ComputeSphericalBasis( 0, 0 );

        int idx_sys;
        double result = 0.;

        // // test in polar coordinates
201
        bool errorWithinBounds = true;
202
203
204
205
206
207
208
209
210
211
212
213
214
215
        for( double my = -1.0; my < 1.0; my += 0.1 ) {
            for( double phi = 0.0; phi < 2 * M_PI; phi += 0.1 ) {
                moment2 = testBase.ComputeSphericalBasis( my, phi );
                moment1 = testBase.ComputeSphericalBasis( -my, M_PI + phi );

                for( int l_idx = 0; l_idx <= int( maxMomentDegree ); l_idx++ ) {
                    for( int k_idx = -l_idx; k_idx <= l_idx; k_idx++ ) {
                        idx_sys = testBase.GlobalIdxBasis( l_idx, k_idx );

                        if( l_idx % 2 == 0 )
                            result = moment2[idx_sys] - moment1[idx_sys];
                        else
                            result = moment2[idx_sys] + moment1[idx_sys];

216
                        if( std::fabs( result ) > 1e2 * std::numeric_limits<double>::epsilon() ) errorWithinBounds = false;
217
218
219
220
                    }
                }
            }
        }
221
        REQUIRE( errorWithinBounds );
222
    }
223
}