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

11
// externals
12
#include <mpi.h>
steffen.schotthoefer's avatar
steffen.schotthoefer committed
13

14
15
#include <fstream>

16
17
//#include <chrono>

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

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

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

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

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

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

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

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

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

    // Solver output
    _outputFields = std::vector( _nTotalEntries, std::vector( _nCells, 0.0 ) );
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

    // Delete that
    /*
        std::string filename = "obFunc.csv";
        std::ofstream myfile;
        myfile.open( filename );

        NewtonOptimizer* optimizer = new NewtonOptimizer( settings );

        Vector alpha( 4, 0.0 );
        Vector u( 4, 0.0 );
        u[0] = -1.11;
        u[1] = 0;
        u[2] = 0;
        u[4] = 0;

        for( double a1 = -1000; a1 <= 100; a1 += 10.0 / 100.0 ) {
            for( double a2 = -60; a2 <= 100; a2 += 60.0 / 40.0 ) {
                alpha[0] = a1;
                alpha[1] = a2;
                myfile << a1 << ", " << a2 << ", " << optimizer->ComputeObjFunc( alpha, u, _moments ) << "\n";
            }
        }
        myfile.close();
     */
steffen.schotthoefer's avatar
steffen.schotthoefer committed
87
88
}

89
90
91
MNSolver::~MNSolver() {
    delete _entropy;
    delete _optimizer;
92
    delete _basis;
93
}
94

steffen.schotthoefer's avatar
steffen.schotthoefer committed
95
96
97
98
99
100
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;
}

101
void MNSolver::ComputeMoments() {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
102
    double my, phi;
103
104

    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
105
106
        my  = _quadPointsSphere[idx_quad][0];
        phi = _quadPointsSphere[idx_quad][1];
107

108
        _moments[idx_quad] = _basis->ComputeSphericalBasis( my, phi );
109
    }
110
111
112
113
114
}

Vector MNSolver::ConstructFlux( unsigned idx_cell ) {

    // ---- Integration of Moment of flux ----
115
    double entropyL, entropyR, entropyFlux;
116

117
    Vector flux( _nTotalEntries, 0.0 );
118

119
    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
120
121
122

        entropyFlux = 0.0;    // Reset temorary flux

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

125
126
127
128
129
        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 {
130
                entropyR = _entropy->EntropyPrimeDual( blaze::dot( _alpha[_neighbors[idx_cell][idx_neigh]], _moments[idx_quad] ) );
131
            }
132
            entropyFlux += _g->Flux( _quadPoints[idx_quad], entropyL, entropyR, _normals[idx_cell][idx_neigh] );
133
        }
134
        flux += _moments[idx_quad] * ( _weights[idx_quad] * entropyFlux );
135
136

        // ------- Relizablity Reconstruction Step ----
137
    }
138
    return flux;
139
140
}

141
142
void MNSolver::ComputeRealizableSolution( unsigned idx_cell ) {
    double entropyReconstruction = 0.0;
143
    _sol[idx_cell]               = 0;
144
    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
145
        // Make entropyReconstruction a member vector, s.t. it does not have to be re-evaluated in ConstructFlux
146
        entropyReconstruction = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _moments[idx_quad] ) );
147
        _sol[idx_cell] += _moments[idx_quad] * ( _weights[idx_quad] * entropyReconstruction );
148
149
150
    }
}

steffen.schotthoefer's avatar
steffen.schotthoefer committed
151
152
153
154
155
156
157
158
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?)
159
    VectorVector psiNew = _sol;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
160
161
162
    double dFlux        = 1e10;
    Vector fluxNew( _nCells, 0.0 );
    Vector fluxOld( _nCells, 0.0 );
163
    VectorVector solTimesArea = _sol;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
164
165

    double mass1 = 0;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
166
    for( unsigned i = 0; i < _nCells; ++i ) {
167
168
        _solverOutput[i] = _sol[i][0];
        mass1 += _sol[i][0] * _areas[i];
steffen.schotthoefer's avatar
steffen.schotthoefer committed
169
170
    }

171
172
173
174
175
    dFlux   = blaze::l2Norm( fluxNew - fluxOld );
    fluxOld = fluxNew;

    Save( -1 );

steffen.schotthoefer's avatar
steffen.schotthoefer committed
176
177
178
179
    if( rank == 0 ) log->info( "{:10}   {:10}", "t", "dFlux" );
    if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", -1.0, dFlux, mass1 );

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

182
        // Loop over the grid cells
183
        //        for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
184

185
186
        // solTimesArea[idx_cell] = _sol[idx_cell] * _areas[idx_cell];    // reconstrucor need moments, not control volume averaged moments!
        // }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
187

188
        // ------- Reconstruction Step ----------------
189

190
        _optimizer->SolveMultiCell( _alpha, _sol, _moments );
191

192
193
        // Loop over the grid cells
        for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
194

195
196
            // ------- Relizablity Reconstruction Step ----

197
            ComputeRealizableSolution( idx_cell );
198
199

            // ------- Flux Computation Step --------------
200

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

204
            psiNew[idx_cell] = ConstructFlux( idx_cell );
205

206
            // ------ Finite Volume Update Step ------
207

208
209
210
            // NEED TO VECTORIZE
            for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {

211
                psiNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
212
                                               ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
213
                                               - _dE * _sol[idx_cell][idx_system] *
214
215
                                                     ( _sigmaA[idx_energy][idx_cell]                                    /* absorbtion influence */
                                                       + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
216
            }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
217
        }
218

219
        // pseudo time iteration output
220
        double mass = 0.0;
221
222
223
224
225
226
227
        for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
            for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                fluxNew[idx_cell]       = _sol[idx_cell][0];    // zeroth moment is raditation densitiy we are interested in
                _solverOutput[idx_cell] = _sol[idx_cell][0];
                mass += _sol[idx_cell][0] * _areas[idx_cell];
                _outputFields[idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
            }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
228
        }
229

steffen.schotthoefer's avatar
steffen.schotthoefer committed
230
231
        dFlux   = blaze::l2Norm( fluxNew - fluxOld );
        fluxOld = fluxNew;
232
        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
233
        Save( idx_energy );
234

235
        // WriteNNTrainingData( idx_energy );
236
237
238

        // Update Solution
        _sol = psiNew;
steffen.schotthoefer's avatar
steffen.schotthoefer committed
239
240
    }
}
241
242
243
244
245
246
247

void MNSolver::Save() const {
    std::vector<std::string> fieldNames{ "flux" };
    std::vector<double> flux;
    flux.resize( _nCells );

    for( unsigned i = 0; i < _nCells; ++i ) {
248
        flux[i] = _sol[i][0];
249
250
251
252
253
254
255
256
257
258
259
260
    }
    std::vector<std::vector<double>> scalarField( 1, flux );
    std::vector<std::vector<std::vector<double>>> results{ scalarField };
    ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
}

void MNSolver::Save( int currEnergy ) const {
    std::vector<std::string> fieldNames{ "flux" };
    std::vector<std::vector<double>> scalarField( 1, _solverOutput );
    std::vector<std::vector<std::vector<double>>> results{ scalarField };
    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
}
261
262
263
264
265
266
267

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