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

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

16
#include <fstream>
17
18
//#include <chrono>

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

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

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

steffen.schotthoefer's avatar
steffen.schotthoefer committed
32
33
34
35
36
37
38
39
40
41
42
43
44
45
    // 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.

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

steffen.schotthoefer's avatar
steffen.schotthoefer committed
71
72
73
74
75
76
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;
}

77
void MNSolver::ComputeMoments() {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
78
    double my, phi;
79
80

    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
81
82
        my  = _quadPointsSphere[idx_quad][0];
        phi = _quadPointsSphere[idx_quad][1];
83

84
        _moments[idx_quad] = _basis->ComputeSphericalBasis( my, phi );
85
    }
86
87
88
89
90
}

Vector MNSolver::ConstructFlux( unsigned idx_cell ) {

    // ---- Integration of Moment of flux ----
91
    double entropyL, entropyR, entropyFlux;
92

93
    Vector flux( _nTotalEntries, 0.0 );
94

95
    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
96
97
98

        entropyFlux = 0.0;    // Reset temorary flux

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

101
102
103
104
105
        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 {
106
                entropyR = _entropy->EntropyPrimeDual( blaze::dot( _alpha[_neighbors[idx_cell][idx_neigh]], _moments[idx_quad] ) );
107
            }
108
            entropyFlux += _g->Flux( _quadPoints[idx_quad], entropyL, entropyR, _normals[idx_cell][idx_neigh] );
109
        }
110
        flux += _moments[idx_quad] * ( _weights[idx_quad] * entropyFlux );
111
112

        // ------- Relizablity Reconstruction Step ----
113
    }
114
    return flux;
115
116
}

117
118
void MNSolver::ComputeRealizableSolution( unsigned idx_cell ) {
    double entropyReconstruction = 0.0;
119
    _sol[idx_cell]               = 0;
120
    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
121
        // Make entropyReconstruction a member vector, s.t. it does not have to be re-evaluated in ConstructFlux
122
        entropyReconstruction = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _moments[idx_quad] ) );
123
        _sol[idx_cell] += _moments[idx_quad] * ( _weights[idx_quad] * entropyReconstruction );
124
125
126
    }
}

steffen.schotthoefer's avatar
steffen.schotthoefer committed
127
128
129
130
131
132
133
134
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?)
135
    VectorVector psiNew = _sol;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
136

137
    double mass = 0;
138
139
140

    Save( -1 );

141
142
143
    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
144
145

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

148
        // ------- Reconstruction Step ----------------
149

150
        _optimizer->SolveMultiCell( _alpha, _sol, _moments );
151

152
153
        // Loop over the grid cells
        for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
154

155
156
            // ------- Relizablity Reconstruction Step ----

157
            ComputeRealizableSolution( idx_cell );
158
159

            // ------- Flux Computation Step --------------
160

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

164
            psiNew[idx_cell] = ConstructFlux( idx_cell );
165

166
            // ------ Finite Volume Update Step ------
167

168
169
170
            // NEED TO VECTORIZE
            for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {

171
                psiNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
172
                                               ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
173
                                               - _dE * _sol[idx_cell][idx_system] *
174
175
                                                     ( _sigmaA[idx_energy][idx_cell]                                    /* absorbtion influence */
                                                       + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
176
            }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
177
        }
178

179
180
        // Update Solution
        _sol = psiNew;
181

182
183
        // --- VTK and CSV Output ---
        mass = WriteOutputFields();
184
        Save( idx_energy );
185
        WriteNNTrainingData( idx_energy );
186

187
        // --- Screen Output ---
188
        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e}", _energies[idx_energy], mass );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
189
190
    }
}
191

192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
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 );

                _outputFields[0][0].resize( _nCells );
                _outputFieldNames[0][0] = "radiation flux density";
                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 MN Solver!", CURRENT_FUNCTION ); break;
        }
    }
}

232
double MNSolver::WriteOutputFields() {
233
234
    double mass      = 0.0;
    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
235

236
237
238
    // 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
239
240
    }

241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
    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;
            default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
        }
    }
258
259
    return mass;
}
260

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

263
void MNSolver::Save( int currEnergy ) const {
264
    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), _outputFields, _outputFieldNames, _mesh );
265
}
266
267
268
269
270
271
272

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++ ) {
273
        myfile << 0 << ", " << _nTotalEntries << "," << idx_pseudoTime;
274
        for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
275
            myfile << "," << _sol[idx_cell][idx_sys];
276
        }
277
        myfile << " \n" << 1 << ", " << _nTotalEntries << "," << idx_pseudoTime;
278
        for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
279
            myfile << "," << _alpha[idx_cell][idx_sys];
280
        }
281
        myfile << "\n";
282
283
284
    }
    myfile.close();
}