test_quadrature.cpp 7.32 KB
Newer Older
1
2
3
4
5
6
#include "catch.hpp"
#include "quadratures/quadraturebase.h"
#include "settings/globalconstants.h"

#include <vector>

steffen.schotthoefer's avatar
steffen.schotthoefer committed
7
8
std::vector<QUAD_NAME> quadraturenames = {
    QUAD_MonteCarlo, QUAD_GaussLegendreTensorized, QUAD_GaussLegendre1D, QUAD_LevelSymmetric, QUAD_Lebedev, QUAD_LDFESA };
9
10
std::vector<std::vector<int>> quadratureorders = {
    { 4, 5, 6, 7 },                            // Monte Carlo
11
    { 4, 6, 8, 10 },                           // Gauss Legendre
steffen.schotthoefer's avatar
steffen.schotthoefer committed
12
    { 4, 6, 8, 10 },                           // Gauss Legendre 1D
13
14
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 ) {
    double tol = 1e-15;              // For computed quadrature weights
    if( lowAccuracy ) tol = 5e-5;    // Mainly for Lookup Quadratures
steffen.schotthoefer's avatar
steffen.schotthoefer committed
22
    return std::abs( a - b ) < tol;
23
24
25
26
27
}

TEST_CASE( "Quadrature weights sum to 4*pi.", "[quadrature]" ) {
    bool lowAccuracyTesting = false;
    for( auto quadraturename : quadraturenames ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
28

29
        lowAccuracyTesting = false;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
30
31
        if( quadraturename == QUAD_GaussLegendreTensorized || quadraturename == QUAD_GaussLegendre1D || quadraturename == QUAD_LevelSymmetric ||
            quadraturename == QUAD_Lebedev || quadraturename == QUAD_LDFESA )
32
            lowAccuracyTesting = true;
33
34
35

        for( auto quadratureorder : quadratureorders[quadraturename] ) {
            QuadratureBase* Q = QuadratureBase::CreateQuadrature( quadraturename, quadratureorder );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

            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) ",
                            quadraturename,
                            quadratureorder,
                            std::abs( Q->SumUpWeights() - 2 ),
                            lowAccuracyTesting );
                    printf( "Computed result %.15f \n", Q->SumUpWeights() );
                }
                REQUIRE( approxequal( Q->SumUpWeights(), 2, lowAccuracyTesting ) );    // 1D special case
            }
            else {
                if( !approxequal( Q->SumUpWeights(), 4 * M_PI, lowAccuracyTesting ) ) {
                    printf( "Quadrature %d at order %d . Error : %.15f  (low accuracy testing was set to %d) ",
                            quadraturename,
                            quadratureorder,
                            std::abs( Q->SumUpWeights() - 4 * M_PI ),
                            lowAccuracyTesting );
                    printf( "Computed result %.15f \n", Q->SumUpWeights() );
                }
                REQUIRE( approxequal( Q->SumUpWeights(), 4 * M_PI, lowAccuracyTesting ) );
58
59
60
61
62
63
64
65
66
            }
        }
    }
}

TEST_CASE( "Quadrature points are on the unit sphere.", "[quadrature]" ) {
    bool lowAccuracyTesting = false;
    for( auto quadraturename : quadraturenames ) {
        lowAccuracyTesting = false;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
67
68
        if( quadraturename == QUAD_GaussLegendreTensorized || quadraturename == QUAD_GaussLegendre1D || quadraturename == QUAD_LevelSymmetric ||
            quadraturename == QUAD_Lebedev || quadraturename == QUAD_LDFESA )
69
            lowAccuracyTesting = true;
70

steffen.schotthoefer's avatar
steffen.schotthoefer committed
71
72
        if( quadraturename == QUAD_GaussLegendre1D ) continue;    // 1D test case not meaningful here

73
74
75
76
77
78
79
80
81
        for( auto quadratureorder : quadratureorders[quadraturename] ) {
            QuadratureBase* Q   = QuadratureBase::CreateQuadrature( quadraturename, quadratureorder );
            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",
                            quadraturename,
                            quadratureorder,
                            i,
steffen.schotthoefer's avatar
steffen.schotthoefer committed
82
                            std::abs( norm( points[i] ) - 1.0 ),
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
                            lowAccuracyTesting );
                    printf( "Computed result %.15f", norm( points[i] ) );
                }
                REQUIRE( approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) );
            }
        }
    }
}

TEST_CASE( "Nq is actually equal to the number of weights.", "[quadrature]" ) {
    for( auto quadraturename : quadraturenames ) {
        for( auto quadratureorder : quadratureorders[quadraturename] ) {
            QuadratureBase* Q = QuadratureBase::CreateQuadrature( quadraturename, quadratureorder );
            REQUIRE( Q->GetNq() == size( Q->GetWeights() ) );
        }
    }
}

TEST_CASE( "Nq is actually equal to the number of points.", "[quadrature]" ) {
    for( auto quadraturename : quadraturenames ) {
        for( auto quadratureorder : quadratureorders[quadraturename] ) {
            QuadratureBase* Q = QuadratureBase::CreateQuadrature( quadraturename, quadratureorder );
            REQUIRE( Q->GetNq() == size( Q->GetPoints() ) );
        }
    }
}

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

steffen.schotthoefer's avatar
steffen.schotthoefer committed
114
115
double sin( double x, double y, double z ) { return sin( x ); }

116
117
118
119
TEST_CASE( "Integrate a constant function.", "[quadrature]" ) {
    bool lowAccuracyTesting = false;
    for( auto quadraturename : quadraturenames ) {
        lowAccuracyTesting = false;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
120
121
        if( quadraturename == QUAD_GaussLegendreTensorized || quadraturename == QUAD_GaussLegendre1D || quadraturename == QUAD_LevelSymmetric ||
            quadraturename == QUAD_Lebedev || quadraturename == QUAD_LDFESA )
122
            lowAccuracyTesting = true;
123
124
125

        for( auto quadratureorder : quadratureorders[quadraturename] ) {
            QuadratureBase* Q = QuadratureBase::CreateQuadrature( quadraturename, quadratureorder );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147

            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",
                            quadraturename,
                            quadratureorder,
                            std::abs( Q->Integrate( sin ) ),
                            lowAccuracyTesting );
                    printf( "Computed result %.15f", Q->Integrate( sin ) );
                }
                REQUIRE( approxequal( Q->Integrate( sin ), 0, lowAccuracyTesting ) );    // 1D special case
            }
            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",
                            quadraturename,
                            quadratureorder,
                            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 ) );
148
149
150
151
            }
        }
    }
}