Homec4science
Diffusion PHPC-graded_hw_1 (master) Browse

R12662/poisson_wip/simulation.cc
ca0d989eb263master

Last Change
Raw File

simulation.cc
View Options

/* -------------------------------------------------------------------------- */
#include "simulation.hh"
/* -------------------------------------------------------------------------- */
#include <cmath>
#include <iostream>
#include <mpi.h>
#include <vector>
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
Simulation::Simulation(int m, int n, int n_iter, int prank, int psize)
: m_global_m(m), m_global_n(n), m_n_iter(n_iter), m_h_m(1. / m),
m_h_n(1. / n), m_grids(m, n), m_f(m, n),
m_dumper(new DumperASCII(m_grids.old())), m_prank(prank), m_psize(psize) {}
/* -------------------------------------------------------------------------- */
void Simulation::set_initial_conditions() {
for (int i = 0; i < m_global_m; i++) {
for (int j = 0; j < m_global_n; j++) {
m_f(i, j) = -2. * 100. * M_PI * M_PI * std::sin(10. * M_PI * ( i + (int)(m_prank*m_global_n/m_psize)) * m_h_m) *
std::sin(10. * M_PI * j * m_h_n);
} // assuming square global problem size
}
}
/* -------------------------------------------------------------------------- */
std::tuple<float, int> Simulation::compute() {
int s = 0;
float l2 = 0;
int next = m_prank+1;
int prev = m_prank-1;
do {
/// Compute step
l2 = compute_step();
Grid & u = m_grids.current();
// Prepare buffers
float *send_a = 0, *send_b = 0, *recv_a = 0, *recv_b = 0;
send_a = new float [m_global_n];
send_b = new float [m_global_n];
recv_a = new float [m_global_n];
recv_b = new float [m_global_n];
for (int i = 0; i < m_global_n; i++){
send_a[i] = u(1, i);
send_b[i] = u(m_global_m-2, i);
}
/// Launch communications
/* Strategy: even processes first communicate with the process below them */
for (int s = 0; s < 2; s++) {
// Exchange information with block above (unless you're the first block)
if (m_prank > 0 && m_prank%2 != s){
// Send/recieve buffer content
MPI_Request request_a;
MPI_Isend(send_a, m_global_n, MPI_FLOAT, prev, 13, MPI_COMM_WORLD, &request_a);
MPI_Irecv(recv_a, m_global_n, MPI_FLOAT, prev, 13, MPI_COMM_WORLD, &request_a);
MPI_Wait(&request_a, MPI_STATUS_IGNORE);
//Accept data
for (int i = 0; i < m_global_n; i++) {
u(0, i) = recv_a[i];
}
}
// Exchange information with block below (unless you're the last block)
if (m_prank < m_psize - 1 && m_prank%2 == s){
// Send/recieve buffer content
MPI_Request request_b;
MPI_Isend(send_b, m_global_n, MPI_FLOAT, next, 13, MPI_COMM_WORLD, &request_b);
MPI_Irecv(recv_b, m_global_n, MPI_FLOAT, next, 13, MPI_COMM_WORLD, &request_b);
MPI_Wait(&request_b, MPI_STATUS_IGNORE);
// Accept data
for (int i = 0; i < m_global_n; i++) {
u(m_global_m-1, i) = recv_b[i];
}
}
}
// Deallocate used memory once it is safe to do so
//MPI_Wait(&request_a, MPI_STATUS_IGNORE);
delete[] send_a; delete[] recv_a;
//MPI_Wait(&request_b, MPI_STATUS_IGNORE);
delete[] send_b; delete[] recv_b;
m_grids.swap();
++s;
} while (s < m_n_iter);
return std::make_tuple(l2, s);
}
/* -------------------------------------------------------------------------- */
float Simulation::compute_step() {
float l2 = 0.;
Grid & u = m_grids.current();
Grid & uo = m_grids.old();
for (int i = 1; i < m_global_m - 1; i++) {
for (int j = 1; j < m_global_n - 1; j++) {
// computation of the new step
u(i, j) = 0.25 * (uo(i - 1, j) + uo(i + 1, j) + uo(i, j - 1) +
uo(i, j + 1) - m_f(i, j) * m_h_m * m_h_n);
// L2 norm
l2 += (uo(i, j) - u(i, j)) * (uo(i, j) - u(i, j));
}
}
return l2;
}
c4science · Help