pnsolver.cpp 23.2 KB
Newer Older
steffen.schotthoefer's avatar
steffen.schotthoefer committed
1
#include "solvers/pnsolver.h"
2
3
#include "common/config.h"
#include "common/io.h"
4
#include "common/mesh.h"
5
#include "fluxes/numericalflux.h"
6
7
#include "toolboxes/errormessages.h"
#include "toolboxes/textprocessingtoolbox.h"
8

9
10
// externals
#include "spdlog/spdlog.h"
11
12
#include <mpi.h>

13
PNSolver::PNSolver( Config* settings ) : SolverBase( settings ) {
14
15
    _polyDegreeBasis = settings->GetMaxMomentDegree();
    _nSystem         = GlobalIndex( int( _polyDegreeBasis ), int( _polyDegreeBasis ) ) + 1;
16
17

    // Initialize System Matrices
18
19
20
21
22
23
24
25
26
27
28
29
30
    _Ax = SymMatrix( _nSystem );
    _Ay = SymMatrix( _nSystem );
    _Az = SymMatrix( _nSystem );

    _AxPlus  = Matrix( _nSystem, _nSystem, 0 );
    _AxMinus = Matrix( _nSystem, _nSystem, 0 );
    _AxAbs   = Matrix( _nSystem, _nSystem, 0 );
    _AyPlus  = Matrix( _nSystem, _nSystem, 0 );
    _AyMinus = Matrix( _nSystem, _nSystem, 0 );
    _AyAbs   = Matrix( _nSystem, _nSystem, 0 );
    _AzPlus  = Matrix( _nSystem, _nSystem, 0 );
    _AzMinus = Matrix( _nSystem, _nSystem, 0 );
    _AzAbs   = Matrix( _nSystem, _nSystem, 0 );
31

32
    // Initialize Scatter Matrix
33
    _scatterMatDiag = Vector( _nSystem, 0 );
34

35
    // Initialize temporary storages of solution derivatives
36
37
    _solDx = VectorVector( _nCells, Vector( _nSystem, 0.0 ) );
    _solDy = VectorVector( _nCells, Vector( _nSystem, 0.0 ) );
38

39
    // Fill System Matrices
40
    ComputeSystemMatrices();
41

42
    // Compute Decomposition in positive and negative (eigenvalue) parts of flux jacobians
43
    ComputeFluxComponents();
44
45
46
47

    // Compute diagonal of the scatter matrix (it's a diagonal matrix)
    ComputeScatterMatrix();

steffen.schotthoefer's avatar
steffen.schotthoefer committed
48
    // AdaptTimeStep();
49
50
51

    if( settings->GetCleanFluxMat() ) CleanFluxMatrices();

52
53
    // Compute moments of initial condition
    // TODO
54
}
55

56
void PNSolver::IterPreprocessing( unsigned /*idx_iter*/ ) {
57
58
    // Nothing to preprocess for PNSolver
}
59

60
void PNSolver::IterPostprocessing( unsigned /*idx_iter*/ ) {
61
62
    // --- Update Solution ---
    _sol = _solNew;
63

64
65
    // --- Compute Flux for solution and Screen Output ---
    ComputeRadFlux();
steffen.schotthoefer's avatar
steffen.schotthoefer committed
66
67
}

68
69
70
71
72
void PNSolver::ComputeRadFlux() {
    double firstMomentScaleFactor = sqrt( 4 * M_PI );
    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
        _fluxNew[idx_cell] = _sol[idx_cell][0] * firstMomentScaleFactor;
    }
73
74
}

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
75
void PNSolver::FluxUpdate() {
76
    if( _reconsOrder > 1 ) {
77
        _mesh->ReconstructSlopesU( _nSystem, _solDx, _solDy, _sol );    // unstructured reconstruction
78
79
80
81
82
83
84
        //_mesh->ComputeSlopes( _nTotalEntries, _solDx, _solDy, _sol );    // unstructured reconstruction
    }
    // Vector solL( _nTotalEntries );
    // Vector solR( _nTotalEntries );
    auto solL = _sol[2];
    auto solR = _sol[2];

85
86
    // Loop over all spatial cells
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
87
88
89

        // Dirichlet cells stay at IC, farfield assumption
        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
90
91

        // Reset temporary variable psiNew
92
        for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
93
            _solNew[idx_cell][idx_sys] = 0.0;
94
95
96
97
98
99
100
        }

        // Loop over all neighbor cells (edges) of cell j and compute numerical fluxes
        for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[idx_cell].size(); idx_neighbor++ ) {

            // Compute flux contribution and store in psiNew to save memory
            if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neighbor] == _nCells )
