mnsolver.cpp 8.29 KB
Newer Older
1
#include "solvers/mnsolver.h"
2
#include "toolboxes/textprocessingtoolbox.h"
3
//#include <chrono>
4
#include <mpi.h>
steffen.schotthoefer's avatar
steffen.schotthoefer committed
5

6
MNSolver::MNSolver( Config* settings ) : Solver( settings ), _nMaxMomentsOrder( settings->GetMaxMomentDegree() ), _basis( _nMaxMomentsOrder ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
7
8
    // Is this good (fast) code using a constructor list?

9
    _nTotalEntries = GlobalIndex( _nMaxMomentsOrder, int( _nMaxMomentsOrder ) ) + 1;
10
    _quadrature    = QuadratureBase::CreateQuadrature( _settings->GetQuadName(), settings->GetQuadOrder() );
11

steffen.schotthoefer's avatar
steffen.schotthoefer committed
12
13
14
15
16
17
18
19
20
21
22
23
24
25
    // transform sigmaT and sigmaS in sigmaA.
    _sigmaA = VectorVector( _nEnergies, Vector( _nCells, 0 ) );    // Get rid of this extra vektor!

    for( unsigned n = 0; n < _nEnergies; n++ ) {
        for( unsigned j = 0; j < _nCells; j++ ) {
            _sigmaA[n][j] = 0;    //_sigmaT[n][j] - _sigmaS[n][j];
            _sigmaS[n][j] = 1;
        }
    }

    // Initialize Scatter Matrix
    _scatterMatDiag    = Vector( _nTotalEntries, 1.0 );
    _scatterMatDiag[0] = 0.0;    // First entry is zero by construction.

26
27
28
29
30
    // Initialize Entropy
    _entropy = EntropyBase::Create( _settings );

    // Initialize Optimizer
    _optimizer = OptimizerBase::Create( _settings );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
31

32
33
    // Initialize lagrange Multiplier
    _alpha = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
34
35

    // Initialize and Pre-Compute Moments at quadrature points
36
    _moments = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) );
37
    ComputeMoments();
steffen.schotthoefer's avatar
steffen.schotthoefer committed
38
39
}

40
41
42
43
44
MNSolver::~MNSolver() {
    delete _quadrature;
    delete _entropy;
    delete _optimizer;
}
45

steffen.schotthoefer's avatar
steffen.schotthoefer committed
46
47
48
49
50
51
int MNSolver::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;
}

52
void MNSolver::ComputeMoments() {
53
    double my, phi, w;
54
55

    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
56
57
58
59
        phi = _quadrature->GetPointsSphere()[idx_quad][0];
        my  = _quadrature->GetPointsSphere()[idx_quad][1];

        _moments[idx_quad] = _basis.ComputeSphericalBasis( my, phi );
60
61
    }

62
    Vector normalizationFactor( _moments[0].size(), 0.0 );
63
64
65

    // Normalize Moments

66
67
68
69
70
71
72
73
74
75
76
77
78
79
    // for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
    //    w = _quadrature->GetWeights()[idx_quad];
    //
    //    for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
    //        normalizationFactor[idx_sys] += w * _moments[idx_quad][idx_sys] * _moments[idx_quad][idx_sys];
    //    }
    //}
    // std::cout << "norm" << normalizationFactor << "\n";
    //
    // for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
    //    for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
    //        _moments[idx_quad][idx_sys] /= sqrt( normalizationFactor[idx_sys] );
    //    }
    //}
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
}

Vector MNSolver::ConstructFlux( unsigned idx_cell ) {

    // ---- Integration of Moment of flux ----
    double w, entropyL, entropyR, entropyFlux;
    // double x, y, z;
    // std::cout << "--------\n";
    // std::cout << " alphas curr cell" << _alpha[idx_cell] << " alphas neigh cell:\n";

    bool good = _alpha[idx_cell] == 0.0;

    int idx_good = 0;
    for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
        // std::cout << _alpha[idx_neigh] << "\n";
        if( _alpha[idx_neigh] == 0.0 ) idx_good++;
    }
    if( idx_good == 2 && good ) {
        std::cout << "Candidate\n";
    }
100

101
    Vector flux( _nTotalEntries, 0.0 );
102
    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

        w = _quadrature->GetWeights()[idx_quad];

        entropyFlux = 0.0;    // Reset temorary flux

        entropyL = dot( _alpha[idx_cell], _moments[idx_quad] );
        // _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _moments[idx_quad] ) );

        for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
            // Store fluxes in psiNew, to save memory
            if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neigh] == _nCells )
                entropyR = entropyL;
            else {
                entropyR = dot( _alpha[idx_neigh], _moments[idx_quad] );
                //_entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_neigh], _moments[idx_quad] ) );
            }
            entropyFlux += _g->Flux( _quadrature->GetPoints()[idx_quad], entropyL, entropyR, _normals[idx_cell][idx_neigh] );
