add mpi support

5 files changed, 91 insertions, 50 deletions

diff --git a/new/agent.cpp b/new/agent.cpp
index fc8c1f1..97d2e69 100644
--- a/new/agent.cpp
+++ b/new/agent.cpp
@@ -152,7 +152,6 @@ void Agent::advance_surface(Mesh* mesh, size_t step)
   // cell, it must be advected out of the domain
   if (distance_norm > 10.0 * location->get_maximum_length())
     domain[step] = OUTSIDE_FLAG;
-
   
   } else {
     // infiltration detected
diff --git a/new/header.hpp b/new/header.hpp
index 9acc489..cff75d2 100644
--- a/new/header.hpp
+++ b/new/header.hpp
@@ -10,6 +10,7 @@
 #include <algorithm> // min
 #include <float.h>   // DBL_MAX
 #include <math.h>    // fabs, sqrt
+#include <mpi.h>
 
 #include <H5Cpp.h>
 #include <H5File.h>
diff --git a/new/main.cpp b/new/main.cpp
index 5363f89..1e45aac 100644
--- a/new/main.cpp
+++ b/new/main.cpp
@@ -6,41 +6,59 @@
 int main(int argc, char* argv[])
 {
 
+  int rank;
+  int size;
+
+  MPI_Init(&argc, &argv);
+  MPI_Comm_size(MPI_COMM_WORLD, &size);
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
   Mesh mesh = Mesh();
   
   H5::H5File fp_subsurface = H5::H5File(SUBSURFACE_MESH_FILE, H5F_ACC_RDONLY);
   H5::H5File dp_subsurface = H5::H5File(SUBSURFACE_DATA_FILE, H5F_ACC_RDONLY);
   H5::H5File fp_surface = H5::H5File(SURFACE_MESH_FILE, H5F_ACC_RDONLY);
   H5::H5File dp_surface = H5::H5File(SURFACE_DATA_FILE, H5F_ACC_RDONLY);
-  
+
+  if (rank == 0)
+    std::cout << "[~>] Building the subsurface mesh" << std::endl;
+
   mesh.build_subsurface_mesh(fp_subsurface);
   fp_subsurface.close();
 
+  if (rank == 0) {
+    std::cout << "[~~>] Number of cells: " << mesh.get_n_cells(SUBSURFACE_FLAG) << std::endl;
+    std::cout << "[~>] Building the surface mesh (includes coupling with subsurface and may take a while)" << std::endl;
+  }
+  
   mesh.build_surface_mesh(fp_surface, dp_surface);
   fp_surface.close();
 
-  FILE *mesh_fp = fopen("mesh.csv", "w");
-  fprintf(mesh_fp, "x (m), y (m), z (m), index\n");
+  if (rank == 0) {
+
+    int len_surf = mesh.get_n_cells(SURFACE_FLAG);
+
+    std::cout << "[~~>] Number of cells: " << len_surf << std::endl;
+    
+    FILE *mesh_fp = fopen("mesh.csv", "w");
+    fprintf(mesh_fp, "x (m), y (m), z (m), index\n");
+    Cell *c;
 
-  Cell *c;
-  int len_surf = mesh.get_n_cells(SURFACE_FLAG);
-  for (int n = 0; n < len_surf; n ++) {
-    std::cout << n << "/" << len_surf << " ";
+    
+    for (int n = 0; n < len_surf; n ++) {
     c = mesh.get_cell(n, SURFACE_FLAG);
-    std::cout << "read in c ";
     fprintf(mesh_fp, "%f, %f, %f, %d\n",
 	    c->get_coordinate(0),
 	    c->get_coordinate(1),
 	    c->get_coordinate(2), n);
-    std::cout << "done reading." << std::endl;
-  }
-
-  std::cout << "closing file..." << std::endl;
+    }
   
-  fclose(mesh_fp);
+    fclose(mesh_fp);
 
-  std::cout << "reading the environment" << std::endl;
- 
+    std::cout << "[~>] Reading the subsurface environment" << std::endl;
+  
+  }
+  
   mesh.build_environment(dp_subsurface, SUBSURFACE_FLAG);
   dp_subsurface.close();
 
@@ -48,29 +66,36 @@ int main(int argc, char* argv[])
   dp_subsurface.close();
 
   // Place agents in the domain and advance N_ITERATION steps
-  std::vector<Agent> agents(N_AGENT);
 
   // Bootstrap 
-  srand(SEED);
+  srand(SEED + rank);
   int len = mesh.get_n_cells(SUBSURFACE_FLAG);
   int index;
 
   std::cout << "[~>] Releasing " << N_AGENT << " agents" << std::endl;
+
+  size_t partition = N_AGENT / size;
+  size_t low_bound = rank * partition;
+  size_t high_bound = (rank + 1) * partition;
+  if (rank == size - 1) // top processor picks up extras
+    high_bound = N_AGENT;
   
-  for (size_t i = 0; i < N_AGENT; i ++) {
+  for (size_t i = low_bound; i < high_bound; i ++) {
 
+    Agent agent;
+    
     index = rand()%len;
-    agents[i].set_cell_index(index);
-    agents[i].set_path(mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(0),
-		       mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(1),
-		       mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(2), 0);
+    agent.set_cell_index(index);
+    agent.set_path(mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(0),
+		   mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(1),
+		   mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(2), 0);
   
-    agents[i].set_domain(SUBSURFACE_FLAG, 0);
-    agents[i].set_domain(SUBSURFACE_FLAG, 1);
+    agent.set_domain(SUBSURFACE_FLAG, 0);
+    agent.set_domain(SUBSURFACE_FLAG, 1);
 
-    agents[i].set_path(mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(0),
-		       mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(1),
-		       mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(2), 1);
+    agent.set_path(mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(0),
+		   mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(1),
+		   mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(2), 1);
 
     size_t step;
 
@@ -79,9 +104,9 @@ int main(int argc, char* argv[])
     // move agent
     for (step = 1; step < N_ITERATION; step ++) {
       if (step % (N_ITERATION / 5) == 0 && i % (N_AGENT / 5) == 0) {
-	std::cout << "[~~~>] agent[" << i << "] step = " << step << " dt = " << agents[i].get_timestepsize() << " s" << std::endl;
+	std::cout << "[~~~>] agent[" << i << "] step = " << step << " dt = " << agent.get_timestepsize() << " s" << std::endl;
       }
-      agents[i].advance(&mesh, step);
+      agent.advance(&mesh, step);
     }
 
     // write out results
@@ -91,17 +116,18 @@ int main(int argc, char* argv[])
     fprintf(result_fp, "time (s), x (m), y (m), z (m), domain\n");
     for (step = 0; step < N_ITERATION; step ++) {
       fprintf(result_fp, "%f, %f, %f, %f, %d\n",
-	      agents[i].get_time(step),
-	      agents[i].get_path(step, 0),
-	      agents[i].get_path(step, 1),
-	      agents[i].get_path(step, 2),
-	      agents[i].get_domain(step));
+	      agent.get_time(step),
+	      agent.get_path(step, 0),
+	      agent.get_path(step, 1),
+	      agent.get_path(step, 2),
+	      agent.get_domain(step));
     }
-
     fclose(result_fp);
       
   }
+
+  MPI_Finalize();
   
-  return EXIT_SUCCESS;
+  return MPI_SUCCESS;
 
 }
diff --git a/new/makefile b/new/makefile
index 3d0ec9d..4a91e28 100644
--- a/new/makefile
+++ b/new/makefile
@@ -1,6 +1,7 @@
 PRG=agent-provocateur
 
-CC=g++-11
+#CC=g++-11
+CC=mpic++
 CXXFLAGS=-Wall -Wextra -pedantic -O0 -fopenmp
 LFLAGS=-I/opt/hdf5/include -L/opt/hdf5/lib
 SRC=cell.cpp mesh.cpp main.cpp agent.cpp
diff --git a/new/mesh.cpp b/new/mesh.cpp
index 6ae6c22..8112ede 100644
--- a/new/mesh.cpp
+++ b/new/mesh.cpp
@@ -9,7 +9,6 @@
 
 void Mesh::build_surface_mesh(H5::H5File fp, H5::H5File ep)
 {
-  std::cout << "[~>] Building the surface mesh" << std::endl;
 
   get_mesh_metadata(fp, SURFACE_NCOL);
 
@@ -21,8 +20,6 @@ void Mesh::build_surface_mesh(H5::H5File fp, H5::H5File ep)
 
 void Mesh::build_subsurface_mesh(H5::H5File fp)
 {
-  
-  std::cout << "[~>] Building the subsurface mesh" << std::endl;
 
   get_mesh_metadata(fp, SUBSURFACE_NCOL);
 
@@ -41,11 +38,9 @@ void Mesh::build_environment(H5::H5File ep, int flag)
   // Read in velocities in x direction
     
   if (flag == SUBSURFACE_FLAG) {
-    std::cout << "[~>] Building the subsurface environment" << std::endl;
     groupname = "darcy_velocity.cell.0";
     datasetname = SUBSURFACE_DATA_GROUP_NAME_H5;
   } else if (flag == SURFACE_FLAG) {
-    std::cout << "[~>] Building the surface environment" << std::endl;
     groupname = "surface-velocity.cell.0";
     datasetname = SURFACE_DATA_GROUP_NAME_H5;
   }
@@ -85,17 +80,21 @@ void Mesh::build_environment(H5::H5File ep, int flag)
     // Rescale the exchange flux from mole to m/s
     for (int i = 0; i < n_cells_surface; i ++) {
       surface_cells[i].set_velocity(2, surface_cells[i].get_velocity(2) * 1.8e-5 * AREADIV);
+#if VERBOSE
       if (i % (n_cells_surface / MAXPRINT) == 0)
   	printf("[~~~~>] Surface velocity in cell[%d] = { %e, %e, %e } m/s\n", i,
 	       surface_cells[i].get_velocity(0),
 	       surface_cells[i].get_velocity(1),
 	       surface_cells[i].get_velocity(2));
+#endif
     }
 
     // Read in water depth
     read_water_depth(ep);
     
-  } else if (flag == SUBSURFACE_FLAG) {
+  }
+#if VERBOSE
+  else if (flag == SUBSURFACE_FLAG) {
     
     for (int i = 0; i < n_cells_subsurface; i ++) {
       if (i % (n_cells_subsurface / MAXPRINT) == 0)
@@ -106,6 +105,7 @@ void Mesh::build_environment(H5::H5File ep, int flag)
     }
     
   }
+#endif
   
 }
 
@@ -180,8 +180,10 @@ void Mesh::handle_surface_node_data(H5::H5File fp)
   hsize_t *dims_n = new hsize_t[rank];
   dataspace_n.getSimpleExtentDims(dims_n);
 
+#if VERBOSE
   for (int i = 0; i < rank; i ++)
     std::cout << "[~~~>] Dimension in " << i << " of nodes: " << dims_n[i] << std::endl;
+#endif
 
   double *data_n = new double[dims_n[0]*dims_n[1]];
   dataset_n.read(data_n, H5::PredType::IEEE_F64LE);
@@ -208,8 +210,6 @@ void Mesh::build_surface_cells(H5::H5File ep)
 
   for (int i = 0; i < n_cells_surface; i ++) {
 
-    //std::cout << "build cell " << i << std::endl;
-
     int n1 = element_node_indices[i][1];
     int n2 = element_node_indices[i][2];
     int n3 = element_node_indices[i][3];
@@ -282,11 +282,13 @@ void Mesh::build_surface_cells(H5::H5File ep)
     surface_cells[i].set_coupled_neighbor(coupled_neighbor);
     subsurface_cells[coupled_neighbor].set_coupled_neighbor(i);
 
+#if VERBOSE
     if (i%(n_cells_surface/5)==0) {
       std::cout << "[~~~>] Coupling surface and subsurface meshes" << std::endl;
       std::cout << "[~~~~>] Subsurface neighbor of cell " << i << " is cell " << surface_cells[i].get_coupled_neighbor() << std::endl;
       std::cout << "[~~~~>] Surface neighbor of cell " << coupled_neighbor << " is cell " << subsurface_cells[coupled_neighbor].get_coupled_neighbor() << std::endl;
     }
+#endif
     
   }
 
@@ -304,9 +306,11 @@ void Mesh::handle_subsurface_node_data(H5::H5File fp)
   hsize_t *dims_n = new hsize_t[rank];
   dataspace_n.getSimpleExtentDims(dims_n);
 
+#if VERBOSE
   for (int i = 0; i < rank; i ++)
     std::cout << "[~~~>] Dimension in " << i << " of nodes: " << dims_n[i] << std::endl;
-
+#endif
+  
   double *data_n = new double[dims_n[0]*dims_n[1]];
   dataset_n.read(data_n, H5::PredType::IEEE_F64LE);
 
@@ -415,8 +419,10 @@ void Mesh::get_mesh_metadata(H5::H5File fp, int flag)
   hsize_t *dims_e = new hsize_t[rank];
   dataspace_e.getSimpleExtentDims(dims_e);
 
+#if VERBOSE
   for (int i = 0; i < rank; i ++)
       std::cout << "[~~~>] Dimension in " << i << " of element node map: " << dims_e[i] << std::endl;
+#endif
 
   int *data_e = new int[dims_e[0]*dims_e[1]];
   dataset_e.read(data_e, H5::PredType::STD_I32LE);
@@ -439,8 +445,6 @@ void Mesh::get_mesh_metadata(H5::H5File fp, int flag)
     n_cells_surface = ny_e;
     surface_cells.resize(ny_e);
   }
