added checkerboard test case to mn and pn

Former-commit-id: 54c8a262

code/include/common/globalconstants.h

@@ -65,10 +65,10 @@ inline std::map<std::string, SOLVER_NAME> Solver_Map{
     { "SN_SOLVER", SN_SOLVER }, { "CSD_SN_SOLVER", CSD_SN_SOLVER }, { "PN_SOLVER", PN_SOLVER }, { "MN_SOLVER", MN_SOLVER } };
 
 // Entropy functional
-enum ENTROPY_NAME { QUADRATIC, MAXWELL_BOLZMANN, BOSE_EINSTEIN, FERMI_DIRAC };
+enum ENTROPY_NAME { QUADRATIC, MAXWELL_BOLTZMANN, BOSE_EINSTEIN, FERMI_DIRAC };
 
 inline std::map<std::string, ENTROPY_NAME> Entropy_Map{
-    { "QUADRATIC", QUADRATIC }, { "MAXWELL_BOLZMANN", MAXWELL_BOLZMANN }, { "BOSE_EINSTEIN", BOSE_EINSTEIN }, { "FERMI_DIRAC", FERMI_DIRAC } };
+    { "QUADRATIC", QUADRATIC }, { "MAXWELL_BOLTZMANN", MAXWELL_BOLTZMANN }, { "BOSE_EINSTEIN", BOSE_EINSTEIN }, { "FERMI_DIRAC", FERMI_DIRAC } };
 
 // Optimizer
 enum OPTIMIZER_NAME { NEWTON, ML };

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/src/common/config.cpp

@@ -250,6 +250,7 @@ void Config::SetConfigOptions() {
     // Entropy related options
     /*! @brief Entropy Functional \n DESCRIPTION: Entropy functional used for the MN_Solver \n DEFAULT QUADRTATIC @ingroup Config. */
     AddEnumOption( "ENTROPY_FUNCTIONAL", _entropyName, Entropy_Map, QUADRATIC );
+    AddEnumOption( "ENTROPY_FUNCTIONAL", _entropyName, Entropy_Map, MAXWELL_BOLTZMANN );
     /*! @brief Optimizer Name \n DESCRIPTION:  Optimizer used to determine the minimal Entropy reconstruction \n DEFAULT NEWTON \ingroup Config */
     AddEnumOption( "ENTROPY_OPTIMIZER", _entropyOptimizerName, Optimizer_Map, NEWTON );
 

code/src/entropies/entropybase.cpp

@@ -7,7 +7,7 @@ EntropyBase* EntropyBase::Create( Config* settings ) {
 
     switch( settings->GetEntropyName() ) {
         case QUADRATIC: return new QuadraticEntropy();
-        case MAXWELL_BOLZMANN: return new MaxwellBoltzmannEntropy();
+        case MAXWELL_BOLTZMANN: return new MaxwellBoltzmannEntropy();
         default: return new QuadraticEntropy();
     }
 }

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] = 1.0 / 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

@@ -175,6 +175,8 @@ void MNSolver::Solve() {
                                                      ( _sigmaT[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];
         }
 
         // Update Solution

code/src/solvers/pnsolver.cpp

@@ -117,6 +117,7 @@ void PNSolver::Solve() {
                                                            + _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
                 }
             }
+            psiNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
         }
 
         // Update Solution