Hot fix mpi deploy

CMakeLists.txt

@@ -9,7 +9,6 @@ option( BUILD_UNITY "enables unity build for faster compile times" ON )
 option( BUILD_CODE_COV "enables compiler option required for code coverage analysis" OFF )
 option( BUILD_ML "enables build with tensorflow backend access" OFF )
 option( BUILD_MPI "enables build with MPI access" OFF )
-
 #################################################
 
 
@@ -26,15 +25,19 @@ if( BUILD_UNITY AND NOT BUILD_CODE_COV )
 endif()
 #################################################
 
-
 ### LIBRARIES ###################################
 set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake ${CMAKE_MODULE_PATH})
 
 find_package( OpenMP REQUIRED )
 
+message(STATUS "MPI build flag: ${BUILD_MPI}")
 if( BUILD_MPI )
+    add_definitions(-DIMPORT_MPI)
     find_package( MPI REQUIRED )
     include_directories( ${MPI_INCLUDE_PATH} )
+    message(STATUS "C++ compiler: ${CMAKE_CXX_COMPILER}")
+    message(STATUS "MPI C++ compiler: ${MPI_CXX_COMPILER}")
+
 endif()
 
 find_package( LAPACK REQUIRED )
@@ -63,6 +66,8 @@ include_directories( ${CMAKE_SOURCE_DIR}/ext/spdlog/include )
 
 if( BUILD_MPI )
     set( CORE_LIBRARIES ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} ${MPI_LIBRARIES} ${VTK_LIBRARIES} OpenMP::OpenMP_CXX -lstdc++fs )
+    message( STATUS "MPI: Libraries loaded" )
+
 else()
     set( CORE_LIBRARIES ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}  ${VTK_LIBRARIES} OpenMP::OpenMP_CXX -lstdc++fs )
 endif()

include/common/config.hpp

@@ -48,9 +48,11 @@ class Config
     // std::vector<double> _1dIntegrationBounds; /*!< @brief Quadrature Order*/
 
     // Mesh
-    unsigned _nCells;                      /*!< @brief Number of cells in the mesh */
-    unsigned short _dim;                   /*!< @brief spatial dimensionality of the mesh/test case */
-    bool _forcedConnectivityWrite; /*!< @brief If true, the meshconnectivity is always computed and written to .con file */
+    unsigned _nCells;                            /*!< @brief Number of cells in the mesh */
+    unsigned short _dim;                         /*!< @brief spatial dimensionality of the mesh/test case */
+    bool _forcedConnectivityWrite;               /*!< @brief If true, the meshconnectivity is always computed and written to .con file */
+    bool _loadrestartSolution;                   /*!< @brief If true, the simulation loads a restart solution from file */
+    unsigned long _saveRestartSolutionFrequency; /*!< @brief Frequency for saving restart solution to file */
 
     // Boundary Conditions
     /*!< @brief List of all Pairs (marker, BOUNDARY_TYPE), e.g. (farfield,DIRICHLET).
@@ -312,6 +314,8 @@ class Config
     unsigned GetNCells() const { return _nCells; }
     unsigned short GetDim() const { return _dim; }
     bool inline GetForcedConnectivity() const { return _forcedConnectivityWrite; }
+    bool inline GetLoadRestartSolution() const { return _loadrestartSolution; }
+    unsigned long inline GetSaveRestartSolutionFrequency() const { return _saveRestartSolutionFrequency; }
 
     // Solver Structure
     bool inline GetHPC() const { return _HPC; }

include/common/globalconstants.hpp

@@ -62,6 +62,7 @@ enum QUAD_NAME {
     QUAD_LevelSymmetric,
     QUAD_Lebedev,
     QUAD_LDFESA,
+    QUAD_SphericalTessalation2D,
     QUAD_Product,
     QUAD_Rectangular1D,
     QUAD_Rectangular2D,
@@ -82,6 +83,7 @@ inline std::map<std::string, QUAD_NAME> Quadrature_Map{ { "MONTE_CARLO", QUAD_Mo
                                                         { "LEVEL_SYMMETRIC", QUAD_LevelSymmetric },
                                                         { "LEBEDEV", QUAD_Lebedev },
                                                         { "LDFESA", QUAD_LDFESA },
+                                                        { "TESSALATION", QUAD_SphericalTessalation2D },
                                                         { "MIDPOINT_1D", QUAD_Midpoint1D },
                                                         { "MIDPOINT_2D", QUAD_Midpoint2D },
                                                         { "MIDPOINT_3D", QUAD_Midpoint3D },
@@ -195,7 +197,8 @@ enum SCALAR_OUTPUT {
     TOTAL_PARTICLE_ABSORPTION_HORIZONTAL,
     PROBE_MOMENT_TIME_TRACE,
     VAR_ABSORPTION_GREEN,
-    VAR_ABSORPTION_GREEN_LINE
+    ABSORPTION_GREEN_BLOCK,
+    ABSORPTION_GREEN_LINE
 };
 
 inline std::map<std::string, SCALAR_OUTPUT> ScalarOutput_Map{ { "ITER", ITER },
@@ -219,7 +222,8 @@ inline std::map<std::string, SCALAR_OUTPUT> ScalarOutput_Map{ { "ITER", ITER },
                                                               { "TOTAL_PARTICLE_ABSORPTION_HORIZONTAL", TOTAL_PARTICLE_ABSORPTION_HORIZONTAL },
                                                               { "PROBE_MOMENT_TIME_TRACE", PROBE_MOMENT_TIME_TRACE },
                                                               { "VAR_ABSORPTION_GREEN", VAR_ABSORPTION_GREEN },
-                                                              { "VAR_ABSORPTION_GREEN_LINE", VAR_ABSORPTION_GREEN_LINE } };
+                                                              { "ABSORPTION_GREEN_BLOCK", ABSORPTION_GREEN_BLOCK },
+                                                              { "ABSORPTION_GREEN_LINE", ABSORPTION_GREEN_LINE } };
 
 // Spherical Basis Name
 enum SPHERICAL_BASIS_NAME { SPHERICAL_HARMONICS, SPHERICAL_MONOMIALS, SPHERICAL_MONOMIALS_ROTATED };

include/common/io.hpp

@@ -43,6 +43,26 @@ std::string ParseArguments( int argc, char* argv[] );
 
 void PrintLogHeader( std::string inputFile );
 
+void WriteRestartSolution( const std::string& baseOutputFile,
+                           const std::vector<double>& solution,
+                           const std::vector<double>& scalarFlux,
+                           int rank,
+                           int idx_iter,
+                           double totalAbsorptionHohlraumCenter,
+                           double totalAbsorptionHohlraumVertical,
+                           double totalAbsorptionHohlraumHorizontal,
+                           double totalAbsorption );
+
+int LoadRestartSolution( const std::string& baseInputFile,
+                         std::vector<double>& solution,
+                         std::vector<double>& scalarFlux,
+                         int rank,
+                         unsigned long nCells,
+                         double& totalAbsorptionHohlraumCenter,
+                         double& totalAbsorptionHohlraumVertical,
+                         double& totalAbsorptionHohlraumHorizontal,
+                         double& totalAbsorption );
+
 // Matrix createSU2MeshFromImage( std::string imageName, std::string SU2Filename ); Deprecated
 
 #endif    // IO_H

include/common/mesh.hpp

@@ -129,10 +129,14 @@ class Mesh
      *  @return cell_idx: unsigned */
     unsigned GetCellOfKoordinate( const double x, const double y ) const;
 
