added new solver output to sn and csdsn

code/include/common/globalconstants.h

@@ -76,9 +76,9 @@ enum OPTIMIZER_NAME { NEWTON, ML };
 inline std::map<std::string, OPTIMIZER_NAME> Optimizer_Map{ { "NEWTON", NEWTON }, { "ML", ML } };
 
 // Volume output
-enum VOLUME_OUTPUT { ANALYTIC, MINIMAL, MOMENTS, DUAL_MOMENTS };
+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 } };
+    { "ANALYTIC", ANALYTIC }, { "MINIMAL", MINIMAL }, { "MOMENTS", MOMENTS }, { "DUAL_MOMENTS", DUAL_MOMENTS }, { "DOSE", DOSE } };
 
 #endif    // GLOBAL_CONSTANTS_H

code/include/solvers/csdsnsolver.h

@@ -25,12 +25,11 @@ class CSDSNSolver : public SNSolver
     /**
      * @brief Solve functions runs main time loop
      */
-    virtual void Solve();
-    /**
-     * @brief Output solution to VTK file
-     */
-    virtual void Save() const;
-    virtual void Save( int currEnergy ) const;
+    void Solve() override;
+
+  private:
+    void PrepareOutputFields() override;
+    double WriteOutputFields( unsigned idx_pseudoTime ) override;
 };
 
 #endif    // CSDSNSOLVER_H

code/include/solvers/mnsolver.h

@@ -19,9 +19,7 @@ class MNSolver : public Solver
     /*! @brief MNSolver destructor */
     ~MNSolver();
 
-    void Solve() override;                      /*! @brief Solve functions runs main time loop */
-    void Save() const override;                 /*! @brief Save Output solution to VTK file */
-    void Save( int currEnergy ) const override; /*! @brief Save Output solution at given energy (pseudo time) to VTK file */
+    void Solve() override; /*! @brief Solve functions runs main time loop */
 
   private:
     unsigned _nTotalEntries;    /*! @brief: Total number of equations in the system */
@@ -46,15 +44,10 @@ class MNSolver : public Solver
                                            Layout: _nCells x _nTotalEntries*/
     OptimizerBase* _optimizer; /*! @brief: Class to solve minimal entropy problem */
 
-    // 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::string>> _outputFieldNames;     /*! @brief: Names of the outputFields: dimensions (GroupID,FieldID) */
-
     // ---- Member functions ---
 
     /*! @brief Initializes the output groups and fields of this solver and names the fields */
-    void PrepareOutputFields();
+    void PrepareOutputFields() override;
 
     /*! @brief Function that writes NN Training Data in a .csv file */
     void WriteNNTrainingData( unsigned idx_pseudoTime );
@@ -62,7 +55,7 @@ class MNSolver : public Solver
     /*! @brief Function that prepares VTK export and csv export of the current solver iteration
         @returns: Mass of current iteration
     */
-    double WriteOutputFields( unsigned idx_pseudoTime );
+    double WriteOutputFields( unsigned idx_pseudoTime ) override;
 
     /*! @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

@@ -11,9 +11,7 @@ class PNSolver : public Solver
      */
     PNSolver( Config* settings );
 
-    virtual void Solve() override;              /*! @brief Solve functions runs main time loop. (Run Solver) */
-    void Save() const override;                 /*! @brief Save Output solution to VTK file */
-    void Save( int currEnergy ) const override; /*! @brief Save Output solution at given energy (pseudo time) to VTK file */
+    virtual void Solve() override; /*! @brief Solve functions runs main time loop. (Run Solver) */
 
   protected:
     unsigned _nTotalEntries; /*! @brief: total number of equations in the system */
@@ -42,21 +40,16 @@ class PNSolver : public Solver
     Vector _scatterMatDiag; /*! @brief: diagonal of the scattering matrix (its a diagonal matrix by construction). Contains eigenvalues of the
                                scattering kernel.  */
 
-    // 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::string>> _outputFieldNames;     /*! @brief: Names of the outputFields: dimensions (GroupID,FieldID) */
-
     // ---- Member functions ----
 
     // IO
     /*! @brief Initializes the output groups and fields of this solver and names the fields */
-    void PrepareOutputFields();
+    void PrepareOutputFields() override;
 
     /*! @brief Function that prepares VTK export and csv export of the current solver iteration
         @returns: Mass of current iteration
     */
