abstracted solver routine to base class, created option for screen output,...

abstracted solver routine to base class, created option for screen output, introduced solver postprocessing. Done for SN Pn and MN


Former-commit-id: c8ec6724

code/include/common/config.h

@@ -93,7 +93,11 @@ class Config
     // Output Options
     unsigned short _nVolumeOutput;            /*!< @brief Number of volume outputs */
     std::vector<VOLUME_OUTPUT> _volumeOutput; /*!< @brief Output groups for volume output*/
-    unsigned short _outputFrequency;          /*!< @brief Frequency of vtk write of volume output*/
+    unsigned short _volumeOutputFrequency;    /*!< @brief Frequency of vtk write of volume output*/
+
+    unsigned short _nScreenOutput;            /*!< @brief Number of screen outputs */
+    std::vector<SCREEN_OUTPUT> _screenOutput; /*!< @brief Output groups for screen output*/
+    unsigned short _screenOutputFrequency;    /*!< @brief Frequency of screen output*/
 
     // --- Parsing Functionality and Initializing of Options ---
     /*!
@@ -263,8 +267,11 @@ class Config
     // Output Structure
     std::vector<VOLUME_OUTPUT> inline GetVolumeOutput() { return _volumeOutput; }
     unsigned short inline GetNVolumeOutput() { return _nVolumeOutput; }
-    unsigned short inline GetOutputFrequency() { return _outputFrequency; }
+    unsigned short inline GetVolumeOutputFrequency() { return _volumeOutputFrequency; }
 
+    std::vector<SCREEN_OUTPUT> inline GetScreenOutput() { return _screenOutput; }
+    unsigned short inline GetNScreenOutput() { return _nScreenOutput; }
+    unsigned short inline GetScreenOutputFrequency() { return _screenOutputFrequency; }
     // ---- Setters for option structure
 
     // Quadrature Structure

code/include/common/globalconstants.h

@@ -81,4 +81,9 @@ enum VOLUME_OUTPUT { ANALYTIC, MINIMAL, MOMENTS, DUAL_MOMENTS, DOSE };
 inline std::map<std::string, VOLUME_OUTPUT> VolOutput_Map{
     { "ANALYTIC", ANALYTIC }, { "MINIMAL", MINIMAL }, { "MOMENTS", MOMENTS }, { "DUAL_MOMENTS", DUAL_MOMENTS }, { "DOSE", DOSE } };
 
+// Screen output
+enum SCREEN_OUTPUT { ITER, MASS, RMS_FLUX };
+
+inline std::map<std::string, SCREEN_OUTPUT> ScreenOutput_Map{ { "ITER", ITER }, { "MASS", MASS }, { "RMS_FLUX", RMS_FLUX } };
+
 #endif    // GLOBAL_CONSTANTS_H

code/include/solvers/mnsolver.h

@@ -19,8 +19,6 @@ class MNSolver : public Solver
     /*! @brief MNSolver destructor */
     ~MNSolver();
 
-    void Solve() override; /*! @brief Solve functions runs main time loop */
-
   private:
     unsigned _nTotalEntries;    /*! @brief: Total number of equations in the system */
     unsigned short _LMaxDegree; /*! @brief: Max Order of Moments */
@@ -84,5 +82,9 @@ class MNSolver : public Solver
     void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) override;
     void FluxUpdate( VectorVector& psiNew ) override;
     void IterPreprocessing() override;
+    void IterPostprocessing();
+
+    // Helper
+    void ComputeRadFlux();
 };
 #endif    // MNSOLVER_H

code/include/solvers/pnsolver.h

@@ -11,8 +11,6 @@ class PNSolver : public Solver
      */
     PNSolver( Config* settings );
 
-    virtual void Solve() override; /*! @brief Solve functions runs main time loop. (Run Solver) */
-
   protected:
     unsigned _nTotalEntries; /*! @brief: total number of equations in the system */
     unsigned _LMaxDegree;    /*! @brief: maximal degree of the spherical harmonics basis*/
@@ -47,14 +45,17 @@ class PNSolver : public Solver
     void PrepareOutputFields() override;
 
     /*! @brief Function that prepares VTK export and csv export of the current solver iteration
-        @returns: Mass of current iteration
-    */
+        @returns: Mass of current iteration     */
     double WriteOutputFields( unsigned idx_pseudoTime ) override;
 
     // Solver
     void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) override;
     void FluxUpdate( VectorVector& psiNew ) override;
     void IterPreprocessing() override;
