added hw2 skeleton code

hw2/skeleton/p1/Makefile

+CXX=h5pcc
+CXXFLAGS = -Wpedantic -Wall -Wextra -std=c++11 -lstdc++ -O3 -g 
+
+all: main
+
+%.o: %.cpp
+	$(CXX) -c -o $@ $< $(CXXFLAGS)
+clean:
+	rm -rf *.o *.xmf *h5 *txt main
+
+main: main.o auxiliary.o wave.o
+	$(CXX) $(CXXFLAGS)-I. -o main $^
+
+run:
+	make clean ; make ; mpirun -n 8 ./main 100 2 1.0	
+run1:
+	make clean ; make ; mpirun -n 1 ./main 100 1 1.0	
+
+.PHONY: all clean run run1

hw2/skeleton/p1/README

+HOW TO RUN
+----------
+To run the code on Euler remove all previously loaded nodules (purge) and load the following modules:
+
+module purge; 
+module load new; 
+module load gcc; 
+module load open_mpi/3.0.0; 
+module load hdf5
+
+You can request 8 cores and an interative shell to work with:
+bsub -n 8 -W  04:00 -Is bash
+
+Then compile and run by typing
+make run
+
+
+WHAT TO MODIFY
+--------------
+You only need to change the code in wave.cpp
+When the code runs, it prints a checksum. This number should be the same for the solution you provide.
+
+
+CODE OUTPUT 
+-----------
+You are not asked to visualize any results in this question, but we thought I'd be nice if you can see what you are actually solving.
+Other outputs are .h5 files which can be used to visualize the solution. Each .h5 file corresponds to an .xmf file.
+You can use Paraview to open the .xmf files and visualize the results.
+Alternatively, you may comment out the line
+#define USE_HDF5
+with '//' in wave.h. In this case, the code will output simple .txt files with the solution, saved as
+x y z u
+...

hw2/skeleton/p1/auxiliary.cpp

