linesource analytical solution avail for no scattering. Option is called...

linesource analytical solution avail for no scattering. Option is called ANALYTIC. linesource analytical with scattering started. Analytic output implemented for MN solver

code/include/common/globalconstants.h

@@ -76,8 +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 { MINIMAL, MOMENTS, DUAL_MOMENTS };
+enum VOLUME_OUTPUT { ANALYTIC, MINIMAL, MOMENTS, DUAL_MOMENTS };
 
-inline std::map<std::string, VOLUME_OUTPUT> VolOutput_Map{ { "MINIMAL", MINIMAL }, { "MOMENTS", MOMENTS }, { "DUAL_MOMENTS", DUAL_MOMENTS } };
+inline std::map<std::string, VOLUME_OUTPUT> VolOutput_Map{
+    { "ANALYTIC", ANALYTIC }, { "MINIMAL", MINIMAL }, { "MOMENTS", MOMENTS }, { "DUAL_MOMENTS", DUAL_MOMENTS } };
 
 #endif    // GLOBAL_CONSTANTS_H

code/include/problems/linesource.h

@@ -3,14 +3,40 @@
 
 #include "problembase.h"
 
-class LineSource_SN : public ProblemBase
+class LineSource : public ProblemBase
+{
+  private:
+    LineSource() = delete;
+
+  public:
+    LineSource( Config* settings, Mesh* mesh );
+
+    ~LineSource();
+
+    /*!< @brief: Exact analytical solution for the Line Source Test Case at
+         @param: x: x coordinate of exact solution
+                 y: y coordinate of exact solution
+                 t: time of the exact solution
+                 sigma_s: scattering cross section of the exact solution
+         @return: exact solution at x,y,t,scatteringXS
+    */
+    double GetAnalyticalSolution( double x, double y, double t, double sigma_s ) override;
+
+  private:
+    double ComputeHelperIntegral( double R, double t );
+    double HelperRho_ptc( double R, double t );
+    double HelperRho_ptc1( double R, double t );
+    double HelperRho_ptc2( double R, double t );
+};
+
+class LineSource_SN : public LineSource
 {
   private:
     LineSource_SN() = delete;
 
   public:
     LineSource_SN( Config* settings, Mesh* mesh );
-    virtual ~LineSource_SN();
+    ~LineSource_SN();
 
     virtual VectorVector GetScatteringXS( const Vector& energies );
     virtual VectorVector GetTotalXS( const Vector& energies );
@@ -41,7 +67,7 @@ class LineSource_SN_Pseudo1D_Physics : public LineSource_SN_Pseudo1D
     Vector GetTotalXSE( const Vector& energies ) override;
 };
 