101
                _solNew[idx_cell] += _g->Flux(
102
                    _AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, _sol[idx_cell], _sol[idx_cell], _normals[idx_cell][idx_neighbor] );
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
            else {
                switch( _reconsOrder ) {
                    // first order solver
                    case 1:
                        _solNew[idx_cell] += _g->Flux( _AxPlus,
                                                       _AxMinus,
                                                       _AyPlus,
                                                       _AyMinus,
                                                       _AzPlus,
                                                       _AzMinus,
                                                       _sol[idx_cell],
                                                       _sol[_neighbors[idx_cell][idx_neighbor]],
                                                       _normals[idx_cell][idx_neighbor] );
                        break;
                    // second order solver
                    case 2:
                        // left status of interface
                        solL = _sol[idx_cell] + _solDx[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neighbor][0] - _cellMidPoints[idx_cell][0] ) +
                               _solDy[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neighbor][1] - _cellMidPoints[idx_cell][1] );
                        // right status of interface
                        solR = _sol[_neighbors[idx_cell][idx_neighbor]] +
                               _solDx[_neighbors[idx_cell][idx_neighbor]] *
                                   ( _interfaceMidPoints[idx_cell][idx_neighbor][0] - _cellMidPoints[_neighbors[idx_cell][idx_neighbor]][0] ) +
                               _solDy[_neighbors[idx_cell][idx_neighbor]] *
                                   ( _interfaceMidPoints[idx_cell][idx_neighbor][1] - _cellMidPoints[_neighbors[idx_cell][idx_neighbor]][1] );
                        // positivity checker (if not satisfied, deduce to first order)
                        // I manually turned it off here since Pn produces negative solutions essentially
                        // if( min(solL) < 0.0 || min(solR) < 0.0 ) {
                        //    solL = _sol[idx_cell];
                        //    solR = _sol[_neighbors[idx_cell][idx_neighbor]];
                        //}

                        // flux evaluation
                        _solNew[idx_cell] +=
                            _g->Flux( _AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, solL, solR, _normals[idx_cell][idx_neighbor] );
                        break;
                    // default: first order solver
                    default:
                        _solNew[idx_cell] += _g->Flux( _AxPlus,
                                                       _AxMinus,
                                                       _AyPlus,
                                                       _AyMinus,
                                                       _AzPlus,
                                                       _AzMinus,
                                                       _sol[idx_cell],
                                                       _sol[_neighbors[idx_cell][idx_neighbor]],
                                                       _normals[idx_cell][idx_neighbor] );
                }
            }
152
153
154
155
        }
    }
}

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
156
void PNSolver::FVMUpdate( unsigned idx_energy ) {
157
    // Loop over all spatial cells
158
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
159
160
161
        // Dirichlet cells stay at IC, farfield assumption
        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
        // Flux update
162
        for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
163
164
            _solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] - ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
                                         - _dE * _sol[idx_cell][idx_sys] *
Steffen Schotthöfer's avatar
update    
Steffen Schotthöfer committed
165
166
167
                                               ( _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] /* scattering influence */
                                                 + _sigmaT[idx_energy][idx_cell] );
            /* total xs influence  */    // Vorzeichenfehler!
168
        }
169
170
        // Source Term
        _solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
171
172
173
    }
}