120
        }
121
122
123
124

        //   std::cout << "\n entropy Flux at cell [" << idx_cell << "] and quad [" << idx_quad << "]: " << entropyFlux << "\n";
        // integrate
        flux += _moments[idx_quad] * ( w * entropyFlux );
125
    }
126
127
128
129

    // std::cout << "\n integrate entropy Flux at cell [" << idx_cell << "] and : " << flux << "\n";

    return flux;
130
131
}

steffen.schotthoefer's avatar
steffen.schotthoefer committed
132
133
134
135
136
137
138
139
140
141
142
143
void MNSolver::Solve() {

    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );

    auto log = spdlog::get( "event" );

    // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
    VectorVector psiNew = _psi;
    double dFlux        = 1e10;
    Vector fluxNew( _nCells, 0.0 );
    Vector fluxOld( _nCells, 0.0 );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
144
145

    double mass1 = 0;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
146
147
    for( unsigned i = 0; i < _nCells; ++i ) {
        _solverOutput[i] = _psi[i][0];
steffen.schotthoefer's avatar
steffen.schotthoefer committed
148
        mass1 += _psi[i][0] * _areas[i];
steffen.schotthoefer's avatar
steffen.schotthoefer committed
149
150
    }

151
152
153
154
155
    dFlux   = blaze::l2Norm( fluxNew - fluxOld );
    fluxOld = fluxNew;

    Save( -1 );

steffen.schotthoefer's avatar
steffen.schotthoefer committed
156
157
158
    if( rank == 0 ) log->info( "{:10}   {:10}", "t", "dFlux" );
    if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", -1.0, dFlux, mass1 );

159
160
161
162
    // Time measurement
    // auto start = chrono::steady_clock::now();
    // auto end   = chrono::steady_clock::now();

steffen.schotthoefer's avatar
steffen.schotthoefer committed
163
164
    // Loop over energies (pseudo-time of continuous slowing down approach)

165
    for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
166

167
168
        // Loop over the grid cells
        for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
169

170
171
            // ------- Reconstruction Step -------

172
            _alpha[idx_cell] = _psi[idx_cell];    // _optimizer->Solve( _psi[idx_cell] );
173

174
            // ------- Flux Computation Step ---------
175

176
177
            // Dirichlet Boundaries are finished now
            if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
178

179
            psiNew[idx_cell] = ConstructFlux( idx_cell );
180

181
            // ------ FVM Step ------
182

183
184
185
            // NEED TO VECTORIZE
            for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {

186
187
188
189
190
191
                psiNew[idx_cell][idx_system] =
                    _psi[idx_cell][idx_system] - ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system]; /* cell averaged flux */

                //- _dE * _psi[idx_cell][idx_system] *
                //      ( _sigmaA[idx_energy][idx_cell]                                    /* absorbtion influence */
                //        + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
192
            }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
193
        }
194
195
196
        _psi = psiNew;

        // pseudo time iteration output
197
198
199
200
201
        double mass = 0.0;
        for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
            fluxNew[idx_cell]       = _psi[idx_cell][0];    // zeroth moment is raditation densitiy we are interested in
            _solverOutput[idx_cell] = _psi[idx_cell][0];
            mass += _psi[idx_cell][0] * _areas[idx_cell];
steffen.schotthoefer's avatar
steffen.schotthoefer committed
202
        }
203

steffen.schotthoefer's avatar
steffen.schotthoefer committed
204
205
        dFlux   = blaze::l2Norm( fluxNew - fluxOld );
        fluxOld = fluxNew;
206
        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
207
        Save( idx_energy );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
208
209
    }
}
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229

void MNSolver::Save() const {
    std::vector<std::string> fieldNames{ "flux" };
    std::vector<double> flux;
    flux.resize( _nCells );

    for( unsigned i = 0; i < _nCells; ++i ) {
        flux[i] = _psi[i][0];
    }
    std::vector<std::vector<double>> scalarField( 1, flux );
    std::vector<std::vector<std::vector<double>>> results{ scalarField };
    ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
}

void MNSolver::Save( int currEnergy ) const {
    std::vector<std::string> fieldNames{ "flux" };
    std::vector<std::vector<double>> scalarField( 1, _solverOutput );
    std::vector<std::vector<std::vector<double>>> results{ scalarField };
    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
}