enabled multifield output for pn and mn solver. Unit tests not yet adapted!....

enabled multifield output for pn and mn solver. Unit tests not yet adapted!. PNSolver.solve still messy, needs cleanup!

code/include/solvers/mnsolver.h

@@ -50,11 +50,16 @@ class MNSolver : public Solver
     OptimizerBase* _optimizer; /*! @brief: Class to solve minimal entropy problem */
 
     // Output related members
-    std::vector<std::vector<double>> _outputFields; /*! @brief: Protoype output */
+    std::vector<std::vector<double>> _outputFields; /*! @brief: Protoype output for multiple output fields. Will replace _solverOutput */
 
     /*! @brief Function that writes NN Training Data in a .csv file */
     void WriteNNTrainingData( unsigned idx_pseudoTime );
 
+    /*! @brief Function that prepares VTK export and csv export of the current solver iteration
+        @returns: Mass of current iteration
+    */
+    double WriteOutputFields();
+
     // Member functions
     /*! @brief : computes the global index of the moment corresponding to basis function (l,k)
      *  @param : degree l, it must hold: 0 <= l <=_nq

code/include/solvers/pnsolver.h

@@ -41,6 +41,9 @@ class PNSolver : public Solver
 
     Vector _scatterMatDiag; /*! @brief: diagonal of the scattering matrix (its a diagonal matrix by construction) */
 
+    // Output related members
+    std::vector<std::vector<double>> _outputFields; /*! @brief: Protoype output for multiple output fields. Will replace _solverOutput */
+
     /*! @brief: parameter functions for setting up system matrix
      *  @param: degree l, it must hold: 0 <= l <=_nq
      *  @param : order k, it must hold: -l <=k <= l

code/include/solvers/solverbase.h

@@ -50,7 +50,7 @@ class Solver
 
     // Solution related members
     VectorVector _sol;                 /*! @brief solution of the PDE, e.g. angular flux or moments */
-    std::vector<double> _solverOutput; /*! @brief PROTOTYPE: Outputfield for solver */
+    std::vector<double> _solverOutput; /*! @brief PROTOTYPE: Outputfield for solver ==> Will be replaced by _outputFields in the near future */
 
     // we will have to add a further dimension for quadPoints and weights once we start with multilevel SN
 

code/input/exampleMN.cfg

@@ -28,11 +28,11 @@ CFL_NUMBER = 0.7
 % Final time for simulation
 TIME_FINAL = 0.3
 % Maximal Moment degree
-MAX_MOMENT_SOLVER = 0
+MAX_MOMENT_SOLVER = 1
 %
 %% Entropy settings
 ENTROPY_FUNCTIONAL = MAXWELL_BOLZMANN
-ENTROPY_OPTIMIZER = ML
+ENTROPY_OPTIMIZER = NEWTON
 %
 % ----- Newton Solver Specifications ----
 %

code/input/examplePN.cfg

@@ -14,8 +14,8 @@ OUTPUT_FILE = examplePN
 % Log directory
 LOG_DIR = ../result/logs
 % Mesh File
-MESH_FILE = linesource.su2
-%MESH_FILE = linesource_debug.su2
+%MESH_FILE = linesource.su2
+MESH_FILE = linesource_debug.su2
 
 %
 PROBLEM = LINESOURCE

code/src/common/io.cpp

@@ -96,21 +96,41 @@ void ExportVTK( const std::string fileName,
             grid->SetCells( VTK_QUAD, cellArray );
         }
 
-        for( unsigned i = 0; i < results.size(); i++ ) {
-            auto cellData = vtkDoubleArraySP::New();
-            cellData->SetName( fieldNames[i].c_str() );
-            switch( results[i].size() ) {
-                case 1:
-                    for( unsigned j = 0; j < numCells; j++ ) {
-                        cellData->InsertNextValue( results[i][0][j] );
-                    }
-                    break;
-                default:
-                    auto log = spdlog::get( "event" );
-                    ErrorMessages::Error( "Please implement output for results of size " + std::to_string( results[i].size() ) + "!",
-                                          CURRENT_FUNCTION );
+        for( unsigned idx_group = 0; idx_group < results.size(); idx_group++ ) {
+            // why is there a group dimension for results, but not for names?
+            // Right now, this is ok, since we only have one group... TODO!
+
+            for( unsigned idx_field = 0; idx_field < results[idx_group].size(); idx_field++ ) {    // Loop over all output fields
+
+                auto cellData = vtkDoubleArraySP::New();
+                cellData->SetName( fieldNames[idx_field].c_str() );
+
+                for( unsigned idx_cell = 0; idx_cell < numCells; idx_cell++ ) {
+                    cellData->InsertNextValue( results[idx_group][idx_field][idx_cell] );
+                }
+
+                /*
+                switch( results[idx_group].size() ) {
+                    case 1:
+                        for( unsigned j = 0; j < numCells; j++ ) {
+                            cellData->InsertNextValue( results[idx_group][0][j] );
+                        }
+                        break;
+                    default:
+                        // for( unsigned l = 0; l < results[i].size(); l++ ) {    // Loop over all output fields
+                        //    for( unsigned j = 0; j < numCells; j++ ) {
+                        //        cellData->InsertNextValue( results[i][l][j] );
+                        //    }
+                        // }
+
+                        auto log = spdlog::get( "event" );
+                        ErrorMessages::Error( "Please implement output for results of size " + std::to_string( results[idx_group].size() ) + "!",
+                                              CURRENT_FUNCTION );
+                        break;
+                }*/
+
+                grid->GetCellData()->AddArray( cellData );
             }
-            grid->GetCellData()->AddArray( cellData );
         }
 
         grid->SetPoints( pts );

code/src/solvers/mnsolver.cpp

@@ -132,24 +132,14 @@ void MNSolver::Solve() {
 
     // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
     VectorVector psiNew = _sol;
-    double dFlux        = 1e10;
-    Vector fluxNew( _nCells, 0.0 );
-    Vector fluxOld( _nCells, 0.0 );
-    VectorVector solTimesArea = _sol;
-
-    double mass1 = 0;
-    for( unsigned i = 0; i < _nCells; ++i ) {
-        _solverOutput[i] = _sol[i][0];
-        mass1 += _sol[i][0] * _areas[i];
-    }
 
-    dFlux   = blaze::l2Norm( fluxNew - fluxOld );
-    fluxOld = fluxNew;
+    double mass = 0;
 
     Save( -1 );
 
-    if( rank == 0 ) log->info( "{:10}   {:10}", "t", "dFlux" );
-    if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", -1.0, dFlux, mass1 );
+    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
 
     // Loop over energies (pseudo-time of continuous slowing down approach)
     for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
@@ -185,55 +175,42 @@ void MNSolver::Solve() {
             }
         }
 
-        // pseudo time iteration output
-        double mass = 0.0;
-        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];
-            }
-        }
+        // Update Solution
+        _sol = psiNew;
 
