test_quadrature.cpp 13 KB
Newer Older
1
#include "catch.hpp"
2
#include "common/config.h"
3
#include "common/globalconstants.h"
4
5
6
7
#include "quadratures/quadraturebase.h"

#include <vector>

steffen.schotthoefer's avatar
steffen.schotthoefer committed
8
9
std::vector<QUAD_NAME> quadraturenames = {
    QUAD_MonteCarlo, QUAD_GaussLegendreTensorized, QUAD_GaussLegendre1D, QUAD_LevelSymmetric, QUAD_Lebedev, QUAD_LDFESA };
10

11
12
std::vector<std::vector<int>> quadratureorders = {
    { 4, 5, 6, 7 },                            // Monte Carlo
13
    { 4, 6, 8, 10 },                           // Gauss Legendre
steffen.schotthoefer's avatar
steffen.schotthoefer committed
14
    { 4, 6, 8, 10 },                           // Gauss Legendre 1D
15
16
17
18
19
20
21
    { 2, 4, 6, 8, 10, 12, 14, 16, 18, 20 },    // Available Orders for LevelSymmetric
    { 3,  5,  7,  9,  11, 13, 15, 17, 19, 21, 23,  25,  27,  29,  31,  35,
      41, 47, 53, 59, 65, 71, 77, 83, 89, 95, 101, 107, 113, 119, 125, 131 },    // Available orders for Lebedev
    { 1, 2, 3 }                                                                  // Available Orders for LDFESA
};

bool approxequal( double a, double b, bool lowAccuracy = false ) {
22
    double tol = 1e-12;              // For computed quadrature weights
23
    if( lowAccuracy ) tol = 5e-5;    // Mainly for Lookup Quadratures
steffen.schotthoefer's avatar
steffen.schotthoefer committed
24
    return std::abs( a - b ) < tol;
25
26
}

27
28
29
30
double f( double x, double y, double z ) {
    return x * x + y * y + z * z;    // == 1
}

31
double sin( double x, double /*y*/, double /*z*/ ) { return sin( x ); }
32
33

TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
34
    std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/quadratures/unit_quadrature.cfg";
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

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

    SECTION( "Quadrature weights sum to 4*pi.", "[quadrature]" ) {

        bool lowAccuracyTesting = false;
        for( auto quadraturename : quadraturenames ) {
            // Set quadName
            config->SetQuadName( quadraturename );

            lowAccuracyTesting = false;
            if( quadraturename == QUAD_GaussLegendreTensorized || quadraturename == QUAD_GaussLegendre1D || quadraturename == QUAD_LevelSymmetric ||
                quadraturename == QUAD_Lebedev || quadraturename == QUAD_LDFESA )
                lowAccuracyTesting = true;

            for( auto quadratureorder : quadratureorders[quadraturename] ) {
                // Set quadOrder
                config->SetQuadOrder( quadratureorder );

55
                QuadratureBase* Q = QuadratureBase::Create( config );
56
57
58
59
60
61
62
63
64
65
66

                if( quadraturename == QUAD_GaussLegendre1D ) {
                    if( !approxequal( Q->SumUpWeights(), 2, lowAccuracyTesting ) ) {
                        printf( "Quadrature %d at order %d . Error : %.15f  (low accuracy testing was set to %d) ",
                                config->GetQuadName(),
                                config->GetQuadOrder(),
                                std::abs( Q->SumUpWeights() - 2 ),
                                lowAccuracyTesting );
                        printf( "Computed result %.15f \n", Q->SumUpWeights() );
                    }
                    REQUIRE( approxequal( Q->SumUpWeights(), 2, lowAccuracyTesting ) );    // 1D special case
steffen.schotthoefer's avatar
steffen.schotthoefer committed
67
                }
68
69
70
71
72
73
74
75
76
77
78
79
80
81
                else {
                    if( !approxequal( Q->SumUpWeights(), 4 * M_PI, lowAccuracyTesting ) ) {
                        printf( "Quadrature %d at order %d . Error : %.15f  (low accuracy testing was set to %d) ",
                                config->GetQuadName(),
                                config->GetQuadOrder(),
                                std::abs( Q->SumUpWeights() - 4 * M_PI ),
                                lowAccuracyTesting );
                        printf( "Computed result %.15f \n", Q->SumUpWeights() );
                    }
                    REQUIRE( approxequal( Q->SumUpWeights(), 4 * M_PI, lowAccuracyTesting ) );
                }
                // Special case for Gauss Legendre with half weights
                if( quadraturename == QUAD_GaussLegendreTensorized ) {
                    config->SetSNAllGaussPts( false );
82
                    QuadratureBase* Q = QuadratureBase::Create( config );
83
84
85
86
87
88
89
90
91
92
93
                    if( !approxequal( Q->SumUpWeights(), 4 * M_PI, lowAccuracyTesting ) ) {
                        printf( "Quadrature %d at order %d . Error : %.15f  (low accuracy testing was set to %d). Reduced number of quadrature "
                                "points used. \n",
                                config->GetQuadName(),
                                config->GetQuadOrder(),
                                std::abs( Q->SumUpWeights() - 4 * M_PI ),
                                lowAccuracyTesting );
                        printf( "Computed result %.15f \n", Q->SumUpWeights() );
                    }
                    REQUIRE( approxequal( Q->SumUpWeights(), 4 * M_PI, lowAccuracyTesting ) );
                    config->SetSNAllGaussPts( true );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
94
                }
95
96
97
98
            }
        }
    }

