pnsolver.cpp 23.8 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"
Tianbai Xiao's avatar
Tianbai Xiao committed
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 ) : Solver( settings ) {
14
15
    _LMaxDegree    = settings->GetMaxMomentDegree();
    _nTotalEntries = GlobalIndex( int( _LMaxDegree ), int( _LMaxDegree ) ) + 1;
16
17
18
19
20
21
22
23

    // Initialize System Matrices
    _Ax = SymMatrix( _nTotalEntries );
    _Ay = SymMatrix( _nTotalEntries );
    _Az = SymMatrix( _nTotalEntries );

    _AxPlus  = Matrix( _nTotalEntries, _nTotalEntries, 0 );
    _AxMinus = Matrix( _nTotalEntries, _nTotalEntries, 0 );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
24
    _AxAbs   = Matrix( _nTotalEntries, _nTotalEntries, 0 );
25
26
    _AyPlus  = Matrix( _nTotalEntries, _nTotalEntries, 0 );
    _AyMinus = Matrix( _nTotalEntries, _nTotalEntries, 0 );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
27
    _AyAbs   = Matrix( _nTotalEntries, _nTotalEntries, 0 );
28
29
    _AzPlus  = Matrix( _nTotalEntries, _nTotalEntries, 0 );
    _AzMinus = Matrix( _nTotalEntries, _nTotalEntries, 0 );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
30
    _AzAbs   = Matrix( _nTotalEntries, _nTotalEntries, 0 );
31

32
    // Initialize Scatter Matrix
steffen.schotthoefer's avatar
steffen.schotthoefer committed
33
    _scatterMatDiag = Vector( _nTotalEntries, 0 );
34

Tianbai Xiao's avatar
Tianbai Xiao committed
35
    // Initialize temporary storages of solution derivatives
Tianbai Xiao's avatar
Tianbai Xiao committed
36
37
38
    _solDx = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
    _solDy = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );

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

    // Solver output
56
    PrepareVolumeOutput();
57
}
58

59
60
61
void PNSolver::IterPreprocessing() {
    // Nothing to preprocess for PNSolver
}
62

63
64
65
void PNSolver::IterPostprocessing() {
    // --- Update Solution ---
    _sol = _solNew;
66

67
68
    // --- Compute Flux for solution and Screen Output ---
    ComputeRadFlux();
steffen.schotthoefer's avatar
steffen.schotthoefer committed
69
70
}

71
72
73
74
75
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;
    }
76
77
}

78
void PNSolver::FluxUpdate() {
Tianbai Xiao's avatar
Tianbai Xiao committed
79
80
    if( _reconsOrder > 1 ) {
        _mesh->ReconstructSlopesU( _nTotalEntries, _solDx, _solDy, _sol );    // unstructured reconstruction
Tianbai Xiao's avatar
Tianbai Xiao committed
81
        //_mesh->ComputeSlopes( _nTotalEntries, _solDx, _solDy, _sol );    // unstructured reconstruction
Tianbai Xiao's avatar
Tianbai Xiao committed
82
    }
Tianbai Xiao's avatar
Tianbai Xiao committed
83
84
85
86
    //Vector solL( _nTotalEntries );
    //Vector solR( _nTotalEntries );
    auto solL = _sol[2];
    auto solR = _sol[2];
Tianbai Xiao's avatar
Tianbai Xiao committed
87

88
89
    // Loop over all spatial cells
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
90
91
92

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

        // Reset temporary variable psiNew
        for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
96
            _solNew[idx_cell][idx_sys] = 0.0;
97
98
99
100
101
102
103
        }

        // 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 )
104
                _solNew[idx_cell] += _g->Flux(
105
                    _AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, _sol[idx_cell], _sol[idx_cell], _normals[idx_cell][idx_neighbor] );
Tianbai Xiao's avatar
Tianbai Xiao committed
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
            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
Tianbai Xiao's avatar
Tianbai Xiao committed
127
                        solR = _sol[_neighbors[idx_cell][idx_neighbor]] +
Tianbai Xiao's avatar
Tianbai Xiao committed
128
129
                               _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] );
130
131
                        // positivity checker (if not satisfied, deduce to first order)
                        // I manually turned it off here since Pn produces negative solutions essentially
Tianbai Xiao's avatar
Tianbai Xiao committed
132
133
134
135
                        //if( min(solL) < 0.0 || min(solR) < 0.0 ) {
                        //    solL = _sol[idx_cell];
                        //    solR = _sol[_neighbors[idx_cell][idx_neighbor]];
                        //}
Tianbai Xiao's avatar
Tianbai Xiao committed
136
137
138
139
140
141
142
143
144
145
                        // flux evaluation
                        _solNew[idx_cell] += _g->Flux( _AxPlus,
                                                       _AxMinus,
                                                       _AyPlus,
                                                       _AyMinus,
                                                       _AzPlus,
                                                       _AzMinus,
                                                       solL,
                                                       solR,
                                                       _normals[idx_cell][idx_neighbor] );
