checkerboard IC supports monomial basis now. Renamed to checkerboard_moment...

checkerboard IC supports monomial basis now. Renamed to checkerboard_moment for clarity. mn solver works with 2d quadrature now


Former-commit-id: 45f3488d

code/include/problems/checkerboard.h

@@ -3,6 +3,8 @@
 
 #include "problembase.h"
 
+class SphericalBase;
+
 class Checkerboard_SN : public ProblemBase
 {
   private:
@@ -24,13 +26,14 @@ class Checkerboard_SN : public ProblemBase
     virtual VectorVector SetupIC() override;
 };
 
-class Checkerboard_PN : public ProblemBase
+class Checkerboard_Moment : public ProblemBase
 {
   private:
     Vector _scatteringXS; /*!< @brief Vector of scattering crosssections len: numCells  */
     Vector _totalXS;      /*!< @brief Vector of total crosssections. len: numCells*/
+    SphericalBase* _basis; /*!< @brief Class to compute and store current spherical harmonics basis */
 
-    Checkerboard_PN() = delete;
+    Checkerboard_Moment() = delete;
 
     bool isAbsorption( const Vector& pos ) const; /*!< @return True if pos is in absorption region, False otherwise */
     bool isSource( const Vector& pos ) const;     /*!< @return True if pos is in source region, False otherwise */
@@ -46,8 +49,8 @@ class Checkerboard_PN : public ProblemBase
     int GlobalIndex( int l, int k ) const;
 
   public:
-    Checkerboard_PN( Config* settings, Mesh* mesh );
-    virtual ~Checkerboard_PN();
+    Checkerboard_Moment( Config* settings, Mesh* mesh );
+    virtual ~Checkerboard_Moment();
 
     virtual VectorVector GetScatteringXS( const Vector& energies ) override;
     virtual VectorVector GetTotalXS( const Vector& energies ) override;

code/input/checkerboard_MN.cfg

@@ -18,6 +18,7 @@ MESH_FILE = meshes/checkerboard.su2
 % ---- Problem specifications ----
 %
 PROBLEM = CHECKERBOARD
+SPATIAL_DIM = 2
 %
 % ---- Solver specifications ----
 %

code/src/fluxes/upwindflux.cpp

@@ -3,7 +3,7 @@
 UpwindFlux::UpwindFlux() : NumericalFlux() {}
 
 double UpwindFlux::Flux( const Vector& Omega, double psiL, double psiR, const Vector& n ) const {
-    double inner = Omega[0] * n[0] + Omega[1] * n[1];
+    double inner = Omega[0] * n[0] + Omega[1] * n[1]; // Only use x and y axis in 2d case
     if( inner > 0 ) {
         return inner * psiL;
     }

code/src/problems/checkerboard.cpp

 #include "problems/checkerboard.h"
 #include "common/config.h"
 #include "common/mesh.h"
+#include "toolboxes/sphericalbase.h"
 
 // ---- Checkerboard Sn ----
 // Constructor for Ckeckerboard case with Sn
@@ -66,9 +67,12 @@ bool Checkerboard_SN::isSource( const Vector& pos ) const {
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 
-// ---- Checkerboard Pn ----
+// ---- Checkerboard Moments ----
+
 // Constructor for checkerboard case with Pn
-Checkerboard_PN::Checkerboard_PN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
+Checkerboard_Moment::Checkerboard_Moment( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
+
+    _basis = SphericalBase::Create(_settings);
     _physics = nullptr;
 
     // Initialise crosssections = 1 (scattering)
@@ -85,13 +89,15 @@ Checkerboard_PN::Checkerboard_PN( Config* settings, Mesh* mesh ) : ProblemBase(
     }
 }
 
-Checkerboard_PN::~Checkerboard_PN() {}
+Checkerboard_Moment::~Checkerboard_Moment() {
+    delete _basis;
+}
 
-VectorVector Checkerboard_PN::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
+VectorVector Checkerboard_Moment::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
 
-VectorVector Checkerboard_PN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
+VectorVector Checkerboard_Moment::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
 
-std::vector<VectorVector> Checkerboard_PN::GetExternalSource( const Vector& /*energies*/ ) {
+std::vector<VectorVector> Checkerboard_Moment::GetExternalSource( const Vector& /*energies*/ ) {
     VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
     auto cellMids = _mesh->GetCellMidPoints();
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
@@ -100,13 +106,13 @@ std::vector<VectorVector> Checkerboard_PN::GetExternalSource( const Vector& /*en
     return std::vector<VectorVector>( 1u, Q );
 }
 
-VectorVector Checkerboard_PN::SetupIC() {
-    int ntotalEquations = GlobalIndex( _settings->GetMaxMomentDegree(), _settings->GetMaxMomentDegree() ) + 1;
+VectorVector Checkerboard_Moment::SetupIC() {
+    int ntotalEquations = _basis->GetBasisSize();
     VectorVector psi( _mesh->GetNumCells(), Vector( ntotalEquations, 1e-10 ) );
     return psi;
 }
 
-bool Checkerboard_PN::isAbsorption( const Vector& pos ) const {
+bool Checkerboard_Moment::isAbsorption( const Vector& pos ) const {
     // Check whether pos is in absorption region
     std::vector<double> lbounds{ 1, 2, 3, 4, 5 };
     std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
@@ -121,16 +127,10 @@ bool Checkerboard_PN::isAbsorption( const Vector& pos ) const {
     return false;
 }
 
-bool Checkerboard_PN::isSource( const Vector& pos ) const {
+bool Checkerboard_Moment::isSource( const Vector& pos ) const {
     // Check whether pos is in source region
     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

@@ -30,7 +30,7 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
         }
         case PROBLEM_Checkerboard: {
             if( settings->GetSolverName() == PN_SOLVER || settings->GetSolverName() == MN_SOLVER )
-                return new Checkerboard_PN( settings, mesh );
+                return new Checkerboard_Moment( settings, mesh );
             else
                 return new Checkerboard_SN( settings, mesh );    // default
         }