174
void PNSolver::ComputeSystemMatrices() {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
175
176
    int idx_col      = 0;
    unsigned idx_row = 0;
177
178

    // loop over columns of A
179
    for( int idx_lOrder = 0; idx_lOrder <= int( _polyDegreeBasis ); idx_lOrder++ ) {     // index of legendre polynom
180
        for( int idx_kOrder = -idx_lOrder; idx_kOrder <= idx_lOrder; idx_kOrder++ ) {    // second index of legendre function
steffen.schotthoefer's avatar
steffen.schotthoefer committed
181
            idx_row = unsigned( GlobalIndex( idx_lOrder, idx_kOrder ) );
182
183

            // flux matrix in direction x
steffen.schotthoefer's avatar
steffen.schotthoefer committed
184
185
186
187
188
189
190
191
192
193
194
195
            {
                if( idx_kOrder != -1 ) {

                    if( CheckIndex( idx_lOrder - 1, kMinus( idx_kOrder ) ) ) {
                        idx_col                             = GlobalIndex( idx_lOrder - 1, kMinus( idx_kOrder ) );
                        _Ax( idx_row, unsigned( idx_col ) ) = 0.5 * CTilde( idx_lOrder - 1, std::abs( idx_kOrder ) - 1 );
                    }

                    if( CheckIndex( idx_lOrder + 1, kMinus( idx_kOrder ) ) ) {
                        idx_col                             = GlobalIndex( idx_lOrder + 1, kMinus( idx_kOrder ) );
                        _Ax( idx_row, unsigned( idx_col ) ) = -0.5 * DTilde( idx_lOrder + 1, std::abs( idx_kOrder ) - 1 );
                    }
196
197
                }

steffen.schotthoefer's avatar
steffen.schotthoefer committed
198
199
200
                if( CheckIndex( idx_lOrder - 1, kPlus( idx_kOrder ) ) ) {
                    idx_col                             = GlobalIndex( idx_lOrder - 1, kPlus( idx_kOrder ) );
                    _Ax( idx_row, unsigned( idx_col ) ) = -0.5 * ETilde( idx_lOrder - 1, std::abs( idx_kOrder ) + 1 );
201
202
                }

steffen.schotthoefer's avatar
steffen.schotthoefer committed
203
204
205
206
                if( CheckIndex( idx_lOrder + 1, kPlus( idx_kOrder ) ) ) {
                    idx_col                             = GlobalIndex( idx_lOrder + 1, kPlus( idx_kOrder ) );
                    _Ax( idx_row, unsigned( idx_col ) ) = 0.5 * FTilde( idx_lOrder + 1, std::abs( idx_kOrder ) + 1 );
                }
207
208
209
            }

            // flux matrix in direction y
steffen.schotthoefer's avatar
steffen.schotthoefer committed
210
211
212
213
214
215
216
217
218
219
220
            {
                if( idx_kOrder != 1 ) {
                    if( CheckIndex( idx_lOrder - 1, -kMinus( idx_kOrder ) ) ) {
                        idx_col                             = GlobalIndex( idx_lOrder - 1, -kMinus( idx_kOrder ) );
                        _Ay( idx_row, unsigned( idx_col ) ) = -0.5 * Sgn( idx_kOrder ) * CTilde( idx_lOrder - 1, std::abs( idx_kOrder ) - 1 );
                    }

                    if( CheckIndex( idx_lOrder + 1, -kMinus( idx_kOrder ) ) ) {
                        idx_col                             = GlobalIndex( idx_lOrder + 1, -kMinus( idx_kOrder ) );
                        _Ay( idx_row, unsigned( idx_col ) ) = 0.5 * Sgn( idx_kOrder ) * DTilde( idx_lOrder + 1, std::abs( idx_kOrder ) - 1 );
                    }
221
222
                }

steffen.schotthoefer's avatar
steffen.schotthoefer committed
223
224
225
                if( CheckIndex( idx_lOrder - 1, -kPlus( idx_kOrder ) ) ) {
                    idx_col                             = GlobalIndex( idx_lOrder - 1, -kPlus( idx_kOrder ) );
                    _Ay( idx_row, unsigned( idx_col ) ) = -0.5 * Sgn( idx_kOrder ) * ETilde( idx_lOrder - 1, std::abs( idx_kOrder ) + 1 );
226
227
                }

steffen.schotthoefer's avatar
steffen.schotthoefer committed
228
229
230
231
                if( CheckIndex( idx_lOrder + 1, -kPlus( idx_kOrder ) ) ) {
                    idx_col                             = GlobalIndex( idx_lOrder + 1, -kPlus( idx_kOrder ) );
                    _Ay( idx_row, unsigned( idx_col ) ) = 0.5 * Sgn( idx_kOrder ) * FTilde( idx_lOrder + 1, std::abs( idx_kOrder ) + 1 );
                }
232
233
234
            }

            // flux matrix in direction z
steffen.schotthoefer's avatar
steffen.schotthoefer committed
235
236
237
238
239
            {
                if( CheckIndex( idx_lOrder - 1, idx_kOrder ) ) {
                    idx_col                             = GlobalIndex( idx_lOrder - 1, idx_kOrder );
                    _Az( idx_row, unsigned( idx_col ) ) = AParam( idx_lOrder - 1, idx_kOrder );
                }
240

steffen.schotthoefer's avatar
steffen.schotthoefer committed
241
242
243
244
                if( CheckIndex( idx_lOrder + 1, idx_kOrder ) ) {
                    idx_col                             = GlobalIndex( idx_lOrder + 1, idx_kOrder );
                    _Az( idx_row, unsigned( idx_col ) ) = BParam( idx_lOrder + 1, idx_kOrder );
                }
245
246
247
248
249
250
            }
        }
    }
}