99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
    SECTION( "Quadrature points are on the unit sphere.", "[quadrature]" ) {

        bool lowAccuracyTesting = false;
        for( auto quadraturename : quadraturenames ) {
            // Set quadName
            config->SetQuadName( quadraturename );

            lowAccuracyTesting = false;
            if( quadraturename == QUAD_GaussLegendreTensorized || quadraturename == QUAD_GaussLegendre1D || quadraturename == QUAD_LevelSymmetric ||
                quadraturename == QUAD_Lebedev || quadraturename == QUAD_LDFESA )
                lowAccuracyTesting = true;

            if( quadraturename == QUAD_GaussLegendre1D ) continue;    // 1D test case not meaningful here

            for( auto quadratureorder : quadratureorders[quadraturename] ) {
                // Set quadOrder
                config->SetQuadOrder( quadratureorder );

117
118
                bool errorWithinBounds = true;

119
                QuadratureBase* Q   = QuadratureBase::Create( config );
120
121
122
123
124
125
126
127
128
129
130
                VectorVector points = Q->GetPoints();
                for( unsigned i = 0; i < Q->GetNq(); i++ ) {
                    if( !approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) ) {
                        printf( "Quadrature %d at order %d . Errorous index: %d | Error : %.15f  (low accuracy testing was set to %d). \n",
                                config->GetQuadName(),
                                config->GetQuadOrder(),
                                i,
                                std::abs( norm( points[i] ) - 1.0 ),
                                lowAccuracyTesting );
                        printf( "Computed result %.15f", norm( points[i] ) );
                    }
131
                    if( !approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) ) errorWithinBounds = false;
132
133
134
135
136
                }

                // Special case for Gauss Legendre with half weights
                if( quadraturename == QUAD_GaussLegendreTensorized ) {
                    config->SetSNAllGaussPts( false );
137
                    QuadratureBase* Q = QuadratureBase::Create( config );
138
139
140
141
142
143
144
145
146
147
148
149
150
151

                    VectorVector points = Q->GetPoints();
                    for( unsigned i = 0; i < Q->GetNq(); i++ ) {
                        if( !approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) ) {
                            printf( "Quadrature %d at order %d . Errorous index: %d | Error : %.15f  (low accuracy testing was set to %d). Reduced "
                                    "number of quadrature "
                                    "points used. \n",
                                    config->GetQuadName(),
                                    config->GetQuadOrder(),
                                    i,
                                    std::abs( norm( points[i] ) - 1.0 ),
                                    lowAccuracyTesting );
                            printf( "Computed result %.15f", norm( points[i] ) );
                        }
152
                        if( !approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) ) errorWithinBounds = false;
153
154
155
                    }

                    config->SetSNAllGaussPts( true );
156
                }
157
                REQUIRE( errorWithinBounds );
158
159
160
161
            }
        }
    }

162
163
164
165
166
167
168
169
170
171
    SECTION( "Nq is actually equal to the number of weights.", "[quadrature]" ) {

        for( auto quadraturename : quadraturenames ) {
            // Set quadName
            config->SetQuadName( quadraturename );

            for( auto quadratureorder : quadratureorders[quadraturename] ) {
                // Set quadOrder
                config->SetQuadOrder( quadratureorder );

172
                QuadratureBase* Q = QuadratureBase::Create( config );
173
174
175
176
177
178
                REQUIRE( Q->GetNq() == size( Q->GetWeights() ) );

                // Special case for Gauss Legendre with half weights

                if( quadraturename == QUAD_GaussLegendreTensorized ) {
                    config->SetSNAllGaussPts( false );
179
                    QuadratureBase* Q = QuadratureBase::Create( config );
180
181
182
183
                    REQUIRE( Q->GetNq() == size( Q->GetWeights() ) );
                    config->SetSNAllGaussPts( true );
                }
            }
184
185
186
        }
    }