+    void IterPostprocessing();
+
+    // Helper
+    void ComputeRadFlux();
 
     // Initialization
     /*! @brief: parameter functions for setting up system matrix

code/include/solvers/snsolver.h

@@ -14,15 +14,10 @@ class SNSolver : public Solver
     Vector _weights;          /*!  @brief quadrature weights, dim(_weights) = (_nq) */
 
   public:
-    /**
-     * @brief SNSolver constructor
+    /*! @brief SNSolver constructor
      * @param settings stores all needed information
      */
     SNSolver( Config* settings );
-    /**
-     * @brief Solve functions runs main time loop
-     */
-    void Solve() override;
 
   private:
     void PrepareOutputFields() override;
@@ -32,6 +27,10 @@ class SNSolver : public Solver
     void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) override;
     void FluxUpdate( VectorVector& psiNew ) override;
     void IterPreprocessing() override;
+    void IterPostprocessing() override;
+
+    // Helper
+    void ComputeRadFlux();
 
     // --- Member variables ---
 };

code/include/solvers/solverbase.h

@@ -5,6 +5,9 @@
 #include "common/globalconstants.h"
 #include "common/typedef.h"
 
+// externals
+#include "spdlog/spdlog.h"
+
 // Forward Declarations
 class NumericalFlux;
 class Mesh;
@@ -53,12 +56,19 @@ class Solver
     std::vector<double> _solverOutput; /*! @brief LEGACY: Outputfield for solver ==> Will be replaced by _outputFields in the near future */
 
     // Output related members
-    std::vector<std::vector<std::vector<double>>> _outputFields; /*! @brief: Solver Output: dimensions (GroupID,FieldID,CellID). !Protoype output for
-                                                                    multiple output fields. Will replace _solverOutput */
+    std::vector<std::vector<std::vector<double>>> _outputFields; /*! @brief: Solver Output: dimensions (GroupID,FieldID,CellID).*/
     std::vector<std::vector<std::string>> _outputFieldNames;     /*! @brief: Names of the outputFields: dimensions (GroupID,FieldID) */
-
     // we will have to add a further dimension for quadPoints and weights once we start with multilevel SN
 
+    // Output related members
+    std::vector<double> _screenOutputFields;          /*! @brief: Solver Output: dimensions (FieldID). */
+    std::vector<std::string> _screenOutputFieldNames; /*! @brief: Names of the outputFields: dimensions (FieldID) */
+
+    // Internal Members
+    VectorVector _solNew; /*! @brief: VectorVector to store the new flux and later the new solution per iteration */    // REPLACES psiNEW
+    Vector _fluxNew; /*! @brief: Vector to store the new Flux */
+    Vector _flux;    /*! @brief: Vector to store the old Flux */
+
     // ---- Member functions ----
 
     // IO
@@ -66,27 +76,38 @@ class Solver
     virtual void PrepareOutputFields() = 0;
 
     /*! @brief Function that prepares VTK export and csv export of the current solver iteration
-        @returns: Mass of current iteration
-    */
+        @returns: Mass of current iteration    */
     virtual double WriteOutputFields( unsigned idx_pseudoTime ) = 0;
 
-    // Solver
-    ///*! @brief Computes the finite volume update for all cells and stores it in psiNew    */
-    // virtual void FVMUpdate( VectorVector& psiNew ) = 0;
+    /*! @brief: Initialized the output fields and their Names for the Screenoutput */
+    void PrepareScreenOutputFields();
+
+    /*! @brief Function that Screen Output and prints to Screen/Logger */
+    void WriteScreenOutputFields( unsigned idx_pseudoTime );
+
+    /*! @brief Prints ScreenOutputFields to Screen and to logger */
+    void PrintScreen( std::shared_ptr<spdlog::logger> log );
 
+    // Solver
     /*! @brief Performs preprocessing for the current solver iteration */
     virtual void IterPreprocessing() = 0;
+    /*! @brief Performs postprocessing for the current solver iteration */
+    virtual void IterPostprocessing() = 0;
     /*! @brief Constructs  the flux update for the current iteration and stores it in psiNew*/
     virtual void FluxUpdate( VectorVector& psiNew ) = 0;
     /*! @brief Computes the finite Volume update step for the current iteration */
     virtual void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) = 0;
 
+    // Helper
     /**
      * @brief ComputeTimeStep calculates the maximal stable time step
      * @param cfl is cfl number
      */
     double ComputeTimeStep( double cfl ) const;
 
+    /*! @brief: Computes the flux of the solution to check conservation properties */
+    virtual void ComputeRadFlux() = 0;
+
   public:
     /**
      * @brief Solver constructor
@@ -106,7 +127,7 @@ class Solver
     /**
      * @brief Solve functions runs main time loop
      */
