test_numericalflux.cpp 4.52 KB
Newer Older
Jonas Kusch's avatar
Jonas Kusch committed
1
2
3
#include <numeric>

#include "catch.hpp"
4
#include "common/typedef.h"
5
#include "fluxes/upwindflux.h"
Jonas Kusch's avatar
Jonas Kusch committed
6

Jonas Kusch's avatar
Jonas Kusch committed
7
TEST_CASE( "unit numericalflux tests", "[numericalflux]" ) {
Jonas Kusch's avatar
Jonas Kusch committed
8

9
    // Test Flux Jacobians
10
11
12
13
    Matrix AxP( 4, 4, 0 );
    Matrix AxM( 4, 4, 0 );
    Matrix AyP( 4, 4, 0 );
    Matrix AyM( 4, 4, 0 );
14
15

    // Fill Matrices with used values
16
17
18
19
    AxP( 0, 0 ) = 1;    // 0.5;
    AxP( 1, 1 ) = 0;    // 0.25;
    AyP( 0, 0 ) = 0;    // 0.25;
    AyP( 1, 1 ) = 1;    // 0.125;
20

21
22
23
24
    AxM( 2, 2 ) = -1;    // -0.25;
    AxM( 3, 3 ) = 0;     // -0.125;
    AyM( 2, 2 ) = 0;     // -0.5;
    AyM( 3, 3 ) = -1;    // -0.25;
Jonas Kusch's avatar
Jonas Kusch committed
25
26

    SECTION( "test symmetry" ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
27
        UpwindFlux g;
28
29
30
31
32
        Vector resultFluxPlus( 4, 0.0 );
        Vector resultFluxMinus( 4, 0.0 );
        Vector psiLVect{ 1.0, 2.0, 3.0, 4.0 };
        Vector psiRVect{ 3.0, 4.0, 1.0, 2.0 };

33
34
35
        bool fluxCancelingScalar = true;
        bool fluxCancelingMatrix = true;

Jonas Kusch's avatar
Jonas Kusch committed
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
        for( unsigned l = 0; l < 10; ++l ) {
            Vector n( 2, 0.0 );
            std::default_random_engine generator;
            std::normal_distribution<double> distribution( -10.0, 10.0 );
            // sample normal
            n[0] = distribution( generator );
            n[1] = distribution( generator );
            // sample edge length
            double length = distribution( generator );

            // scale normalized normal
            n = length * n / sqrt( n[0] * n[0] + n[1] * n[1] );

            // sample solution at neighboring cells
            double psiL = distribution( generator ) + 10.0;
            double psiR = distribution( generator ) + 10.0;

            // sample omega (note that omega now does not need to lie on unit sphere)
            Vector omega( 2, 0.0 );
            std::normal_distribution<double> uni( -1.0, 1.0 );
            omega[0] = uni( generator );
            omega[1] = uni( generator );

59
60
61
            // ---- Test all fluxes ----

            // a) Upwinding methods - scalar
62
63
            if( std::fabs( g.Flux( omega, psiL, psiR, n ) - ( -g.Flux( omega, psiR, psiL, -n ) ) ) > 1e2 * std::numeric_limits<double>::epsilon() )
                fluxCancelingScalar = false;
64

65
66
67
68
            // b) Upwinding methods - Systems  MatrixFlux, 4x4 Matrices, i.e. P1 case;

            resultFluxPlus  = g.Flux( AxP, AxM, AyP, AyM, AyP, AyM, psiLVect, psiRVect, n );
            resultFluxMinus = g.Flux( AxP, AxM, AyP, AyM, AyP, AyM, psiRVect, psiLVect, -n );
69

70
            if( blaze::norm( resultFluxPlus + resultFluxMinus ) > 1e2 * std::numeric_limits<double>::epsilon() ) fluxCancelingMatrix = false;
71
        }
72
73
        REQUIRE( fluxCancelingScalar );
        REQUIRE( fluxCancelingMatrix );
Jonas Kusch's avatar
Jonas Kusch committed
74
75
76
    }

    SECTION( "test consistency" ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
77
        UpwindFlux g;
78
79
80
81
82
        Vector n( 2, 0.0 );
        Vector resultFlux( 4, 0 );
        Vector resultFluxOrig( 4, 0 );
        Vector psiLVect{ 1.0, 2.0, 3.0, 4.0 };

83
84
85
        bool consitencyScalar = true;
        bool consitencyMatrix = true;

Jonas Kusch's avatar
Jonas Kusch committed
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
        for( unsigned l = 0; l < 10; ++l ) {
            std::default_random_engine generator;
            std::normal_distribution<double> distribution( -10.0, 10.0 );
            // sample normal
            n[0] = distribution( generator );
            n[1] = distribution( generator );
            // sample edge length
            double length = distribution( generator );

            // scale normalized normal
            n = length * n / sqrt( n[0] * n[0] + n[1] * n[1] );

            // sample solution at neighboring cells
            double psi = distribution( generator ) + 10.0;

            // sample omega (note that omega now does not need to lie on unit sphere)
            Vector omega( 2, 0.0 );
            std::normal_distribution<double> uni( -1.0, 1.0 );
            omega[0] = uni( generator );
            omega[1] = uni( generator );

            double inner = omega[0] * n[0] + omega[1] * n[1];

109
110
111
112
            // ---- Test all fluxes ----

            // a) Upwinding methods - scalar

113
            if( std::fabs( g.Flux( omega, psi, psi, n ) - inner * psi ) > 1e2 * std::numeric_limits<double>::epsilon() ) consitencyScalar = false;
114

115
            // b) Upwinding methods -  MatrixFlux, 4x4 Matrices, i.e. P1 case;
116

117
118
            resultFluxOrig = ( n[0] * ( AxP + AxM ) * psiLVect + +n[1] * ( AyP + AyM ) * psiLVect );
            resultFlux     = g.Flux( AxP, AxM, AyP, AyM, AyP, AyM, psiLVect, psiLVect, n );
119

120
            if( blaze::norm( resultFlux - resultFluxOrig ) > 1e2 * std::numeric_limits<double>::epsilon() ) consitencyMatrix = false;
121
        }
122
123
        REQUIRE( consitencyScalar );
        REQUIRE( consitencyMatrix );
Jonas Kusch's avatar
Jonas Kusch committed
124
125
    }
}