mnsolver.cpp 13.8 KB
Newer Older
1
#include "solvers/mnsolver.h"
2
3
#include "common/config.h"
#include "common/io.h"
4
#include "common/mesh.h"
5
6
7
#include "entropies/entropybase.h"
#include "fluxes/numericalflux.h"
#include "optimizers/optimizerbase.h"
8
#include "problems/problembase.h"
9
#include "quadratures/quadraturebase.h"
10
#include "toolboxes/errormessages.h"
11
#include "toolboxes/sphericalbase.h"
12
#include "toolboxes/textprocessingtoolbox.h"
13

14
// externals
15
#include "spdlog/spdlog.h"
16
#include <mpi.h>
steffen.schotthoefer's avatar
steffen.schotthoefer committed
17

18
#include <fstream>
19

20
MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
21

22
    _LMaxDegree    = settings->GetMaxMomentDegree();
23
24
    _basis         = SphericalBase::Create( _settings );
    _nTotalEntries = _basis->GetBasisSize();
25
26

    // build quadrature object and store quadrature points and weights
27
28
29
    _quadPoints = _quadrature->GetPoints();
    _weights    = _quadrature->GetWeights();
    //_nq               = _quadrature->GetNq();
30
    _quadPointsSphere = _quadrature->GetPointsSphere();
31
32
33
34
35
    //_settings->SetNQuadPoints( _nq );

    // Initialize temporary storages of alpha derivatives
    _solDx = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
    _solDy = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
36

37
38
39
    // Initialize Scatter Matrix --
    _scatterMatDiag = Vector( _nTotalEntries, 0.0 );
    ComputeScatterMatrix();
steffen.schotthoefer's avatar
steffen.schotthoefer committed
40

41
42
43
44
45
    // Initialize Entropy
    _entropy = EntropyBase::Create( _settings );

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

47
48
    // Initialize lagrange Multiplier
    _alpha = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
49
50

    // Initialize and Pre-Compute Moments at quadrature points
51
    _moments = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) );
52

53
    ComputeMoments();
steffen.schotthoefer's avatar
steffen.schotthoefer committed
54
55
}

56
57
58
MNSolver::~MNSolver() {
    delete _entropy;
    delete _optimizer;
59
    delete _basis;
60
}
61

62
63
64
65
66
67
68
69
70
void MNSolver::ComputeScatterMatrix() {

    // --- Isotropic ---
    _scatterMatDiag[0] = -1.0;
    for( unsigned idx_diag = 1; idx_diag < _nTotalEntries; idx_diag++ ) {
        _scatterMatDiag[idx_diag] = 0.0;
    }
}

71
void MNSolver::ComputeMoments() {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
72
    double my, phi;
73
74

    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
75
76
        my  = _quadPointsSphere[idx_quad][0];
        phi = _quadPointsSphere[idx_quad][1];
77

78
        _moments[idx_quad] = _basis->ComputeSphericalBasis( my, phi );
79
    }
80
81
82
83
}

Vector MNSolver::ConstructFlux( unsigned idx_cell ) {

84
    //--- Integration of moments of flux ---
85
    double entropyL, entropyR, entropyFlux;
86

87
    Vector flux( _nTotalEntries, 0.0 );
88

89
90
91
    //--- Temporary storages of reconstructed alpha ---
    Vector alphaL( _nTotalEntries, 0.0 );
    Vector alphaR( _nTotalEntries, 0.0 );
92

93
94
    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
        entropyFlux = 0.0;    // reset temorary flux
95

96
        for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
97
98
99
100
101
102
103
104
105
106
107
            // Left side reconstruction
            if( _reconsOrder > 1 ) {
                alphaL = _alpha[idx_cell] + _solDx[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neigh][0] - _cellMidPoints[idx_cell][0] ) +
                         _solDy[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neigh][1] - _cellMidPoints[idx_cell][1] );
            }
            else {
                alphaL = _alpha[idx_cell];
            }
            entropyL = _entropy->EntropyPrimeDual( blaze::dot( alphaL, _moments[idx_quad] ) );

            // Right side reconstruction
108
109
110
            if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neigh] == _nCells )
                entropyR = entropyL;
            else {
111
112
113
114
115
116
117
118
119
120
121
                if( _reconsOrder > 1 ) {
                    alphaR = _alpha[_neighbors[idx_cell][idx_neigh]] +
                             _solDx[_neighbors[idx_cell][idx_neigh]] *
                                 ( _interfaceMidPoints[idx_cell][idx_neigh][0] - _cellMidPoints[_neighbors[idx_cell][idx_neigh]][0] ) +
                             _solDy[_neighbors[idx_cell][idx_neigh]] *
                                 ( _interfaceMidPoints[idx_cell][idx_neigh][1] - _cellMidPoints[_neighbors[idx_cell][idx_neigh]][1] );
                }
                else {
                    alphaR = _alpha[_neighbors[idx_cell][idx_neigh]];
                }
                entropyR = _entropy->EntropyPrimeDual( blaze::dot( alphaR, _moments[idx_quad] ) );
122
            }
123
124

            // Entropy flux