-    virtual void Solve() = 0;
+    virtual void Solve();
 
     void Save() const;                 /*! @brief Save Output solution to VTK file */
     void Save( int currEnergy ) const; /*! @brief Save Output solution at given energy (pseudo time) to VTK file */

code/src/common/config.cpp

@@ -279,8 +279,12 @@ void Config::SetConfigOptions() {
     // Output related options
     /*! @brief Volume output \n DESCRIPTION: Describes output groups to write to vtk \ingroup Config */
     AddEnumListOption( "VOLUME_OUTPUT", _nVolumeOutput, _volumeOutput, VolOutput_Map );
-    /*! @brief Output Frequency \n DESCRIPTION: Describes output write frequency \n DEFAULT 0 ,i.e. only last value \ingroup Config */
-    AddUnsignedShortOption( "OUTPUT_FREQUENCY", _outputFrequency, 0 );
+    /*! @brief Volume Output Frequency \n DESCRIPTION: Describes output write frequency \n DEFAULT 0 ,i.e. only last value \ingroup Config */
+    AddUnsignedShortOption( "OUTPUT_FREQUENCY", _volumeOutputFrequency, 0 );
+    /*! @brief Screen output \n DESCRIPTION: Describes scren output groups \ingroup Config */
+    AddEnumListOption( "SCREEN_OUTPUT", _nScreenOutput, _screenOutput, ScreenOutput_Map );
+    /*! @brief Screen Output Frequency \n DESCRIPTION: Describes screen output write frequency \n DEFAULT 0 ,i.e. only last value \ingroup Config */
+    AddUnsignedShortOption( "SCREEN_OUTPUT_FREQUENCY", _screenOutputFrequency, 0 );
 }
 
 void Config::SetConfigParsing( string case_filename ) {
@@ -497,6 +501,14 @@ void Config::SetPostprocessing() {
         _nVolumeOutput = 1;
         _volumeOutput.push_back( MINIMAL );
     }
+
+    // Set default screen output
+    if( _nScreenOutput == 0 ) {
+        _nScreenOutput = 3;
+        _screenOutput.push_back( ITER );
+        _screenOutput.push_back( RMS_FLUX );
+        _screenOutput.push_back( MASS );
+    }
 }
 
 void Config::SetOutput() {

code/src/solvers/csdsnsolver.cpp

@@ -198,7 +198,7 @@ double CSDSNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
         mass += flux[idx_cell] * _areas[idx_cell];
     }
 
-    if( ( _settings->GetOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetOutputFrequency() == 0 ) ||
+    if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
         ( idx_pseudoTime == _nEnergies - 1 ) /* need sol at last iteration */ ) {
 
         for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {

code/src/solvers/mnsolver.cpp

@@ -125,44 +125,6 @@ void MNSolver::ComputeRealizableSolution( unsigned idx_cell ) {
     }
 }
 
-void MNSolver::Solve() {
-
-    int rank;
-    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
-
-    auto log = spdlog::get( "event" );
-
-    VectorVector psiNew = _sol;
-
-    double mass = 0;
-
-    if( rank == 0 ) log->info( "{:10}  {:10}", "t", "mass" );
-
-    // Loop over energies (pseudo-time of continuous slowing down approach)
-    for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
-
-        // ----- Iteration Preprocessing -----
-        IterPreprocessing();
-
-        // ------- Flux Computation Step --------------
-        FluxUpdate( psiNew );
-
-        // ------ Finite Volume Update Step ------
-        FVMUpdate( psiNew, idx_energy );
-
-        // Update Solution
-        _sol = psiNew;
-
-        // --- VTK and CSV Output ---
-        mass = WriteOutputFields( idx_energy );
-        Save( idx_energy );
-        // WriteNNTrainingData( idx_energy );
-
-        // --- Screen Output ---
-        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e}", _energies[idx_energy], mass );
-    }
-}
-
 void MNSolver::IterPreprocessing() {
 
     // ------- Reconstruction Step ----------------
@@ -175,6 +137,21 @@ void MNSolver::IterPreprocessing() {
     }
 }
 
+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;
+    }
+}
+
 void MNSolver::FluxUpdate( VectorVector& psiNew ) {
     // Loop over the grid cells
     for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
@@ -278,14 +255,14 @@ double MNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
     }
     mass *= firstMomentScaleFactor;
 