+#include "wave.h"
+
+WaveEquation::WaveEquation(int a_N, int a_procs_per_dim) {
+  MPI_Comm_size(MPI_COMM_WORLD, &size);
+  procs_per_dim = a_procs_per_dim;
+  Ntot = a_N;
+  h = L / a_N;
+  N = a_N / procs_per_dim;
+
+  // the chosen numerical method is stable if  dt <= h/sqrt(3)
+  dt = h / sqrt(3.0);
+
+  FindCoordinates();
+  
+  origin[0] = coords[0] * N * h;
+  origin[1] = coords[1] * N * h;
+  origin[2] = coords[2] * N * h;
+  u = new double[(N + 2) * (N + 2) * (N + 2)];
+  for (int i0 = 0; i0 < N; i0++) {
+    double x0 = origin[0] + i0 * h + 0.5 * h;
+    for (int i1 = 0; i1 < N; i1++) {
+      double x1 = origin[1] + i1 * h + 0.5 * h;
+      for (int i2 = 0; i2 < N; i2++) {
+        double x2 = origin[2] + i2 * h + 0.5 * h;
+        u[(i0 + 1) * (N + 2) * (N + 2) + (i1 + 1) * (N + 2) + (i2 + 1)] =
+            Initial_Condition(x0, x1, x2);
+      }
+    }
+  }
+
+  u_old = new double[(N + 2) * (N + 2) * (N + 2)];
+  u_new = new double[(N + 2) * (N + 2) * (N + 2)];
+
+  for (int i0 = 1; i0 <= N; i0++)
+    for (int i1 = 1; i1 <= N; i1++)
+      for (int i2 = 1; i2 <= N; i2++) {
+        int m = i2 + i1 * (N + 2) + i0 * (N + 2) * (N + 2);
+        u_new[m] = u[m];
+        u_old[m] = u[m];
+        // assuming that u_old = u is equivalent to du/dt(t=0) = 0
+      }
+
+  aux = dt * dt / h / h;
+}
+
+WaveEquation::~WaveEquation() {
+  delete[] u;
+  delete[] u_old;
+  delete[] u_new;
+}
+
+double WaveEquation::Initial_Condition(double x0, double x1, double x2) {
+  double r =
+      (x0 - 0.5) * (x0 - 0.5) + (x1 - 0.5) * (x1 - 0.5); //+ (x2-0.5)*(x2-0.5);
+  return exp(-r / 0.1);
+}
+
+// do not change this function
+void WaveEquation::Print(int kt = 0) {
+#ifdef USE_HDF5
+
+  std::string name = "Wave__" + std::to_string(kt);
+  if (rank == 0) {
+    FILE *xmf = 0;
+    xmf = fopen((name + ".xmf").c_str(), "w");
+
+    std::stringstream s_head;
+
+    s_head << "<?xml version=\"1.0\" ?>\n";
+    s_head << "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []>\n";
+    s_head << "<Xdmf Version=\"2.0\">\n";
+    s_head << " <Domain>\n";
+
+    s_head << "     <Time Value=\"" << std::setprecision(10) << std::setw(10)
+           << t << "\"/>\n\n";
+
+    fprintf(xmf, (s_head.str()).c_str());
+
+    std::stringstream s;
+    s << "   <Grid GridType=\"Uniform\">\n";
+    s << "     <Topology TopologyType=\"3DCoRectMesh\" Dimensions=\" "
+      << Ntot + 1 << " " << Ntot + 1 << " " << Ntot + 1 << "\"/>\n\n";
+    s << "       <Geometry GeometryType=\"ORIGIN_DXDYDZ\">\n";
+    s << "          <DataItem Dimensions=\"3\" NumberType=\"Double\" "
+         "Precision=\"8\" Format=\"XML\">\n";
+    s << "                 " << 0 << " " << 0 << " " << 0 << "\n";
+    s << "          </DataItem>\n";
+    s << "          <DataItem Dimensions=\"3\" NumberType=\"Double\" "
+         "Precision=\"8\" Format=\"XML\">\n";
+    s << "                 " << std::setprecision(10) << std::setw(15) << h
+      << " " << std::setw(15) << h << " " << std::setw(15) << h << "\n";
+    s << "          </DataItem>\n";
+    s << "     </Geometry>\n\n";
+
+    s << "     <Attribute Name=\"data\" AttributeType=\" "
+      << "Scalar"
+      << "\" Center=\"Cell\">\n";
+    s << "       <DataItem Dimensions=\" " << Ntot << " " << Ntot << " " << Ntot
+      << " " << std::setw(10) << 1 << "\" NumberType=\"Float\" Precision=\" "
+      << sizeof(double) << "\" Format=\"HDF\">\n";
+    std::stringstream name_ss;
+    name_ss << "data"; // << std::setfill('0') << std::setw(10) << mm ;
+    std::string tmp = name_ss.str();
+    s << "        " << (name + ".h5").c_str() << ":/" << tmp << "\n";
+    s << "       </DataItem>\n";
+    s << "     </Attribute>\n";
+
+    s << "   </Grid>\n";
+
+    std::string st = s.str();
+
+    fprintf(xmf, st.c_str());
+
+    std::stringstream s_tail;
+    s_tail << " </Domain>\n";
+    s_tail << "</Xdmf>\n";
+
+    fprintf(xmf, (s_tail.str()).c_str());
+    fclose(xmf);
+  }
+
+  hid_t file_id, dataset_id, fspace_id, fapl_id, mspace_id;
+
+  H5open();
+  fapl_id = H5Pcreate(H5P_FILE_ACCESS);
+
+  H5Pset_fapl_mpio(fapl_id, cart_comm, MPI_INFO_NULL);
+
+  file_id =
+      H5Fcreate((name + ".h5").c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, fapl_id);
+
+  H5Pclose(fapl_id);
+
+  double *array_all = new double[N * N * N];
+
+  hsize_t count[4] = {static_cast<unsigned int>(N),
+                      static_cast<unsigned int>(N),
+                      static_cast<unsigned int>(N), 1};
+  hsize_t dims[4] = {static_cast<unsigned int>(Ntot),
+                     static_cast<unsigned int>(Ntot),
+                     static_cast<unsigned int>(Ntot), 1};
+  hsize_t offset[4] = {static_cast<unsigned int>(coords[0] * N),
+                       static_cast<unsigned int>(coords[1] * N),
+                       static_cast<unsigned int>(coords[2] * N), 0};
+
+  for (int i0 = 1; i0 <= N; i0++) {
+    for (int i1 = 1; i1 <= N; i1++) {
+      for (int i2 = 1; i2 <= N; i2++) {
+        array_all[(i2 - 1) + (i1 - 1) * N + (i0 - 1) * N * N] =
+            u[i2 + i1 * (N + 2) + i0 * (N + 2) * (N + 2)];
+      }
+    }
+  }
+
+  fapl_id = H5Pcreate(H5P_DATASET_XFER);
+  H5Pset_dxpl_mpio(fapl_id, H5FD_MPIO_COLLECTIVE);
+  fspace_id = H5Screate_simple(4, dims, NULL);
+
+  dataset_id = H5Dcreate(file_id, "data", H5T_NATIVE_DOUBLE, fspace_id,
+                         H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+
+  fspace_id = H5Dget_space(dataset_id);
+
+  H5Sselect_hyperslab(fspace_id, H5S_SELECT_SET, offset, NULL, count, NULL);
+
+  mspace_id = H5Screate_simple(4, count, NULL);
+
+  H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, mspace_id, fspace_id, fapl_id,
+           array_all);
+
+  H5Sclose(mspace_id);
+  H5Sclose(fspace_id);
+  H5Dclose(dataset_id);
+  H5Pclose(fapl_id);
+  H5Fclose(file_id);
+  H5close();
+
+  delete[] array_all;
+
+#else
+
+  string name = "Wave__" + std::to_string(kt) + ".txt";
+  MPI_File file;
+  MPI_File_open(cart_comm, name.c_str(), MPI_MODE_RDWR | MPI_MODE_CREATE,
+                MPI_INFO_NULL, &file);
+
+  stringstream ss;
+
+  for (int i0 = 0; i0 < N; i0++) {
+    double x0 = origin[0] + i0 * h + 0.5 * h;
+    for (int i1 = 0; i1 < N; i1++) {
+      double x1 = origin[1] + i1 * h + 0.5 * h;
+      for (int i2 = 0; i2 < N; i2++) {
+        double x2 = origin[2] + i2 * h + 0.5 * h;
+        ss << x0 << " " << x1 << " " << x2 << " "
+           << u[(i2 + 1) + (i1 + 1) * (N + 2) + (i0 + 1) * (N + 2) * (N + 2)]
+           << "\n";
+      }
+    }
+  }
+
+  std::string asciiData = ss.str();
+
+  MPI_Offset my_length = asciiData.size() * sizeof(char);
+
+  MPI_Offset offset;
+  MPI_Exscan(&my_length, &offset, 1, MPI_OFFSET, MPI_SUM, cart_comm);
+
+  MPI_File_write_at_all(file, offset, asciiData.data(), asciiData.size(),
+                        MPI_CHAR, MPI_STATUS_IGNORE);
+  MPI_File_close(&file);
+
+#endif
+}
+
+// do not change this function
+void WaveEquation::UpdateGridPoint(int i0, int i1, int i2) {
+  int m = i0 * (N + 2) * (N + 2) + i1 * (N + 2) + i2;
+  int ip1 = (i0 + 1) * (N + 2) * (N + 2) + (i1) * (N + 2) + (i2);
+  int im1 = (i0 - 1) * (N + 2) * (N + 2) + (i1) * (N + 2) + (i2);
+  int jp1 = (i0) * (N + 2) * (N + 2) + (i1 + 1) * (N + 2) + (i2);
+  int jm1 = (i0) * (N + 2) * (N + 2) + (i1 - 1) * (N + 2) + (i2);
+  int kp1 = (i0) * (N + 2) * (N + 2) + (i1) * (N + 2) + (i2 + 1);
+  int km1 = (i0) * (N + 2) * (N + 2) + (i1) * (N + 2) + (i2 - 1);
+  u_new[m] =
+      2.0 * u[m] - u_old[m] +
+      aux * (u[ip1] + u[im1] + u[jp1] + u[jm1] + u[kp1] + u[km1] - 6.0 * u[m]);
+}