125
            entropyFlux += _g->Flux( _quadPoints[idx_quad], entropyL, entropyR, _normals[idx_cell][idx_neigh] );
126
        }
127
        // Solution flux
128
        flux += _moments[idx_quad] * ( _weights[idx_quad] * entropyFlux );
129
  }
130
    return flux;
131
132
}

133
134
void MNSolver::ComputeRealizableSolution( unsigned idx_cell ) {
    double entropyReconstruction = 0.0;
135
    _sol[idx_cell]               = 0;
136
    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
137
        // Make entropyReconstruction a member vector, s.t. it does not have to be re-evaluated in ConstructFlux
138
        entropyReconstruction = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _moments[idx_quad] ) );
139
        _sol[idx_cell] += _moments[idx_quad] * ( _weights[idx_quad] * entropyReconstruction );
140
141
142
    }
}

143
void MNSolver::IterPreprocessing( unsigned idx_pseudotime ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
144

145
    // ------- Reconstruction Step ----------------
146

147
    _optimizer->SolveMultiCell( _alpha, _sol, _moments );
148

149
    // ------- Relizablity Preservation Step ----
150
151
152
153
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
        ComputeRealizableSolution( idx_cell );
    }
}
154

155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
void MNSolver::IterPostprocessing() {
    // --- Update Solution ---
    _sol = _solNew;

    // --- Compute Flux for solution and Screen Output ---
    ComputeRadFlux();
}

void MNSolver::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;
    }
}

170
void MNSolver::FluxUpdate() {
171
172
173
174
    if( _reconsOrder > 1 ) {
        _mesh->ReconstructSlopesU( _nTotalEntries, _solDx, _solDy, _alpha );    // unstructured reconstruction
    }

175
176
177
178
    // Loop over the grid cells
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
        // Dirichlet Boundaries stay
        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
179
        _solNew[idx_cell] = ConstructFlux( idx_cell );
180
181
    }
}
182

183
void MNSolver::FVMUpdate( unsigned idx_energy ) {
184
185
186
187
    // Loop over the grid cells
    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
        // Dirichlet Boundaries stay
        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
188

189
        for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
190

191
192
193
194
195
            _solNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
                                            ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_system] /* cell averaged flux */
                                            - _dE * _sol[idx_cell][idx_system] *
                                                  ( _sigmaT[idx_energy][idx_cell]                                    /* absorbtion influence */
                                                    + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
196
        }
197

198
        _solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
steffen.schotthoefer's avatar
steffen.schotthoefer committed
199
200
    }
}
201

202
void MNSolver::PrepareVolumeOutput() {
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();

    _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 );

219
220
                _outputFields[idx_group][0].resize( _nCells );
                _outputFieldNames[idx_group][0] = "radiation flux density";
221
222
223
224
225
226
227
                break;

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

228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
                if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
                    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][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
                            _outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
                                std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                        }
                    }
                }
                else {
                    for( unsigned idx_l = 0; idx_l <= _LMaxDegree; idx_l++ ) {
                        unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
                        for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {
                            _outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
                            _outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
                                std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                        }
245
246
247
                    }
                }
                break;
248
249
250
251
252
253

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

254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
                if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
                    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][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
                            _outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
                                std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                        }
                    }
                }
                else {
                    for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
                        unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
                        for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {
                            _outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
                            _outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
                                std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                        }
271
272
273
274
275
276
277
278
279
280
281
282
                    }
                }
                break;

            case ANALYTIC:
                // one entry per cell
                _outputFields[idx_group].resize( 1 );
                _outputFieldNames[idx_group].resize( 1 );
                _outputFields[idx_group][0].resize( _nCells );
                _outputFieldNames[idx_group][0] = std::string( "analytic radiation flux density" );
                break;

283
284
285
286
287
            default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
        }
    }
}

288
void MNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
289
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
290

291
    // Check if volume output fields are written to file this iteration
292
    if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
293
        ( idx_pseudoTime == _nEnergies - 1 ) /* need sol at last iteration */ ) {
294

295
296
297
        for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
            switch( _settings->GetVolumeOutput()[idx_group] ) {
                case MINIMAL:
298
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
299
                        _outputFields[idx_group][0][idx_cell] = _fluxNew[idx_cell];
300
                    }
301
302
303
304
305
306
                    break;
                case MOMENTS:
                    for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
                        for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                            _outputFields[idx_group][idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
                        }
307
                    }
308
309
310
311
312
313
314
315
316
317
318
                    break;
                case DUAL_MOMENTS:
                    for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
                        for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                            _outputFields[idx_group][idx_sys][idx_cell] = _alpha[idx_cell][idx_sys];
                        }
                    }
                    break;
                case ANALYTIC:
                    // Compute total "mass" of the system ==> to check conservation properties
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
319

320
                        double time = idx_pseudoTime * _dE;
321

322
323
324
325
                        _outputFields[idx_group][0][idx_cell] = _problem->GetAnalyticalSolution(
                            _mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, _sigmaS[idx_pseudoTime][idx_cell] );
                    }
                    break;
326

327
328
                default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
            }
329
330
        }
    }
331
}