added checkerboard; scaling might need work

code/include/problems/checkerboard.h

@@ -3,20 +3,42 @@
 
 #include "problembase.h"
 
-class Checkerboard : public ProblemBase
+class Checkerboard_SN : public ProblemBase
 {
   private:
     Vector _scatteringXS;
     Vector _totalXS;
 
-    Checkerboard() = delete;
+    Checkerboard_SN() = delete;
 
     bool isAbsorption( const Vector& pos ) const;
     bool isSource( const Vector& pos ) const;
 
   public:
-    Checkerboard( Config* settings, Mesh* mesh );
-    virtual ~Checkerboard();
+    Checkerboard_SN( Config* settings, Mesh* mesh );
+    virtual ~Checkerboard_SN();
+
+    virtual VectorVector GetScatteringXS( const Vector& energies );
+    virtual VectorVector GetTotalXS( const Vector& energies );
+    virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
+    virtual VectorVector SetupIC();
+};
+
+class Checkerboard_PN : public ProblemBase
+{
+  private:
+    Vector _scatteringXS;
+    Vector _totalXS;
+
+    Checkerboard_PN() = delete;
+
+    bool isAbsorption( const Vector& pos ) const;
+    bool isSource( const Vector& pos ) const;
+    int GlobalIndex( int l, int k ) const;
+
+  public:
+    Checkerboard_PN( Config* settings, Mesh* mesh );
+    virtual ~Checkerboard_PN();
 
     virtual VectorVector GetScatteringXS( const Vector& energies );
     virtual VectorVector GetTotalXS( const Vector& energies );

code/include/solvers/pnsolver.h

@@ -19,7 +19,6 @@ class PNSolver : public Solver
     unsigned _nTotalEntries; /*! @brief: total number of equations in the system */
     unsigned _LMaxDegree;    /*! @brief: maximal degree of the spherical harmonics basis*/
 
-    VectorVector _sigmaA; /*!  @brief: Absorption coefficient for all energies*/
     // System Matrix for x, y and z flux
     //    ==> not needed after computation of A+ and A- ==> maybe safe only temporarly and remove as member?
     SymMatrix _Ax; /*! @brief:  Flux Jacbioan in x direction */

code/src/problems/checkerboard.cpp

@@ -2,7 +2,7 @@
 #include "common/config.h"
 #include "common/mesh.h"
 
-Checkerboard::Checkerboard( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
+Checkerboard_SN::Checkerboard_SN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
     _physics      = nullptr;
     _scatteringXS = Vector( _mesh->GetNumCells(), 1.0 );
     _totalXS      = Vector( _mesh->GetNumCells(), 1.0 );
@@ -15,13 +15,13 @@ Checkerboard::Checkerboard( Config* settings, Mesh* mesh ) : ProblemBase( settin
     }
 }
 
-Checkerboard::~Checkerboard() {}
+Checkerboard_SN::~Checkerboard_SN() {}
 
-VectorVector Checkerboard::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
+VectorVector Checkerboard_SN::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
 
-VectorVector Checkerboard::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
+VectorVector Checkerboard_SN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
 
-std::vector<VectorVector> Checkerboard::GetExternalSource( const Vector& energies ) {
+std::vector<VectorVector> Checkerboard_SN::GetExternalSource( const Vector& energies ) {
     VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
     auto cellMids = _mesh->GetCellMidPoints();
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
@@ -30,12 +30,12 @@ std::vector<VectorVector> Checkerboard::GetExternalSource( const Vector& energie
     return std::vector<VectorVector>( 1u, Q );
 }
 
-VectorVector Checkerboard::SetupIC() {
+VectorVector Checkerboard_SN::SetupIC() {
     VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
     return psi;
 }
 
-bool Checkerboard::isAbsorption( const Vector& pos ) const {
+bool Checkerboard_SN::isAbsorption( const Vector& pos ) const {
     std::vector<double> lbounds{ 1, 2, 3, 4, 5 };
     std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
     for( unsigned k = 0; k < lbounds.size(); ++k ) {
@@ -49,9 +49,72 @@ bool Checkerboard::isAbsorption( const Vector& pos ) const {
     return false;
 }
 
-bool Checkerboard::isSource( const Vector& pos ) const {
+bool Checkerboard_SN::isSource( const Vector& pos ) const {
     if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
         return true;
     else
         return false;
 }
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+
+Checkerboard_PN::Checkerboard_PN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
+    _physics      = nullptr;
+    _scatteringXS = Vector( _mesh->GetNumCells(), 1.0 );
+    _totalXS      = Vector( _mesh->GetNumCells(), 1.0 );
+    auto cellMids = _mesh->GetCellMidPoints();
+    for( unsigned j = 0; j < cellMids.size(); ++j ) {
+        if( isAbsorption( cellMids[j] ) ) {
+            _scatteringXS[j] = 0.0;
+            _totalXS[j]      = 10.0;
+        }
+    }
+}
+
+Checkerboard_PN::~Checkerboard_PN() {}
+
+VectorVector Checkerboard_PN::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
+
+VectorVector Checkerboard_PN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
+
+std::vector<VectorVector> Checkerboard_PN::GetExternalSource( const Vector& energies ) {
+    VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
+    auto cellMids = _mesh->GetCellMidPoints();
+    for( unsigned j = 0; j < cellMids.size(); ++j ) {
+        if( isSource( cellMids[j] ) ) Q[j] = std::sqrt( 4 * M_PI );    // isotropic source
+    }
+    return std::vector<VectorVector>( 1u, Q );
+}
+
+VectorVector Checkerboard_PN::SetupIC() {
+    int ntotalEquations = GlobalIndex( _settings->GetMaxMomentDegree(), _settings->GetMaxMomentDegree() ) + 1;
+    VectorVector psi( _mesh->GetNumCells(), Vector( ntotalEquations, 1e-10 ) );
+    return psi;
+}
+
+bool Checkerboard_PN::isAbsorption( const Vector& pos ) const {
+    std::vector<double> lbounds{ 1, 2, 3, 4, 5 };
+    std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
+    for( unsigned k = 0; k < lbounds.size(); ++k ) {
+        for( unsigned l = 0; l < lbounds.size(); ++l ) {
+            if( ( l + k ) % 2 == 1 || ( k == 2 && l == 2 ) || ( k == 2 && l == 4 ) ) continue;
+            if( pos[0] >= lbounds[k] && pos[0] <= ubounds[k] && pos[1] >= lbounds[l] && pos[1] <= ubounds[l] ) {
+                return true;
+            }
+        }
+    }
+    return false;
+}
+
+bool Checkerboard_PN::isSource( const Vector& pos ) const {
+    if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
+        return true;
+    else
+        return false;
+}
+
+int Checkerboard_PN::GlobalIndex( int l, int k ) const {
+    int numIndicesPrevLevel  = l * l;    // number of previous indices untill level l-1
+    int prevIndicesThisLevel = k + l;    // number of previous indices in current level
+    return numIndicesPrevLevel + prevIndicesThisLevel;
+}

code/src/problems/problembase.cpp

@@ -22,7 +22,12 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
             else
                 return new LineSource_SN( settings, mesh );    // default
         }
-        case PROBLEM_Checkerboard: return new Checkerboard( settings, mesh );
+        case PROBLEM_Checkerboard: {
+            if( settings->GetSolverName() == PN_SOLVER || settings->GetSolverName() == MN_SOLVER )
+                return new Checkerboard_PN( settings, mesh );
+            else
+                return new Checkerboard_SN( settings, mesh );    // default
+        }
         case PROBLEM_ElectronRT: return new ElectronRT( settings, mesh );
         case PROBLEM_WaterPhantom: return new WaterPhantom( settings, mesh );
         case PROBLEM_LineSource_Pseudo_1D: return new LineSource_SN_Pseudo1D( settings, mesh );

code/src/solvers/mnsolver.cpp

@@ -142,7 +142,8 @@ void MNSolver::Solve() {
 
     for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
 
-        // Loop over the grid cells
+// Loop over the grid cells
+#pragma omp parallel for
         for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
 
             // ------- Reconstruction Step -------
@@ -165,9 +166,10 @@ void MNSolver::Solve() {
                 psiNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
                                                ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
                                                - _dE * _sol[idx_cell][idx_system] *
-                                                     ( _sigmaA[idx_energy][idx_cell]                                    /* absorbtion influence */
-                                                       + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
+                                                     ( ( _sigmaT[idx_energy][idx_cell] - _sigmaS[idx_energy][idx_cell] ) /* absorbtion influence */
+                                                       + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] );  /* scattering influence */
             }
+            psiNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
         }
         _sol = psiNew;
 
@@ -175,26 +177,20 @@ void MNSolver::Solve() {
         double mass = 0.0;
         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];
+            _solverOutput[idx_cell] = std::sqrt( 4 * PI ) * _sol[idx_cell][0];
             mass += _sol[idx_cell][0] * _areas[idx_cell];
         }
 
         dFlux   = blaze::l2Norm( fluxNew - fluxOld );
         fluxOld = fluxNew;
         if( rank == 0 ) log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
-        Save( idx_energy );
+        // Save( idx_energy );
     }
 }
 
 void MNSolver::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<double>> scalarField( 1, _solverOutput );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };
     ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
 }

code/src/solvers/pnsolver.cpp

@@ -13,16 +13,6 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
     _LMaxDegree    = settings->GetMaxMomentDegree();
     _nTotalEntries = GlobalIndex( int( _LMaxDegree ), int( _LMaxDegree ) ) + 1;
 
-    // transform sigmaT and sigmaS in sigmaA.
-    _sigmaA = VectorVector( _nEnergies, Vector( _nCells, 0 ) );    // Get rid of this extra vektor!
-
-    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;
-        }
-    }
-
     // Initialize System Matrices
     _Ax = SymMatrix( _nTotalEntries );
     _Ay = SymMatrix( _nTotalEntries );