hw2/skeleton/p1/main.cpp

+#include "wave.h"
+
+int main(int argc, char **argv) {
+  MPI_Init(&argc, &argv);
+
+  int rank, size;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  MPI_Comm_size(MPI_COMM_WORLD, &size);
+
+  if (argc != 4) {
+    if (rank == 0) {
+      std::cout << " Incorrect number of inputs. Run as:\n";
+      std::cout << " mpirun -n x ./main N y t \n";
+      std::cout << " where \n";
+      std::cout << " N= number of grid points per direction\n";
+      std::cout << " y= number of ranks per direction (=(x)^(1/3))\n";
+      std::cout << " t= final time of simulation \n";
+      std::cout << " Aborting...\n";
+    }
+    int err = 1;
+    MPI_Abort(MPI_COMM_WORLD, err);
+  }
+
+  int points = atoi(argv[1]);
+  int procs_per_dim = atoi(argv[2]);
+  double t_end = atof(argv[3]);
+
+  if (size != procs_per_dim * procs_per_dim * procs_per_dim) {
+    if (rank == 0)
+      std::cout << " Incorrect number of ranks per direction. Aborting... \n";
+    int err = 2;
+    MPI_Abort(MPI_COMM_WORLD, err);
+  }
+
+  if (points % procs_per_dim != 0) {
+    if (rank == 0)
+      std::cout << " Grid points per direction must be divisible by the number "
+                   "of ranks per direction. Aborting... \n";
+    int err = 3;
+    MPI_Abort(MPI_COMM_WORLD, err);
+  }
+
+  WaveEquation *Simulation = new WaveEquation(points, procs_per_dim);
+
+  double start = MPI_Wtime();
+  Simulation->run(t_end);
+  double finish = MPI_Wtime();
+
+  if (rank == 0)
+    std::cout << "Total time = " << finish - start << "\n";
+
+  delete Simulation;
+
+  MPI_Finalize();
+  return 0;
+}
\ No newline at end of file

