solverbase.cpp 18.7 KB
Newer Older
1
#include "solvers/solverbase.h"
2
3
4
5
#include "common/config.h"
#include "common/globalconstants.h"
#include "common/io.h"
#include "common/mesh.h"
6
7
#include "fluxes/numericalflux.h"
#include "problems/problembase.h"
steffen.schotthoefer's avatar
steffen.schotthoefer committed
8
#include "quadratures/quadraturebase.h"
jannick.wolters's avatar
jannick.wolters committed
9
#include "solvers/csdpnsolver.h"
Jonas Kusch's avatar
Jonas Kusch committed
10
#include "solvers/csdsnsolver.h"
11
12
13
#include "solvers/csdsolvertrafofp.h"
#include "solvers/csdsolvertrafofp2d.h"
#include "solvers/csdsolvertrafofpsh2d.h"
14
#include "toolboxes/textprocessingtoolbox.h"
15

16
#include "solvers/mnsolver.h"
steffen.schotthoefer's avatar
steffen.schotthoefer committed
17
18
#include "solvers/pnsolver.h"
#include "solvers/snsolver.h"
Jonas Kusch's avatar
Jonas Kusch committed
19

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
20
21
#include <mpi.h>

22
SolverBase::SolverBase( Config* settings ) {
23
24
    _settings = settings;

Jonas Kusch's avatar
Jonas Kusch committed
25
26
    // @TODO save parameters from settings class

27
    // build mesh and store  and store frequently used params
28
29
30
31
32
    _mesh      = LoadSU2MeshFromFile( settings );
    _areas     = _mesh->GetCellAreas();
    _neighbors = _mesh->GetNeighbours();
    _normals   = _mesh->GetNormals();
    _nCells    = _mesh->GetNumCells();
33
    _settings->SetNCells( _nCells );
Jonas Kusch's avatar
Jonas Kusch committed
34

35
    // build quadrature object and store frequently used params
36
    _quadrature = QuadratureBase::Create( settings );
37
    _nq         = _quadrature->GetNq();
38
    _settings->SetNQuadPoints( _nq );
Jonas Kusch's avatar
Jonas Kusch committed
39

40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
    // build slope related params
    _reconstructor = new Reconstructor( settings );
    _reconsOrder   = _reconstructor->GetReconsOrder();

    auto nodes = _mesh->GetNodes();
    auto cells = _mesh->GetCells();
    std::vector<std::vector<Vector>> interfaceMidPoints( _nCells, std::vector<Vector>( _mesh->GetNumNodesPerCell(), Vector( 2, 1e-10 ) ) );
    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
        for( unsigned k = 0; k < _mesh->GetDim(); ++k ) {
            for( unsigned j = 0; j < _neighbors[idx_cell].size() - 1; ++j ) {
                interfaceMidPoints[idx_cell][j][k] = 0.5 * ( nodes[cells[idx_cell][j]][k] + nodes[cells[idx_cell][j + 1]][k] );
            }
            interfaceMidPoints[idx_cell][_neighbors[idx_cell].size() - 1][k] =
                0.5 * ( nodes[cells[idx_cell][_neighbors[idx_cell].size() - 1]][k] + nodes[cells[idx_cell][0]][k] );
        }
    }
    _interfaceMidPoints = interfaceMidPoints;
    _cellMidPoints      = _mesh->GetCellMidPoints();

    _psiDx = VectorVector( _nCells, Vector( _nq, 0.0 ) );
    _psiDy = VectorVector( _nCells, Vector( _nq, 0.0 ) );

62
    // set time step
63
64
    _dE        = ComputeTimeStep( settings->GetCFL() );
    _nEnergies = unsigned( settings->GetTEnd() / _dE );
jannick.wolters's avatar
jannick.wolters committed
65
66
    _energies.resize( _nEnergies );
    for( unsigned i = 0; i < _nEnergies; ++i ) _energies[i] = ( i + 1 ) * _dE;
Jonas Kusch's avatar
Jonas Kusch committed
67

68
    // setup problem  and store frequently used params
69
    _problem = ProblemBase::Create( _settings, _mesh );
70
    _sol     = _problem->SetupIC();
71
72
    _solNew  = _sol;    // setup temporary sol variable

Jonas Kusch's avatar
Jonas Kusch committed
73
74
75
    _sigmaT = _problem->GetTotalXS( _energies );
    _sigmaS = _problem->GetScatteringXS( _energies );
    _Q      = _problem->GetExternalSource( _energies );
76
77

    // setup numerical flux
