150 lines
4.1 KiB
C++
150 lines
4.1 KiB
C++
// Exercise 5 – Parallelization of Exercise 3
|
||
//HOW TO RUN IT ON MY MAC:
|
||
//cd ~/Desktop/jacobi_oo_stl/jacobi_oo_stl
|
||
//export CPPFLAGS="-I/opt/homebrew/opt/libomp/include"
|
||
//export LDFLAGS="-L/opt/homebrew/opt/libomp/lib"
|
||
//export PATH="/opt/homebrew/opt/llvm/bin:$PATH"
|
||
|
||
//clang++ -std=c++17 -O3 \
|
||
-Xpreprocessor -fopenmp \
|
||
$CPPFLAGS \
|
||
main.cpp \
|
||
jacsolve.cpp \
|
||
getmatrix.cpp \
|
||
geom.cpp \
|
||
vdop.cpp \
|
||
userset.cpp \
|
||
$LDFLAGS -lomp \
|
||
-o Ex5
|
||
|
||
//./Ex5
|
||
|
||
#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;
|
||
}
|
||
|