void PNSolver::ComputeFluxComponents() {
251
252
253
    Vector eigenValues( _nSystem, 0 );
    Vector eigenValuesX( _nSystem, 0 );
    Vector eigenValuesY( _nSystem, 0 );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
254

255
    MatrixCol eigenVectors( _nSystem, _nSystem, 0 );    // ColumnMatrix for _AxPlus * eigenVectors Multiplication via SIMD
256
257
258
259
260
    // --- For x Direction ---
    {
        blaze::eigen( _Ax, eigenValues, eigenVectors );    // Compute Eigenvalues and Eigenvectors

        // Compute Flux Matrices A+ and A-
261
        for( unsigned idx_ij = 0; idx_ij < _nSystem; idx_ij++ ) {
262
263
            if( eigenValues[idx_ij] >= 0 ) {
                _AxPlus( idx_ij, idx_ij ) = eigenValues[idx_ij];    // positive part of Diagonal Matrix stored in _AxPlus
steffen.schotthoefer's avatar
steffen.schotthoefer committed
264
                _AxAbs( idx_ij, idx_ij )  = eigenValues[idx_ij];
265
266
267
            }
            else {
                _AxMinus( idx_ij, idx_ij ) = eigenValues[idx_ij];    // negative part of Diagonal Matrix stored in _AxMinus
steffen.schotthoefer's avatar
steffen.schotthoefer committed
268
                _AxAbs( idx_ij, idx_ij )   = -eigenValues[idx_ij];
269
270
271
272
273
            }
        }

        _AxPlus  = eigenVectors * _AxPlus;    // col*row minimum performance
        _AxMinus = eigenVectors * _AxMinus;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
274
        _AxAbs   = eigenVectors * _AxAbs;
275
276
277
        blaze::transpose( eigenVectors );
        _AxPlus  = _AxPlus * eigenVectors;    // row*col maximum performance
        _AxMinus = _AxMinus * eigenVectors;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
278
279
        _AxAbs   = _AxAbs * eigenVectors;

280
        // eigenValuesX = eigenValues;
281
282
283
284
285
286
    }
    // --- For y Direction -------
    {
        blaze::eigen( _Ay, eigenValues, eigenVectors );    // Compute Eigenvalues and Eigenvectors

        // Compute Flux Matrices A+ and A-
287
        for( unsigned idx_ij = 0; idx_ij < _nSystem; idx_ij++ ) {
288
289
            if( eigenValues[idx_ij] >= 0 ) {
                _AyPlus( idx_ij, idx_ij ) = eigenValues[idx_ij];    // positive part of Diagonal Matrix stored in _AxPlus
steffen.schotthoefer's avatar
steffen.schotthoefer committed
290
                _AyAbs( idx_ij, idx_ij )  = eigenValues[idx_ij];
291
292
293
            }
            else {
                _AyMinus( idx_ij, idx_ij ) = eigenValues[idx_ij];    // negative part of Diagonal Matrix stored in _AxMinus
steffen.schotthoefer's avatar
steffen.schotthoefer committed
294
                _AyAbs( idx_ij, idx_ij )   = -eigenValues[idx_ij];
295
296
297
298
299
            }
        }

        _AyPlus  = eigenVectors * _AyPlus;
        _AyMinus = eigenVectors * _AyMinus;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
300
        _AyAbs   = eigenVectors * _AyAbs;
301
302
303
        blaze::transpose( eigenVectors );
        _AyPlus  = _AyPlus * eigenVectors;
        _AyMinus = _AyMinus * eigenVectors;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
304
305
        _AyAbs   = _AyAbs * eigenVectors;

306
        // eigenValuesY = eigenValues;
307
308
309
310
311
312
    }
    // --- For z Direction -------
    {
        blaze::eigen( _Az, eigenValues, eigenVectors );    // Compute Eigenvalues and Eigenvectors

        // Compute Flux Matrices A+ and A-
313
        for( unsigned idx_ij = 0; idx_ij < _nSystem; idx_ij++ ) {
314
315
            if( eigenValues[idx_ij] >= 0 ) {
                _AzPlus( idx_ij, idx_ij ) = eigenValues[idx_ij];    // positive part of Diagonal Matrix stored in _AxPlus
steffen.schotthoefer's avatar
steffen.schotthoefer committed
316
                _AzAbs( idx_ij, idx_ij )  = eigenValues[idx_ij];
317
318
319
            }
            else {
                _AzMinus( idx_ij, idx_ij ) = eigenValues[idx_ij];    // negative part of Diagonal Matrix stored in _AxMinus
steffen.schotthoefer's avatar
steffen.schotthoefer committed
320
                _AzAbs( idx_ij, idx_ij )   = -eigenValues[idx_ij];
321
322
323
324
325
            }
        }

        _AzPlus  = eigenVectors * _AzPlus;
        _AzMinus = eigenVectors * _AzMinus;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
326
        _AzAbs   = eigenVectors * _AzAbs;
327
328
329
        blaze::transpose( eigenVectors );
        _AzPlus  = _AzPlus * eigenVectors;
        _AzMinus = _AzMinus * eigenVectors;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
330
        _AzAbs   = _AzAbs * eigenVectors;
331
    }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
332
333

    // Compute Spectral Radius
334
335
336
337
338
339
340
    // std::cout << "Eigenvalues x direction " << eigenValuesX << "\n";
    // std::cout << "Eigenvalues y direction " << eigenValuesY << "\n";
    // std::cout << "Eigenvalues z direction " << eigenValues << "\n";
    //
    // std::cout << "Spectral Radius X " << blaze::max( blaze::abs( eigenValuesX ) ) << "\n";
    // std::cout << "Spectral Radius Y " << blaze::max( blaze::abs( eigenValuesY ) ) << "\n";
    // std::cout << "Spectral Radius Z " << blaze::max( blaze::abs( eigenValues ) ) << "\n";
steffen.schotthoefer's avatar
steffen.schotthoefer committed
341

342
    //_combinedSpectralRadius = blaze::max( blaze::abs( eigenValues + eigenValuesX + eigenValuesY ) );
343
    // std::cout << "Spectral Radius combined " << _combinedSpectralRadius << "\n";
344
345
}