@@ -106,6 +96,7 @@ void PNSolver::Solve() {
 
     for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
         // Loop over all spatial cells
+#pragma omp parallel for
         for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
             if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;    // Dirichlet cells stay at IC, farfield assumption
 
@@ -144,9 +135,12 @@ void PNSolver::Solve() {
                     psiNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
                                                    ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
                                                    - _dE * _sol[idx_cell][idx_system] *
-                                                         ( _sigmaA[idx_energy][idx_cell]                                    /* absorbtion influence */
-                                                           + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
+                                                         ( ( _sigmaT[idx_energy][idx_cell] - _sigmaS[idx_energy][idx_cell] ) +
+                                                           _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
                 }
+
+                // add isotropic source to 0-th moment
+                psiNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
             }
         }
         _sol = psiNew;
@@ -154,7 +148,7 @@ void PNSolver::Solve() {
         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];
+            _solverOutput[i] = std::sqrt( 4 * PI ) * _sol[i][0];
             mass += _sol[i][0] * _areas[i];
         }
 
@@ -163,7 +157,7 @@ void PNSolver::Solve() {
         if( rank == 0 ) {
             log->info( "{:03.8f}   {:01.5e} {:01.5e}", _energies[idx_energy], dFlux, mass );
         }
-        Save( idx_energy );
+        // Save( idx_energy );
     }
 }
 
@@ -369,13 +363,7 @@ 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<double>> scalarField( 1, _solverOutput );
     std::vector<std::vector<std::vector<double>>> results{ scalarField };
     ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
 }