-class LineSource_PN : public ProblemBase
+class LineSource_PN : public LineSource
 {
   private:
     LineSource_PN() = delete;
@@ -58,7 +84,7 @@ class LineSource_PN : public ProblemBase
 
   public:
     LineSource_PN( Config* settings, Mesh* mesh );
-    virtual ~LineSource_PN();
+    ~LineSource_PN();
 
     virtual VectorVector GetScatteringXS( const Vector& energies ) override;
     virtual VectorVector GetTotalXS( const Vector& energies ) override;

code/include/problems/problembase.h

@@ -79,6 +79,15 @@ class ProblemBase
      */
     virtual VectorVector SetupIC() = 0;
 
+    /*!< @brief: Exact analytical solution for the Line Source Test Case at
+         @param: x: x coordinate of exact solution
+                 y: y coordinate of exact solution
+                 t: time of the exact solution
+                 scatteringXS: scattering cross section of the exact solution
+         @return: exact solution at x,y,t,scatteringXS
+    */ // Default is set to 0. ~> if no analytical solution is available.
+    double virtual GetAnalyticalSolution( double x, double y, double t, double scatteringXS ) { return 0.0; }
+
     /**
      * @brief Physics constructor
      * @param settings stores all needed user information

code/include/solvers/mnsolver.h

@@ -65,7 +65,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();
+    double WriteOutputFields( unsigned idx_pseudoTime );
 
     /*! @brief : computes the global index of the moment corresponding to basis function (l,k)
      *  @param : degree l, it must hold: 0 <= l <=_nq

code/input/exampleMN.cfg

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%  Linesource Benchmarking File PN     %
+%  Linesource Benchmarking File MN     %
 %  Author <Steffen Schotthöfer>		   %
 %  Date   01.07.2020                   %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -10,12 +10,12 @@
 % Output directory
 OUTPUT_DIR = ../result
 % Output file
-OUTPUT_FILE = exampleMN_feature_output
+OUTPUT_FILE = exampleMN_dualM_Analytic
 % 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
 %
@@ -60,4 +60,4 @@ QUAD_ORDER = 12
 %
 % ----- Output ---- 
 %
-VOLUME_OUTPUT = (MINIMAL, MOMENTS)
+VOLUME_OUTPUT = (ANALYTIC, MOMENTS, DUAL_MOMENTS)

code/src/common/config.cpp

@@ -432,21 +432,39 @@ void Config::SetPostprocessing() {
     for( unsigned short idx_volOutput = 0; idx_volOutput < _nVolumeOutput; idx_volOutput++ ) {
         switch( _solverName ) {
             case SN_SOLVER:
-                if( _volumeOutput[idx_volOutput] != MINIMAL ) {
-                    ErrorMessages::Error( "SN_SOLVER only supports volume output MINIMAL.\nPlease check your .cfg file.", CURRENT_FUNCTION );
+                if( _volumeOutput[idx_volOutput] != MINIMAL && _volumeOutput[idx_volOutput] != ANALYTIC ) {
+                    ErrorMessages::Error( "SN_SOLVER only supports volume output MINIMAL and ANALYTIC.\nPlease check your .cfg file.",
+                                          CURRENT_FUNCTION );
+                }
+                if( _volumeOutput[idx_volOutput] == ANALYTIC && _problemName != PROBLEM_LineSource ) {
+                    ErrorMessages::Error(
+                        "Analytical solution (VOLUME_OUTPUT=ANALYTIC) is only available for the PROBLEM=LINESOURCE.\nPlease check your .cfg file.",
+                        CURRENT_FUNCTION );
                 }
                 break;
             case MN_SOLVER:
-                if( _volumeOutput[idx_volOutput] != MINIMAL || _volumeOutput[idx_volOutput] != MOMENTS ) {
-                    ErrorMessages::Error( "MN_SOLVER only supports volume output MINIMAL and MOMENTS.\nPlease check your .cfg file.",
-                                          CURRENT_FUNCTION );
+                if( _volumeOutput[idx_volOutput] != MINIMAL && _volumeOutput[idx_volOutput] != MOMENTS &&
+                    _volumeOutput[idx_volOutput] != DUAL_MOMENTS && _volumeOutput[idx_volOutput] != ANALYTIC ) {
+                    ErrorMessages::Error(
+                        "MN_SOLVER only supports volume output ANALYTIC, MINIMAL, MOMENTS and DUAL_MOMENTS.\nPlease check your .cfg file.",
+                        CURRENT_FUNCTION );
+                }
+                if( _volumeOutput[idx_volOutput] == ANALYTIC && _problemName != PROBLEM_LineSource ) {
+                    ErrorMessages::Error(
+                        "Analytical solution (VOLUME_OUTPUT=ANALYTIC) is only available for the PROBLEM=LINESOURCE.\nPlease check your .cfg file.",
+                        CURRENT_FUNCTION );
                 }
                 break;
             case PN_SOLVER:
-                if( _volumeOutput[idx_volOutput] != MINIMAL || _volumeOutput[idx_volOutput] != MOMENTS ) {
-                    ErrorMessages::Error( "PN_SOLVER only supports volume output MINIMAL and MOMENTS.\nPlease check your .cfg file.",
+                if( _volumeOutput[idx_volOutput] != MINIMAL && _volumeOutput[idx_volOutput] != MOMENTS && _volumeOutput[idx_volOutput] != ANALYTIC ) {
+                    ErrorMessages::Error( "PN_SOLVER only supports volume output ANALYTIC, MINIMAL and MOMENTS.\nPlease check your .cfg file.",
                                           CURRENT_FUNCTION );
                 }
+                if( _volumeOutput[idx_volOutput] == ANALYTIC && _problemName != PROBLEM_LineSource ) {
+                    ErrorMessages::Error(
+                        "Analytical solution (VOLUME_OUTPUT=ANALYTIC) is only available for the PROBLEM=LINESOURCE.\nPlease check your .cfg file.",
+                        CURRENT_FUNCTION );
+                }
                 break;
             case CSD_SN_SOLVER:
                 if( _volumeOutput[idx_volOutput] != MINIMAL ) {

code/src/problems/linesource.cpp

@@ -3,9 +3,67 @@
 #include "common/mesh.h"
 #include "physics.h"
 
+// ---- Linesource ----
+
+LineSource::LineSource( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }
+LineSource::~LineSource() {}
+
+double LineSource::GetAnalyticalSolution( double x, double y, double t, double sigma_s ) {
+
+    double solution = 0.0;
+    double R        = sqrt( x * x + y * y );
+
+    if( sigma_s == 0.0 ) {
+        if( ( t - R ) > 0 ) {
+            solution = 1 / ( 2 * M_PI * t * sqrt( t * t - R * R ) );
+        }
+    }
+    else {
+        double gamma = R / t;
+
+        if( ( 1 - gamma ) > 0 ) {
+            solution = exp( -t ) / ( 2 * M_PI * t * t * sqrt( 1 - gamma * gamma ) );
+
+            double integral = ComputeHelperIntegral( R, t );
+            solution += 2 * t * integral;
+        }
+    }
+    return solution;
+}
+
+double LineSource::ComputeHelperIntegral( double R, double t ) {
+    double result = 0;
+
+    return result;
+}
+
+double LineSource::HelperRho_ptc( double R, double t ) {
+    double result = HelperRho_ptc1( R, t ) + HelperRho_ptc2( R, t );
+
+    return result;
+}
+
+double LineSource::HelperRho_ptc1( double R, double t ) {
+    double gamma  = R / t;
+    double result = exp( -t ) / ( 4.0 * M_PI * R * t ) * log( ( 1 + gamma ) / ( 1 - gamma ) );
+
+    return result;
+}
+
+double LineSource::HelperRho_ptc2( double R, double t ) {
+    double gamma = R / t;
+
+    double result = 0;
+    if( 1 - gamma > 0 ) {
+        result = exp( -t ) / ( 32 * M_PI * M_PI * R ) * ( 1 - gamma * gamma );
+    }
+
+    return result;
+}
+
 // ---- LineSource_SN ----
 
-LineSource_SN::LineSource_SN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }
+LineSource_SN::LineSource_SN( Config* settings, Mesh* mesh ) : LineSource( settings, mesh ) {}
 
 LineSource_SN::~LineSource_SN() {}
 
@@ -64,7 +122,7 @@ int LineSource_PN::GlobalIndex( int l, int k ) const {
     return numIndicesPrevLevel + prevIndicesThisLevel;
 }
 
-LineSource_PN::LineSource_PN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }
+LineSource_PN::LineSource_PN( Config* settings, Mesh* mesh ) : LineSource( settings, mesh ) {}
 
 LineSource_PN::~LineSource_PN() {}
 

code/src/solvers/csdsnsolver.cpp

@@ -10,6 +10,7 @@
 #include <mpi.h>
 
 CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
+
     _dose = std::vector<double>( _settings->GetNCells(), 0.0 );
 
     // Set angle and energies

code/src/solvers/mnsolver.cpp

 #include "solvers/mnsolver.h"
 #include "common/config.h"
 #include "common/io.h"
+#include "common/mesh.h"
 #include "entropies/entropybase.h"
 #include "fluxes/numericalflux.h"
 #include "optimizers/optimizerbase.h"
+#include "problems/problembase.h"
 #include "quadratures/quadraturebase.h"
 #include "solvers/sphericalharmonics.h"
 #include "toolboxes/errormessages.h"
@@ -35,7 +37,7 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
     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;
+            _sigmaS[n][j] = 0;
         }
     }
 
@@ -180,7 +182,7 @@ void MNSolver::Solve() {
         _sol = psiNew;
 
         // --- VTK and CSV Output ---
-        mass = WriteOutputFields();
+        mass = WriteOutputFields( idx_energy );
         Save( idx_energy );
         WriteNNTrainingData( idx_energy );
 
@@ -224,12 +226,36 @@ void MNSolver::PrepareOutputFields() {
                     }
                 }
                 break;
+
+            case DUAL_MOMENTS:
+                // As many entries as there are moments in the system
+                _outputFields[idx_group].resize( _nTotalEntries );
+                _outputFieldNames[idx_group].resize( _nTotalEntries );
+
+                for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
+                    for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
+                        _outputFields[idx_group][GlobalIndex( idx_l, idx_k )].resize( _nCells );
+
+                        _outputFieldNames[idx_group][GlobalIndex( idx_l, idx_k )] =
+                            std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
+                    }
+                }
+                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 MN Solver!", CURRENT_FUNCTION ); break;
         }
     }
 }
 
-double MNSolver::WriteOutputFields() {
+double MNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
     double mass      = 0.0;
     unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
 
@@ -252,6 +278,25 @@ double MNSolver::WriteOutputFields() {
                     }
                 }
                 break;
+            case DUAL_MOMENTS:
+                for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
+                    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+                        _outputFields[idx_group][idx_sys][idx_cell] = _alpha[idx_cell][idx_sys];
+                    }
+                }
+                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;
+                    double sigma = 0;
+
+                    _outputFields[idx_group][0][idx_cell] = _problem->GetAnalyticalSolution(
+                        _mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, sigma );
+                }
+                break;
+
             default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
         }
     }