Tianbai Xiao's avatar
Tianbai Xiao committed
146
                        break;
Tianbai Xiao's avatar
Tianbai Xiao committed
147
148
149
150
151
152
153
154
155
156
157
158
159
                    // 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] );
                }
            }
160
161
162
163
        }
    }
}

164
void PNSolver::FVMUpdate( unsigned idx_energy ) {
165
    // Loop over all spatial cells
166
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
167
168
169
        // Dirichlet cells stay at IC, farfield assumption
        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
        // Flux update
170
        for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
171
172
173
174
            _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] *
                                               ( _sigmaT[idx_energy][idx_cell]                                 /* absorbtion influence */
                                                 + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] ); /* scattering influence */
175
        }
176
177
        // Source Term
        _solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
178
179
180
    }
}

181
void PNSolver::ComputeSystemMatrices() {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
182
183
    int idx_col      = 0;
    unsigned idx_row = 0;
184
185

    // loop over columns of A
186
    for( int idx_lOrder = 0; idx_lOrder <= int( _LMaxDegree ); idx_lOrder++ ) {          // index of legendre polynom
187
        for( int idx_kOrder = -idx_lOrder; idx_kOrder <= idx_lOrder; idx_kOrder++ ) {    // second index of legendre function
steffen.schotthoefer's avatar
steffen.schotthoefer committed
188
            idx_row = unsigned( GlobalIndex( idx_lOrder, idx_kOrder ) );
189
190

            // flux matrix in direction x
steffen.schotthoefer's avatar
steffen.schotthoefer committed
191
192
193
194
195
196
197
198
199
200
201
202
            {
                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 );
                    }
203
204
                }

steffen.schotthoefer's avatar
steffen.schotthoefer committed
205
206
207
                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 );
208
209
                }

steffen.schotthoefer's avatar
steffen.schotthoefer committed
210
211
212
213
                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 );
                }
214
215
216
            }

            // flux matrix in direction y
steffen.schotthoefer's avatar
steffen.schotthoefer committed
217
218
219
220
221
222
223
224
225
226
227
            {
                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 );
                    }
228
229
                }

steffen.schotthoefer's avatar
steffen.schotthoefer committed
230
231
232
                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 );
233
234
                }

steffen.schotthoefer's avatar
steffen.schotthoefer committed
235
236
237
238
                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 );
                }
239
240
241
            }

            // flux matrix in direction z
steffen.schotthoefer's avatar
steffen.schotthoefer committed
242
243
244
245
246
            {
                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 );
                }
247

steffen.schotthoefer's avatar
steffen.schotthoefer committed
248
249
250
251
                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 );
                }
252
253
254
255
256
257
            }
        }
    }
}

void PNSolver::ComputeFluxComponents() {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
258
259
260
261
    Vector eigenValues( _nTotalEntries, 0 );
    Vector eigenValuesX( _nTotalEntries, 0 );
    Vector eigenValuesY( _nTotalEntries, 0 );

262
263
264
265
266
267
268
269
270
    MatrixCol eigenVectors( _nTotalEntries, _nTotalEntries, 0 );    // ColumnMatrix for _AxPlus * eigenVectors Multiplication via SIMD
    // --- For x Direction ---
    {
        blaze::eigen( _Ax, eigenValues, eigenVectors );    // Compute Eigenvalues and Eigenvectors

        // Compute Flux Matrices A+ and A-
        for( unsigned idx_ij = 0; idx_ij < _nTotalEntries; idx_ij++ ) {
            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
271
                _AxAbs( idx_ij, idx_ij )  = eigenValues[idx_ij];
272
273
274
            }
            else {
                _AxMinus( idx_ij, idx_ij ) = eigenValues[idx_ij];    // negative part of Diagonal Matrix stored in _AxMinus
steffen.schotthoefer's avatar
steffen.schotthoefer committed
275
                _AxAbs( idx_ij, idx_ij )   = -eigenValues[idx_ij];
276
277
278
279
280
            }
        }

        _AxPlus  = eigenVectors * _AxPlus;    // col*row minimum performance
        _AxMinus = eigenVectors * _AxMinus;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
281
        _AxAbs   = eigenVectors * _AxAbs;
282
283
284
        blaze::transpose( eigenVectors );
        _AxPlus  = _AxPlus * eigenVectors;    // row*col maximum performance
        _AxMinus = _AxMinus * eigenVectors;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
285
286
        _AxAbs   = _AxAbs * eigenVectors;

287
        // eigenValuesX = eigenValues;
288
289
290
291
292
293
294
295
296
    }
    // --- For y Direction -------
    {
        blaze::eigen( _Ay, eigenValues, eigenVectors );    // Compute Eigenvalues and Eigenvectors

        // Compute Flux Matrices A+ and A-
        for( unsigned idx_ij = 0; idx_ij < _nTotalEntries; idx_ij++ ) {
            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
297
                _AyAbs( idx_ij, idx_ij )  = eigenValues[idx_ij];
298
299
300
            }
            else {
                _AyMinus( idx_ij, idx_ij ) = eigenValues[idx_ij];    // negative part of Diagonal Matrix stored in _AxMinus
steffen.schotthoefer's avatar
steffen.schotthoefer committed
301
                _AyAbs( idx_ij, idx_ij )   = -eigenValues[idx_ij];
302
303
304
305
306
            }
        }

        _AyPlus  = eigenVectors * _AyPlus;
        _AyMinus = eigenVectors * _AyMinus;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