-    if( ( _settings->GetOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetOutputFrequency() == 0 ) ||
+    if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
         ( idx_pseudoTime == _nEnergies - 1 ) /* need sol at last iteration */ ) {
 
         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] = firstMomentScaleFactor * _sol[idx_cell][0];
+                        _outputFields[idx_group][0][idx_cell] = _fluxNew[idx_cell];
                     }
                     break;
                 case MOMENTS:

code/src/solvers/pnsolver.cpp

@@ -51,45 +51,23 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
     PrepareOutputFields();
 }
 
-void PNSolver::Solve() {
-
-    int rank;
-    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
-
-    auto log = spdlog::get( "event" );
-
-    VectorVector psiNew = _sol;
-
-    double mass = 0;
-
-    if( rank == 0 ) log->info( "{:10}   {:10}", "t", "mass" );
-
-    // Loop over energies (pseudo-time of continuous slowing down approach)
-    for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
-
-        // --- Prepare Boundaries and temp variables
-        IterPreprocessing();
-
-        // --- Compute Fluxes ---
-        FluxUpdate( psiNew );
-
-        // --- Finite Volume Update ---
-        FVMUpdate( psiNew, idx_energy );
-
-        // --- Update Solution ---
-        _sol = psiNew;
+void PNSolver::IterPreprocessing() {
+    // Nothing to preprocess for PNSolver
+}
 
-        // --- VTK and CSV Output ---
-        mass = WriteOutputFields( idx_energy );
-        Save( idx_energy );
+void PNSolver::IterPostprocessing() {
+    // --- Update Solution ---
+    _sol = _solNew;
 
-        // --- Screen Output ---
-        if( rank == 0 ) log->info( "{:03.8f}  {:01.5e}", _energies[idx_energy], mass );
-    }
+    // --- Compute Flux for solution and Screen Output ---
+    ComputeRadFlux();
 }
 
-void PNSolver::IterPreprocessing() {
-    // Nothing to preprocess for PNSolver
+void PNSolver::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;
+    }
 }
 
 void PNSolver::FluxUpdate( VectorVector& psiNew ) {
@@ -392,13 +370,13 @@ double PNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
     }
     mass *= firstMomentScaleFactor;
 
-    if( ( _settings->GetOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetOutputFrequency() == 0 ) ||
+    if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
         ( idx_pseudoTime == _nEnergies - 1 ) /* need sol at last iteration */ ) {
         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] = firstMomentScaleFactor * _sol[idx_cell][0];
+                        _outputFields[idx_group][0][idx_cell] = _fluxNew[idx_cell];
                     }
                     break;
                 case MOMENTS:

code/src/solvers/snsolver.cpp

@@ -24,54 +24,22 @@ SNSolver::SNSolver( Config* settings ) : Solver( settings ) {
     PrepareOutputFields();
 }
 
-void SNSolver::Solve() {
-
-    int rank;
-    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
-
-    auto log = spdlog::get( "event" );
-
-    VectorVector psiNew = _sol;
-
-    double mass  = 0.0;
-    double dFlux = 1e10;
-    Vector fluxNew( _nCells, 0.0 );
-    Vector fluxOld( _nCells, 0.0 );
-
-    if( rank == 0 ) log->info( "{:10}   {:10}", "t", "dFlux", "mass" );
-
-    // loop over energies (pseudo-time)
-    for( unsigned idx_energy = 0; idx_energy < _nEnergies; ++idx_energy ) {
-
-        // --- Prepare Boundaries and temp variables
-        IterPreprocessing();
-
-        // --- Compute Fluxes ---
-        FluxUpdate( psiNew );
-
-        // --- Compute Finite Volume Step ---
-        FVMUpdate( psiNew, idx_energy );
-
-        // --- Update Solution ---
-        _sol = psiNew;
-
-        // --- VTK and CSV Output ---
-        mass = WriteOutputFields( idx_energy );
-        Save( idx_energy );
+void SNSolver::IterPreprocessing() {
+    // Nothing to do for SNSolver
+}
 
-        // --- Screen Output ---
-        for( unsigned i = 0; i < _nCells; ++i ) {
-            fluxNew[i] = dot( _sol[i], _weights );
-        }
+void SNSolver::IterPostprocessing() {
+    // --- Update Solution ---
+    _sol = _solNew;
 
-        dFlux   = blaze::l2Norm( fluxNew - fluxOld );
-        fluxOld = fluxNew;
-        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
-    }
+    // --- Compute Flux for solution and Screen Output ---
+    ComputeRadFlux();
 }
 
-void SNSolver::IterPreprocessing() {
-    // Nothing to do for SNSolver
+void SNSolver::ComputeRadFlux() {
+    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+        _fluxNew[idx_cell] = blaze::dot( _sol[idx_cell], _weights );
+    }
 }
 
 void SNSolver::FluxUpdate( VectorVector& psiNew ) {
@@ -263,14 +231,14 @@ double SNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
         mass += flux[idx_cell] * _areas[idx_cell];
     }
 
-    if( ( _settings->GetOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetOutputFrequency() == 0 ) ||
+    if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
         ( idx_pseudoTime == _nEnergies - 1 ) /* need sol at last iteration */ ) {
 
         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] = flux[idx_cell];
+                        _outputFields[idx_group][0][idx_cell] = _fluxNew[idx_cell];
                     }
                     break;
 

