linesource.cpp 5.86 KB
Newer Older
1
#include "problems/linesource.h"
steffen.schotthoefer's avatar
steffen.schotthoefer committed
2
#include "mesh.h"
3
#include "settings/config.h"
4

5
// ---- LineSource_SN ----
6

7
LineSource_SN::LineSource_SN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }
8

9
10
11
LineSource_SN::~LineSource_SN(){};

VectorVector LineSource_SN::GetScatteringXS( const std::vector<double>& energies ) {
jannick.wolters's avatar
jannick.wolters committed
12
    return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) );
13
14
}

15
VectorVector LineSource_SN::GetTotalXS( const std::vector<double>& energies ) {
jannick.wolters's avatar
jannick.wolters committed
16
    return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) );
17
}
18

19
std::vector<VectorVector> LineSource_SN::GetExternalSource( const std::vector<double>& energies ) {
20
    return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
jannick.wolters's avatar
jannick.wolters committed
21
22
}

23
std::vector<double> LineSource_SN::GetStoppingPower( const std::vector<double>& energies ) {
24
    // @TODO
25
    return std::vector<double>( energies.size(), 1.0 );
26
27
}

28
VectorVector LineSource_SN::SetupIC() {
29
    VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
jannick.wolters's avatar
jannick.wolters committed
30
    auto cellMids = _mesh->GetCellMidPoints();
jannick.wolters's avatar
jannick.wolters committed
31
    double t      = 3.2e-4;    // pseudo time for gaussian smoothing
32
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
jannick.wolters's avatar
jannick.wolters committed
33
34
        double x = cellMids[j][0];
        double y = cellMids[j][1];
steffen.schotthoefer's avatar
steffen.schotthoefer committed
35
        psi[j]   = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) );
36
37
38
    }
    return psi;
}
39

40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
// ---- LineSource_SN_peudo1D ----

LineSource_SN_Pseudo1D::LineSource_SN_Pseudo1D( Config* settings, Mesh* mesh ) : LineSource_SN( settings, mesh ) {}

VectorVector LineSource_SN_Pseudo1D::SetupIC() {
    VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
    auto cellMids = _mesh->GetCellMidPoints();
    double t      = 3.2e-4;    // pseudo time for gaussian smoothing
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
        double x = cellMids[j][0];
        psi[j]   = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x ) / ( 4 * t ) );
    }
    return psi;
}

55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
// ---- LineSource_PN ----

int LineSource_PN::GlobalIndex( int l, int k ) const {
    int numIndicesPrevLevel  = l * l;    // number of previous indices untill level l-1
    int prevIndicesThisLevel = k + l;    // number of previous indices in current level
    return numIndicesPrevLevel + prevIndicesThisLevel;
}

LineSource_PN::LineSource_PN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }

LineSource_PN::~LineSource_PN(){};

VectorVector LineSource_PN::GetScatteringXS( const std::vector<double>& energies ) {
    return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) );
}

VectorVector LineSource_PN::GetTotalXS( const std::vector<double>& energies ) {
    return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) );
}

std::vector<VectorVector> LineSource_PN::GetExternalSource( const std::vector<double>& energies ) {
    return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
}

std::vector<double> LineSource_PN::GetStoppingPower( const std::vector<double>& energies ) {
    // @TODO
    return std::vector<double>( energies.size(), 0.0 );
}

VectorVector LineSource_PN::SetupIC() {
    // Compute number of equations in the system
86
    int ntotalEquations = GlobalIndex( _settings->GetMaxMomentDegree(), _settings->GetMaxMomentDegree() ) + 1;
87
88
89
90
91

    VectorVector psi( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) );    // zero could lead to problems?
    VectorVector cellMids = _mesh->GetCellMidPoints();

    // Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments are zero.
92
93
    double t      = 3.2e-4;    // pseudo time for gaussian smoothing (Approx to dirac git impulse)
    double offset = 0.0;
94
    for( unsigned j = 0; j < cellMids.size(); ++j ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
95
96
97
        double x = cellMids[j][0];
        double y = cellMids[j][1];    // (x- 0.5) * (x- 0.5)

98
99
        psi[j][0] =
            1.0 / ( 4.0 * M_PI * t ) *
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
            std::exp(
                -( ( x - offset ) * ( x - offset ) + ( y - offset ) * ( y - offset ) ) /
                ( 4 * t ) );    // +
                                // 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x - offset ) * ( x - offset ) + ( y - offset ) * ( y - offset ) ) / ( 4 *
                                // t ) ) + 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + offset ) * ( x + offset ) + ( y - offset ) * ( y - offset ) )
                                // / ( 4 * t ) ) + 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x - offset ) * ( x - offset ) + ( y + offset ) * ( y +
                                // offset ) ) / ( 4 * t ) ) + 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + offset ) * ( x + offset ) + ( y + offset )
                                // * ( y + offset ) ) / ( 4 * t ) ) + 1.0 / ( 4.0 * M_PI * t ) *
                                //    std::exp( -( ( x - 2 * offset ) * ( x - 2 * offset ) + ( y - 2 * offset ) * ( y - 2 * offset ) ) / ( 4 * t ) ) +
                                // 1.0 / ( 4.0 * M_PI * t ) *
                                //    std::exp( -( ( x + 2 * offset ) * ( x + 2 * offset ) + ( y - 2 * offset ) * ( y - 2 * offset ) ) / ( 4 * t ) ) +
                                // 1.0 / ( 4.0 * M_PI * t ) *
                                //    std::exp( -( ( x - 2 * offset ) * ( x - 2 * offset ) + ( y + 2 * offset ) * ( y + 2 * offset ) ) / ( 4 * t ) ) +
                                // 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + 2 * offset ) * ( x + 2 * offset ) + ( y + 2 * offset ) * ( y + 2 *
                                // offset ) ) / ( 4 * t ) );
115
    }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
116

117
    // Debugging jump test case
steffen.schotthoefer's avatar
steffen.schotthoefer committed
118
119
120
121
122
123
    // for( unsigned j = 0; j < cellMids.size(); ++j ) {
    //    if( cellMids[j][0] > -0.2 && cellMids[j][1] > -0.2 && cellMids[j][1] < 0.2 && cellMids[j][0] < 0.2 )
    //        psi[j][0] = 1.0;
    //    else
    //        psi[j][0] = 0.0;
    //}
steffen.schotthoefer's avatar
steffen.schotthoefer committed
124
    return psi;
125
}