steffen.schotthoefer's avatar
steffen.schotthoefer committed
78
    _g = NumericalFlux::Create();
79
80

    // boundary type
81
    _boundaryCells = _mesh->GetBoundaryTypes();
steffen.schotthoefer's avatar
steffen.schotthoefer committed
82
83

    // Solver Output
84
    _solverOutput.resize( _nCells );    // LEGACY! Only used for CSD SN
jonas.kusch's avatar
jonas.kusch committed
85

86
87
    PrepareScreenOutput();     // Screen Output
    PrepareHistoryOutput();    // History Output
88
89
90
91

    // initialize Helper Variables
    _fluxNew = Vector( _nCells, 0 );
    _flux    = Vector( _nCells, 0 );
92
93
94
95

    // write density
    _density = _problem->GetDensity( _mesh->GetCellMidPoints() );
    //_density = std::vector( _mesh->GetCellMidPoints().size(), 0.0 );
Jonas Kusch's avatar
Jonas Kusch committed
96
97
}

98
SolverBase::~SolverBase() {
99
    delete _quadrature;
100
101
102
103
    delete _mesh;
    delete _problem;
}

104
SolverBase* SolverBase::Create( Config* settings ) {
105
106
    switch( settings->GetSolverName() ) {
        case SN_SOLVER: return new SNSolver( settings );
107

108
        case PN_SOLVER: return new PNSolver( settings );
109
        case MN_SOLVER: return new MNSolver( settings );
110
111
112
113
        case CSD_SN_SOLVER: return new CSDSNSolver( settings );
        case CSD_SN_FOKKERPLANCK_TRAFO_SOLVER: return new CSDSolverTrafoFP( settings );
        case CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D: return new CSDSolverTrafoFP2D( settings );
        case CSD_SN_FOKKERPLANCK_TRAFO_SH_SOLVER_2D: return new CSDSolverTrafoFPSH2D( settings );
jannick.wolters's avatar
jannick.wolters committed
114
        case CSD_PN_SOLVER: return new CSDPNSolver( settings );
115
        default: ErrorMessages::Error( "Creator for the chosen solver does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
116
    }
117
118
    ErrorMessages::Error( "Creator for the chosen solver does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
    return nullptr;    // This code is never reached. Just to disable compiler warnings.
119
}
steffen.schotthoefer's avatar
steffen.schotthoefer committed
120

121
void SolverBase::Solve() {
122

123
124
    // --- Preprocessing ---
    PrepareVolumeOutput();
125

126
    DrawPreSolverOutput();
127

128
129
130
131
    // Adjust maxIter, depending if we have a normal run or a csd Run
    _maxIter = _nEnergies;
    if( _settings->GetIsCSD() ) {
        _maxIter = _nEnergies - 1;    // Since CSD does not go the last energy step
132
133
    }

134
135
136
    // Preprocessing before first pseudo time step
    SolverPreprocessing();

137
    // Loop over energies (pseudo-time of continuous slowing down approach)
138
    for( unsigned iter = 0; iter < _maxIter; iter++ ) {
139
140

        // --- Prepare Boundaries and temp variables
141
        IterPreprocessing( iter );
142
143
144
145
146
147
148

        // --- Compute Fluxes ---
        FluxUpdate();

        // --- Finite Volume Update ---
        FVMUpdate( iter );

149
        // --- Iter Postprocessing ---
jonas.kusch's avatar
jonas.kusch committed
150
        IterPostprocessing( iter );
151

152
        // --- Solver Output ---
153
        WriteVolumeOutput( iter );
154
155
156
157
        WriteScalarOutput( iter );
        PrintScreenOutput( iter );
        PrintHistoryOutput( iter );
        PrintVolumeOutput( iter );
158
    }
159
160
161
162

    // --- Postprocessing ---

    DrawPostSolverOutput();
163
164
}

165
void SolverBase::PrintVolumeOutput() const { ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); }
166

167
void SolverBase::PrintVolumeOutput( int currEnergy ) const {
168
    if( _settings->GetVolumeOutputFrequency() != 0 && currEnergy % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) {
169
        ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), _outputFields, _outputFieldNames, _mesh );
steffen.schotthoefer's avatar
steffen.schotthoefer committed
170
    }
171
    if( currEnergy == (int)_maxIter - 1 ) {    // Last iteration write without suffix.
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
172
173
        ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh );
    }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
174
}
175