-    double WriteOutputFields( unsigned idx_pseudoTime );
+    double WriteOutputFields( unsigned idx_pseudoTime ) override;
 
     // Solver
     /*! @brief: parameter functions for setting up system matrix

code/include/solvers/snsolver.h

@@ -23,11 +23,10 @@ class SNSolver : public Solver
      * @brief Solve functions runs main time loop
      */
     void Solve() override;
-    /**
-     * @brief Output solution to VTK file
-     */
-    void Save() const override;
-    void Save( int currEnergy ) const override;
+
+  private:
+    void PrepareOutputFields() override;
+    double WriteOutputFields( unsigned idx_pseudoTime ) override;
 };
 
 #endif    // SNSOLVER_H

code/include/solvers/solverbase.h

@@ -52,8 +52,24 @@ class Solver
     VectorVector _sol;                 /*! @brief solution of the PDE, e.g. angular flux or moments */
     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::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
 
+    // ---- Member functions ----
+
+    // IO
+    /*! @brief Initializes the output groups and fields of this solver and names the fields */
+    virtual void PrepareOutputFields() = 0;
+
+    /*! @brief Function that prepares VTK export and csv export of the current solver iteration
+        @returns: Mass of current iteration
+    */
+    virtual double WriteOutputFields( unsigned idx_pseudoTime ) = 0;
+
     /**
      * @brief ComputeTimeStep calculates the maximal stable time step
      * @param cfl is cfl number
@@ -81,12 +97,8 @@ class Solver
      */
     virtual void Solve() = 0;
 
-    /**
-     * @brief Output solution to VTK file
-     */
-    virtual void Save() const = 0;
-
-    virtual void Save( int currEnergy ) const = 0;
+    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 */
 };
 
 #endif    // SOLVER_H

code/input/exampleSN.cfg

@@ -39,3 +39,8 @@ BC_DIRICHLET = ( void )
 QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
 % Quadrature Order
 QUAD_ORDER = 8
+%
+% ----- Output ---- 
+%
+VOLUME_OUTPUT = (ANALYTIC,  MINIMAL)
+OUTPUT_FREQUENCY = 1

code/src/solvers/csdsnsolver.cpp

@@ -4,6 +4,7 @@
 #include "fluxes/numericalflux.h"
 #include "kernels/scatteringkernelbase.h"
 #include "problems/problembase.h"
+#include "toolboxes/errormessages.h"
 
 // externals
 #include "spdlog/spdlog.h"
@@ -38,6 +39,9 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
 
     // Get patient density
     _density = Vector( _nCells, 1.0 );
+
+    // Solver output
+    PrepareOutputFields();
 }
 
 void CSDSNSolver::Solve() {
@@ -138,7 +142,11 @@ void CSDSNSolver::Solve() {
             _solverOutput[j] = fluxNew[j];
         }
 
-        // Save( n );
+        // --- VTK and CSV Output ---
+        WriteOutputFields( n );
+        Save( n );
+
+        // --- Screen Output ---
         dFlux   = blaze::l2Norm( fluxNew - fluxOld );
         fluxOld = fluxNew;
         if( rank == 0 ) log->info( "{:03.8f}  {:01.5e}  {:01.5e}", _energies[n], _dE / densityMin, dFlux );
@@ -146,18 +154,70 @@ void CSDSNSolver::Solve() {
     }
 }
 
-void CSDSNSolver::Save() const {
-    std::vector<std::string> fieldNames{ "dose" };
-    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
-    std::vector<std::vector<double>> scalarField( 1, _dose );
-    std::vector<std::vector<std::vector<double>>> results{ scalarField };
-    ExportVTK( _settings->GetOutputFile(), results, fieldNamesWrapper, _mesh );
+void CSDSNSolver::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[idx_group][0].resize( _nCells );
+                _outputFieldNames[idx_group][0] = "radiation flux density";
+                break;
+
+            case DOSE:
+                // one entry per cell
+                _outputFields[idx_group].resize( 1 );
+                _outputFieldNames[idx_group].resize( 1 );
+                _outputFields[idx_group][0].resize( _nCells );
+                _outputFieldNames[idx_group][0] = std::string( "raditation dose" );
+                break;
+
+            default: ErrorMessages::Error( "Volume Output Group not defined for CSD SN Solver!", CURRENT_FUNCTION ); break;
+        }
+    }
 }
 