346
void PNSolver::ComputeScatterMatrix() {
347

348
    // --- Isotropic ---
Steffen Schotthöfer's avatar
update    
Steffen Schotthöfer committed
349
    _scatterMatDiag[0] = -1;
350
    for( unsigned idx_diag = 1; idx_diag < _nSystem; idx_diag++ ) {
351
        _scatterMatDiag[idx_diag] = 0.0;
352
    }
353
354
}

355
double PNSolver::LegendrePoly( double x, int l ) {    // Legacy. TO BE DELETED
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
    // Pre computed low order polynomials for faster computation
    switch( l ) {
        case 0: return 1;
        case 1: return x;
        case 2:    // 0.5*(3*x*x - 1)
            return 1.5 * x * x - 0.5;
        case 3:    // 0.5* (5*x*x*x -3 *x)
            return 2.5 * x * x * x - 1.5 * x;
        case 4:    // 1/8*(35x^4-30x^2 + 3)
            return 4.375 * x * x * x * x - 3.75 * x * x + 0.375;
        case 5:    // 1/8(63x^5-70x^3 + 15*x )
            return 7.875 * x * x * x * x * x - 8.75 * x * x * x + 1.875 * x;
        case 6:    // 1/16(231x^6-315x^4+105x^2-5)
            return 14.4375 * x * x * x * x * x * x - 19.6875 * x * x * x * x + 6.5625 * x * x - 3.125;
        default: ErrorMessages::Error( "Legendre Polynomials only implemented up to order 6", CURRENT_FUNCTION ); return 0;
    }
}
373

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
374
void PNSolver::PrepareVolumeOutput() {
375
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
steffen.schotthoefer's avatar
steffen.schotthoefer committed
376

377
378
379
380
381
382
383
384
385
386
387
388
389
390
    _outputFieldNames.resize( nGroups );
    _outputFields.resize( nGroups );

    // Prepare all OutputGroups ==> Specified in option VOLUME_OUTPUT
    for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
        // Prepare all Output Fields per group

        // Different procedure, depending on the Group...
        switch( _settings->GetVolumeOutput()[idx_group] ) {
            case MINIMAL:
                // Currently only one entry ==> rad flux
                _outputFields[idx_group].resize( 1 );
                _outputFieldNames[idx_group].resize( 1 );

391
392
                _outputFields[idx_group][0].resize( _nCells );
                _outputFieldNames[idx_group][0] = "radiation flux density";
393
394
395
396
                break;

            case MOMENTS:
                // As many entries as there are moments in the system
397
398
                _outputFields[idx_group].resize( _nSystem );
                _outputFieldNames[idx_group].resize( _nSystem );
399

400
                for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
401
402
403
404
405
406
407
408
409
410
411
412
413
                    for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
                        _outputFields[idx_group][GlobalIndex( idx_l, idx_k )].resize( _nCells );

                        _outputFieldNames[idx_group][GlobalIndex( idx_l, idx_k )] =
                            std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                    }
                }
                break;
            default: ErrorMessages::Error( "Volume Output Group not defined for PN Solver!", CURRENT_FUNCTION ); break;
        }
    }
}

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
414
void PNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
415
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
416