-        dFlux   = blaze::l2Norm( fluxNew - fluxOld );
-        fluxOld = fluxNew;
-        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
+        // --- VTK and CSV Output ---
+        mass = WriteOutputFields();
         Save( idx_energy );
-
         WriteNNTrainingData( idx_energy );
 
-        // Update Solution
-        _sol = psiNew;
+        // --- Screen Output ---
+        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[idx_energy], mass );
     }
 }
 
-void MNSolver::Save() const {
-    std::vector<std::string> fieldNames{ "flux" };
-    std::vector<double> flux;
-    flux.resize( _nCells );
+double MNSolver::WriteOutputFields() {
+    double mass = 0.0;
 
-    for( unsigned i = 0; i < _nCells; ++i ) {
-        flux[i] = _sol[i][0];
+    for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
+        for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+            mass += _sol[idx_cell][0] * _areas[idx_cell];
+            _outputFields[idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
+        }
     }
-    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" };
+    return mass;
+}
 
-    // Multifield ooutput
-    // std::vector<std::string> fieldNames( _nTotalEntries, "flux" );
-    // for( unsigned idx_sys = 1; idx_sys < _nTotalEntries; idx_sys++ ) {
-    //    fieldNames[idx_sys] = std::string( "flux_moment_" + std::to_string( idx_sys ) );
-    //}
+void MNSolver::Save() const { Save( _nEnergies - 1 ); }
 
-    std::vector<std::vector<double>> scalarField( 1, _solverOutput );
-    std::vector<std::vector<std::vector<double>>> results{ scalarField };
+void MNSolver::Save( int currEnergy ) const {
+    std::vector<std::string> fieldNames( _nTotalEntries, "flux" );
+    for( int idx_l = 0; idx_l <= (int)_nMaxMomentsOrder; idx_l++ ) {
+        for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
+            fieldNames[GlobalIndex( idx_l, idx_k )] = std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
+        }
+    }
+    std::vector<std::vector<std::vector<double>>> results{ _outputFields };
     ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
 }
 

code/src/solvers/pnsolver.cpp

@@ -71,6 +71,9 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
     // std::cout << "_AzMinus :\n" << _AzMinus << "\n ";
     //
     // std::cout << "_nCells: " << _nCells << "\n";
+
+    // Solver output
+    _outputFields = std::vector( _nTotalEntries, std::vector( _nCells, 0.0 ) );    // WILL REPLACE _solverOutput
 }
 
 void PNSolver::Solve() {
@@ -149,20 +152,26 @@ void PNSolver::Solve() {
                 }
             }
         }
+        // Update Solution
         _sol = psiNew;
 
+        // pseudo time iteration output
         double mass = 0.0;
-        for( unsigned i = 0; i < _nCells; ++i ) {
-            fluxNew[i]       = _sol[i][0];    // zeroth moment is raditation densitiy we are interested in
-            _solverOutput[i] = _sol[i][0];
-            mass += _sol[i][0] * _areas[i];
+        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];
+            }
         }
 
+        // Screen Output
         dFlux   = blaze::l2Norm( fluxNew - fluxOld );
         fluxOld = fluxNew;
-        if( rank == 0 ) {
-            log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
-        }
+        if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
         Save( idx_energy );
     }
 }
@@ -367,23 +376,16 @@ double PNSolver::LegendrePoly( double x, int l ) {    // Legacy. TO BE DELETED
     }
 }
 
-void PNSolver::Save() const {
-    std::vector<std::string> fieldNames{ "flux" };
-    std::vector<double> flux;
-    flux.resize( _nCells );
-
-    for( unsigned i = 0; i < _nCells; ++i ) {
-        flux[i] = _sol[i][0];
-    }
-    std::vector<std::vector<double>> scalarField( 1, flux );
-    std::vector<std::vector<std::vector<double>>> results{ scalarField };
-    ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
-}
+void PNSolver::Save() const { Save( _nEnergies - 1 ); }
 
 void PNSolver::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 };
+    std::vector<std::string> fieldNames( _nTotalEntries, "flux" );
+    for( int idx_l = 0; idx_l <= _LMaxDegree; idx_l++ ) {
+        for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
+            fieldNames[GlobalIndex( idx_l, idx_k )] = std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
+        }
+    }
+    std::vector<std::vector<std::vector<double>>> results{ _outputFields };
     ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
 }