hw2/skeleton/p1/wave.cpp

+#include "wave.h"
+
+/********************************************************************/ 
+/* Subquestion a: change the following function and use a Cartesian */ 
+/* topology to find coords[3], rank_plus[3] and rank_minus[3]       */
+/********************************************************************/
+void WaveEquation::FindCoordinates() {
+
+  int p = procs_per_dim;
+
+  cart_comm = MPI_COMM_WORLD;
+  MPI_Comm_rank(cart_comm, &rank);
+
+  coords[0] = rank / (p * p);
+
+  coords[1] = (rank - coords[0] * (p * p)) / p;
+
+  coords[2] = (rank - coords[0] * (p * p) - coords[1] * p) % p;
+
+  int coor_0_plus = (coords[0] + 1 + p) % p;
+  int coor_1_plus = (coords[1] + 1 + p) % p;
+  int coor_2_plus = (coords[2] + 1 + p) % p;
+
+  int coor_0_minus = (coords[0] - 1 + p) % p;
+  int coor_1_minus = (coords[1] - 1 + p) % p;
+  int coor_2_minus = (coords[2] - 1 + p) % p;
+
+  rank_plus[0] = (p * p) * coor_0_plus + coords[1] * p + coords[2];
+  rank_plus[1] = coords[0] * p * p + p * coor_1_plus + coords[2];
+  rank_plus[2] = coords[0] * p * p + coords[1] * p + coor_2_plus;
+
+  rank_minus[0] = (p * p) * coor_0_minus + coords[1] * p + coords[2];
+  rank_minus[1] = coords[0] * p * p + p * coor_1_minus + coords[2];
+  rank_minus[2] = coords[0] * p * p + coords[1] * p + coor_2_minus;
+}
+
+// This function is not needed when custom datatypes are used
+void WaveEquation::pack_face(double *pack, int array_of_sizes[3],
+                             int array_of_subsizes[3], int array_of_starts[3]) {
+  int p0 = array_of_subsizes[0];
+  int p1 = array_of_subsizes[1];
+  int p2 = array_of_subsizes[2];
+
+  int n1 = array_of_sizes[1];
+  int n2 = array_of_sizes[2];
+
+  for (int i0 = array_of_starts[0]; i0 < array_of_starts[0] + p0; i0++)
+    for (int i1 = array_of_starts[1]; i1 < array_of_starts[1] + p1; i1++)
+      for (int i2 = array_of_starts[2]; i2 < array_of_starts[2] + p2; i2++) {
+        int i = (i0 - array_of_starts[0]) * p1 * p2 +
+                (i1 - array_of_starts[1]) * p2 + (i2 - array_of_starts[2]);
+        pack[i] = *(u + i0 * n1 * n2 + i1 * n2 + i2);
+      }
+}
+
+// This function is not needed when custom datatypes are used
+void WaveEquation::unpack_face(double *pack, int array_of_sizes[3],
+                               int array_of_subsizes[3],
+                               int array_of_starts[3]) {
+  int p0 = array_of_subsizes[0];
+  int p1 = array_of_subsizes[1];
+  int p2 = array_of_subsizes[2];
+
+  int n1 = array_of_sizes[1];
+  int n2 = array_of_sizes[2];
+
+  for (int i0 = array_of_starts[0]; i0 < array_of_starts[0] + p0; i0++)
+    for (int i1 = array_of_starts[1]; i1 < array_of_starts[1] + p1; i1++)
+      for (int i2 = array_of_starts[2]; i2 < array_of_starts[2] + p2; i2++) {
+        int i = (i0 - array_of_starts[0]) * p1 * p2 +
+                (i1 - array_of_starts[1]) * p2 + (i2 - array_of_starts[2]);
+        *(u + i0 * n1 * n2 + i1 * n2 + i2) = pack[i];
+      }
+}
+
+
+
+
+
+/********************************************************************/ 
+/* Subquestion b: you should no longer need the functions pack_face */
+/* and unpack_face nor should you need to allocate memory by using  */
+/* double *pack[6] and double *unpack[6].                           */
+/********************************************************************/
+void WaveEquation::run(double t_end) {
+
+  t = 0;
+
+  /********************************************************************/ 
+  /* Subquestion b: you need to define 12 custom datatypes.           */
+  /* For sending data, six datatypes (one per face) are required.     */
+  /* For receiving data, six more datatypes are required.             */
+  /* You should use MPI_Type_create_subarray for those datatypes.     */
+  /********************************************************************/
+
+
+  /* Subquestion b: You can delete this part when the subquestion is completed. */
+  /************************************************************************************************/
+  double *pack[6];
+  pack[0] = new double[N * N];
+  pack[1] = new double[N * N];
+  pack[2] = new double[N * N];
+  pack[3] = new double[N * N];
+  pack[4] = new double[N * N];
+  pack[5] = new double[N * N];
+  double *unpack[6];
+  unpack[0] = new double[N * N];
+  unpack[1] = new double[N * N];
+  unpack[2] = new double[N * N];
+  unpack[3] = new double[N * N];
+  unpack[4] = new double[N * N];
+  unpack[5] = new double[N * N];
+  /************************************************************************************************/
+
+
+
+
+
+
+  /* Subquestion b: Create and commit custom datatypes here */
+  /************************************************************************************************/
+  MPI_Datatype SEND_FACE_PLUS [3];
+  MPI_Datatype SEND_FACE_MINUS[3]; 
+
+  MPI_Datatype RECV_FACE_PLUS [3];
+  MPI_Datatype RECV_FACE_MINUS[3]; 
+  /************************************************************************************************/
+
+
+
+
+  int count = 0;
+  do {
+    if (count % 100 == 0) {
+      if (rank == 0)
+        std::cout << count << " t=" << t << "\n";
+      Print(count);
+    }
+
+
+  /* Subquestion b: You can delete this part when the subquestion is completed. */
+  /************************************************************************************************/
+    int array_of_sizes[3] = {N + 2, N + 2, N + 2};
+
+    {
+      int array_of_subsizes[3] = {1, N, N};
+      int array_of_starts[3] = {N, 1, 1};
+      pack_face(pack[0], array_of_sizes, array_of_subsizes, array_of_starts);
+    }
+
+    {
+      int array_of_subsizes[3] = {1, N, N};
+      int array_of_starts[3] = {1, 1, 1};
+      pack_face(pack[1], array_of_sizes, array_of_subsizes, array_of_starts);
+    }
+
+    {
+      int array_of_subsizes[3] = {N, 1, N};
+      int array_of_starts[3] = {1, N, 1};
+      pack_face(pack[2], array_of_sizes, array_of_subsizes, array_of_starts);
+    }
+
+    {
+      int array_of_subsizes[3] = {N, 1, N};
+      int array_of_starts[3] = {1, 1, 1};
+      pack_face(pack[3], array_of_sizes, array_of_subsizes, array_of_starts);
+    }
+
+    {
+      int array_of_subsizes[3] = {N, N, 1};
+      int array_of_starts[3] = {1, 1, N};
+      pack_face(pack[4], array_of_sizes, array_of_subsizes, array_of_starts);
+    }
+
+    {
+      int array_of_subsizes[3] = {N, N, 1};
+      int array_of_starts[3] = {1, 1, 1};
+      pack_face(pack[5], array_of_sizes, array_of_subsizes, array_of_starts);
+    }
+  /************************************************************************************************/
+
+
+
+
+
+    MPI_Request request[12];
+
+    /* Subquestion b: Replace the sends and receives with ones that correspond to custom datatypes*/
+    /**********************************************************************************************/
+    MPI_Irecv(unpack[0], N * N, MPI_DOUBLE, rank_minus[0], 100, cart_comm,&request[0]);
+    MPI_Isend(pack  [0], N * N, MPI_DOUBLE, rank_plus [0], 100, cart_comm,&request[1]);
+
+    MPI_Irecv(unpack[1], N * N, MPI_DOUBLE, rank_plus [0], 101, cart_comm, &request[2]);
+    MPI_Isend(pack  [1], N * N, MPI_DOUBLE, rank_minus[0], 101, cart_comm, &request[3]);
+
+    MPI_Irecv(unpack[2], N * N, MPI_DOUBLE, rank_minus[1], 200, cart_comm, &request[4]);
+    MPI_Isend(pack  [2], N * N, MPI_DOUBLE, rank_plus [1], 200, cart_comm, &request[5]);
+
+    MPI_Irecv(unpack[3], N * N, MPI_DOUBLE, rank_plus [1], 201, cart_comm, &request[6]);
+    MPI_Isend(pack  [3], N * N, MPI_DOUBLE, rank_minus[1], 201, cart_comm, &request[7]);
+
+    MPI_Irecv(unpack[4], N * N, MPI_DOUBLE, rank_minus[2], 300, cart_comm, &request[8]);
+    MPI_Isend(pack  [4], N * N, MPI_DOUBLE, rank_plus [2], 300, cart_comm, &request[9]);
+
+    MPI_Irecv(unpack[5], N * N, MPI_DOUBLE, rank_plus [2], 301, cart_comm, &request[10]);
+    MPI_Isend(pack  [5], N * N, MPI_DOUBLE, rank_minus[2], 301, cart_comm, &request[11]);
+    /**********************************************************************************************/
+
+    // Wait for communication to finish
+    MPI_Waitall(12, &request[0], MPI_STATUSES_IGNORE);
+
+
+
+
+  /* Subquestion b: You can delete this part when the subquestion is completed. */
+  /************************************************************************************************/
+    {
+      int array_of_subsizes[3] = {1, N, N};
+      int array_of_starts[3] = {0, 1, 1};
+      unpack_face(unpack[0], array_of_sizes, array_of_subsizes,
+                  array_of_starts);
+    }
+    {
+      int array_of_subsizes[3] = {1, N, N};
+      int array_of_starts[3] = {N + 1, 1, 1};
+      unpack_face(unpack[1], array_of_sizes, array_of_subsizes,
+                  array_of_starts);
+    }
+    {
+      int array_of_subsizes[3] = {N, 1, N};
+      int array_of_starts[3] = {1, 0, 1};
+      unpack_face(unpack[2], array_of_sizes, array_of_subsizes,
+                  array_of_starts);
+    }
+    {
+      int array_of_subsizes[3] = {N, 1, N};
+      int array_of_starts[3] = {1, N + 1, 1};
+      unpack_face(unpack[3], array_of_sizes, array_of_subsizes,
+                  array_of_starts);
+    }
+    {
+      int array_of_subsizes[3] = {N, N, 1};
+      int array_of_starts[3] = {1, 1, 0};
+      unpack_face(unpack[4], array_of_sizes, array_of_subsizes,
+                  array_of_starts);
+    }