From 2833f2e67ca27ae6c0b6d1062a1de1be36701e75 Mon Sep 17 00:00:00 2001 From: Ilhan Özgen Xian Date: Fri, 21 May 2021 20:30:42 -0700 Subject: add mpi support --- new/agent.cpp | 1 - new/header.hpp | 1 + new/main.cpp | 98 +++++++++++++++++++++++++++++++++++++--------------------- new/makefile | 3 +- new/mesh.cpp | 38 ++++++++++++++++------- 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 // min #include // DBL_MAX #include // fabs, sqrt +#include #include #include 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 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 } } -- cgit v1.2.3-13-gbd6f