417
    if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
418
        ( idx_pseudoTime == _nEnergies - 1 ) /* need sol at last iteration */ ) {
419
420
421
        for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
            switch( _settings->GetVolumeOutput()[idx_group] ) {
                case MINIMAL:
422
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
423
                        _outputFields[idx_group][0][idx_cell] = _fluxNew[idx_cell];
424
                    }
425
426
                    break;
                case MOMENTS:
427
                    for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
428
                        for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
429
                            _outputFields[idx_group][idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
430
431
432
433
434
                        }
                    }
                    break;
                default: ErrorMessages::Error( "Volume Output Group not defined for PN Solver!", CURRENT_FUNCTION ); break;
            }
435
436
        }
    }
437
438
}

439
void PNSolver::CleanFluxMatrices() {
440
441
    for( unsigned idx_row = 0; idx_row < _nSystem; idx_row++ ) {
        for( unsigned idx_col = 0; idx_col < _nSystem; idx_col++ ) {
442
443
444
            if( std::abs( _AxAbs( idx_row, idx_col ) ) < 0.00000000001 ) _AxAbs( idx_row, idx_col ) = 0.0;
            if( std::abs( _AxPlus( idx_row, idx_col ) ) < 0.00000000001 ) _AxPlus( idx_row, idx_col ) = 0.0;
            if( std::abs( _AxMinus( idx_row, idx_col ) ) < 0.00000000001 ) _AxMinus( idx_row, idx_col ) = 0.0;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
445

446
447
448
            if( std::abs( _AyAbs( idx_row, idx_col ) ) < 0.00000000001 ) _AyAbs( idx_row, idx_col ) = 0.0;
            if( std::abs( _AyPlus( idx_row, idx_col ) ) < 0.00000000001 ) _AyPlus( idx_row, idx_col ) = 0.0;
            if( std::abs( _AyMinus( idx_row, idx_col ) ) < 0.00000000001 ) _AyMinus( idx_row, idx_col ) = 0.0;
449

450
451
452
453
454
455
            if( std::abs( _AzAbs( idx_row, idx_col ) ) < 0.00000000001 ) _AzAbs( idx_row, idx_col ) = 0.0;
            if( std::abs( _AzPlus( idx_row, idx_col ) ) < 0.00000000001 ) _AzPlus( idx_row, idx_col ) = 0.0;
            if( std::abs( _AzMinus( idx_row, idx_col ) ) < 0.00000000001 ) _AzMinus( idx_row, idx_col ) = 0.0;
        }
    }
}
456

