Merge pull request #33 from motrieu/feature/finalMpi/simon

Feature/final mpi/simon

Author: simongrether

src/CMakeLists.txt

@@ -1,23 +1,29 @@
 cmake_minimum_required(VERSION 3.8)
 
 # Define the project name.
-project(numsim)
+project(numsim_parallel)
 
 # Specify the name of the executable (${PROJECT_NAME} which is equal to what was set in the project() command).
 # Also specify the source files.
 add_executable(${PROJECT_NAME}
    main.cpp
    computation/computation.cpp
+   computation/computationParallel.cpp
    discretization/centralDifferences.cpp
    discretization/discretization.cpp
    discretization/donorCell.cpp
    discretization/staggeredGrid.cpp
-   output_writer/output_writer.cpp
+   output_writer/output_writer_paraview_parallel.cpp
    output_writer/output_writer_paraview.cpp
+   output_writer/output_writer_text_parallel.cpp
    output_writer/output_writer_text.cpp
+   output_writer/output_writer.cpp
+   partitioning/partitioning.cpp
    pressure_solver/gaussSeidel.cpp
    pressure_solver/pressureSolver.cpp
+   pressure_solver/pressureSolverParallel.cpp
    pressure_solver/sor.cpp
+   pressure_solver/sorParallel.cpp
    settings/settings.cpp
    storage/array2d.cpp
    storage/fieldVariable.cpp
@@ -53,3 +59,19 @@ set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 message("CMAKE_BUILD_TYPE: ${CMAKE_BUILD_TYPE}")
+
+find_package(MPI REQUIRED)
+
+include_directories(${MPI_INCLUDE_PATH})
+target_link_libraries(${PROJECT_NAME} ${MPI_LIBRARIES})
+
+if(MPI_COMPILE_FLAGS)
+  set_target_properties(${PROJECT_NAME} PROPERTIES
+    COMPILE_FLAGS "${MPI_COMPILE_FLAGS}")
+endif()
+
+if(MPI_LINK_FLAGS)
+  set_target_properties(${PROJECT_NAME} PROPERTIES
+    LINK_FLAGS "${MPI_LINK_FLAGS}")
+endif()
+

src/computation/computation.cpp

@@ -32,7 +32,7 @@ void Computation::runSimulation()
 
         time += dt_;
 
-        (*outputWriterParaview_).writeFile(time);
+        // (*outputWriterParaview_).writeFile(time);
         // (*outputWriterText_).writeFile(time);
     }
 }
@@ -48,9 +48,11 @@ void Computation::initialize(int argc, char *argv[])
     settings_.loadFromFile(filename);
 
     // calculates mesh width in x- and y-direction based on given parameters
-    const double meshWidthX = settings_.physicalSize[0] / settings_.nCells[0];
-    const double meshWidthY = settings_.physicalSize[1] / settings_.nCells[1];
-    meshWidth_ = {meshWidthX, meshWidthY};
+    dx_ = settings_.physicalSize[0] / settings_.nCells[0];
+    dy_ = settings_.physicalSize[1] / settings_.nCells[1];
+    meshWidth_ = {dx_, dy_};
+    dxSquared_ = dx_ * dx_;
+    dySquared_ = dy_ * dy_;
 
     // either the Central Differences or the Donor cell scheme is used
     if (settings_.useDonorCell)
@@ -66,7 +68,7 @@ void Computation::initialize(int argc, char *argv[])
     else
         throw std::invalid_argument("Only SOR and GaussSeidel are supported as pressure solvers.");
 
-    outputWriterParaview_ = std::make_unique<OutputWriterParaview>(discretization_);
+    // outputWriterParaview_ = std::make_unique<OutputWriterParaview>(discretization_);
     // outputWriterText_ = std::make_unique<OutputWriterText>(discretization_);
 
     // boundary conditions for u and v on the boundary faces only need to be set once in the beginning of the computation
@@ -137,9 +139,7 @@ void Computation::applyPreliminaryBCOnBoundary()
 
 void Computation::computeTimeStepWidth()
 {
-    const double dx = meshWidth_[0];
-    const double dy = meshWidth_[1];
-    const double dtDiffusive = (settings_.re/2.0) * (dx*dx * dy*dy) / (dx*dx + dy*dy);
+    const double dtDiffusive = (settings_.re/2.0) * (dxSquared_ * dySquared_) / (dxSquared_ + dySquared_);
 
     double uAbsMax = 0.0;
     for (int i=(*discretization_).uIBegin()-1; i < (*discretization_).uIEnd()+1; i++)
@@ -162,12 +162,17 @@ void Computation::computeTimeStepWidth()
         }
     }
 
