129 lines
4.2 KiB
C++
129 lines
4.2 KiB
C++
// MPI code in C++.
|
|
// See [Gropp/Lusk/Skjellum, "Using MPI", p.33/41 etc.]
|
|
// and /opt/mpich/include/mpi2c++/comm.h for details
|
|
|
|
#include "geom.h"
|
|
#include "getmatrix.h"
|
|
#include "jacsolve.h"
|
|
#include "userset.h"
|
|
#include "vdop.h"
|
|
|
|
#include <chrono> // timing
|
|
#include <cmath>
|
|
#include <iostream>
|
|
using namespace std;
|
|
using namespace std::chrono; // timing
|
|
|
|
|
|
int main(int, char ** )
|
|
{
|
|
const int numprocs = 1;
|
|
const int myrank = 0;
|
|
|
|
if (myrank == 0)
|
|
{
|
|
cout << "\n There are " << numprocs << " processes running.\n \n";
|
|
}
|
|
|
|
const auto procx = static_cast<int>(sqrt(numprocs + 0.0));
|
|
const int procy = procx;
|
|
|
|
if (procy * procx != numprocs)
|
|
{
|
|
cout << "\n Wrong number of processors !\n \n";
|
|
}
|
|
else
|
|
{
|
|
// #####################################################################
|
|
// Here starts the real code
|
|
// #####################################################################
|
|
//bool ScaleUp = !true;
|
|
int nx, ny, NXglob, NYglob; /* number of local intervals on (xl,xr)=:nx, (yb,yt)=:ny */
|
|
//nx = 1024;
|
|
//ny = 1024;
|
|
nx = 100;
|
|
ny = 100;
|
|
NXglob = nx * procx;
|
|
NYglob = ny * procy;
|
|
cout << "Intervalls: " << NXglob << " x " << NYglob << endl;
|
|
|
|
// ##################### STL ###########################################
|
|
{
|
|
Mesh_2d_3_square const mesh(nx, ny);
|
|
//mesh.Debug();
|
|
|
|
CRS_Matrix SK(mesh); // CRS matrix
|
|
//SK.Debug();
|
|
vector<double> uv(SK.Nrows(), 0.0); // temperature
|
|
vector<double> fv(SK.Nrows(), 0.0); // r.h.s.
|
|
|
|
SK.CalculateLaplace(fv);
|
|
//SK.Debug();
|
|
|
|
//mesh.SetU(uv); // deprecated
|
|
//mesh.SetF(fv); // deprecated
|
|
// Two ways to initialize the vector
|
|
//mesh.SetValues(uv,f_zero); // functional
|
|
mesh.SetValues(uv, [](double x, double y) -> double {return 0.0 * x *y;} ); // lambda function
|
|
|
|
SK.ApplyDirichletBC(uv, fv);
|
|
//SK.Compare2Old(nnode, id, ik, sk);
|
|
//SK.Debug();
|
|
|
|
auto tstart = system_clock::now(); // start timer
|
|
|
|
JacobiSolve(SK, fv, uv ); // solve the system of equations
|
|
|
|
auto tend = system_clock::now(); // end timer
|
|
auto duration = duration_cast<microseconds>(tend - tstart);
|
|
auto t1 = static_cast<double>(duration.count()) / 1e6 ; // t1 in seconds
|
|
cout << "JacobiSolve: timing in sec. : " << t1 << endl;
|
|
|
|
//CompareVectors(uv, nnode, u, 1e-6); // Check correctness
|
|
|
|
//mesh.SaveVectorP("t.dat", uv);
|
|
//mesh.Visualize(uv);
|
|
|
|
}
|
|
// ##################### STL ###########################################
|
|
{
|
|
//Mesh_2d_3_matlab const mesh("square_tiny.txt");
|
|
Mesh_2d_3_matlab const mesh("square_100.txt");
|
|
//Mesh_2d_3_matlab const mesh("L_shape.txt");
|
|
//mesh.Debug();
|
|
|
|
CRS_Matrix SK(mesh); // CRS matrix
|
|
//SK.Debug();
|
|
|
|
vector<double> uv(SK.Nrows(), 0.0); // temperature
|
|
vector<double> fv(SK.Nrows(), 0.0); // r.h.s.
|
|
|
|
SK.CalculateLaplace(fv);
|
|
//SK.Debug();
|
|
|
|
//mesh.SetU(uv); // deprecated
|
|
// Two ways to initialize the vector
|
|
//mesh.SetValues(uv,f_zero); // user function
|
|
mesh.SetValues(uv, [](double x, double y) -> double {return 0.0 * x *y;} ); // lambda function
|
|
|
|
SK.ApplyDirichletBC(uv, fv);
|
|
//SK.Compare2Old(nnode, id, ik, sk);
|
|
//SK.Debug();
|
|
|
|
auto tstart = system_clock::now(); // start timer
|
|
|
|
JacobiSolve(SK, fv, uv ); // solve the system of equations
|
|
|
|
auto tend = system_clock::now(); // end timer
|
|
auto duration = duration_cast<microseconds>(tend - tstart);
|
|
auto t1 = static_cast<double>(duration.count()) / 1e6 ; // t1 in seconds
|
|
cout << "JacobiSolve: timing in sec. : " << t1 << endl;
|
|
|
|
//mesh.Write_ascii_matlab("uv.txt", uv);
|
|
//mesh.Visualize(uv);
|
|
}
|
|
return 0;
|
|
}
|
|
}
|
|
|
|
|