-void CSDSNSolver::Save( int currEnergy ) const {
-    std::vector<std::string> fieldNames{ "dose" };
-    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
-    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, fieldNamesWrapper, _mesh );
+double CSDSNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
+    double mass      = 0.0;
+    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
+
+    // Compute total "mass" of the system ==> to check conservation properties
+    std::vector<double> flux( _nCells, 0.0 );
+    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+        flux[idx_cell] = dot( _sol[idx_cell], _weights );
+        mass += flux[idx_cell] * _areas[idx_cell];
+    }
+
+    if( ( _settings->GetOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetOutputFrequency() == 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];
+                    }
+                    break;
+
+                case DOSE:
+                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+                        _outputFields[idx_group][0][idx_cell] = _dose[idx_cell];    // Remove DOSE
+                    }
+                    break;
+
+                default: ErrorMessages::Error( "Volume Output Group not defined for CSD SN Solver!", CURRENT_FUNCTION ); break;
+            }
+        }
+    }
+    return mass;
 }

code/src/solvers/mnsolver.cpp

@@ -309,14 +309,6 @@ double MNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
     return mass;
 }
 
-void MNSolver::Save() const { ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); }
-
-void MNSolver::Save( int currEnergy ) const {
-    if( _settings->GetOutputFrequency() != 0 && currEnergy % (unsigned)_settings->GetOutputFrequency() == 0 ) {
-        ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), _outputFields, _outputFieldNames, _mesh );
-    }
-}
-
 void MNSolver::WriteNNTrainingData( unsigned idx_pseudoTime ) {
     std::string filename = "trainNN.csv";
     std::ofstream myfile;

code/src/solvers/pnsolver.cpp

@@ -406,14 +406,6 @@ double PNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
     return mass;
 }
 
-void PNSolver::Save() const { ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); }
-
-void PNSolver::Save( int currEnergy ) const {
-    if( _settings->GetOutputFrequency() != 0 && currEnergy % (unsigned)_settings->GetOutputFrequency() == 0 ) {
-        ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), _outputFields, _outputFieldNames, _mesh );
-    }
-}
-
 void PNSolver::CleanFluxMatrices() {
     for( unsigned idx_row = 0; idx_row < _nTotalEntries; idx_row++ ) {
         for( unsigned idx_col = 0; idx_col < _nTotalEntries; idx_col++ ) {

code/src/solvers/snsolver.cpp

@@ -4,6 +4,7 @@
 #include "common/mesh.h"
 #include "fluxes/numericalflux.h"
 #include "kernels/scatteringkernelbase.h"
+#include "problems/problembase.h"
 #include "quadratures/quadraturebase.h"
 
 // externals
@@ -18,6 +19,9 @@ SNSolver::SNSolver( Config* settings ) : Solver( settings ) {
     ScatteringKernel* k = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), _quadrature );
     _scatteringKernel   = k->GetScatteringKernel();
     delete k;
+
+    // Solver output
+    PrepareOutputFields();
 }
 
 void SNSolver::Solve() {
@@ -32,7 +36,7 @@ void SNSolver::Solve() {
     // left and right angular flux of interface, used in numerical flux evaluation
     double psiL;
     double psiR;
-
+    double mass;
     // derivatives of angular flux in x and y directions
     VectorVector psiDx( _nCells, Vector( _nq, 0.0 ) );
     VectorVector psiDy( _nCells, Vector( _nq, 0.0 ) );
@@ -70,16 +74,7 @@ void SNSolver::Solve() {
     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
     if( rank == 0 ) log->info( "{:10}   {:10}", "t", "dFlux", "mass" );
 
-    // Remove
-    double mass1 = 0;
-    for( unsigned i = 0; i < _nCells; ++i ) {
-        _solverOutput[i] = dot( _sol[i], _weights );
-    }
-    // Compute total "mass" of the system ==> to check conservation properties
-    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
-        mass1 += _solverOutput[idx_cell] * _areas[idx_cell];    // Should probably go to postprocessing
-    }
-    if( rank == 0 ) log->info( "{:01.5e}   {:01.5e} {:01.5e}", 0.0, dFlux, mass1 );
+    // if( rank == 0 ) log->info( "{:01.5e}   {:01.5e} {:01.5e}", 0.0, dFlux, mass1 );
 
     // loop over energies (pseudo-time)
     for( unsigned n = 0; n < _nEnergies; ++n ) {
@@ -152,42 +147,90 @@ void SNSolver::Solve() {
         }
         _sol = psiNew;
         for( unsigned i = 0; i < _nCells; ++i ) {
-            fluxNew[i]       = dot( _sol[i], _weights );
-            _solverOutput[i] = fluxNew[i];
+            fluxNew[i] = dot( _sol[i], _weights );
         }
 
-        double mass = 0;
-
-        // Compute total "mass" of the system ==> to check conservation properties
-        for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
-            mass += _solverOutput[idx_cell] * _areas[idx_cell];    // Should probably go to postprocessing
-        }
+        // --- VTK and CSV Output ---
+        mass = WriteOutputFields( n );
+        Save( n );
 
-        // Save( n );
+        // --- Screen Output ---
         dFlux   = blaze::l2Norm( fluxNew - fluxOld );
         fluxOld = fluxNew;
         if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[n], dFlux, mass );
     }
 }
 