-    const double dtConvectiveU = dx / uAbsMax;
-    const double dtConvectiveV = dy / vAbsMax;
+    double dtConvectiveU = settings_.maximumDt;
+    double dtConvectiveV = settings_.maximumDt;
+    if (uAbsMax > 0.0)
+        dtConvectiveU = dx_ / uAbsMax;
+    if (vAbsMax > 0.0)
+        dtConvectiveV = dy_ / vAbsMax;
 
     // makes sure that all stability conditions (the convective conditions and the diffusive condition) are fulfilled
     // and that the demanded maximal time step is not exceeded
-    dt_ = settings_.tau * std::min({dtDiffusive, dtConvectiveU, dtConvectiveV, settings_.maximumDt});
+    dt_ = settings_.tau * std::min({dtDiffusive, dtConvectiveU, dtConvectiveV});
+    dt_ = std::min(dt_, settings_.maximumDt);
 }
 
 void Computation::computePreliminaryVelocities()
@@ -204,8 +209,8 @@ void Computation::computeRightHandSide()
     {
         for (int j=(*discretization_).pJBegin(); j < (*discretization_).pJEnd(); j++)
         {
-            const double fDiffQuotient = ((*discretization_).f(i,j) - (*discretization_).f(i-1,j)) / meshWidth_[0];
-            const double gDiffQuotient = ((*discretization_).g(i,j) - (*discretization_).g(i,j-1)) / meshWidth_[1];
+            const double fDiffQuotient = ((*discretization_).f(i,j) - (*discretization_).f(i-1,j)) / dx_;
+            const double gDiffQuotient = ((*discretization_).g(i,j) - (*discretization_).g(i,j-1)) / dy_;
             (*discretization_).rhs(i,j) = (1.0/dt_) * (fDiffQuotient + gDiffQuotient);
         }
     }
@@ -222,15 +227,15 @@ void Computation::computeVelocities()
     {
         for (int j=(*discretization_).uJBegin(); j < (*discretization_).uJEnd(); j++)
         {
-            (*discretization_).u(i,j) = (*discretization_).f(i,j) - (dt_/meshWidth_[0])
+            (*discretization_).u(i,j) = (*discretization_).f(i,j) - (dt_/dx_)
                                             * ((*discretization_).p(i+1,j) - (*discretization_).p(i,j));
         }
     }
     for (int i=(*discretization_).vIBegin(); i < (*discretization_).vIEnd(); i++)
     {
         for (int j=(*discretization_).vJBegin(); j < (*discretization_).vJEnd(); j++)
         {
-            (*discretization_).v(i,j) = (*discretization_).g(i,j) - (dt_/meshWidth_[1])
+            (*discretization_).v(i,j) = (*discretization_).g(i,j) - (dt_/dy_)
                                             * ((*discretization_).p(i,j+1) - (*discretization_).p(i,j));
         }
     }

src/computation/computation.h

@@ -32,7 +32,7 @@ class Computation
     ///        pressure p, preliminary velocities F and G and the right hand side of the pressure Poisson problem (rhs)
     void runSimulation();
 
-private:
+protected:
 
     /// @brief applies boundary conditions of u and v on the actual boundary faces
     ///        manipulates the in the data structure stored values such that they are overwritten
@@ -70,6 +70,7 @@ class Computation
     ///        manipulates the in the data structure stored values such that they are overwritten
     void computeVelocities();
 
+
     /// @brief contains parameters used for the current scenario
     Settings settings_;
 
@@ -84,11 +85,23 @@ class Computation
 
     /// @brief unique pointer to the output writer which produces readable txt-files
     std::unique_ptr<OutputWriterText> outputWriterText_;
+
 
     /// @brief array of doubles that contains the mesh width in x- and y-direction
     std::array<double,2> meshWidth_;
 
     /// @brief current time step width, is set in each time iteration
     double dt_;
 
+    /// @brief mesh width in x-direction
+    double dx_;
+
+    /// @brief squared mesh width in x-direction
+    double dxSquared_;
+
+    /// @brief mesh width in y-direction
+    double dy_;
+
+    /// @brief squared mesh width in y-direction
+    double dySquared_;
 };