changed indices to unsigned long to prevent overflow for large meshes

include/solvers/snsolver_hpc.hpp

@@ -20,26 +20,26 @@ class SNSolverHPC
   private:
     int _rank;
     int _numProcs;
-    unsigned _localNSys;
-    unsigned _startSysIdx;
-    unsigned _endSysIdx;
+    unsigned long _localNSys;
+    unsigned long _startSysIdx;
+    unsigned long _endSysIdx;
 
     double _currTime;  /*!< @brief wall-time after current iteration */
     Config* _settings; /*!< @brief config class for global information */
     Mesh* _mesh;
     ProblemBase* _problem;
 
     // Time
-    unsigned _nIter; /*!< @brief number of time steps, for non CSD, this equals _nEnergies, for _csd, _maxIter =_nEnergies-1*/
-    double _dT;      /*!< @brief energy/time step size */
+    unsigned long _nIter; /*!< @brief number of time steps, for non CSD, this equals _nEnergies, for _csd, _maxIter =_nEnergies-1*/
+    double _dT;           /*!< @brief energy/time step size */
 
     // Mesh related members, memory optimized
-    unsigned _nCells; /*!< @brief number of spatial cells */
-    unsigned _nSys;   /*!< @brief number of equations in the transport system, i.e. num quad pts */
-    unsigned _nq;     /*!< @brief number of quadrature points */
-    unsigned _nDim;
-    unsigned _nNbr;
-    unsigned _nNodes;
+    unsigned long _nCells; /*!< @brief number of spatial cells */
+    unsigned long _nSys;   /*!< @brief number of equations in the transport system, i.e. num quad pts */
+    unsigned long _nq;     /*!< @brief number of quadrature points */
+    unsigned long _nDim;
+    unsigned long _nNbr;
+    unsigned long _nNodes;
 
     std::vector<double> _areas;         /*!< @brief surface area of all spatial cells,
                                            dim(_areas) = _NCells */
@@ -111,9 +111,9 @@ class SNSolverHPC
     double _curAbsorptionHohlraumVertical;
     double _curAbsorptionHohlraumHorizontal;
     double _varAbsorptionHohlraumGreen;
-    std::vector<  std::vector<unsigned>> _probingCellsHohlraum; /*!< @brief Indices of cells that contain a probing sensor */
-    std::vector<double> _probingMoments;         /*!< @brief Solution Momnets at the probing cells that contain a probing sensor */
-    unsigned _probingMomentsTimeIntervals;       /*!< @brief Solution Momnets at the probing cells that contain a probing sensor */
+    std::vector<std::vector<unsigned>> _probingCellsHohlraum; /*!< @brief Indices of cells that contain a probing sensor */
+    std::vector<double> _probingMoments;                      /*!< @brief Solution Momnets at the probing cells that contain a probing sensor */
+    unsigned _probingMomentsTimeIntervals;                    /*!< @brief Solution Momnets at the probing cells that contain a probing sensor */
 
     unsigned _nProbingCellsLineGreen;                 /*!< @brief Number of sampling cells that contain a probing sensor for the sliding window */
     std::vector<unsigned> _probingCellsLineGreen;     /*!< @brief Indices of cells that contain a probing sensor for the sliding window */
@@ -198,8 +198,8 @@ class SNSolverHPC
     void DrawPostSolverOutput();
 
     // Helper
-    unsigned Idx2D( unsigned idx1, unsigned idx2, unsigned len2 );
-    unsigned Idx3D( unsigned idx1, unsigned idx2, unsigned idx3, unsigned len2, unsigned len3 );
+    unsigned long Idx2D( unsigned long idx1, unsigned long idx2, unsigned long len2 );
+    unsigned long Idx3D( unsigned long idx1, unsigned long idx2, unsigned long idx3, unsigned long len2, unsigned long len3 );
     bool IsAbsorptionLattice( double x, double y ) const;
     void ComputeCellsPerimeterLattice();
 

src/common/mesh.cpp

@@ -604,8 +604,8 @@ unsigned Mesh::GetCellOfKoordinate( const double x, const double y ) const {
 std::vector<unsigned> Mesh::GetCellsofBall( const double x, const double y, const double r ) const {
     std::vector<unsigned> cells_in_ball;
 
-    // Experimental parallel implementation
-    // #pragma omp parallel for
+// Experimental parallel implementation
+#pragma omp parallel for
     for( unsigned idx_cell = 0; idx_cell < _numCells; idx_cell++ ) {
         // Assume GetCellCenter returns the center coordinates of the cell
         double cell_x = _cellMidPoints[idx_cell][0];
@@ -615,7 +615,7 @@ std::vector<unsigned> Mesh::GetCellsofBall( const double x, const double y, cons
         double distance = std::sqrt( ( cell_x - x ) * ( cell_x - x ) + ( cell_y - y ) * ( cell_y - y ) );
 
         if( distance <= r ) {
-            // #pragma omp critical
+#pragma omp critical
             { cells_in_ball.push_back( idx_cell ); }
         }
     }

src/problems/symmetrichohlraum.cpp

@@ -130,19 +130,23 @@ void SymmetricHohlraum::SetGhostCells() {
         if( vq[idx_q][0] < 0.0 ) right_inflow[idx_q] = 1.0;
     }
 
+#pragma omp parallel for
     for( unsigned idx_cell = 0; idx_cell < _mesh->GetNumCells(); idx_cell++ ) {
         double x = _mesh->GetCellMidPoints()[idx_cell][0];
         double y = _mesh->GetCellMidPoints()[idx_cell][1];
 
         if( cellBoundaries[idx_cell] == BOUNDARY_TYPE::NEUMANN || cellBoundaries[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) {
-            if( y < -0.6 )
-                ghostCellMap.insert( { idx_cell, vertical_flow } );
-            else if( y > 0.6 )
-                ghostCellMap.insert( { idx_cell, vertical_flow } );
-            else if( x < -0.6 )
-                ghostCellMap.insert( { idx_cell, left_inflow } );
-            else if( x > 0.6 )
-                ghostCellMap.insert( { idx_cell, right_inflow } );
+#pragma omp critical
+            {
+                if( y < -0.6 )
+                    ghostCellMap.insert( { idx_cell, vertical_flow } );
+                else if( y > 0.6 )
+                    ghostCellMap.insert( { idx_cell, vertical_flow } );
+                else if( x < -0.6 )
+                    ghostCellMap.insert( { idx_cell, left_inflow } );
+                else if( x > 0.6 )
+                    ghostCellMap.insert( { idx_cell, right_inflow } );
+            }
         }
     }
     _ghostCells = ghostCellMap;
@@ -318,6 +322,7 @@ std::vector<unsigned> SymmetricHohlraum::linspace2D( const std::vector<double>&
     double stepX = ( end[0] - start[0] ) / ( num_points - 1 );
     double stepY = ( end[1] - start[1] ) / ( num_points - 1 );
 
+#pragma omp parallel for
     for( unsigned i = 0; i < num_points; ++i ) {
         double x = start[0] + i * stepX;
         double y = start[1] + i * stepY;