code/src/solvers/solverbase.cpp

@@ -36,6 +36,8 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
     // setup problem  and store frequently used params
     _problem = ProblemBase::Create( _settings, _mesh );
     _sol     = _problem->SetupIC();
+    _solNew  = _sol;    // setup temporary sol variable
+
     //_s       = _problem->GetStoppingPower( _energies );
     _sigmaT = _problem->GetTotalXS( _energies );
     _sigmaS = _problem->GetScatteringXS( _energies );
@@ -48,7 +50,13 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
     _boundaryCells = _mesh->GetBoundaryTypes();
 
     // Solver Output
-    _solverOutput.resize( _nCells );    // Currently only Flux
+    _solverOutput.resize( _nCells );    // LEGACY! Only used for CSD SN
+    // Screen Output
+    PrepareScreenOutputFields();
+
+    // initialize Helper Variables
+    _fluxNew = Vector( _nCells, 0 );
+    _flux    = Vector( _nCells, 0 );
 }
 
 Solver::~Solver() {
@@ -81,7 +89,97 @@ Solver* Solver::Create( Config* settings ) {
 void Solver::Save() const { ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); }
 
 void Solver::Save( int currEnergy ) const {
-    if( _settings->GetOutputFrequency() != 0 && currEnergy % (unsigned)_settings->GetOutputFrequency() == 0 ) {
+    if( _settings->GetVolumeOutputFrequency() != 0 && currEnergy % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) {
         ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), _outputFields, _outputFieldNames, _mesh );
     }
 }
+
+void Solver::PrepareScreenOutputFields() {
+    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] ) {
+            case MASS: _screenOutputFieldNames[idx_field] = "mass"; break;
+
+            case ITER: _screenOutputFieldNames[idx_field] = "iter"; break;
+
+            case RMS_FLUX: _screenOutputFieldNames[idx_field] = "RMS_flux"; break;
+
+            default: ErrorMessages::Error( "Screen Output Group not defined!", CURRENT_FUNCTION ); break;
+        }
+    }
+}
+
+void Solver::WriteScreenOutputFields( unsigned idx_pseudoTime ) {
+
+    unsigned nFields = (unsigned)_settings->GetNScreenOutput();
+    double mass      = 0.0;
+
+    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;
+
+            case ITER: _screenOutputFields[idx_field] = idx_pseudoTime; break;
+
+            case RMS_FLUX:
+                _screenOutputFields[idx_field] = blaze::l2Norm( _fluxNew - _flux );
+                _flux                          = _fluxNew;
+                break;
+
+            default: ErrorMessages::Error( "Screen Output Group not defined!", CURRENT_FUNCTION ); break;
+        }
+    }
+}
+
+void Solver::PrintScreen( std::shared_ptr<spdlog::logger> log ) {
+    log->info( "{:03.8f}   {:01.5e} {:01.5e}", _screenOutputFields[0], _screenOutputFields[1], _screenOutputFields[2] );
+}
+
+void Solver::Solve() {
+    int rank;
+    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
+
+    auto log = spdlog::get( "event" );
+
+    double mass = 0;
+
+    if( rank == 0 ) log->info( "{:10}   {:10}", "t", "mass" );
+
+    // Loop over energies (pseudo-time of continuous slowing down approach)
+    for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
+
+        // --- Prepare Boundaries and temp variables
+        IterPreprocessing();
+
+        // --- Compute Fluxes ---
+        FluxUpdate( _solNew );
+
+        // --- Finite Volume Update ---
+        FVMUpdate( _solNew, idx_energy );
+
+        // --- Postprocessing ---
+        IterPostprocessing();
+
+        // --- VTK and CSV Output ---
+        mass = WriteOutputFields( idx_energy );
+        Save( idx_energy );
+
+        // --- Screen Output ---
+        WriteScreenOutputFields( idx_energy );
+        PrintScreen( log );
+    }
+}