mnsolver.cpp 12.6 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
10
#include "quadratures/quadraturebase.h"
#include "solvers/sphericalharmonics.h"
11
#include "toolboxes/errormessages.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
//#include <chrono>

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

23
    // Is this good (fast) code using a constructor list?
24
25
    _LMaxDegree    = settings->GetMaxMomentDegree();
    _nTotalEntries = GlobalIndex( _LMaxDegree, int( _LMaxDegree ) ) + 1;
26
27
28
29
30
31
32

    // build quadrature object and store quadrature points and weights
    _quadPoints       = _quadrature->GetPoints();
    _weights          = _quadrature->GetWeights();
    _nq               = _quadrature->GetNq();
    _quadPointsSphere = _quadrature->GetPointsSphere();
    _settings->SetNQuadPoints( _nq );
33

34
    // This must be shifted to Problem Class
steffen.schotthoefer's avatar
steffen.schotthoefer committed
35
36
    for( unsigned n = 0; n < _nEnergies; n++ ) {
        for( unsigned j = 0; j < _nCells; j++ ) {
37
38
            _sigmaT[n][j] = 1;    //_sigmaT[n][j] - _sigmaS[n][j];
            _sigmaS[n][j] = 1;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
39
40
41
        }
    }

42
43
44
    // Initialize Scatter Matrix --
    _scatterMatDiag = Vector( _nTotalEntries, 0.0 );
    ComputeScatterMatrix();
steffen.schotthoefer's avatar
steffen.schotthoefer committed
45

46
47
48
49
50
    // Initialize Entropy
    _entropy = EntropyBase::Create( _settings );

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

52
53
    // Initialize lagrange Multiplier
    _alpha = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
54
55

    // Initialize and Pre-Compute Moments at quadrature points
56
    _basis = new SphericalHarmonics( _LMaxDegree );
57

58
    _moments = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) );
59
    ComputeMoments();
60
61

    // Solver output
62
    PrepareOutputFields();
steffen.schotthoefer's avatar
steffen.schotthoefer committed
63
64
}

65
66
67
MNSolver::~MNSolver() {
    delete _entropy;
    delete _optimizer;
68
    delete _basis;
69
}
70

71
72
73
74
75
76
77
78
79
void MNSolver::ComputeScatterMatrix() {

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

steffen.schotthoefer's avatar
steffen.schotthoefer committed
80
81
82
83
84
85
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;
}

86
void MNSolver::ComputeMoments() {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
87
    double my, phi;
88
89

    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
90
91
        my  = _quadPointsSphere[idx_quad][0];
        phi = _quadPointsSphere[idx_quad][1];
92

93
        _moments[idx_quad] = _basis->ComputeSphericalBasis( my, phi );
94
    }
95
96
97
98
99
}

Vector MNSolver::ConstructFlux( unsigned idx_cell ) {

    // ---- Integration of Moment of flux ----
100
    double entropyL, entropyR, entropyFlux;
101

102
    Vector flux( _nTotalEntries, 0.0 );
103

104
    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
105
106
107

        entropyFlux = 0.0;    // Reset temorary flux

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

110
111
112
113
114
        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 {
115
                entropyR = _entropy->EntropyPrimeDual( blaze::dot( _alpha[_neighbors[idx_cell][idx_neigh]], _moments[idx_quad] ) );
116
            }
117
            entropyFlux += _g->Flux( _quadPoints[idx_quad], entropyL, entropyR, _normals[idx_cell][idx_neigh] );
118
        }
119
        flux += _moments[idx_quad] * ( _weights[idx_quad] * entropyFlux );
120
121

        // ------- Relizablity Reconstruction Step ----
122
    }
123
    return flux;
124
125
}

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

steffen.schotthoefer's avatar
steffen.schotthoefer committed
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?)
144
    VectorVector psiNew = _sol;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
145

146
    double mass = 0;
147
148
149

    Save( -1 );

150
151
152
    if( rank == 0 ) log->info( "{:10}  {:10}", "t", "mass" );

    // if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", -1.0, dFlux, mass1 ); Should be deleted
steffen.schotthoefer's avatar
steffen.schotthoefer committed
153
154

    // Loop over energies (pseudo-time of continuous slowing down approach)
155
    for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
156

157
        // ------- Reconstruction Step ----------------
158

159
        _optimizer->SolveMultiCell( _alpha, _sol, _moments );
160

161
162
        // Loop over the grid cells
        for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
163

164
165
            // ------- Relizablity Reconstruction Step ----

166
            ComputeRealizableSolution( idx_cell );
167
168

            // ------- Flux Computation Step --------------
169

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

173
            psiNew[idx_cell] = ConstructFlux( idx_cell );
174

175
            // ------ Finite Volume Update Step ------
176

177
178
179
            // NEED TO VECTORIZE
            for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {

180
                psiNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
181
                                               ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
182
                                               - _dE * _sol[idx_cell][idx_system] *
183
                                                     ( _sigmaT[idx_energy][idx_cell]                                    /* absorbtion influence */
184
                                                       + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
185
            }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
186
        }
187

188
189
        // Update Solution
        _sol = psiNew;
190

191
        // --- VTK and CSV Output ---
192
        mass = WriteOutputFields( idx_energy );
193
        Save( idx_energy );
194
        WriteNNTrainingData( idx_energy );
195

196
        // --- Screen Output ---
197
        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e}", _energies[idx_energy], mass );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
198
199
    }
}
200

201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
void MNSolver::PrepareOutputFields() {
    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 );

218
219
                _outputFields[idx_group][0].resize( _nCells );
                _outputFieldNames[idx_group][0] = "radiation flux density";
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
                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;
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259

            case DUAL_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( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
                    }
                }
                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;

260
261
262
263
264
            default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
        }
    }
}

265
double MNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
266
267
    double mass      = 0.0;
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
268

269
270
271
    // Compute total "mass" of the system ==> to check conservation properties
    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
        mass += _sol[idx_cell][0] * _areas[idx_cell];    // Should probably go to postprocessing
272
273
    }

274
275
276
277
278
279
280
281
282
283
284
285
286
287
    for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
        switch( _settings->GetVolumeOutput()[idx_group] ) {
            case MINIMAL:
                for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                    _outputFields[idx_group][0][idx_cell] = _sol[idx_cell][0];
                }
                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];
                    }
                }
                break;
288
289
290
291
292
293
294
295
296
297
298
299
            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 ) {

                    double time  = idx_pseudoTime * _dE;
300
                    double sigma = 1;
301
302
303
304
305
306

                    _outputFields[idx_group][0][idx_cell] = _problem->GetAnalyticalSolution(
                        _mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, sigma );
                }
                break;

307
308
309
            default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
        }
    }
310
311
    return mass;
}
312

313
void MNSolver::Save() const { ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); }
314

315
void MNSolver::Save( int currEnergy ) const {
316
    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), _outputFields, _outputFieldNames, _mesh );
317
}
318
319
320
321
322
323
324

void MNSolver::WriteNNTrainingData( unsigned idx_pseudoTime ) {
    std::string filename = "trainNN.csv";
    std::ofstream myfile;
    myfile.open( filename, std::ofstream::app );

    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
325
        myfile << 0 << ", " << _nTotalEntries << "," << idx_pseudoTime;
326
        for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
327
            myfile << "," << _sol[idx_cell][idx_sys];
328
        }
329
        myfile << " \n" << 1 << ", " << _nTotalEntries << "," << idx_pseudoTime;
330
        for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
331
            myfile << "," << _alpha[idx_cell][idx_sys];
332
        }
333
        myfile << "\n";
334
335
336
    }
    myfile.close();
}