-void SNSolver::Save() const {
-    std::vector<std::string> fieldNames{ "flux" };
-    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
+void SNSolver::PrepareOutputFields() {
+    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
 
-    std::vector<double> flux( _nCells, 0.0 );
-    for( unsigned i = 0; i < _nCells; ++i ) {
-        flux[i] = dot( _sol[i], _weights );
+    _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[idx_group][0].resize( _nCells );
+                _outputFieldNames[idx_group][0] = "radiation flux density";
+                break;
+
+            case ANALYTIC:
+                // one entry per cell
+                _outputFields[idx_group].resize( 1 );
+                _outputFieldNames[idx_group].resize( 1 );
+                _outputFields[idx_group][0].resize( _nCells );
+                _outputFieldNames[idx_group][0] = std::string( "analytic radiation flux density" );
+                break;
+
+            default: ErrorMessages::Error( "Volume Output Group not defined for SN Solver!", CURRENT_FUNCTION ); break;
+        }
     }
-    std::vector<std::vector<double>> scalarField( 1, flux );
-    std::vector<std::vector<std::vector<double>>> results{ scalarField };
-    ExportVTK( _settings->GetOutputFile(), results, fieldNamesWrapper, _mesh );
 }
 
-void SNSolver::Save( int currEnergy ) const {
-    std::vector<std::string> fieldNames{ "flux" };
-    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
+double SNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
+    double mass      = 0.0;
+    unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
+
+    // Compute total "mass" of the system ==> to check conservation properties
+    std::vector<double> flux( _nCells, 0.0 );
+    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+        flux[idx_cell] = dot( _sol[idx_cell], _weights );
+        mass += flux[idx_cell] * _areas[idx_cell];
+    }
+
+    if( ( _settings->GetOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetOutputFrequency() == 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];
+                    }
+                    break;
+
+                case ANALYTIC:
+                    // Compute total "mass" of the system ==> to check conservation properties
+                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+
+                        double time = idx_pseudoTime * _dE;
+
+                        _outputFields[idx_group][0][idx_cell] = _problem->GetAnalyticalSolution(
+                            _mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, _sigmaS[idx_pseudoTime][idx_cell] );
+                    }
+                    break;
 
-    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, fieldNamesWrapper, _mesh );
+                default: ErrorMessages::Error( "Volume Output Group not defined for SN Solver!", CURRENT_FUNCTION ); break;
+            }
+        }
+    }
+    return mass;
 }

code/src/solvers/solverbase.cpp

@@ -78,17 +78,10 @@ Solver* Solver::Create( Config* settings ) {
     }
 }
 
-void Solver::Save() const {
-    std::vector<std::string> fieldNames{ "flux" };
-    std::vector<std::vector<std::string>> fieldNamesWrapper{ fieldNames };
+void Solver::Save() const { ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); }
 
-    std::vector<double> flux;
-    flux.resize( _nCells );
-
-    for( unsigned i = 0; i < _nCells; ++i ) {
-        flux[i] = _sol[i][0];
+void Solver::Save( int currEnergy ) const {
+    if( _settings->GetOutputFrequency() != 0 && currEnergy % (unsigned)_settings->GetOutputFrequency() == 0 ) {
+        ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), _outputFields, _outputFieldNames, _mesh );
     }
-    std::vector<std::vector<double>> scalarField( 1, flux );
-    std::vector<std::vector<std::vector<double>>> results{ scalarField };
-    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( _nEnergies ), results, fieldNamesWrapper, _mesh );
 }