307
        _AyAbs   = eigenVectors * _AyAbs;
308
309
310
        blaze::transpose( eigenVectors );
        _AyPlus  = _AyPlus * eigenVectors;
        _AyMinus = _AyMinus * eigenVectors;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
311
312
        _AyAbs   = _AyAbs * eigenVectors;

313
        // eigenValuesY = eigenValues;
314
315
316
317
318
319
320
321
322
    }
    // --- For z Direction -------
    {
        blaze::eigen( _Az, eigenValues, eigenVectors );    // Compute Eigenvalues and Eigenvectors

        // Compute Flux Matrices A+ and A-
        for( unsigned idx_ij = 0; idx_ij < _nTotalEntries; idx_ij++ ) {
            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
323
                _AzAbs( idx_ij, idx_ij )  = eigenValues[idx_ij];
324
325
326
            }
            else {
                _AzMinus( idx_ij, idx_ij ) = eigenValues[idx_ij];    // negative part of Diagonal Matrix stored in _AxMinus
steffen.schotthoefer's avatar
steffen.schotthoefer committed
327
                _AzAbs( idx_ij, idx_ij )   = -eigenValues[idx_ij];
328
329
330
331
332
            }
        }

        _AzPlus  = eigenVectors * _AzPlus;
        _AzMinus = eigenVectors * _AzMinus;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
333
        _AzAbs   = eigenVectors * _AzAbs;
334
335
336
        blaze::transpose( eigenVectors );
        _AzPlus  = _AzPlus * eigenVectors;
        _AzMinus = _AzMinus * eigenVectors;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
337
        _AzAbs   = _AzAbs * eigenVectors;
338
    }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
339
340

    // Compute Spectral Radius
341
342
343
344
345
346
347
    // 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
348

349
    // _combinedSpectralRadius = blaze::max( blaze::abs( eigenValues + eigenValuesX + eigenValuesY ) );
350
    // std::cout << "Spectral Radius combined " << _combinedSpectralRadius << "\n";
351
352
}

353
void PNSolver::ComputeScatterMatrix() {
354

355
    // --- Isotropic ---
356
    _scatterMatDiag[0] = -1.0;
357
    for( unsigned idx_diag = 1; idx_diag < _nTotalEntries; idx_diag++ ) {
358
        _scatterMatDiag[idx_diag] = 0.0;
359
    }
360
361
}

362
double PNSolver::LegendrePoly( double x, int l ) {    // Legacy. TO BE DELETED
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
    // 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;
    }
}
380

381
void PNSolver::PrepareVolumeOutput() {
382
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
steffen.schotthoefer's avatar
steffen.schotthoefer committed
383

384
385
386
387
388
389
390
391
392
393
394
395
396
397
    _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 );

398
399
                _outputFields[idx_group][0].resize( _nCells );
                _outputFieldNames[idx_group][0] = "radiation flux density";
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
                break;

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

                for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
                    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;
        }
    }
}

421
void PNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
422
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
423

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

446
447
448
449
450
451
void PNSolver::CleanFluxMatrices() {
    for( unsigned idx_row = 0; idx_row < _nTotalEntries; idx_row++ ) {
        for( unsigned idx_col = 0; idx_col < _nTotalEntries; idx_col++ ) {
            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
452

453
454
455
            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;
456

457
458
459
460
461
462
            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;
        }
    }
}
463

464
465
466
467
468
469
470
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 );
}
471

472
473
474
475
476
477
478
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
479

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

487
488
489
490
491
492
double PNSolver::FTilde( int l, int k ) const {
    if( k == 1 )
        return std::sqrt( 2 ) * FParam( l, k );
    else
        return FParam( l, k );
}
493

494
495
496
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 ) ) );
}
497

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

500
501
502
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 ) ) ) );
}
503

504
505
506
double PNSolver::DParam( int l, int k ) const {
    return std::sqrt( double( ( l - k ) * ( l - k - 1 ) ) / double( ( ( 2 * l + 1 ) * ( 2 * l - 1 ) ) ) );
}
507

508
509
510
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 ) ) ) );
}
511

512
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
513

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

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

518
519
520
521
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;
522
523
}

524
bool PNSolver::CheckIndex( int l, int k ) const {
525
    if( l >= 0 && l <= int( _LMaxDegree ) ) {
526
        if( k >= -l && k <= l ) return true;
527
    }
528
    return false;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
529
530
}

531
532
533
534
535
int PNSolver::Sgn( int k ) const {
    if( k >= 0 )
        return 1;
    else
        return -1;
536
}