176
// --- Helper ---
177
double SolverBase::ComputeTimeStep( double cfl ) const {
178
179
180
    double maxEdge = -1.0;
    for( unsigned j = 0; j < _nCells; j++ ) {
        for( unsigned l = 0; l < _normals[j].size(); l++ ) {
181
            double currentEdge = _areas[j] / norm( _normals[j][l] );
182
            if( currentEdge > maxEdge ) maxEdge = currentEdge;
183
184
185
186
187
        }
    }
    return cfl * maxEdge;
}

188
// --- IO ----
189
void SolverBase::PrepareScreenOutput() {
190
191
192
193
194
195
196
197
198
199
200
    unsigned nFields = (unsigned)_settings->GetNScreenOutput();

    _screenOutputFieldNames.resize( nFields );
    _screenOutputFields.resize( nFields );

    // Prepare all output Fields ==> Specified in option SCREEN_OUTPUT
    for( unsigned idx_field = 0; idx_field < nFields; idx_field++ ) {
        // Prepare all Output Fields per group

        // Different procedure, depending on the Group...
        switch( _settings->GetScreenOutput()[idx_field] ) {
201
            case MASS: _screenOutputFieldNames[idx_field] = "Mass"; break;
202

203
            case ITER: _screenOutputFieldNames[idx_field] = "Iter"; break;
204

205
            case RMS_FLUX: _screenOutputFieldNames[idx_field] = "RMS flux"; break;
206

207
208
209
210
211
            case VTK_OUTPUT: _screenOutputFieldNames[idx_field] = "VTK out"; break;

            case CSV_OUTPUT: _screenOutputFieldNames[idx_field] = "CSV out"; break;

            default: ErrorMessages::Error( "Screen output field not defined!", CURRENT_FUNCTION ); break;
212
        }
213
214
    }
}
steffen.schotthoefer's avatar
steffen.schotthoefer committed
215

216
void SolverBase::WriteScalarOutput( unsigned iteration ) {
steffen.schotthoefer's avatar
steffen.schotthoefer committed
217

218
219
220
    unsigned nFields = (unsigned)_settings->GetNScreenOutput();
    double mass      = 0.0;

221
    // -- Screen Output
222
223
224
225
226
227
228
229
230
231
232
    for( unsigned idx_field = 0; idx_field < nFields; idx_field++ ) {
        // Prepare all Output Fields per group
        // Different procedure, depending on the Group...
        switch( _settings->GetScreenOutput()[idx_field] ) {
            case MASS:
                for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                    mass += _fluxNew[idx_cell] * _areas[idx_cell];
                }
                _screenOutputFields[idx_field] = mass;
                break;

233
            case ITER: _screenOutputFields[idx_field] = iteration; break;
234

235
            case RMS_FLUX: _screenOutputFields[idx_field] = blaze::l2Norm( _fluxNew - _flux ); break;
236

237
238
239
            case VTK_OUTPUT:
                _screenOutputFields[idx_field] = 0;
                if( ( _settings->GetVolumeOutputFrequency() != 0 && iteration % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
240
                    ( iteration == _maxIter - 1 ) /* need sol at last iteration */ ) {
241
242
243
244
                    _screenOutputFields[idx_field] = 1;
                }
                break;

245
246
247
            case CSV_OUTPUT:
                _screenOutputFields[idx_field] = 0;
                if( ( _settings->GetHistoryOutputFrequency() != 0 && iteration % (unsigned)_settings->GetHistoryOutputFrequency() == 0 ) ||
248
                    ( iteration == _maxIter - 1 ) /* need sol at last iteration */ ) {
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
                    _screenOutputFields[idx_field] = 1;
                }
                break;
            default: ErrorMessages::Error( "Screen output group not defined!", CURRENT_FUNCTION ); break;
        }
    }

    // --- History output ---
    nFields = (unsigned)_settings->GetNHistoryOutput();

    std::vector<SCALAR_OUTPUT> screenOutputFields = _settings->GetScreenOutput();
    for( unsigned idx_field = 0; idx_field < nFields; idx_field++ ) {

        // Check first, if the field was already filled by screenoutput writer!
        std::vector<SCALAR_OUTPUT>::iterator itScreenOutput =
            std::find( screenOutputFields.begin(), screenOutputFields.end(), _settings->GetHistoryOutput()[idx_field] );

        // Prepare all Output Fields per group
        // Different procedure, depending on the Group...
        switch( _settings->GetHistoryOutput()[idx_field] ) {
            case MASS:
                if( screenOutputFields.end() == itScreenOutput ) {
                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
                        mass += _fluxNew[idx_cell] * _areas[idx_cell];
                    }
                    _historyOutputFields[idx_field] = mass;
                }
                else {
                    _historyOutputFields[idx_field] = *itScreenOutput;
                }
                break;

            case ITER: _historyOutputFields[idx_field] = iteration; break;

            case RMS_FLUX:
                if( screenOutputFields.end() == itScreenOutput ) {
                    _screenOutputFields[idx_field] = blaze::l2Norm( _fluxNew - _flux );
                }
                else {
                    _historyOutputFields[idx_field] = *itScreenOutput;
                }
                break;

            case VTK_OUTPUT:
                _historyOutputFields[idx_field] = 0;
                if( ( _settings->GetVolumeOutputFrequency() != 0 && iteration % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
295
                    ( iteration == _maxIter - 1 ) /* need sol at last iteration */ ) {
296
297
298
299
300
301
302
                    _historyOutputFields[idx_field] = 1;
                }
                break;

            case CSV_OUTPUT:
                _historyOutputFields[idx_field] = 0;
                if( ( _settings->GetHistoryOutputFrequency() != 0 && iteration % (unsigned)_settings->GetHistoryOutputFrequency() == 0 ) ||
303
                    ( iteration == _maxIter - 1 ) /* need sol at last iteration */ ) {
304
305
306
307
308
                    _historyOutputFields[idx_field] = 1;
                }
                break;

            default: ErrorMessages::Error( "History output group not defined!", CURRENT_FUNCTION ); break;
309
310
        }
    }
311
    _flux = _fluxNew;
312
313
}