-
-  std::cout << "[~~~>] Number of cells: " << ny_e << std::endl;
   
 }
 
@@ -457,8 +461,10 @@ void Mesh::read_velocity(H5::H5File fp, std::string groupname, std::string datas
   hsize_t *dims = new hsize_t[rank];
   dataspace.getSimpleExtentDims(dims);
 
+#if VERBOSE
   for (int i = 0; i < rank; i ++)
     std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl;
+#endif
   
   double *data = new double[dims[0]*dims[1]];
   dataset.read(data, H5::PredType::IEEE_F64LE);
@@ -486,18 +492,22 @@ void Mesh::read_water_depth(H5::H5File fp)
   hsize_t *dims = new hsize_t[rank];
   dataspace.getSimpleExtentDims(dims);
 
+#if VERBOSE
   for (int i = 0; i < rank; i ++)
     std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl;
+#endif
   
   double *data = new double[dims[0]*dims[1]];
   dataset.read(data, H5::PredType::IEEE_F64LE);
 
   for (int i = 0; i < n_cells_surface; i ++) {
     surface_cells[i].set_water_depth(data[i]);
+#if VERBOSE
     if (i%(n_cells_surface/5)==0) {
       printf("[~~~~>] Surface ponded depth in cell[%d] =  %e m\n", i,
 	     surface_cells[i].get_water_depth());
     }
+#endif
   }
   
 }
@@ -515,18 +525,22 @@ void Mesh::read_elevation(H5::H5File fp)
   hsize_t *dims = new hsize_t[rank];
   dataspace.getSimpleExtentDims(dims);
 
+#if VERBOSE  
   for (int i = 0; i < rank; i ++)
     std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl;
+#endif
   
   double *data = new double[dims[0]*dims[1]];
   dataset.read(data, H5::PredType::IEEE_F64LE);
-
+  
   for (int i = 0; i < n_cells_surface; i ++) {
     surface_cells[i].set_coordinate(2, data[i]);
+#if VERBOSE
     if (i%(n_cells_surface/5)==0) {
       printf("[~~~~>] Surface elevation in cell[%d] =  %e m\n", i,
 	     surface_cells[i].get_coordinate(2));
     }
+#endif
   }
   
 }