187
188
189
190
191
192
193
194
195
196
    SECTION( "Nq is actually equal to the number of points.", "[quadrature]" ) {

        for( auto quadraturename : quadraturenames ) {
            // Set quadName
            config->SetQuadName( quadraturename );

            for( auto quadratureorder : quadratureorders[quadraturename] ) {
                // Set quadOrder
                config->SetQuadOrder( quadratureorder );

197
                QuadratureBase* Q = QuadratureBase::Create( config );
198
199
200
201
202
                REQUIRE( Q->GetNq() == size( Q->GetPoints() ) );

                // Special case for Gauss Legendre with half weights
                if( quadraturename == QUAD_GaussLegendreTensorized ) {
                    config->SetSNAllGaussPts( false );
203
                    QuadratureBase* Q = QuadratureBase::Create( config );
204
205
206
207
                    REQUIRE( Q->GetNq() == size( Q->GetPoints() ) );
                    config->SetSNAllGaussPts( true );
                }
            }
208
209
210
        }
    }

211
    SECTION( "Integrate a constant function.", "[quadrature]" ) {
212

213
214
215
216
        bool lowAccuracyTesting = false;
        for( auto quadraturename : quadraturenames ) {
            // Set quadName
            config->SetQuadName( quadraturename );
217

218
219
220
221
222
223
224
225
226
            lowAccuracyTesting = false;
            if( quadraturename == QUAD_GaussLegendreTensorized || quadraturename == QUAD_GaussLegendre1D || quadraturename == QUAD_LevelSymmetric ||
                quadraturename == QUAD_Lebedev || quadraturename == QUAD_LDFESA )
                lowAccuracyTesting = true;

            for( auto quadratureorder : quadratureorders[quadraturename] ) {
                // Set quadOrder
                config->SetQuadOrder( quadratureorder );

227
                QuadratureBase* Q = QuadratureBase::Create( config );
228
229
230
231
232
233
234
235
236
237
238

                if( quadraturename == QUAD_GaussLegendre1D ) {
                    if( !approxequal( Q->Integrate( sin ), 0, lowAccuracyTesting ) ) {
                        printf( "Quadrature %d at order %d :  Error : %.15f (low accuracy testing was set to %d)\n",
                                config->GetQuadName(),
                                config->GetQuadOrder(),
                                std::abs( Q->Integrate( sin ) ),
                                lowAccuracyTesting );
                        printf( "Computed result %.15f", Q->Integrate( sin ) );
                    }
                    REQUIRE( approxequal( Q->Integrate( sin ), 0, lowAccuracyTesting ) );    // 1D special case
steffen.schotthoefer's avatar
steffen.schotthoefer committed
239
                }
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
                else {
                    if( !approxequal( Q->Integrate( f ), 4.0 * M_PI, lowAccuracyTesting ) ) {
                        printf( "Quadrature %d at order %d :  Error : %.15f (low accuracy testing was set to %d)\n",
                                config->GetQuadName(),
                                config->GetQuadOrder(),
                                std::abs( Q->Integrate( f ) - 4.0 * M_PI ),
                                lowAccuracyTesting );
                        printf( "Computed result %.15f", Q->Integrate( f ) );
                    }
                    REQUIRE( approxequal( Q->Integrate( f ), 4.0 * M_PI, lowAccuracyTesting ) );
                }

                // Special case for Gauss Legendre with half weights
                if( quadraturename == QUAD_GaussLegendreTensorized ) {
                    config->SetSNAllGaussPts( false );
255
                    QuadratureBase* Q = QuadratureBase::Create( config );
256
257
258
259
260
261

                    if( !approxequal( Q->Integrate( f ), 4.0 * M_PI, lowAccuracyTesting ) ) {
                        printf(
                            "Quadrature %d at order %d :  Error : %.15f (low accuracy testing was set to %d). Reduced number of quad points used\n",
                            config->GetQuadName(),
                            config->GetQuadOrder(),
steffen.schotthoefer's avatar
steffen.schotthoefer committed
262
263
                            std::abs( Q->Integrate( f ) - 4.0 * M_PI ),
                            lowAccuracyTesting );
264
265
266
267
268
                        printf( "Computed result %.15f", Q->Integrate( f ) );
                    }
                    REQUIRE( approxequal( Q->Integrate( f ), 4.0 * M_PI, lowAccuracyTesting ) );

                    config->SetSNAllGaussPts( true );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
269
                }
270
271
272
273
            }
        }
    }
}