314
void SolverBase::PrintScreenOutput( unsigned iteration ) {
315
316
    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
317
318
319
320
321
322
323
    auto log = spdlog::get( "event" );

    unsigned strLen  = 10;    // max width of one column
    char paddingChar = ' ';

    // assemble the line to print
    std::string lineToPrint = "| ";
324
    std::string tmp;
325
    for( unsigned idx_field = 0; idx_field < _settings->GetNScreenOutput(); idx_field++ ) {
326
        tmp = std::to_string( _screenOutputFields[idx_field] );
327
328
329
330

        // Format outputs correctly
        std::vector<SCALAR_OUTPUT> integerFields    = { ITER };
        std::vector<SCALAR_OUTPUT> scientificFields = { RMS_FLUX, MASS };
331
        std::vector<SCALAR_OUTPUT> booleanFields    = { VTK_OUTPUT, CSV_OUTPUT };
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347

        if( !( integerFields.end() == std::find( integerFields.begin(), integerFields.end(), _settings->GetScreenOutput()[idx_field] ) ) ) {
            tmp = std::to_string( (int)_screenOutputFields[idx_field] );
        }
        else if( !( booleanFields.end() == std::find( booleanFields.begin(), booleanFields.end(), _settings->GetScreenOutput()[idx_field] ) ) ) {
            tmp = "no";
            if( (bool)_screenOutputFields[idx_field] ) tmp = "yes";
        }
        else if( !( scientificFields.end() ==
                    std::find( scientificFields.begin(), scientificFields.end(), _settings->GetScreenOutput()[idx_field] ) ) ) {

            std::stringstream ss;
            ss << _screenOutputFields[idx_field];
            tmp = ss.str();
            tmp.erase( std::remove( tmp.begin(), tmp.end(), '+' ), tmp.end() );    // removing the '+' sign
        }
348

349
350
351
352
353
354
355
        if( strLen > tmp.size() )    // Padding
            tmp.insert( 0, strLen - tmp.size(), paddingChar );
        else if( strLen < tmp.size() )    // Cutting
            tmp.resize( strLen );

        lineToPrint += tmp + " |";
    }
356
    if( rank == 0 ) {
357
358
359
        if( _settings->GetScreenOutputFrequency() != 0 && iteration % (unsigned)_settings->GetScreenOutputFrequency() == 0 ) {
            log->info( lineToPrint );
        }
360
        else if( iteration == _maxIter - 1 ) {    // Always print last iteration
361
362
            log->info( lineToPrint );
        }
363
    }
364
365
}

