Solution for hw3 added

hw3/solution/Makefile

+CXX=mpicxx  #h5pcc
+CXXFLAGS = -Wpedantic -Wall -Wextra -std=c++11 -lstdc++ -g -O3 -fopenmp 
+
+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_mpi:
+	make clean ; make ; mpirun -n 64 ./main 512 4 0.25
+
+run_hybrid_v1:
+	export OMP_NUM_THREADS=8 ; make clean ; make ; mpirun -n 8 --map-by ppr:3:node ./main 512 2 0.25
+
+run_hybrid_v2:
+	export OMP_NUM_THREADS=8 ; make clean ; make ; mpirun -n 8 --map-by ppr:2:node ./main 512 2 0.25
+
+develop:
+	export OMP_NUM_THREADS=8 ; make clean ; make ; mpirun -n 8 ./main 256 2 0.1
+
+
+.PHONY: all clean run_mpi run_hybrid_v1 run_hybrid_v2 develop

hw3/solution/main.cpp

+#include "wave.h"
+
+int main(int argc, char **argv) {
+  int provided;
+  MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &provided);
+
+  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 = std::stoi(argv[1]);
+  int procs_per_dim = std::stoi(argv[2]);
+  double t_end = std::stof(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);
+  }
+
+  int nthreads = omp_get_max_threads();
+  int threads_per_dim = pow(nthreads, 1.0 / 3.0) + 0.5;
+  if (nthreads != threads_per_dim * threads_per_dim * threads_per_dim) {
+    if (rank == 0)
+      std::cout
+          << " Number of OPENMP threads must be a cubic number. Aborting... \n";
+    int err = 4;
+    MPI_Abort(MPI_COMM_WORLD, err);
+  }
+
+  WaveEquation Simulation(points, procs_per_dim);
+
+  Simulation.run(t_end);
+
+  MPI_Finalize();
+  return 0;
+}

hw3/solution/wave.h

+#pragma once
+
+#include <cassert>
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+#include <mpi.h>
+#include <omp.h>
+#include <sstream>
+#include <string>
+#include <vector>
+
+//#define USE_HDF5 // comment this line if you do not want hdf5 file output
+                 // (suitable for visualization with Paraview)
+
+#ifdef USE_HDF5
+#include <hdf5.h>
+#endif
+
+#define L 1.0 // domain size (unit cube) with Ntot x Ntot x Ntot grid points
+
+struct WaveEquation {
+  int N;     // grid points per direction for this rank
+  double h;  // grid spacing (dx = dy = dz = h)
+  double dt; // timestep
+  double t;  // current time
+  double *u; // solution vector
+  double *u_old;
+  double *u_new;
+
+  int Ntot;          // total grid points per direction
+  int procs_per_dim; // ranks per direction
+
+  int nthreads;
+  int threads_per_dim;
+
+  int size; // =total ranks = procs_per_dim*procs_per_dim*procs_per_dim
+  int rank;
+
+  int rank_plus[3];  // neighboring ranks
+  int rank_minus[3]; //
+
+  MPI_Comm cart_comm;
+
+  double origin[3]; // physical space (x0,x1,x2) coordinates of 1st grid point
+                    // of this rank
+  int coords[3];    // index space coordinates of this rank
+
+  double aux;
+
+  ~WaveEquation();
+  WaveEquation(int a_N, int a_procs_per_dim);
+  void FindCoordinates();
+  double Initial_Condition(double x0, double x1, double x2);
+  void UpdateGridPoint(int i0, int i1, int i2);
+  void Print(int kt);
+  void pack_face(double *pack, int array_of_sizes[3], int array_of_subsizes[3],
+                 int array_of_starts[3]);
+  void unpack_face(double *pack, int array_of_sizes[3],
+                   int array_of_subsizes[3], int array_of_starts[3]);
+  void run(double t_end);
+
+  double **pack;
+  double **unpack;
+
+  void pack_all();
+  void unpack_all();
+
+  void thread_coordinates(int tid, int threads_per_dim, int &ti0, int &ti1,
+                          int &ti2);
+};