-      /*! @brief Returns index of cells contained in the ball around the coordinate (x,y) with radius r
+    /*! @brief Returns index of cells contained in the ball around the coordinate (x,y) with radius r
      *  @return cell_idxs:   std::vector<unsigned> */
     std::vector<unsigned> GetCellsofBall( const double x, const double y, const double r ) const;
 
+    /*! @brief Returns index of cells contained in the rectangle with specified corner coordinates/*
+     *  @return cell_idxs:   std::vector<unsigned> */
+    std::vector<unsigned> GetCellsofRectangle( const std::vector<std::vector<double>>& cornercoordinates ) const;
+
     /*! @brief ComputeSlopes calculates the slope in every cell into x and y direction using the divergence theorem.
      *  @param nq is number of quadrature points
      *  @param psiDerX is slope in x direction (gets computed. Slope is stored here)

include/quadratures/qsphericaltessalation.hpp

@@ -0,0 +1,37 @@
+#ifndef QSPHERICALTRIANGLE_H
+#define QSPHERICALTRIANGLE_H
+
+#include "quadraturebase.hpp"
+
+class QSphericalTessalation : public QuadratureBase
+{
+    // Implementation is done accordingly to Kendall Atkinson 1981, Australian Matematical Society.
+
+  public:
+    QSphericalTessalation( Config* settings );
+
+    virtual ~QSphericalTessalation() {}
+
+  protected:
+    virtual bool CheckOrder();
+    virtual inline void SetName() override { _name = "Tesselated Spherical Triangle Quadrature"; }
+    void SetNq() override;
+    void SetPointsAndWeights() override;
+    void SetConnectivity() override;
+
+  private:
+    std::array<double, 3> compute_centroid( const std::array<std::array<double, 3>, 3>& triangle );
+    std::array<double, 3> map_to_unit_sphere( const std::array<double, 3>& point );
+    double dot_product( const std::array<double, 3>& v1, const std::array<double, 3>& v2 );
+    double angle_between_vectors( const std::array<double, 3>& a, const std::array<double, 3>& b, const std::array<double, 3>& c );
+    double spherical_triangle_area( const std::array<double, 3>& a, const std::array<double, 3>& b, const std::array<double, 3>& c );
+    std::array<double, 3> midpoint( const std::array<double, 3>& p1, const std::array<double, 3>& p2 );
+    void reflect_and_permute( const std::vector<std::array<double, 3>>& points,
+                              const std::vector<double>& weights,
+                              std::vector<std::array<double, 3>>& full_points,
+                              std::vector<double>& full_weights );
+
+    std::vector<std::array<std::array<double, 3>, 3>> generate_tessellation( const std::array<std::array<double, 3>, 3>& triangle, int order );
+};
+
+#endif    // QSPHERICALTRIANGLE_H

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 */
+    int _idx_start_iter;  /*!< @brief index of first iteration */
     // 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,14 +111,18 @@ 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 */
 
-    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 */
-    std::vector<double> _absorptionValsIntegrated;    /*!< @brief Avg Absorption value at the sampleing points of lineGreen */
-    std::vector<double> _varAbsorptionValsIntegrated; /*!< @brief Var in Avg Absorption value at the sampleing points of lineGreen */
+    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 */
+    std::vector<double> _absorptionValsLineSegment; /*!< @brief Avg Absorption value at the sampleing points of lineGreen */
+
+    unsigned _nProbingCellsBlocksGreen;
+    std::vector<std::vector<unsigned>> _probingCellsBlocksGreen; /*!< @brief Indices of cells that contain a probing sensor blocks */
+    std::vector<double> _absorptionValsBlocksGreen;              /*!< @brief Avg Absorption value at the sampleing blocks of lineGreen */
 
     // Design parameters
     std::vector<double> _cornerUpperLeftGreen; /*!< @brief Coord of corner of the green area (minus thickness/2 of it) relative to the green center */
@@ -198,8 +202,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();