366
void SolverBase::PrepareHistoryOutput() {
367
    unsigned nFields = (unsigned)_settings->GetNHistoryOutput();
368

369
370
    _historyOutputFieldNames.resize( nFields );
    _historyOutputFields.resize( nFields );
371

372
373
374
    // Prepare all output Fields ==> Specified in option SCREEN_OUTPUT
    for( unsigned idx_field = 0; idx_field < nFields; idx_field++ ) {
        // Prepare all Output Fields per group
375

376
        // Different procedure, depending on the Group...
377
        switch( _settings->GetHistoryOutput()[idx_field] ) {
378
            case MASS: _historyOutputFieldNames[idx_field] = "Mass"; break;
379

380
            case ITER: _historyOutputFieldNames[idx_field] = "Iter"; break;
381

382
            case RMS_FLUX: _historyOutputFieldNames[idx_field] = "RMS_flux"; break;
383

384
            case VTK_OUTPUT: _historyOutputFieldNames[idx_field] = "VTK_out"; break;
385

386
            case CSV_OUTPUT: _historyOutputFieldNames[idx_field] = "CSV_out"; break;
387

388
            default: ErrorMessages::Error( "History output field not defined!", CURRENT_FUNCTION ); break;
389
        }
390
    }
391
}
392

393
void SolverBase::PrintHistoryOutput( unsigned iteration ) {
394
395
396
397
398
399
400
    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
    auto log = spdlog::get( "tabular" );

    // assemble the line to print
    std::string lineToPrint = "";
    std::string tmp;
401
    for( int idx_field = 0; idx_field < _settings->GetNScreenOutput() - 1; idx_field++ ) {
402
403
404
405
406
407
408
409
410
411
        tmp = std::to_string( _screenOutputFields[idx_field] );
        lineToPrint += tmp + ",";
    }
    tmp = std::to_string( _screenOutputFields[_settings->GetNScreenOutput() - 1] );
    lineToPrint += tmp;    // Last element without comma

    if( rank == 0 ) {
        if( _settings->GetHistoryOutputFrequency() != 0 && iteration % (unsigned)_settings->GetHistoryOutputFrequency() == 0 ) {
            log->info( lineToPrint );
        }
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
412
        else if( iteration == _nEnergies - 1 ) {    // Always print last iteration
413
414
            log->info( lineToPrint );
        }
steffen.schotthoefer's avatar
steffen.schotthoefer committed
415
416
    }
}
417

418
void SolverBase::DrawPreSolverOutput() {
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
    // MPI
    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );

    // Logger
    auto log    = spdlog::get( "event" );
    auto logCSV = spdlog::get( "tabular" );

    std::string hLine = "--";

    if( rank == 0 ) {
        unsigned strLen  = 10;    // max width of one column
        char paddingChar = ' ';

        // Assemble Header for Screen Output
        std::string lineToPrint = "| ";
        std::string tmpLine     = "------------";
        for( unsigned idxFields = 0; idxFields < _settings->GetNScreenOutput(); idxFields++ ) {
            std::string tmp = _screenOutputFieldNames[idxFields];

            if( strLen > tmp.size() )    // Padding
                tmp.insert( 0, strLen - tmp.size(), paddingChar );
            else if( strLen < tmp.size() )    // Cutting
                tmp.resize( strLen );

            lineToPrint += tmp + " |";
            hLine += tmpLine;
        }
        log->info( "---------------------------- Solver Starts -----------------------------" );
        log->info( "| The simulation will run for {} iterations.", _nEnergies );
        log->info( hLine );
        log->info( lineToPrint );
        log->info( hLine );

        std::string lineToPrintCSV = "";
        for( int idxFields = 0; idxFields < _settings->GetNHistoryOutput() - 1; idxFields++ ) {
            std::string tmp = _historyOutputFieldNames[idxFields];
            lineToPrintCSV += tmp + ",";
        }
        lineToPrintCSV += _historyOutputFieldNames[_settings->GetNHistoryOutput() - 1];
        logCSV->info( lineToPrintCSV );
    }
}

463
void SolverBase::DrawPostSolverOutput() {
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
    // MPI
    int rank;
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );

    // Logger
    auto log = spdlog::get( "event" );

    std::string hLine = "--";

    if( rank == 0 ) {
        unsigned strLen  = 10;    // max width of one column
        char paddingChar = ' ';

        // Assemble Header for Screen Output
        std::string lineToPrint = "| ";
        std::string tmpLine     = "------------";
        for( unsigned idxFields = 0; idxFields < _settings->GetNScreenOutput(); idxFields++ ) {
            std::string tmp = _screenOutputFieldNames[idxFields];

            if( strLen > tmp.size() )    // Padding
                tmp.insert( 0, strLen - tmp.size(), paddingChar );
            else if( strLen < tmp.size() )    // Cutting
                tmp.resize( strLen );

            lineToPrint += tmp + " |";
            hLine += tmpLine;
        }
        log->info( hLine );
        log->info( "| Postprocessing screen output goes here." );
        log->info( "--------------------------- Solver Finished ----------------------------" );
    }
}

497
void SolverBase::SolverPreprocessing() {}