added IC for monomial solver for grade up to 1

Former-commit-id: c1ecd6e8

code/input/exampleMN.cfg

@@ -26,22 +26,22 @@ SCATTER_COEFF = 1
 %
 % Choice of basis functions
 SPHERICAL_BASIS = SPHERICAL_HARMONICS
-%
+% Maximal Moment degree
+MAX_MOMENT_SOLVER = 1
 % Solver type
 SOLVER = MN_SOLVER
 % CFL number
 CFL_NUMBER = 0.7
 % Final time for simulation
 TIME_FINAL = 0.3
-% Maximal Moment degree
-MAX_MOMENT_SOLVER = 2
+
+% Reconstruction order (spatial flux)
+RECONS_ORDER = 1
 %
 % ---- Entropy settings ----
 %
 ENTROPY_FUNCTIONAL = MAXWELL_BOLTZMANN
 ENTROPY_OPTIMIZER = NEWTON
-% Reconstruction order
-RECONS_ORDER = 2
 %
 % ----- Newton Solver Specifications ----
 %
@@ -66,8 +66,8 @@ QUAD_ORDER = 8
 %
 % ----- Output ---- 
 %
-VOLUME_OUTPUT = (ANALYTIC,  MINIMAL, MOMENTS, DUAL_MOMENTS)
-VOLUME_OUTPUT_FREQUENCY = 0
+VOLUME_OUTPUT = (MINIMAL, MOMENTS)
+VOLUME_OUTPUT_FREQUENCY = 1
 SCREEN_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)
 SCREEN_OUTPUT_FREQUENCY = 1
 HISTORY_OUTPUT = (ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT)

code/src/common/config.cpp

@@ -233,7 +233,7 @@ void Config::SetConfigOptions() {
     AddEnumOption( "PROBLEM", _problemName, Problem_Map, PROBLEM_ElectronRT );
     /*! @brief Solver \n DESCRIPTION: Solver used for problem \n DEFAULT SN_SOLVER @ingroup Config. */
     AddEnumOption( "SOLVER", _solverName, Solver_Map, SN_SOLVER );
-    /*! @brief RECONS_ORDER \n DESCRIPTION: Reconstruction order for solver \n DEFAULT 1 \ingroup Config.*/
+    /*! @brief RECONS_ORDER \n DESCRIPTION: Reconstruction order for solver (spatial flux) \n DEFAULT 1 \ingroup Config.*/
     AddUnsignedShortOption( "RECONS_ORDER", _reconsOrder, 1 );
     /*! @brief CleanFluxMatrices \n DESCRIPTION:  If true, very low entries (10^-10 or smaller) of the flux matrices will be set to zero,
      * to improve floating point accuracy \n DEFAULT false \ingroup Config */
@@ -448,6 +448,10 @@ void Config::SetPostprocessing() {
                               CURRENT_FUNCTION );
     }
 
+    if( GetSolverName() == MN_SOLVER && GetSphericalBasisName() == SPHERICAL_MONOMIALS && GetMaxMomentDegree() > 1 ) {
+        ErrorMessages::Error( "MN Solver only with monomial basis only stable up to degree 1. This is a TODO.", CURRENT_FUNCTION );
+    }
+
     // --- Output Postprocessing ---
 
     // Volume Output Postprocessing

code/src/problems/linesource.cpp

 #include "problems/linesource.h"
 #include "common/config.h"
 #include "common/mesh.h"
+#include "quadratures/quadraturebase.h"
 #include "toolboxes/physics.h"
 #include "toolboxes/sphericalbase.h"
 #include <complex>
@@ -182,17 +183,45 @@ VectorVector LineSource_PN::SetupIC() {
     // In case of PN, spherical basis is per default SPHERICAL_HARMONICS
     SphericalBase* tempBase  = SphericalBase::Create( _settings );
     unsigned ntotalEquations = tempBase->GetBasisSize();
-    delete tempBase;    // Only temporally needed
 
     VectorVector psi( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) );    // zero could lead to problems?
     VectorVector cellMids = _mesh->GetCellMidPoints();
 
+    Vector uIC( ntotalEquations, 0 );
+
+    if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
+        QuadratureBase* quad          = QuadratureBase::Create( _settings );
+        VectorVector quadPointsSphere = quad->GetPointsSphere();
+        Vector w                      = quad->GetWeights();
+
+        double my, phi;
+        VectorVector moments = VectorVector( quad->GetNq(), Vector( tempBase->GetBasisSize(), 0.0 ) );
+
+        for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
+            my                = quadPointsSphere[idx_quad][0];
+            phi               = quadPointsSphere[idx_quad][1];
+            moments[idx_quad] = tempBase->ComputeSphericalBasis( my, phi );
+        }
+        // Integrate <1*m> to get factors for monomial basis in isotropic scattering
+        for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) {
+            uIC += w[idx_quad] * moments[idx_quad];
+        }
+        delete quad;
+    }
+
     // Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments - exept first - are zero.
     double t = 3.2e-4;    // pseudo time for gaussian smoothing (Approx to dirac impulse)
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
         double x  = cellMids[j][0];
         double y  = cellMids[j][1];    // (x- 0.5) * (x- 0.5)
         psi[j][0] = /*sqrt( 4 * M_PI ) * */ 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) );
+
+        double c = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) );
+
+        if( _settings->GetSphericalBasisName() == SPHERICAL_MONOMIALS ) {
+            psi[j] = c * uIC / uIC[0];    // Remember scaling
+        }
     }
+    delete tempBase;    // Only temporally needed
     return psi;
 }

code/src/quadratures/quadraturebase.cpp

@@ -90,7 +90,7 @@ void QuadratureBase::PrintWeights() {
 void QuadratureBase::PrintPoints() {
     auto log = spdlog::get( "event" );
     for( unsigned i = 0; i < _nq; i++ ) {
-        log->info( "{0}, {1}, {2}", _points[i][0], _points[i][1], _points[i][2], );
+        log->info( "{0}, {1}, {2}", _points[i][0], _points[i][1], _points[i][2] );
     }
 }
 void QuadratureBase::PrintPointsAndWeights() {