457
458
459
460
461
462
463
double PNSolver::CTilde( int l, int k ) const {
    if( k < 0 ) return 0.0;
    if( k == 0 )
        return std::sqrt( 2 ) * CParam( l, k );
    else
        return CParam( l, k );
}
464

465
466
467
468
469
470
471
double PNSolver::DTilde( int l, int k ) const {
    if( k < 0 ) return 0.0;
    if( k == 0 )
        return std::sqrt( 2 ) * DParam( l, k );
    else
        return DParam( l, k );
}
steffen.schotthoefer's avatar
steffen.schotthoefer committed
472

473
474
475
476
477
478
double PNSolver::ETilde( int l, int k ) const {
    if( k == 1 )
        return std::sqrt( 2 ) * EParam( l, k );
    else
        return EParam( l, k );
}
479

480
481
482
483
484
485
double PNSolver::FTilde( int l, int k ) const {
    if( k == 1 )
        return std::sqrt( 2 ) * FParam( l, k );
    else
        return FParam( l, k );
}
486

487
488
489
double PNSolver::AParam( int l, int k ) const {
    return std::sqrt( double( ( l - k + 1 ) * ( l + k + 1 ) ) / double( ( 2 * l + 3 ) * ( 2 * l + 1 ) ) );
}
490

491
double PNSolver::BParam( int l, int k ) const { return std::sqrt( double( ( l - k ) * ( l + k ) ) / double( ( ( 2 * l + 1 ) * ( 2 * l - 1 ) ) ) ); }
492

493
494
495
double PNSolver::CParam( int l, int k ) const {
    return std::sqrt( double( ( l + k + 1 ) * ( l + k + 2 ) ) / double( ( ( 2 * l + 3 ) * ( 2 * l + 1 ) ) ) );
}
496

497
498
499
double PNSolver::DParam( int l, int k ) const {
    return std::sqrt( double( ( l - k ) * ( l - k - 1 ) ) / double( ( ( 2 * l + 1 ) * ( 2 * l - 1 ) ) ) );
}
500

501
502
503
double PNSolver::EParam( int l, int k ) const {
    return std::sqrt( double( ( l - k + 1 ) * ( l - k + 2 ) ) / double( ( ( 2 * l + 3 ) * ( 2 * l + 1 ) ) ) );
}
504

505
double PNSolver::FParam( int l, int k ) const { return std::sqrt( double( ( l + k ) * ( l + k - 1 ) ) / double( ( 2 * l + 1 ) * ( 2 * l - 1 ) ) ); }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
506

507
int PNSolver::kPlus( int k ) const { return k + Sgn( k ); }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
508

509
int PNSolver::kMinus( int k ) const { return k - Sgn( k ); }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
510

511
512
513
514
int PNSolver::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;
515
516
}

517
bool PNSolver::CheckIndex( int l, int k ) const {
518
    if( l >= 0 && l <= int( _polyDegreeBasis ) ) {
519
        if( k >= -l && k <= l ) return true;
520
    }
521
    return false;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
522
523
}

524
525
526
527
528
int PNSolver::Sgn( int k ) const {
    if( k >= 0 )
        return 1;
    else
        return -1;
529
}