Main File

Sheet5/Ex5_Second_Attempt/main.cpp

@@ -0,0 +1,130 @@
+// Exercise 5 – Parallelization of Exercise 3
+
+#include "geom.h"
+#include "getmatrix.h" //in here it is defined the parallel function 
+#include "jacsolve.h" //in here it is defined the parallel function 
+#include "userset.h"
+#include "vdop.h"
+
+#include <chrono>
+#include <cmath>
+#include <iostream>
+#include <omp.h>
+
+using namespace std;
+using namespace std::chrono;
+
+int main(int, char **)
+{
+    // ------------------------------------------------------------
+    // OpenMP info
+    // ------------------------------------------------------------
+    #pragma omp parallel
+    {
+        #pragma omp single
+        cout << "Using " << omp_get_num_threads()
+             << " OpenMP threads\n\n";
+    }
+
+    // ------------------------------------------------------------
+    // Domain decomposition (MPI-style layout kept, MPI=1)
+    // ------------------------------------------------------------
+    const int numprocs = 1;
+    const int myrank   = 0;
+
+    if (myrank == 0)
+        cout << "There are " << numprocs << " processes running.\n\n";
+
+    const int procx = static_cast<int>(sqrt(numprocs + 0.0));
+    const int procy = procx;
+
+    if (procx * procy != numprocs)
+    {
+        cout << "Wrong number of processors!\n";
+        return 1;
+    }
+
+    // ------------------------------------------------------------
+    // Problem size
+    // ------------------------------------------------------------
+    int nx = 100;
+    int ny = 100;
+    int NXglob = nx * procx;
+    int NYglob = ny * procy;
+
+    cout << "Intervals: " << NXglob << " x " << NYglob << "\n";
+
+// ##################### STL ###########################################
+        {
+        Mesh_2d_3_square const mesh(nx, ny);
+        //mesh.Debug();
+
+        CRS_Matrix SK(mesh);                  // sparse FEM matrix
+        //SK.Debug();
+        vector<double> u(SK.Nrows(), 0.0);    // solution
+        vector<double> f(SK.Nrows(), 0.0);    // RHS
+
+        // Parallel FEM matrix assembly
+        auto t0 = high_resolution_clock::now();
+
+        SK.CalculateLaplace_omp(f);   // PARALLEL assembly
+        //SK.Debug();
+
+        auto t1 = high_resolution_clock::now();
+        cout << "FEM assembly time (parallel): "
+             << duration<double>(t1 - t0).count()
+             << " s\n";
+
+        // Initial condition
+        mesh.SetValues(u, [](double x, double y) -> double {return 0.0 * x *y;} ); // lambda function
+
+        // Apply Dirichlet BC
+        SK.ApplyDirichletBC(u, f);
+        //SK.Debug();
+
+        // Parallel Jacobi solver
+        t0 = high_resolution_clock::now();
+
+        JacobiSolve_omp(SK, f, u);    // PARALLEL Jacobi
+
+        t1 = high_resolution_clock::now();
+        cout << "JacobiSolve time (parallel): "
+             << duration<double>(t1 - t0).count()
+             << " s\n";
+    }
+
+    // ============================================================
+    {
+        Mesh_2d_3_matlab const mesh("square_100.txt");
+
+        CRS_Matrix SK(mesh);
+        //SK.Debug();
+        vector<double> u(SK.Nrows(), 0.0);
+        vector<double> f(SK.Nrows(), 0.0);
+
+        auto t0 = high_resolution_clock::now();
+
+        SK.CalculateLaplace_omp(f);   // PARALLEL assembly
+        //SK.Debug();
+
+        auto t1 = high_resolution_clock::now();
+        cout << "FEM assembly time (parallel, matlab mesh): "
+             << duration<double>(t1 - t0).count()
+             << " s\n";
+
+        mesh.SetValues(u, [](double x, double y) -> double {return 0.0 * x *y;} ); // lambda function
+        SK.ApplyDirichletBC(u, f);
+        //SK.Debug();
+
+        t0 = high_resolution_clock::now();
+
+        JacobiSolve_omp(SK, f, u);    // PARALLEL Jacobi
+
+        t1 = high_resolution_clock::now();
+        cout << "JacobiSolve time (parallel, matlab mesh): "
+             << duration<double>(t1 - t0).count()
+             << " s\n";
+    }
+
+    return 0;
+}