diff options
author | Ilhan Özgen Xian <iozgen@lbl.gov> | 2021-05-25 11:17:08 -0700 |
---|---|---|
committer | Ilhan Özgen Xian <iozgen@lbl.gov> | 2021-05-25 11:17:08 -0700 |
commit | df6adc77e0cbe568832f57212929a9acfb568a4c (patch) | |
tree | 0ab3aa5ed3468e011c344ccf655358b806cf0bd2 | |
parent | 2833f2e67ca27ae6c0b6d1062a1de1be36701e75 (diff) |
clean up
-rw-r--r-- | agent.cpp | 323 | ||||
-rw-r--r-- | agent.hpp | 59 | ||||
-rw-r--r-- | cell.cpp (renamed from new/cell.cpp) | 0 | ||||
-rw-r--r-- | cell.hpp (renamed from new/cell.hpp) | 0 | ||||
-rw-r--r-- | header.hpp | 69 | ||||
-rw-r--r-- | main.cpp | 203 | ||||
-rw-r--r-- | makefile | 10 | ||||
-rwxr-xr-x | merge-results.sh (renamed from new/merge-results.sh) | 0 | ||||
-rw-r--r-- | mesh.cpp | 1334 | ||||
-rw-r--r-- | mesh.hpp | 120 | ||||
-rw-r--r-- | new/agent.cpp | 245 | ||||
-rw-r--r-- | new/agent.hpp | 41 | ||||
-rw-r--r-- | new/header.hpp | 56 | ||||
-rw-r--r-- | new/main.cpp | 133 | ||||
-rw-r--r-- | new/makefile | 33 | ||||
-rw-r--r-- | new/mesh.cpp | 546 | ||||
-rw-r--r-- | new/mesh.hpp | 59 | ||||
-rwxr-xr-x | test.cpp | 140 | ||||
-rw-r--r-- | tools.cpp | 11 | ||||
-rw-r--r-- | tools.hpp | 8 |
20 files changed, 758 insertions, 2632 deletions
@@ -1,196 +1,245 @@ -// -*- mode: c++; c-default-style: "linux" -*- +// -*- mode: c++ -*- -#include "header.hpp" +#include "agent.hpp" -HydroAgent::HydroAgent() +// ------------------------------------------------------------------------ +// -- PUBLIC METHODS +// ------------------------------------------------------------------------ + +void Agent::advance(Mesh* mesh, size_t step) { - Cell *location; + + switch(domain[step-1]) { + case SURFACE_FLAG: + advance_surface(mesh,step); + break; + case SUBSURFACE_FLAG: + advance_subsurface(mesh,step); + break; + default: + set_path(path[step-1][0], path[step-1][1], path[step-1][2], step); + time[step] = time[step-1]; + domain[step] = OUTSIDE_FLAG; + break; + } - probability_i = 0.0; - probability_a = 0.0; } -Cell* HydroAgent::get_location() +// ------------------------------------------------------------------------ +// GETTERS AND SETTERS +// ------------------------------------------------------------------------ + +Agent::Agent() { - return location; + time.resize(N_ITERATION); + domain.resize(N_ITERATION); + path.resize(N_ITERATION); + for (size_t i = 0; i < N_ITERATION; i ++) + path[i].resize(DIM); } -int HydroAgent::set_location(Cell *cell) +void Agent::set_cell_index(int value) { - location = cell; - return 0; + cell_index = value; } -int HydroAgent::set_path(double xcrd, double ycrd, double zcrd, size_t step) +void Agent::set_domain(int value, size_t step) { - x[step] = xcrd; - y[step] = ycrd; - z[step] = zcrd; - return 0; + domain[step] = value; } -int HydroAgent::advance(Mesh *mesh, size_t step) +double Agent::get_time(size_t step) { - int flag = 0; - - if (domain[step-1] == SURFACE_FLAG) { - flag = advance_surface(mesh, step); - } else if (domain[step-1] == SUBSURFACE_FLAG) { - flag = advance_subsurface(mesh, step); - } else if (domain[step-1] == OUTOFBOUNDS_FLAG) { - flag = 1; - x[step] = x[step-1]; - y[step] = y[step-1]; - set_domain(OUTOFBOUNDS_FLAG, step); - #if VERBOSE - printf("[ok] agent is out of bounds. do not advect.\n"); - #endif - } + return time[step]; +} - return flag; +void Agent::set_path(double x, double y, double z, size_t step) +{ + path[step][0] = x; + path[step][1] = y; + path[step][2] = z; } -int HydroAgent::advance_surface(Mesh *mesh, size_t step) +double Agent::get_timestepsize() { + return time_step_size; +} + +double Agent::get_path(size_t step, size_t index) +{ + return path[step][index]; +} + +int Agent::get_domain(size_t step) +{ + return domain[step]; +} + +// ------------------------------------------------------------------------ +// -- PRIVATE METHODS +// ------------------------------------------------------------------------ - double vx = location->get_v(0); - double vy = location->get_v(1); +void Agent::advance_surface(Mesh* mesh, size_t step) +{ + Cell *location = mesh->get_cell(cell_index, SURFACE_FLAG); - double x_agent = x[step-1] + DT_SURFACE * vx; - double y_agent = y[step-1] + DT_SURFACE * vy; + std::vector<double> velocity(DIM); + std::vector<double> coordinates(DIM-1); + std::vector<double> coordinates_neighbor(DIM-1); + std::vector<double> distance(DIM-1); + + domain[step] = SURFACE_FLAG; // set domain for this step, may change later if agent infiltrates - x[step] = x_agent; - y[step] = y_agent; + for (size_t i = 0; i < DIM; i ++) + velocity[i] = location->get_velocity(i); - double x_centroid = location->get_x(); - double y_centroid = location->get_y(); + // CFL criteria + double length = location->get_characteristic_length(); + double velocity_norm = sqrt(velocity[0]*velocity[0] + velocity[1]*velocity[1]); - double dx = x_centroid - x_agent; - double dy = y_centroid - y_agent; + if (velocity_norm > EPS) { + time_step_size = CFL * length / velocity_norm; + } else { + time_step_size = MAXTIMESTEP; + } - double distance = sqrt(dx * dx + dy * dy); + time[step] = time[step-1] + time_step_size; - int cell_index = -1; + // check if agent infiltrates + double h = location->get_water_depth(); + int flag = INFILTRATE_FALSE; - for (int n = 0; n < N_EDGES; n ++) { + if (h > EPS) { // only wet cell can infiltrate + double qe = std::min(velocity[2], 0.0); + double pi = fabs(qe) * time_step_size / h; + double trial = (double) (rand()%100) / 100.0; + if (trial < pi) + flag = INFILTRATE_TRUE; + } - hsize_t neighbor_index = location->get_neighbor(n); + if (flag == INFILTRATE_FALSE) { - double x_neighbor = mesh->get_cells()[neighbor_index].get_x(); - double y_neighbor = mesh->get_cells()[neighbor_index].get_y(); + // advect the agent horizontally, which means the last dimension + // (z-direction) is omitted + for (size_t i = 0; i < DIM-1; i ++) { + coordinates[i] = path[step-1][i] + time_step_size * velocity[i]; + distance[i] = location->get_coordinate(i) - coordinates[i]; + } - dx = x_neighbor - x_agent; - dy = y_neighbor - y_agent; + // horizontal distance to the centroid of current cell + double distance_norm = sqrt(distance[0]*distance[0]+distance[1]*distance[1]); - double _distance = sqrt(dx * dx + dy * dy); + // check if cell has changed + size_t n_cells = mesh->get_n_cells(SURFACE_FLAG); + for (size_t i = 0; i < n_cells; i ++) { - if (distance > _distance) { - cell_index = n; - distance = _distance; + for (size_t j = 0; j < DIM-1; j ++) { + coordinates_neighbor[j] = mesh->get_cell(i, SURFACE_FLAG)->get_coordinate(j); + distance[j] = coordinates_neighbor[j] - coordinates[j]; + } + + double distance_norm_ = sqrt(distance[0]*distance[0]+distance[1]*distance[1]); + + if (distance_norm > distance_norm_) { + cell_index = i; + distance_norm = distance_norm_; } } - if (cell_index >= 0) { - location = &(mesh->get_cells()[location->get_neighbor(cell_index)]); + // check if out of the domain: if the agent travelled thrice the + // maximum length of the cell and still can't find a neighboring + // 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 + domain[step] = SUBSURFACE_FLAG; + coordinates[0] = path[step-1][0]; + coordinates[1] = path[step-1][1]; + cell_index = mesh->get_cell(cell_index, SURFACE_FLAG)->get_coupled_neighbor(); } - if ((x_agent > mesh->get_xmax()) || (x_agent < mesh->get_xmin()) || - (y_agent > mesh->get_ymax()) || (y_agent < mesh->get_ymin())) { - domain[step] = OUTOFBOUNDS_FLAG; - } else { - domain[step] = SURFACE_FLAG; - } + set_path(coordinates[0], coordinates[1], location->get_coordinate(2), step); - return 1; } -int HydroAgent::advance_subsurface(Mesh *mesh, size_t step) -{ - double vx = location->get_v(0); - double vy = location->get_v(1); - double vz = location->get_v(2); +void Agent::advance_subsurface(Mesh* mesh, size_t step) +{ + Cell *location = mesh->get_cell(cell_index, SUBSURFACE_FLAG); - double x_agent = x[step-1] + DT_SUBSURFACE * vx; - double y_agent = y[step-1] + DT_SUBSURFACE * vy; - double z_agent = z[step-1] + DT_SUBSURFACE * vz; + std::vector<double> velocity(DIM); + std::vector<double> coordinates(DIM); + std::vector<double> coordinates_neighbor(DIM); + std::vector<double> distance(DIM); - x[step] = x_agent; - y[step] = y_agent; - z[step] = z_agent; + domain[step] = SUBSURFACE_FLAG; // set domain for this step, may change later if agent infiltrates + + for (size_t i = 0; i < DIM; i ++) + velocity[i] = location->get_velocity(i); - double x_centroid = location->get_x(); - double y_centroid = location->get_y(); - double z_centroid = location->get_z(); + // CFL criteria + double length = location->get_characteristic_length(); + double velocity_norm = sqrt(velocity[0] * velocity[0] + velocity[1] * velocity[1] + velocity[2] * velocity[2]); - double dx = x_centroid - x_agent; - double dy = y_centroid - y_agent; - double dz = z_centroid - z_agent; + if (velocity_norm > EPS) { + time_step_size = CFL * length / velocity_norm; + } else { + time_step_size = MAXTIMESTEP; + } - double distance = sqrt(dx * dx + dy * dy + dz * dz); + time[step] = time[step-1] + time_step_size; - int cell_index = -1; + // advect the agent + for (size_t i = 0; i < DIM; i ++) { + coordinates[i] = path[step-1][i] + time_step_size * velocity[i]; + distance[i] = location->get_coordinate(i) - coordinates[i]; + } + + set_path(coordinates[0], coordinates[1], coordinates[2], step); - for (int n = 0; n < SUBSURFACE_N_NEIGHBORS; n ++) { + // distance to the centroid of current cell + double distance_norm = sqrt(distance[0] * distance[0] + distance[1] * distance[1] + distance[2] * distance[2]); + + // check if cell has changed + size_t n_cells = mesh->get_n_cells(SUBSURFACE_FLAG); + for (size_t i = 0; i < n_cells; i ++) { - hsize_t neighbor_index = location->get_neighbor(n); + for (size_t j = 0; j < DIM; j ++) { + coordinates_neighbor[j] = mesh->get_cell(i, SUBSURFACE_FLAG)->get_coordinate(j); + distance[j] = coordinates_neighbor[j] - coordinates[j]; + } - double x_neighbor = mesh->get_cells()[neighbor_index].get_x(); - double y_neighbor = mesh->get_cells()[neighbor_index].get_y(); - double z_neighbor = mesh->get_cells()[neighbor_index].get_z(); + double distance_norm_ = sqrt(distance[0] * distance[0] + distance[1] * distance[1] + distance[2] * distance[2]); - dx = x_neighbor - x_agent; - dy = y_neighbor - y_agent; - dz = z_neighbor - z_agent; + if (distance_norm > distance_norm_) { + cell_index = i; + distance_norm = distance_norm_; + } + + } - double _distance = sqrt(dx * dx + dy * dy + dz * dz); + // check if out of the domain: if the agent travelled thrice the + // maximum length of the cell and still can't find a neighboring + // cell, it must be advected out of the domain + if (distance_norm > 10.0 * location->get_maximum_length()) { + domain[step] = OUTSIDE_FLAG; + } - if (distance > _distance) { - cell_index = n; - distance = _distance; + // check if exfiltrates + int neighbor_index = mesh->get_cell(cell_index, SUBSURFACE_FLAG)->get_coupled_neighbor(); + double z_surface; + if (neighbor_index > -1) { + // border to surface detected + z_surface = mesh->get_cell(neighbor_index, SURFACE_FLAG)->get_coordinate(2); + if (coordinates[2] >= z_surface) { + // exfiltration detected + cell_index = (int) neighbor_index; + domain[step] = SURFACE_FLAG; + set_path(coordinates[0], coordinates[1], z_surface, step); } - } - if (cell_index >= 0) { - location = &(mesh->get_cells()[location->get_neighbor(cell_index)]); - } - - if ((x_agent > mesh->get_xmax()) || (x_agent < mesh->get_xmin()) || - (y_agent > mesh->get_ymax()) || (y_agent < mesh->get_ymin()) || - (z_agent > mesh->get_zmax()) || (z_agent < mesh->get_zmin())) { - domain[step] = OUTOFBOUNDS_FLAG; - } else { - domain[step] = SURFACE_FLAG; - } - - return 2; -} - -int HydroAgent::set_domain(int domainflag, size_t step) -{ - domain[step] = domainflag; - - return 0; -} - -int HydroAgent::get_domain(size_t step) -{ - return domain[step]; -} - -double* HydroAgent::get_x() -{ - return &x[0]; -} - -double* HydroAgent::get_y() -{ - return &y[0]; -} - -double* HydroAgent::get_z() -{ - return &z[0]; } @@ -1,55 +1,40 @@ -// -*- mode: c++; c-default-style: "linux" -*- +// -*- mode: c++ -*- #ifndef AGENT_HPP #define AGENT_HPP #include "header.hpp" +#include "cell.hpp" +#include "mesh.hpp" -// agent.hpp -// ========= -// HydroAgent to be used in agent-based modeling to track flow paths -// on a given flow field. -// At each time step, the HydroAgent needs to decide whether to -// infiltrate or to get advected following the flow. This is based -// on computed probabilities, see (Reaney, 2008). -// First, check if the agent infiltrates/exfiltrates. If not, check -// if it advects. If yes, advect according to flow velocity. - -class HydroAgent +class Agent { public: - - HydroAgent(); - - Cell* get_location(); - int set_location(Cell *cell); - int set_path(double x, double y, double z, size_t step); - int set_domain(int domainflag, size_t step); - int get_domain(size_t step); - - int advance(Mesh *mesh, size_t step); - int advance_surface(Mesh *mesh, size_t step); - int advance_subsurface(Mesh *mesh, size_t step); - double *get_x(); - double *get_y(); - double *get_z(); + Agent(); -protected: + void set_cell_index(int value); + void set_domain(int value, size_t step); + void set_path(double x, double y, double z, size_t step); + double get_path(size_t step, size_t index); + int get_domain(size_t step); + double get_timestepsize(); + double get_time(size_t step); + + void advance(Mesh* mesh, size_t step); private: - Cell *location; + int cell_index; + double time_step_size; - int domain[N_ITERATION]; - - double x[N_ITERATION]; - double y[N_ITERATION]; - double z[N_ITERATION]; - - double probability_i; // infiltration probability - double probability_a; // advection probability + std::vector<int> domain; + std::vector<double> time; + std::vector<std::vector<double>> path; + + void advance_surface(Mesh *mesh, size_t step); + void advance_subsurface(Mesh *mesh, size_t step); }; @@ -1,49 +1,56 @@ -// -*- mode: c++; c-default-style: "linux" -*- +// -*- mode:c++ -*- -#ifndef HEADER_HPP +#ifndef CONSTANTS_HPP +#define CONSTANTS_HPP #include <iostream> +#include <vector> +#include <string> // string class +//#include <omp.h> +#include <algorithm> // min +#include <float.h> // DBL_MAX +#include <math.h> // fabs, sqrt +#include <mpi.h> #include <H5Cpp.h> #include <H5File.h> -#include "agent.hpp" -#include "mesh.hpp" -#include "tools.hpp" - -#include <float.h> -#include <math.h> -#include <vector> +#define MAXPRINT 5 +#define MAXTIMESTEP 1.0e5 +#define EPS 1.0e-10 -#define VERBOSE 1 -#define VERBOSE_HIGH 0 +#define N_AGENT 100 +#define N_ITERATION 2000 +#define SEED 108 +#define CFL 0.5 -#define N_EDGES 3 -#define SURFACE_N_EDGES 3 -#define SUBSURFACE_N_EDGES 9 -#define SUBSURFACE_N_NODES 6 -#define SUBSURFACE_N_NEIGHBORS 5 // 5 von Neumann neighbors, 11 if Moore neighborhood #define DIM 3 +#define N_NEIGHBORS 5 +#define N_VERTICES 6 +#define N_EDGES 9 -#define RADIUS 200.0 +#define SURFACE_NCOL 4 +#define SUBSURFACE_NCOL 7 -#define N_AGENT 300 -#define N_ITERATION 5000 -#define CFL 0.5 -#define DT_SURFACE 1800.0 -#define DT_SUBSURFACE 1800.0 +#define INFILTRATE_TRUE 1 +#define INFILTRATE_FALSE 2 -#define SUBSURFACE 1 +#define SURFACE_MESH_FILE "../test/visdump_surface_mesh.h5" +#define SURFACE_DATA_FILE "../test/visdump_surface_data.h5" +#define SUBSURFACE_MESH_FILE "../test/visdump_mesh.h5" +#define SUBSURFACE_DATA_FILE "../test/visdump_data.h5" +#define SURFACE_DATA_GROUP_NAME_H5 "11841" +#define SUBSURFACE_DATA_GROUP_NAME_H5 "11841" -#define SURFACE_FLAG 1 -#define SUBSURFACE_FLAG 2 -#define OUTOFBOUNDS_FLAG 3 +#define MESH_MIXED_ELEMENTS "/0/Mesh/MixedElements" +#define MESH_NODES "/0/Mesh/Nodes" -#define MESH_FILE_NAME "test/visdump_surface_mesh.h5" -#define DATA_FILE_NAME "test/visdump_surface_data.h5" -#define SUBSURFACE_MESH_FILE_NAME "test/visdump_mesh.h5" -#define SUBSURFACE_DATA_FILE_NAME "test/visdump_data.h5" +#define AREADIV 0.01 +#define FUZZ 0.1 +#define RADIUS 200.0 -#define DATA_GROUP_NAME_H5 "11841" +#define SUBSURFACE_FLAG 1 +#define SURFACE_FLAG 2 +#define OUTSIDE_FLAG 3 #endif @@ -1,138 +1,133 @@ -// -*- mode: c++; c-default-style: "linux" -*- +// -*- mode: c++ -*- -#include "header.hpp" +#include "mesh.hpp" +#include "agent.hpp" -int main(void) +int main(int argc, char* argv[]) { - printf("\n\n"); - printf(" / /^ agent provocateur\n"); - printf(" o--o =================\n"); - printf("\n\n"); - - printf("[OK] # of agents: %d, # of iterations: %d\n[OK] size of surface time step: %f s\n[OK] size of subsurface time step: %f s\n\n", N_AGENT, N_ITERATION, DT_SURFACE, DT_SUBSURFACE); - - H5::H5File fp = H5::H5File(MESH_FILE_NAME, H5F_ACC_RDONLY); - + 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(); - mesh.build_surface(fp); - - printf("[OK] # of cells: %d\n\n", mesh.get_n_cells()); - fp.close(); + 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 SUBSURFACE - H5::H5File subsurface_fp = H5::H5File(SUBSURFACE_MESH_FILE_NAME, H5F_ACC_RDONLY); + if (rank == 0) + std::cout << "[~>] Building the subsurface mesh" << std::endl; - Mesh subsurface_mesh = Mesh(); - subsurface_mesh.build_subsurface(subsurface_fp); + mesh.build_subsurface_mesh(fp_subsurface); + fp_subsurface.close(); - printf("[OK] # of cells: %d\n\n", subsurface_mesh.get_n_cells()); - - subsurface_fp.close(); -#endif - - H5::H5File dp = H5::H5File(DATA_FILE_NAME, H5F_ACC_RDONLY); - mesh.read_surface_data(dp); - dp.close(); - -#if SUBSURFACE - H5::H5File sdp = H5::H5File(SUBSURFACE_DATA_FILE_NAME, H5F_ACC_RDONLY); - subsurface_mesh.read_subsurface_data(sdp); - sdp.close(); -#endif + 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; + } - // Place agents in the domain and advance N_ITERATION steps - std::vector<HydroAgent> agents(N_AGENT); + mesh.build_surface_mesh(fp_surface, dp_surface); + fp_surface.close(); - srand(108); + if (rank == 0) { - int len = mesh.get_n_cells(); + int len_surf = mesh.get_n_cells(SURFACE_FLAG); - int index; + 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; - printf("[OK] inject agents\n"); - - for (size_t i = 0; i < N_AGENT; i ++) { - index = rand() % len; - - if (fabs(mesh.get_cells()[index].get_v(0)) > 0.0 || - fabs(mesh.get_cells()[index].get_v(1)) > 0.0) { - - agents[i].set_location(&(mesh.get_cells()[index])); - agents[i].set_path(agents[i].get_location()->get_x(), agents[i].get_location()->get_y(), 0.0, 0); - agents[i].set_domain(SURFACE_FLAG, 0); - // bootstrap - agents[i].set_domain(SURFACE_FLAG, 1); - agents[i].set_path(agents[i].get_location()->get_x(), agents[i].get_location()->get_y(), 0.0, 1); - + for (int n = 0; n < len_surf; n ++) { + c = mesh.get_cell(n, SURFACE_FLAG); + fprintf(mesh_fp, "%f, %f, %f, %d\n", + c->get_coordinate(0), + c->get_coordinate(1), + c->get_coordinate(2), n); } - + + fclose(mesh_fp); + + std::cout << "[~>] Reading the subsurface environment" << std::endl; + } + + mesh.build_environment(dp_subsurface, SUBSURFACE_FLAG); + dp_subsurface.close(); - for (size_t step = 1; step < N_ITERATION; step ++) { + mesh.build_environment(dp_surface, SURFACE_FLAG); + dp_subsurface.close(); - if (step % 10 == 0) { - printf("[OK] time step %08lu {\n", step); - } - - for (size_t i = 0; i < N_AGENT; i ++) { - - #if VERBOSE - if (step % 100 == 0) { - if (i < 10) { - printf("[ok] advance agent #%lu\n", i+1); - } else if (i == 10) { - printf("[ok] ...\n"); - } - } - #endif - - agents[i].advance(&mesh, step); - - } - - if (step % 10 == 0) { - printf("[OK] } done\n"); - } - - } + // Place agents in the domain and advance N_ITERATION steps - printf("[OK] simulation completed\n\n"); + // Bootstrap + srand(SEED + rank); + int len = mesh.get_n_cells(SUBSURFACE_FLAG); + int index; - printf("[OK] write out agent data {\n"); + std::cout << "[~>] Releasing " << N_AGENT << " agents" << std::endl; - char output_fname[1024]; + 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; - // write out the results - for (size_t i = 0; i < N_AGENT; i ++) { + for (size_t i = low_bound; i < high_bound; i ++) { - sprintf(output_fname, "agent-%06lu.dat", i); - FILE *result_fp = fopen(output_fname, "w"); + Agent agent; + + index = rand()%len; + 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); + + agent.set_domain(SUBSURFACE_FLAG, 0); + agent.set_domain(SUBSURFACE_FLAG, 1); - #if VERBOSE - if (i < 10) { - printf("[ok] write out %s\n", output_fname); - } else if (i == 10) { - printf("[ok] ...\n"); - } - #endif + 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; - double* x = agents[i].get_x(); - double* y = agents[i].get_y(); + char output_fname[1024]; - for (size_t step = 0; step < N_ITERATION; step ++) { - fprintf(result_fp, "%f %f %f\n", step * DT_SURFACE, x[step], y[step]); + // 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 = " << agent.get_timestepsize() << " s" << std::endl; + } + agent.advance(&mesh, step); } + // write out results + sprintf(output_fname, "results/agent-%06lu.csv", i); + FILE *result_fp = fopen(output_fname, "w"); + + 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", + 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); - + } - printf("[OK] } done\n"); - printf("[OK] finish\n\n"); - - return 0; + MPI_Finalize(); + return MPI_SUCCESS; + } @@ -1,10 +1,10 @@ -PRG=agent +PRG=agent-provocateur -CC=g++-11 -#CXXFLAGS=-Wall -Wextra -pedantic -O0 -fopenmp -CXXFLAGS=-Wall -Wextra -pedantic -O0 +#CC=g++-11 +CC=mpic++ +CXXFLAGS=-Wall -Wextra -pedantic -O0 -fopenmp LFLAGS=-I/opt/hdf5/include -L/opt/hdf5/lib -SRC=mesh.cpp main.cpp tools.cpp agent.cpp +SRC=cell.cpp mesh.cpp main.cpp agent.cpp LIB=-lm -lhdf5 -lhdf5_cpp OBJ=$(SRC:.cpp=.o) diff --git a/new/merge-results.sh b/merge-results.sh index 5459897..5459897 100755 --- a/new/merge-results.sh +++ b/merge-results.sh @@ -1,1128 +1,546 @@ -// -*- mode: c++; c-default-style: "linux" -*- +// -*- mode: c++ -*- -#include "header.hpp" +#include "mesh.hpp" -Cell::Cell() -{ - x = 0.0; - y = 0.0; - z = 0.0; -} -Cell::Cell(double *coordinates) -{ - x = coordinates[0]; - y = coordinates[1]; - z = coordinates[2]; -} +// ------------------------------------------------------------------------ +// -- PUBLIC METHODS +// ------------------------------------------------------------------------ -Cell::Cell(const Cell &cell) +void Mesh::build_surface_mesh(H5::H5File fp, H5::H5File ep) { - x = cell.x; - y = cell.y; - z = cell.z; -} -int Cell::init_neighbors(size_t n_neighbors) -{ - neighbors.resize(n_neighbors); - - return 0; -} + get_mesh_metadata(fp, SURFACE_NCOL); -int Cell::set_coordinates(double* coordinates) -{ - x = coordinates[0]; - y = coordinates[1]; - z = coordinates[2]; + handle_surface_node_data(fp); - return 0; + build_surface_cells(ep); + } -int Cell::set_neighbors(hsize_t *value_array) +void Mesh::build_subsurface_mesh(H5::H5File fp) { - for (size_t i = 0; i < neighbors.size(); i ++) - neighbors[i] = value_array[i]; - - return 0; -} -int Cell::set_neighbor(int index, hsize_t value) -{ - neighbors[index] = value; - return 0; -} + get_mesh_metadata(fp, SUBSURFACE_NCOL); -int Cell::set_h(double value) -{ - h = value; - return 0; -} + handle_subsurface_node_data(fp); -int Cell::set_v(int index, double value) -{ - v[index] = value; - return 0; + build_subsurface_cells(); + } -int Cell::set_cell_id(int value) +void Mesh::build_environment(H5::H5File ep, int flag) { - cell_id = value; - return 0; -} -hsize_t Cell::get_neighbor(int index) -{ - return neighbors[index]; -} + std::string groupname; + std::string datasetname; + + // Read in velocities in x direction + + if (flag == SUBSURFACE_FLAG) { + groupname = "darcy_velocity.cell.0"; + datasetname = SUBSURFACE_DATA_GROUP_NAME_H5; + } else if (flag == SURFACE_FLAG) { + groupname = "surface-velocity.cell.0"; + datasetname = SURFACE_DATA_GROUP_NAME_H5; + } -double Cell::get_x() -{ - return x; -} + read_velocity(ep, groupname, datasetname, flag, 0); -double Cell::get_y() -{ - return y; -} + // Read in velocities in y direction -double Cell::get_z() -{ - return z; -} + if (flag == SUBSURFACE_FLAG) { + groupname = "darcy_velocity.cell.1"; + datasetname = SUBSURFACE_DATA_GROUP_NAME_H5; + } else if (flag == SURFACE_FLAG) { + groupname = "surface-velocity.cell.1"; + datasetname = SURFACE_DATA_GROUP_NAME_H5; + } -double Cell::get_h() -{ - return h; -} + read_velocity(ep, groupname, datasetname, flag, 1); -double Cell::get_v(int index) -{ - return v[index]; -} + // Read velocities in z direction + + if (flag == SUBSURFACE_FLAG) { + groupname = "darcy_velocity.cell.2"; + datasetname = SUBSURFACE_DATA_GROUP_NAME_H5; + } else if (flag == SURFACE_FLAG) { + groupname = "surface-surface_subsurface_flux.cell.0"; + datasetname = SURFACE_DATA_GROUP_NAME_H5; + } -int Cell::get_cell_id() -{ - return cell_id; -} + read_velocity(ep, groupname, datasetname, flag, 2); -Mesh::Mesh() -{ - Mesh::n_cells = 0; -} + + // Surface environment needs adjustments to the vertical velocity + // and needs to read in water depth + + if (flag == SURFACE_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 + } -int Mesh::get_n_cells() -{ - return Mesh::n_cells; + // Read in water depth + read_water_depth(ep); + + } +#if VERBOSE + else if (flag == SUBSURFACE_FLAG) { + + for (int i = 0; i < n_cells_subsurface; i ++) { + if (i % (n_cells_subsurface / MAXPRINT) == 0) + printf("[~~~~>] Darcy velocity in cell[%d] = { %e, %e, %e } m/s\n", i, + subsurface_cells[i].get_velocity(0), + subsurface_cells[i].get_velocity(1), + subsurface_cells[i].get_velocity(2)); + } + + } +#endif + } -Cell* Mesh::get_cells() -{ - return Mesh::cells; -} +// ------------------------------------------------------------------------ +// GETTERS AND SETTERS +// ------------------------------------------------------------------------ -int Mesh::set_min_coordinates(double *coordinates_min) +Mesh::Mesh() { - xmin = coordinates_min[0]; - ymin = coordinates_min[1]; - zmin = coordinates_min[2]; - - return 0; + n_cells_surface = 0; + n_cells_subsurface = 0; + high_coordinates.resize(DIM); + low_coordinates.resize(DIM); } -int Mesh::set_max_coordinates(double *coordinates_max) +int Mesh::get_n_cells(int flag) { - xmax = coordinates_max[0]; - ymax = coordinates_max[1]; - zmax = coordinates_max[2]; - - return 0; + int n_cells; + + if (flag == SUBSURFACE_FLAG) { + n_cells = n_cells_subsurface; + } else { + n_cells = n_cells_surface; + } + + return n_cells; } -double Mesh::get_xmax() +Cell* Mesh::get_cell(int index, int flag) { - return xmax; -} -double Mesh::get_ymax() -{ - return ymax; + if (flag == SUBSURFACE_FLAG) { + return &(subsurface_cells[index]); + } else { + return &(surface_cells[index]); + } } -double Mesh::get_zmax() +void Mesh::set_high_coordinate(int index, double value) { - return zmax; + high_coordinates[index] = value; } -double Mesh::get_xmin() +double Mesh::get_high_coordinate(int index) { - return xmin; + return high_coordinates[index]; } -double Mesh::get_ymin() +void Mesh::set_low_coordinate(int index, double value) { - return ymin; + low_coordinates[index] = value; } -double Mesh::get_zmin() +double Mesh::get_low_coordinate(int index) { - return zmin; + return low_coordinates[index]; } -#if SUBSURFACE -int Mesh::build_subsurface(H5::H5File fp) -{ - - printf("[OK] building the subsurface mesh\n"); - printf("[OK] reading element node map { \n"); - - H5::DataSet dataset_e = fp.openDataSet("/0/Mesh/MixedElements"); - H5::DataType datatype_e = dataset_e.getDataType(); - H5::DataSpace dataspace_e = dataset_e.getSpace(); - - int rank = dataspace_e.getSimpleExtentNdims(); - - printf("[ok] rank of element node map: %d\n", rank); - - hsize_t *dims_e = new hsize_t[rank]; - dataspace_e.getSimpleExtentDims(dims_e); - - for (int i = 0; i < (int) rank; i ++) - printf("[ok] dims[%d] of element node map: %lld\n", i, dims_e[i]); - - int *data_e = new int[dims_e[0]*dims_e[1]]; - - size_t nx_e = 7; - size_t ny_e = dims_e[0] * dims_e[1] / nx_e; - - int **md_e = new int*[ny_e]; - - for (size_t i = 0; i < ny_e; i ++) - md_e[i] = data_e + i * nx_e; - - printf("[ok] read the element map... "); - - dataset_e.read(data_e, H5::PredType::STD_I32LE); +// ------------------------------------------------------------------------ +// -- PRIVATE METHODS +// ------------------------------------------------------------------------ - printf("done\n"); - - #if VERBOSE - printf("[ok] element map:\n"); - - size_t maxprint = 10; - - if (ny_e < maxprint) - maxprint = ny_e; - - for (size_t i = 0; i < maxprint; i ++) - { - printf("[ok]"); - for (size_t j = 0; j < nx_e; j ++) - printf(" %d ", md_e[i][j]); - printf("\n"); - } - - printf("[ok] ...\n"); - #endif - - printf("[OK] } done\n"); - - printf("[OK] reading nodes { \n"); +void Mesh::handle_surface_node_data(H5::H5File fp) +{ - H5::DataSet dataset_n = fp.openDataSet("/0/Mesh/Nodes"); + H5::DataSet dataset_n = fp.openDataSet(MESH_NODES); H5::DataType datatype_n = dataset_n.getDataType(); H5::DataSpace dataspace_n = dataset_n.getSpace(); - rank = dataspace_n.getSimpleExtentNdims(); + int rank = dataspace_n.getSimpleExtentNdims(); - printf("[ok] rank of nodes: %d\n", rank); - hsize_t *dims_n = new hsize_t[rank]; dataspace_n.getSimpleExtentDims(dims_n); - for (int i = 0; i < (int) rank; i ++) - printf("[ok] dims[%d] of nodes: %lld\n", i, dims_n[i]); +#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]]; - double **md_n = new double*[dims_n[0]]; - - for (size_t i = 0; i < dims_n[0]; i ++) - md_n[i] = data_n + i * dims_n[1]; - - printf("[ok] read the nodes... "); - dataset_n.read(data_n, H5::PredType::IEEE_F64LE); - printf("done\n"); - - #if VERBOSE - printf("[ok] nodes:\n"); - - maxprint = 10; - if (dims_n[0] < 10) - maxprint = dims_n[0]; + nodes.resize(dims_n[0]); - for (size_t i = 0; i < maxprint; i ++) - { - printf("[ok]"); - for (size_t j = 0; j < dims_n[1]; j ++) - printf(" %f ", md_n[i][j]); - printf("\n"); - } - - printf("[ok] ...\n"); - #endif + for (size_t i = 0; i < dims_n[0]; i ++) { + nodes[i].resize(dims_n[1]); + for (size_t j = 0; j < dims_n[1]; j ++) + nodes[i][j] = data_n[i * dims_n[1] + j]; + } - printf("[OK] } done\n"); - - printf("[OK] building cells {\n"); - - Mesh::n_cells = ny_e; - Mesh::cells = new Cell[ny_e]; - - double *coords1 = new double[dims_n[1]]; - double *coords2 = new double[dims_n[1]]; - double *coords3 = new double[dims_n[1]]; - double *coords4 = new double[dims_n[1]]; - double *coords5 = new double[dims_n[1]]; - double *coords6 = new double[dims_n[1]]; - - double *coords_avg = new double[dims_n[1]]; - - double *coords_max = new double[DIM]; - double *coords_min = new double[DIM]; - - for (size_t i = 0; i < DIM; i ++) - { - coords_max[i] = DBL_MIN; - coords_min[i] = DBL_MAX; - } - - for (size_t i = 0; i < ny_e; i ++) - { - - cells[i].init_neighbors(SUBSURFACE_N_NEIGHBORS); - - int n1 = md_e[i][1]; - int n2 = md_e[i][2]; - int n3 = md_e[i][3]; - int n4 = md_e[i][4]; - int n5 = md_e[i][5]; - int n6 = md_e[i][6]; - - for (size_t j = 0; j < dims_n[1]; j++) - { - coords1[j] = md_n[n1][j]; - coords2[j] = md_n[n2][j]; - coords3[j] = md_n[n3][j]; - coords4[j] = md_n[n4][j]; - coords5[j] = md_n[n5][j]; - coords6[j] = md_n[n6][j]; - - coords_avg[j] = (coords1[j] + coords2[j] + coords3[j] + coords4[j] + coords5[j] + coords6[j]) / 6.0; - } - - for (size_t k = 0; k < DIM; k ++) - { - coords_max[k] = (coords_avg[k] > coords_max[k]) ? coords_avg[k] : coords_max[k]; - coords_min[k] = (coords_avg[k] < coords_min[k]) ? coords_avg[k] : coords_min[k]; - } - - cells[i].set_cell_id((int) i); - cells[i].set_coordinates(coords_avg); - - } - - set_min_coordinates(coords_min); - set_max_coordinates(coords_max); - - printf("[ok] min: %f %f %f, max: %f %f %f\n", get_xmin(), get_ymin(), get_zmin(), get_xmax(), get_ymax(), get_zmax()); - - #if VERBOSE - - maxprint = 10; - - if (ny_e < maxprint) - maxprint = ny_e; - - for (size_t i = 0; i < maxprint; i ++) - { - printf("[ok] cell[%ld] centroid: %f %f %f\n", i, cells[i].get_x(), cells[i].get_y(), cells[i].get_z()); - } +} - printf("[ok] ...\n"); - - #endif - - delete[] coords_avg; - delete[] coords1; - delete[] coords2; - delete[] coords3; - delete[] coords4; - delete[] coords5; - delete[] coords6; - - printf("[OK] } done\n"); +void Mesh::build_surface_cells(H5::H5File ep) +{ - printf("[OK] initialize edges {\n"); - - int ns[SUBSURFACE_N_NODES]; - - // Szudzik IDs for very large number of edges is not safe - std::vector<std::vector<std::vector<hsize_t>>> edge_ids(ny_e, std::vector<std::vector<hsize_t>>(SUBSURFACE_N_EDGES, std::vector<hsize_t>(2))); - - - for (size_t i = 0; i < ny_e; i ++) { - - for (size_t j = 0; j < SUBSURFACE_N_NODES; j ++) { - ns[j] = md_e[i][j+1]; - } + std::vector<double> coords1(DIM-1); + std::vector<double> coords2(DIM-1); + std::vector<double> coords3(DIM-1); + std::vector<double> coords_avg(DIM-1); - // top of the cell + read_elevation(ep); - if (ns[0] < ns[1]) { - edge_ids[i][0][0] = ns[0]; - edge_ids[i][0][1] = ns[1]; - } else { - edge_ids[i][0][1] = ns[0]; - edge_ids[i][0][0] = ns[1]; - } + for (int i = 0; i < n_cells_surface; i ++) { - if (ns[1] < ns[2]) { - edge_ids[i][1][0] = ns[1]; - edge_ids[i][1][1] = ns[2]; - } else { - edge_ids[i][1][1] = ns[2]; - edge_ids[i][1][0] = ns[1]; - } + int n1 = element_node_indices[i][1]; + int n2 = element_node_indices[i][2]; + int n3 = element_node_indices[i][3]; - if (ns[2] < ns[0]) { - edge_ids[i][2][0] = ns[2]; - edge_ids[i][2][1] = ns[0]; - } else { - edge_ids[i][2][1] = ns[0]; - edge_ids[i][2][0] = ns[2]; - } + std::vector<double> edge_length(DIM); - // bottom of the cell + for (size_t k = 0; k < DIM; k ++) + edge_length[k] = DBL_MAX; - if (ns[3] < ns[4]) { - edge_ids[i][3][0] = ns[3]; - edge_ids[i][3][1] = ns[4]; - } else { - edge_ids[i][3][0] = ns[4]; - edge_ids[i][3][1] = ns[3]; - } - if (ns[4] < ns[5]) { - edge_ids[i][4][0] = ns[4]; - edge_ids[i][4][1] = ns[5]; - } else { - edge_ids[i][4][0] = ns[5]; - edge_ids[i][4][1] = ns[4]; - } + double smallest_distance; + double largest_distance; + + for (size_t j = 0; j < DIM-1; j ++) { - if (ns[5] < ns[3]) { - edge_ids[i][5][0] = ns[5]; - edge_ids[i][5][1] = ns[3]; - } else { - edge_ids[i][5][0] = ns[3]; - edge_ids[i][5][1] = ns[5]; - } + coords1[j] = nodes[n1][j]; + coords2[j] = nodes[n2][j]; + coords3[j] = nodes[n3][j]; - // sides of the cell + coords_avg[j] = (coords1[j] + coords2[j] + coords3[j]) / 3.0; - if (ns[0] < ns[3]) { - edge_ids[i][6][0] = ns[0]; - edge_ids[i][6][1] = ns[3]; - } else { - edge_ids[i][6][0] = ns[3]; - edge_ids[i][6][1] = ns[0]; - } + edge_length[0] = fabs(coords1[j] - coords_avg[j]); + edge_length[1] = fabs(coords2[j] - coords_avg[j]); + edge_length[2] = fabs(coords3[j] - coords_avg[j]); - if (ns[1] < ns[4]) { - edge_ids[i][7][0] = ns[1]; - edge_ids[i][7][1] = ns[4]; - } else { - edge_ids[i][7][0] = ns[4]; - edge_ids[i][7][1] = ns[1]; - } + smallest_distance = edge_length[0]; + largest_distance = edge_length[0]; - if (ns[2] < ns[5]) { - edge_ids[i][8][0] = ns[2]; - edge_ids[i][8][1] = ns[5]; - } else { - edge_ids[i][8][0] = ns[5]; - edge_ids[i][8][1] = ns[2]; + for (size_t k = 1; k < DIM; k ++) { + smallest_distance = std::min(smallest_distance, edge_length[k]); + largest_distance = std::max(largest_distance, edge_length[k]); } } - - printf("[OK] } done\n"); - - printf("[OK] initialize neighbors {\n"); - - hsize_t e00; - hsize_t e01; - hsize_t e10; - hsize_t e11; - hsize_t e20; - hsize_t e21; - hsize_t e30; - hsize_t e31; - hsize_t e40; - hsize_t e41; - hsize_t e50; - hsize_t e51; - hsize_t e60; - hsize_t e61; - hsize_t e70; - hsize_t e71; - hsize_t e80; - hsize_t e81; - - hsize_t _e0; - hsize_t _e1; - for (size_t i = 0; i < ny_e; i ++) { - - double xi = cells[i].get_x(); - double yi = cells[i].get_y(); - double zi = cells[i].get_z(); + surface_cells[i].set_characteristic_length(smallest_distance); + surface_cells[i].set_maximum_length(largest_distance); + + for (size_t k = 0; k < DIM-1; k ++) + surface_cells[i].set_coordinate(k, coords_avg[k]); - for (size_t j = 0; j < SUBSURFACE_N_NEIGHBORS; j ++) - cells[i].set_neighbor(j, -1); - -#if VERBOSE_HIGH - if (i % 100 == 0) - printf("[ok] searching neighbors of cell: %ld/%ld\n", i+1, ny_e); -#endif + surface_cells[i].set_cell_id((int) i); - e00 = edge_ids[i][0][0]; - e01 = edge_ids[i][0][1]; + // find the coupled cell in the subsurface + std::vector<double> distance(DIM); + int coupled_neighbor = -1; // store current coupled neighbor here + double dz = DBL_MAX; - e10 = edge_ids[i][1][0]; - e11 = edge_ids[i][1][1]; + for (int k = 0; k < n_cells_subsurface; k ++) { - e20 = edge_ids[i][2][0]; - e21 = edge_ids[i][2][1]; + for (size_t j = 0; j < DIM; j ++) + distance[j] = subsurface_cells[k].get_coordinate(j) - surface_cells[i].get_coordinate(j); - e30 = edge_ids[i][3][0]; - e31 = edge_ids[i][3][1]; + // Because the mesh is a stacked mesh, the only candidates for the subsurface + // neighbor are the cells located in the same column. These cells will have + // the same centroid in the xy-plane as the surface cell. The problem reduces + // to finding the nearest cell in the vertical direction in the column. - e40 = edge_ids[i][4][0]; - e41 = edge_ids[i][4][1]; - - e50 = edge_ids[i][5][0]; - e51 = edge_ids[i][5][1]; + double xy_distance = sqrt(distance[0]*distance[0] + distance[1]*distance[1]); - e60 = edge_ids[i][6][0]; - e61 = edge_ids[i][6][1]; - - e70 = edge_ids[i][7][0]; - e71 = edge_ids[i][7][1]; - - e80 = edge_ids[i][8][0]; - e81 = edge_ids[i][8][1]; - -#if VERBOSE_HIGH - if (i % 100 == 0) - printf("[ok] indices: [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld]\n", e00, e01, e10, e11, e20, e21, e30, e31, e40, e41, e50, e51, e60, e61, e70, e71, e80, e81); -#endif - - int counter = 0; - - for (size_t k = 0; k < ny_e; k ++) { - - // we will limit our search to a 3D sphere with radius RADIUS around the cell - double xk = cells[k].get_x(); - double yk = cells[k].get_y(); - double zk = cells[k].get_z(); - - double distance = sqrt((xk - xi) * (xk - xi) + (yk - yi) * (yk - yi) + (zk - zi) * (zk - zi)); - - if (i != k && distance <= RADIUS) { - - int _tmp_count = 0; - - for (size_t m = 0; m < SUBSURFACE_N_EDGES; m ++) { - - _e0 = edge_ids[k][m][0]; - _e1 = edge_ids[k][m][1]; - - if ((e00 == _e0 && e01 == _e1) || (e10 == _e0 && e11 == _e1) || - (e20 == _e0 && e21 == _e1) || (e30 == _e0 && e31 == _e1) || - (e40 == _e0 && e41 == _e1) || (e50 == _e0 && e51 == _e1) || - (e60 == _e0 && e61 == _e1) || (e70 == _e0 && e71 == _e1) || - (e80 == _e0 && e81 == _e1)) { - - _tmp_count ++; - - if (_tmp_count == 2) { // a cell k is a neighbor of cell - // i if 2 of their edges match, - // which means that a surface is - // shared between k and i - cells[i].set_neighbor(counter, (hsize_t) k); - counter ++; -#if VERBOSE_HIGH - if (i % 100 == 0) - printf("[ok] found neighbor %d/%d: %ld [%lld, %lld]\n", counter, SUBSURFACE_N_NEIGHBORS, k, _e0, _e1); -#endif - break; - } - - } - } - } - - if (counter == SUBSURFACE_N_NEIGHBORS) + if (xy_distance <= FUZZ) { + dz = std::min(dz, fabs(distance[2])); + coupled_neighbor = k; + if (dz <= 2.0 * subsurface_cells[k].get_characteristic_length()) break; - } - } - - #if VERBOSE - - maxprint = (ny_e < 10) ? ny_e : 10; - - for (size_t i = 0; i < maxprint; i ++) - { - printf("[ok]"); - for (size_t j = 0; j < SUBSURFACE_N_NEIGHBORS; j ++) - { - printf(" %6lld ", cells[i].get_neighbor(j)); - } - printf("\n"); + } - - #endif - - printf("[OK] } done\n"); - - // delete[] szudzik_ids; - delete[] md_e; - delete[] data_e; - delete[] dims_e; - - delete[] md_n; - delete[] data_n; - delete[] dims_n; - - printf("[OK] mesh built\n\n"); - - return 0; - -} - -int Mesh::read_subsurface_data(H5::H5File dp) -{ - printf("[OK] build subsurface environment\n"); + surface_cells[i].set_coupled_neighbor(coupled_neighbor); + subsurface_cells[coupled_neighbor].set_coupled_neighbor(i); - // read darcy velocity x - - H5::Group group = dp.openGroup("darcy_velocity.cell.0"); - - printf("[OK] reading velocity x{ \n"); - printf("[ok] # of objects: %lld\n", group.getNumObjs()); - - H5::DataSet dataset = group.openDataSet(DATA_GROUP_NAME_H5); - H5::DataType datatype = dataset.getDataType(); - H5::DataSpace dataspace = dataset.getSpace(); - - int rank = dataspace.getSimpleExtentNdims(); - - printf("[ok] rank of data: %d\n", rank); - - hsize_t *dims = new hsize_t[rank]; - dataspace.getSimpleExtentDims(dims); - - for (int i = 0; i < (int) rank; i ++) - printf("[ok] dims[%d] of data: %lld\n", i, dims[i]); - - double *data = new double[dims[0]*dims[1]]; - double **md = new double*[dims[0]]; - - for (size_t i = 0; i < dims[0]; i ++) - md[i] = data + i * dims[1]; - - dataset.read(data, H5::PredType::IEEE_F64LE); - #if VERBOSE - size_t maxprint = (dims[0] > 10) ? 10 : dims[0]; - for (size_t i = 0; i < maxprint; i ++) - printf("[ok] %f\n", data[i]); -#endif - - printf("[OK] } done\n"); - printf("[OK] set velocity x { \n"); - - for (int i = 0; i < n_cells; i ++) - cells[i].set_v(0, data[i]); - -#if VERBOSE - maxprint = (dims[0] > 10) ? 10 : dims[0]; - for (size_t i = 0; i < dims[0]; i ++) { - if (cells[i].get_v(0) > 0.0) { - printf("[ok] %f\n", cells[i].get_v(0)); + 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 - - printf("[OK] } done\n"); - return 0; -} - -#endif - -int Mesh::build_surface(H5::H5File fp) -{ - - printf("[OK] building the surface mesh\n"); - printf("[OK] reading element node map { \n"); - H5::DataSet dataset_e = fp.openDataSet("/0/Mesh/MixedElements"); - H5::DataType datatype_e = dataset_e.getDataType(); - H5::DataSpace dataspace_e = dataset_e.getSpace(); - - int rank = dataspace_e.getSimpleExtentNdims(); - - printf("[ok] rank of element node map: %d\n", rank); - - hsize_t *dims_e = new hsize_t[rank]; - dataspace_e.getSimpleExtentDims(dims_e); - - for (int i = 0; i < (int) rank; i ++) - printf("[ok] dims[%d] of element node map: %lld\n", i, dims_e[i]); - - int *data_e = new int[dims_e[0]*dims_e[1]]; - - // structure of hdf5 is <element type> <node 1> <node 2> <node 3> - size_t nx_e = 4; - size_t ny_e = dims_e[0] * dims_e[1] / nx_e; - - int **md_e = new int*[ny_e]; - - for (size_t i = 0; i < ny_e; i ++) - md_e[i] = data_e + i * nx_e; - - printf("[ok] read the element map... "); - - dataset_e.read(data_e, H5::PredType::STD_I32LE); - - printf("done\n"); - - #if VERBOSE - printf("[ok] element map:\n"); - - size_t maxprint = 10; + } - if (ny_e < maxprint) - maxprint = ny_e; - - for (size_t i = 0; i < maxprint; i ++) - { - printf("[ok]"); - for (size_t j = 0; j < nx_e; j ++) - printf(" %d ", md_e[i][j]); - printf("\n"); - } +} - printf("[ok] ...\n"); - #endif +void Mesh::handle_subsurface_node_data(H5::H5File fp) +{ - printf("[OK] } done\n"); - - printf("[OK] reading nodes { \n"); - - H5::DataSet dataset_n = fp.openDataSet("/0/Mesh/Nodes"); + H5::DataSet dataset_n = fp.openDataSet(MESH_NODES); H5::DataType datatype_n = dataset_n.getDataType(); H5::DataSpace dataspace_n = dataset_n.getSpace(); - rank = dataspace_n.getSimpleExtentNdims(); + int rank = dataspace_n.getSimpleExtentNdims(); - printf("[ok] rank of nodes: %d\n", rank); - hsize_t *dims_n = new hsize_t[rank]; dataspace_n.getSimpleExtentDims(dims_n); - for (int i = 0; i < (int) rank; i ++) - printf("[ok] dims[%d] of nodes: %lld\n", i, dims_n[i]); - +#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]]; - double **md_n = new double*[dims_n[0]]; - - for (size_t i = 0; i < dims_n[0]; i ++) - md_n[i] = data_n + i * dims_n[1]; - - printf("[ok] read the nodes... "); - dataset_n.read(data_n, H5::PredType::IEEE_F64LE); - printf("done\n"); - - #if VERBOSE - printf("[ok] nodes:\n"); - - maxprint = 10; - if (dims_n[0] < 10) - maxprint = dims_n[0]; + nodes.resize(dims_n[0]); - for (size_t i = 0; i < maxprint; i ++) - { - printf("[ok]"); - for (size_t j = 0; j < dims_n[1]; j ++) - printf(" %f ", md_n[i][j]); - printf("\n"); - } - - printf("[ok] ...\n"); - #endif + for (size_t i = 0; i < dims_n[0]; i ++) { + nodes[i].resize(dims_n[1]); + for (size_t j = 0; j < dims_n[1]; j ++) + nodes[i][j] = data_n[i * dims_n[1] + j]; + } - printf("[OK] } done\n"); +} - printf("[OK] building cells {\n"); +void Mesh::build_subsurface_cells() +{ - Mesh::n_cells = ny_e; - Mesh::cells = new Cell[Mesh::n_cells]; + std::vector<double> coords1(DIM); + std::vector<double> coords2(DIM); + std::vector<double> coords3(DIM); + std::vector<double> coords4(DIM); + std::vector<double> coords5(DIM); + std::vector<double> coords6(DIM); - double *coords1 = new double[dims_n[1]]; - double *coords2 = new double[dims_n[1]]; - double *coords3 = new double[dims_n[1]]; + std::vector<double> coords_avg(DIM); - double *coords_avg = new double[dims_n[1]]; + for (size_t i = 0; i < DIM; i ++) { + set_high_coordinate(i, DBL_MIN); + set_low_coordinate(i, DBL_MAX); + } - double *coords_max = new double[DIM]; - double *coords_min = new double[DIM]; + double smallest_distance; + double largest_distance; + + for (int i = 0; i < n_cells_subsurface; i ++) { - for (size_t i = 0; i < DIM; i ++) - { - coords_max[i] = DBL_MIN; - coords_min[i] = DBL_MAX; - } + int n1 = element_node_indices[i][1]; + int n2 = element_node_indices[i][2]; + int n3 = element_node_indices[i][3]; + int n4 = element_node_indices[i][4]; + int n5 = element_node_indices[i][5]; + int n6 = element_node_indices[i][6]; - for (size_t i = 0; i < ny_e; i ++) - { + std::vector<double> edge_length(N_VERTICES); - cells[i].init_neighbors(SURFACE_N_EDGES); + for (size_t k = 0; k < N_VERTICES; k ++) + edge_length[k] = DBL_MAX; + + for (size_t j = 0; j < DIM; j++) { - int n1 = md_e[i][1]; - int n2 = md_e[i][2]; - int n3 = md_e[i][3]; - - for (size_t j = 0; j < dims_n[1]; j++) - { - coords1[j] = md_n[n1][j]; - coords2[j] = md_n[n2][j]; - coords3[j] = md_n[n3][j]; - - coords_avg[j] = (coords1[j] + coords2[j] + coords3[j]) / 3.0; - } + coords1[j] = nodes[n1][j]; + coords2[j] = nodes[n2][j]; + coords3[j] = nodes[n3][j]; + coords4[j] = nodes[n4][j]; + coords5[j] = nodes[n5][j]; + coords6[j] = nodes[n6][j]; + + coords_avg[j] = (coords1[j] + coords2[j] + coords3[j] + coords4[j] + coords5[j] + coords6[j]) / 6.0; + + // the characteristic length is the radius of the insphere of + // the wedge + edge_length[0] = fabs(coords1[j] - coords_avg[j]); + edge_length[1] = fabs(coords2[j] - coords_avg[j]); + edge_length[2] = fabs(coords3[j] - coords_avg[j]); + edge_length[3] = fabs(coords4[j] - coords_avg[j]); + edge_length[4] = fabs(coords5[j] - coords_avg[j]); + edge_length[5] = fabs(coords6[j] - coords_avg[j]); + + smallest_distance = edge_length[0]; + largest_distance = edge_length[0]; + + for (size_t k = 1; k < N_VERTICES; k ++) { + smallest_distance = std::min(smallest_distance, edge_length[k]); + largest_distance = std::max(largest_distance, edge_length[k]); + } + + } - for (size_t k = 0; k < DIM; k ++) - { - coords_max[k] = (coords_avg[k] > coords_max[k]) ? coords_avg[k] : coords_max[k]; - coords_min[k] = (coords_avg[k] < coords_min[k]) ? coords_avg[k] : coords_min[k]; - } + subsurface_cells[i].set_characteristic_length(smallest_distance); + subsurface_cells[i].set_maximum_length(largest_distance); + + for (size_t k = 0; k < DIM; k ++) { - cells[i].set_cell_id((int) i); - cells[i].set_coordinates(coords_avg); + // update lowest and highest corner of mesh + high_coordinates[k] = std::max(coords_avg[k], high_coordinates[k]); + low_coordinates[k] = std::min(coords_avg[k], low_coordinates[k]); + subsurface_cells[i].set_coordinate(k, coords_avg[k]); + } + + subsurface_cells[i].set_cell_id((int) i); + + } + +} - set_min_coordinates(coords_min); - set_max_coordinates(coords_max); - - printf("[ok] min: %f %f %f, max: %f %f %f\n", get_xmin(), get_ymin(), get_zmin(), get_xmax(), get_ymax(), get_zmax()); +void Mesh::get_mesh_metadata(H5::H5File fp, int flag) +{ - #if VERBOSE + H5::DataSet dataset_e = fp.openDataSet(MESH_MIXED_ELEMENTS); + H5::DataType datatype_e = dataset_e.getDataType(); + H5::DataSpace dataspace_e = dataset_e.getSpace(); - maxprint = 10; + int rank = dataspace_e.getSimpleExtentNdims(); - if (ny_e < maxprint) - maxprint = ny_e; + hsize_t *dims_e = new hsize_t[rank]; + dataspace_e.getSimpleExtentDims(dims_e); - for (size_t i = 0; i < maxprint; i ++) - { - printf("[ok] cell[%ld] centroid: %f %f %f\n", i, cells[i].get_x(), cells[i].get_y(), cells[i].get_z()); - } +#if VERBOSE + for (int i = 0; i < rank; i ++) + std::cout << "[~~~>] Dimension in " << i << " of element node map: " << dims_e[i] << std::endl; +#endif - printf("[ok] ...\n"); - - #endif + int *data_e = new int[dims_e[0]*dims_e[1]]; + dataset_e.read(data_e, H5::PredType::STD_I32LE); - delete[] coords_avg; - delete[] coords1; - delete[] coords2; - delete[] coords3; + size_t nx_e = flag; + size_t ny_e = dims_e[0] * dims_e[1] / nx_e; - printf("[OK] } done\n"); + element_node_indices.resize(ny_e); - printf("[OK] initialize edges {\n"); - - int ns[N_EDGES]; - hsize_t **szudzik_ids = new hsize_t*[ny_e]; + for (size_t i = 0; i < ny_e; i ++) { + element_node_indices[i].resize(nx_e); + for (size_t j = 0; j < nx_e; j ++) + element_node_indices[i][j] = data_e[i * nx_e + j]; + } - for (size_t i = 0; i < ny_e; i ++) - szudzik_ids[i] = new hsize_t[N_EDGES]; + if (flag == SUBSURFACE_NCOL) { + n_cells_subsurface = ny_e; + subsurface_cells.resize(ny_e); + } else if (flag == SURFACE_NCOL) { + n_cells_surface = ny_e; + surface_cells.resize(ny_e); + } - for (size_t i = 0; i < ny_e; i ++) - { - for (size_t j = 0; j < N_EDGES; j ++) - ns[j] = md_e[i][j+1]; - - for (size_t j = 0; j < N_EDGES - 1; j ++) - szudzik_ids[i][j] = szudzik(ns[j], ns[j+1]); - - szudzik_ids[i][N_EDGES-1] = szudzik(ns[N_EDGES - 1], ns[0]); - } - - #if VERBOSE +} - maxprint = (ny_e < 10) ? ny_e : 10; - - for (size_t i = 0; i < maxprint; i ++) - { - printf("[ok]"); - for (size_t j = 0; j < N_EDGES; j ++) - { - printf(" %6lld ", szudzik_ids[i][j]); - } - printf("\n"); - } +void Mesh::read_velocity(H5::H5File fp, std::string groupname, std::string datasetname, int flag, int index) +{ - #endif + H5::Group group = fp.openGroup(groupname.c_str()); + H5::DataSet dataset = group.openDataSet(datasetname.c_str()); + H5::DataType datatype = dataset.getDataType(); + H5::DataSpace dataspace = dataset.getSpace(); - printf("[OK] } done\n"); - - printf("[OK] initialize neighbors {\n"); - - hsize_t e0; - hsize_t e1; - - for (size_t i = 0; i < ny_e; i ++) - { - for (size_t j = 0; j < N_EDGES; j ++) - { - - e0 = szudzik_ids[i][j]; - cells[i].set_neighbor(j, -1); - - for (size_t k = 0; k < ny_e; k ++) - { - - for (size_t m = 0; m < N_EDGES; m ++) - { - - e1 = szudzik_ids[k][m]; - - if (e0 == e1 && i != k) { - cells[i].set_neighbor(j, (hsize_t) k); - } - - } - } - } - } + int rank = dataspace.getSimpleExtentNdims(); - #if VERBOSE + hsize_t *dims = new hsize_t[rank]; + dataspace.getSimpleExtentDims(dims); - maxprint = (ny_e < 10) ? ny_e : 10; - - for (size_t i = 0; i < maxprint; i ++) - { - printf("[ok]"); - for (size_t j = 0; j < N_EDGES; j ++) - { - printf(" %6lld ", cells[i].get_neighbor(j)); - } - printf("\n"); - } - - #endif +#if VERBOSE + for (int i = 0; i < rank; i ++) + std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl; +#endif - printf("[OK] } done\n"); - - delete[] szudzik_ids; - delete[] md_e; - delete[] data_e; - delete[] dims_e; + double *data = new double[dims[0]*dims[1]]; + dataset.read(data, H5::PredType::IEEE_F64LE); - delete[] md_n; - delete[] data_n; - delete[] dims_n; - - printf("[OK] mesh built\n\n"); + if (flag == SUBSURFACE_FLAG) { + for (int i = 0; i < n_cells_subsurface; i ++) + subsurface_cells[i].set_velocity(index, data[i]); + } else if (flag == SURFACE_FLAG) { + for (int i = 0; i < n_cells_surface; i ++) + surface_cells[i].set_velocity(index, data[i]); + } - return 0; - } -int Mesh::read_surface_data(H5::H5File dp) +void Mesh::read_water_depth(H5::H5File fp) { - - printf("[OK] build surface environment\n"); - - // read velocity x - - H5::Group group = dp.openGroup("surface-velocity.cell.0"); - - printf("[OK] reading velocity x { \n"); - printf("[ok] # of objects: %lld\n", group.getNumObjs()); - - H5::DataSet dataset = group.openDataSet(DATA_GROUP_NAME_H5); - H5::DataType datatype = dataset.getDataType(); - H5::DataSpace dataspace = dataset.getSpace(); + H5::Group group = fp.openGroup("surface-ponded_depth.cell.0"); + H5::DataSet dataset = group.openDataSet(SURFACE_DATA_GROUP_NAME_H5); + H5::DataType datatype = dataset.getDataType(); + H5::DataSpace dataspace = dataset.getSpace(); + int rank = dataspace.getSimpleExtentNdims(); - printf("[ok] rank of data: %d\n", rank); - hsize_t *dims = new hsize_t[rank]; dataspace.getSimpleExtentDims(dims); - for (int i = 0; i < (int) rank; i ++) - printf("[ok] dims[%d] of data: %lld\n", i, dims[i]); - +#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]]; - double **md = new double*[dims[0]]; - - for (size_t i = 0; i < dims[0]; i ++) - md[i] = data + i * dims[1]; - dataset.read(data, H5::PredType::IEEE_F64LE); - - #if VERBOSE - size_t maxprint = (dims[0] > 10) ? 10 : dims[0]; - for (size_t i = 0; i < maxprint; i ++) - printf("[ok] %f\n", data[i]); - #endif - - printf("[OK] } done\n"); - printf("[OK] set velocity x { \n"); - - for (int i = 0; i < n_cells; i ++) - cells[i].set_v(0, data[i]); - - #if VERBOSE - maxprint = (dims[0] > 10) ? 10 : dims[0]; - for (size_t i = 0; i < dims[0]; i ++) { - if (cells[i].get_v(0) > 0.0) { - printf("[ok] %f\n", cells[i].get_v(0)); + + 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 } - #endif - printf("[OK] } done\n"); - - // read velocity y - - H5::Group group_vy = dp.openGroup("surface-velocity.cell.1"); - - printf("[OK] reading velocity y { \n"); - printf("[ok] # of objects: %lld\n", group_vy.getNumObjs()); - - H5::DataSet dataset_vy = group_vy.openDataSet(DATA_GROUP_NAME_H5); - H5::DataType datatype_vy = dataset_vy.getDataType(); - H5::DataSpace dataspace_vy = dataset_vy.getSpace(); - - rank = dataspace_vy.getSimpleExtentNdims(); - - printf("[ok] rank of data: %d\n", rank); - - hsize_t *dims_vy = new hsize_t[rank]; - dataspace_vy.getSimpleExtentDims(dims_vy); - - for (int i = 0; i < (int) rank; i ++) - printf("[ok] dims[%d] of data: %lld\n", i, dims_vy[i]); - - double *data_vy = new double[dims_vy[0]*dims_vy[1]]; - double **md_vy = new double*[dims_vy[0]]; +} - for (size_t i = 0; i < dims_vy[0]; i ++) - md_vy[i] = data_vy + i * dims_vy[1]; +void Mesh::read_elevation(H5::H5File fp) +{ - dataset_vy.read(data_vy, H5::PredType::IEEE_F64LE); + H5::Group group = fp.openGroup("surface-elevation.cell.0"); + H5::DataSet dataset = group.openDataSet(SURFACE_DATA_GROUP_NAME_H5); + H5::DataType datatype = dataset.getDataType(); + H5::DataSpace dataspace = dataset.getSpace(); - #if VERBOSE - maxprint = (dims_vy[0] > 10) ? 10 : dims_vy[0]; - for (size_t i = 0; i < maxprint; i ++) - printf("[ok] %f\n", data_vy[i]); - #endif - - printf("[OK] } done\n"); - printf("[OK] set velocity y { \n"); - - for (int i = 0; i < n_cells; i ++) - cells[i].set_v(1, data_vy[i]); - - #if VERBOSE - maxprint = (dims_vy[0] > 10) ? 10 : dims_vy[0]; - for (size_t i = 0; i < maxprint; i ++) - printf("[ok] %f\n", cells[i].get_v(1)); - #endif - - printf("[OK] } done\n"); - - // read water depth - H5::Group group_h = dp.openGroup("surface-ponded_depth.cell.0"); - - printf("[OK] reading ponded water depth { \n"); - printf("[ok] # of objects: %lld\n", group_vy.getNumObjs()); - - H5::DataSet dataset_h = group_h.openDataSet(DATA_GROUP_NAME_H5); - H5::DataType datatype_h = dataset_h.getDataType(); - H5::DataSpace dataspace_h = dataset_h.getSpace(); - - rank = dataspace_h.getSimpleExtentNdims(); - - printf("[ok] rank of data: %d\n", rank); - - hsize_t *dims_h = new hsize_t[rank]; - dataspace_h.getSimpleExtentDims(dims_h); - - for (int i = 0; i < (int) rank; i ++) - printf("[ok] dims[%d] of data: %lld\n", i, dims_h[i]); - - double *data_h = new double[dims_h[0]*dims_h[1]]; - double **md_h = new double*[dims_h[0]]; + int rank = dataspace.getSimpleExtentNdims(); - for (size_t i = 0; i < dims_h[0]; i ++) - md_h[i] = data_h + i * dims_h[1]; + hsize_t *dims = new hsize_t[rank]; + dataspace.getSimpleExtentDims(dims); - dataset_h.read(data_h, H5::PredType::IEEE_F64LE); - - #if VERBOSE - maxprint = (dims_h[0] > 10) ? 10 : dims_h[0]; - for (size_t i = 0; i < maxprint; i ++) - printf("[ok] %f\n", data_h[i]); - #endif - - printf("[OK] } done\n"); - printf("[OK] set ponded water depth { \n"); - - for (int i = 0; i < n_cells; i ++) - cells[i].set_h(data_h[i]); - - #if VERBOSE - maxprint = (dims_h[0] > 10) ? 10 : dims_h[0]; - for (size_t i = 0; i < maxprint; i ++) - printf("[ok] %f\n", cells[i].get_h()); - #endif - - printf("[OK] } done\n"); - - delete [] md_h; - delete [] data_h; - delete [] md_vy; - delete [] data_vy; - delete [] md; - delete [] data; - delete [] dims; +#if VERBOSE + for (int i = 0; i < rank; i ++) + std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl; +#endif - printf("[OK] data read\n\n"); + double *data = new double[dims[0]*dims[1]]; + dataset.read(data, H5::PredType::IEEE_F64LE); - return 0; + 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 + } } - @@ -1,115 +1,59 @@ -// -*- mode: c++; c-default-style: "linux" -*- +// -*- mode: c++ -*- #ifndef MESH_HPP #define MESH_HPP #include "header.hpp" +#include "cell.hpp" -// mesh.hpp -// ======== -// Mesh topology to get environment variables that inform the -// HydroAgent (agent.hpp). -// Reads in a HDF5 mesh file and builds the topology. The HydroAgent -// stores the cell it was last located in. Then, in the next iteration -// it searches the neighborhood to locate itself. This makes the search -// more efficient. - -class Cell +class Mesh { public: - Cell(); - Cell(double* coordinates); - Cell(const Cell &cell); - - double get_x(); - double get_y(); - double get_z(); - - double get_h(); - double get_v(int index); - - int get_cell_id(); - - hsize_t get_neighbor(int index); + Mesh(); - int set_coordinates(double *coordinates); - - int set_neighbors(hsize_t *value_array); - int set_neighbor(int index, hsize_t value); + void build_surface_mesh(H5::H5File fp, H5::H5File ep); + void build_subsurface_mesh(H5::H5File fp); - int set_h(double value); - int set_v(int index, double value); + void build_environment(H5::H5File ep, int flag); + void build_subsurface_environment(H5::H5File ep); - int set_cell_id(int value); + Cell* get_cell(int index, int flag); - int init_neighbors(size_t n_neighbors); + void set_high_coordinate(int index, double value); + double get_high_coordinate(int index); + void set_low_coordinate(int index, double value); + double get_low_coordinate(int index); -protected: + int get_n_cells(int flag); private: - - int dim; // n_neighbors = 2 * dim - 1 - int cell_id; - - double x; - double y; - double z; - - double h; - double v[3]; - - std::vector<hsize_t> neighbors; - -}; - -class Mesh -{ - -public: - - Mesh(); - int get_n_cells(); + int n_cells_subsurface; + int n_cells_surface; - int build_surface(H5::H5File fp); - int read_surface_data(H5::H5File dp); - -#if SUBSURFACE - int build_subsurface(H5::H5File fp); - int read_subsurface_data(H5::H5File dp); -#endif + std::vector<Cell> subsurface_cells; + std::vector<Cell> surface_cells; - Cell* get_cells(); + std::vector<double> high_coordinates; + std::vector<double> low_coordinates; - int set_min_coordinates(double *coordinates_min); - int set_max_coordinates(double *coordinates_max); - - double get_xmax(); - double get_ymax(); - double get_zmax(); - - double get_xmin(); - double get_ymin(); - double get_zmin(); - -protected: + std::vector<std::vector<double>> nodes; + std::vector<std::vector<int>> element_node_indices; + //std::vector<std::vector<std::vector<hsize_t>>> edge_ids; -private: + void handle_subsurface_node_data(H5::H5File fp); + void build_subsurface_cells(); - int n_cells; - Cell *cells; + void handle_surface_node_data(H5::H5File fp); + void build_surface_cells(H5::H5File fp); - double xmin; - double xmax; - - double ymin; - double ymax; - - double zmin; - double zmax; - + void get_mesh_metadata(H5::H5File fp, int ncol); + void read_velocity(H5::H5File fp, std::string groupname, std::string datasetname, int flag, int index); + void read_water_depth(H5::H5File fp); + void read_elevation(H5::H5File fp); + }; #endif diff --git a/new/agent.cpp b/new/agent.cpp deleted file mode 100644 index 97d2e69..0000000 --- a/new/agent.cpp +++ /dev/null @@ -1,245 +0,0 @@ -// -*- mode: c++ -*- - -#include "agent.hpp" - -// ------------------------------------------------------------------------ -// -- PUBLIC METHODS -// ------------------------------------------------------------------------ - -void Agent::advance(Mesh* mesh, size_t step) -{ - - switch(domain[step-1]) { - case SURFACE_FLAG: - advance_surface(mesh,step); - break; - case SUBSURFACE_FLAG: - advance_subsurface(mesh,step); - break; - default: - set_path(path[step-1][0], path[step-1][1], path[step-1][2], step); - time[step] = time[step-1]; - domain[step] = OUTSIDE_FLAG; - break; - } - -} - -// ------------------------------------------------------------------------ -// GETTERS AND SETTERS -// ------------------------------------------------------------------------ - -Agent::Agent() -{ - time.resize(N_ITERATION); - domain.resize(N_ITERATION); - path.resize(N_ITERATION); - for (size_t i = 0; i < N_ITERATION; i ++) - path[i].resize(DIM); -} - -void Agent::set_cell_index(int value) -{ - cell_index = value; -} - -void Agent::set_domain(int value, size_t step) -{ - domain[step] = value; -} - -double Agent::get_time(size_t step) -{ - return time[step]; -} - -void Agent::set_path(double x, double y, double z, size_t step) -{ - path[step][0] = x; - path[step][1] = y; - path[step][2] = z; -} - -double Agent::get_timestepsize() -{ - return time_step_size; -} - -double Agent::get_path(size_t step, size_t index) -{ - return path[step][index]; -} - -int Agent::get_domain(size_t step) -{ - return domain[step]; -} - -// ------------------------------------------------------------------------ -// -- PRIVATE METHODS -// ------------------------------------------------------------------------ - -void Agent::advance_surface(Mesh* mesh, size_t step) -{ - Cell *location = mesh->get_cell(cell_index, SURFACE_FLAG); - - std::vector<double> velocity(DIM); - std::vector<double> coordinates(DIM-1); - std::vector<double> coordinates_neighbor(DIM-1); - std::vector<double> distance(DIM-1); - - domain[step] = SURFACE_FLAG; // set domain for this step, may change later if agent infiltrates - - for (size_t i = 0; i < DIM; i ++) - velocity[i] = location->get_velocity(i); - - // CFL criteria - double length = location->get_characteristic_length(); - double velocity_norm = sqrt(velocity[0]*velocity[0] + velocity[1]*velocity[1]); - - if (velocity_norm > EPS) { - time_step_size = CFL * length / velocity_norm; - } else { - time_step_size = MAXTIMESTEP; - } - - time[step] = time[step-1] + time_step_size; - - // check if agent infiltrates - double h = location->get_water_depth(); - int flag = INFILTRATE_FALSE; - - if (h > EPS) { // only wet cell can infiltrate - double qe = std::min(velocity[2], 0.0); - double pi = fabs(qe) * time_step_size / h; - double trial = (double) (rand()%100) / 100.0; - if (trial < pi) - flag = INFILTRATE_TRUE; - } - - if (flag == INFILTRATE_FALSE) { - - // advect the agent horizontally, which means the last dimension - // (z-direction) is omitted - for (size_t i = 0; i < DIM-1; i ++) { - coordinates[i] = path[step-1][i] + time_step_size * velocity[i]; - distance[i] = location->get_coordinate(i) - coordinates[i]; - } - - // horizontal distance to the centroid of current cell - double distance_norm = sqrt(distance[0]*distance[0]+distance[1]*distance[1]); - - // check if cell has changed - size_t n_cells = mesh->get_n_cells(SURFACE_FLAG); - for (size_t i = 0; i < n_cells; i ++) { - - for (size_t j = 0; j < DIM-1; j ++) { - coordinates_neighbor[j] = mesh->get_cell(i, SURFACE_FLAG)->get_coordinate(j); - distance[j] = coordinates_neighbor[j] - coordinates[j]; - } - - double distance_norm_ = sqrt(distance[0]*distance[0]+distance[1]*distance[1]); - - if (distance_norm > distance_norm_) { - cell_index = i; - distance_norm = distance_norm_; - } - - } - - // check if out of the domain: if the agent travelled thrice the - // maximum length of the cell and still can't find a neighboring - // 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 - domain[step] = SUBSURFACE_FLAG; - coordinates[0] = path[step-1][0]; - coordinates[1] = path[step-1][1]; - cell_index = mesh->get_cell(cell_index, SURFACE_FLAG)->get_coupled_neighbor(); - } - - set_path(coordinates[0], coordinates[1], location->get_coordinate(2), step); - -} - - -void Agent::advance_subsurface(Mesh* mesh, size_t step) -{ - Cell *location = mesh->get_cell(cell_index, SUBSURFACE_FLAG); - - std::vector<double> velocity(DIM); - std::vector<double> coordinates(DIM); - std::vector<double> coordinates_neighbor(DIM); - std::vector<double> distance(DIM); - - domain[step] = SUBSURFACE_FLAG; // set domain for this step, may change later if agent infiltrates - - for (size_t i = 0; i < DIM; i ++) - velocity[i] = location->get_velocity(i); - - // CFL criteria - double length = location->get_characteristic_length(); - double velocity_norm = sqrt(velocity[0] * velocity[0] + velocity[1] * velocity[1] + velocity[2] * velocity[2]); - - if (velocity_norm > EPS) { - time_step_size = CFL * length / velocity_norm; - } else { - time_step_size = MAXTIMESTEP; - } - - time[step] = time[step-1] + time_step_size; - - // advect the agent - for (size_t i = 0; i < DIM; i ++) { - coordinates[i] = path[step-1][i] + time_step_size * velocity[i]; - distance[i] = location->get_coordinate(i) - coordinates[i]; - } - - set_path(coordinates[0], coordinates[1], coordinates[2], step); - - // distance to the centroid of current cell - double distance_norm = sqrt(distance[0] * distance[0] + distance[1] * distance[1] + distance[2] * distance[2]); - - // check if cell has changed - size_t n_cells = mesh->get_n_cells(SUBSURFACE_FLAG); - for (size_t i = 0; i < n_cells; i ++) { - - for (size_t j = 0; j < DIM; j ++) { - coordinates_neighbor[j] = mesh->get_cell(i, SUBSURFACE_FLAG)->get_coordinate(j); - distance[j] = coordinates_neighbor[j] - coordinates[j]; - } - - double distance_norm_ = sqrt(distance[0] * distance[0] + distance[1] * distance[1] + distance[2] * distance[2]); - - if (distance_norm > distance_norm_) { - cell_index = i; - distance_norm = distance_norm_; - } - - } - - // check if out of the domain: if the agent travelled thrice the - // maximum length of the cell and still can't find a neighboring - // cell, it must be advected out of the domain - if (distance_norm > 10.0 * location->get_maximum_length()) { - domain[step] = OUTSIDE_FLAG; - } - - // check if exfiltrates - int neighbor_index = mesh->get_cell(cell_index, SUBSURFACE_FLAG)->get_coupled_neighbor(); - double z_surface; - if (neighbor_index > -1) { - // border to surface detected - z_surface = mesh->get_cell(neighbor_index, SURFACE_FLAG)->get_coordinate(2); - if (coordinates[2] >= z_surface) { - // exfiltration detected - cell_index = (int) neighbor_index; - domain[step] = SURFACE_FLAG; - set_path(coordinates[0], coordinates[1], z_surface, step); - } - } - -} diff --git a/new/agent.hpp b/new/agent.hpp deleted file mode 100644 index f99c2ab..0000000 --- a/new/agent.hpp +++ /dev/null @@ -1,41 +0,0 @@ -// -*- mode: c++ -*- - -#ifndef AGENT_HPP -#define AGENT_HPP - -#include "header.hpp" -#include "cell.hpp" -#include "mesh.hpp" - -class Agent -{ - -public: - - Agent(); - - void set_cell_index(int value); - void set_domain(int value, size_t step); - void set_path(double x, double y, double z, size_t step); - double get_path(size_t step, size_t index); - int get_domain(size_t step); - double get_timestepsize(); - double get_time(size_t step); - - void advance(Mesh* mesh, size_t step); - -private: - - int cell_index; - double time_step_size; - - std::vector<int> domain; - std::vector<double> time; - std::vector<std::vector<double>> path; - - void advance_surface(Mesh *mesh, size_t step); - void advance_subsurface(Mesh *mesh, size_t step); - -}; - -#endif diff --git a/new/header.hpp b/new/header.hpp deleted file mode 100644 index cff75d2..0000000 --- a/new/header.hpp +++ /dev/null @@ -1,56 +0,0 @@ -// -*- mode:c++ -*- - -#ifndef CONSTANTS_HPP -#define CONSTANTS_HPP - -#include <iostream> -#include <vector> -#include <string> // string class -//#include <omp.h> -#include <algorithm> // min -#include <float.h> // DBL_MAX -#include <math.h> // fabs, sqrt -#include <mpi.h> - -#include <H5Cpp.h> -#include <H5File.h> - -#define MAXPRINT 5 -#define MAXTIMESTEP 1.0e5 -#define EPS 1.0e-10 - -#define N_AGENT 100 -#define N_ITERATION 2000 -#define SEED 108 -#define CFL 0.5 - -#define DIM 3 -#define N_NEIGHBORS 5 -#define N_VERTICES 6 -#define N_EDGES 9 - -#define SURFACE_NCOL 4 -#define SUBSURFACE_NCOL 7 - -#define INFILTRATE_TRUE 1 -#define INFILTRATE_FALSE 2 - -#define SURFACE_MESH_FILE "../test/visdump_surface_mesh.h5" -#define SURFACE_DATA_FILE "../test/visdump_surface_data.h5" -#define SUBSURFACE_MESH_FILE "../test/visdump_mesh.h5" -#define SUBSURFACE_DATA_FILE "../test/visdump_data.h5" -#define SURFACE_DATA_GROUP_NAME_H5 "11841" -#define SUBSURFACE_DATA_GROUP_NAME_H5 "11841" - -#define MESH_MIXED_ELEMENTS "/0/Mesh/MixedElements" -#define MESH_NODES "/0/Mesh/Nodes" - -#define AREADIV 0.01 -#define FUZZ 0.1 -#define RADIUS 200.0 - -#define SUBSURFACE_FLAG 1 -#define SURFACE_FLAG 2 -#define OUTSIDE_FLAG 3 - -#endif diff --git a/new/main.cpp b/new/main.cpp deleted file mode 100644 index 1e45aac..0000000 --- a/new/main.cpp +++ /dev/null @@ -1,133 +0,0 @@ -// -*- mode: c++ -*- - -#include "mesh.hpp" -#include "agent.hpp" - -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(); - - 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; - - - for (int n = 0; n < len_surf; n ++) { - c = mesh.get_cell(n, SURFACE_FLAG); - fprintf(mesh_fp, "%f, %f, %f, %d\n", - c->get_coordinate(0), - c->get_coordinate(1), - c->get_coordinate(2), n); - } - - fclose(mesh_fp); - - std::cout << "[~>] Reading the subsurface environment" << std::endl; - - } - - mesh.build_environment(dp_subsurface, SUBSURFACE_FLAG); - dp_subsurface.close(); - - mesh.build_environment(dp_surface, SURFACE_FLAG); - dp_subsurface.close(); - - // Place agents in the domain and advance N_ITERATION steps - - // Bootstrap - 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 = low_bound; i < high_bound; i ++) { - - Agent agent; - - index = rand()%len; - 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); - - agent.set_domain(SUBSURFACE_FLAG, 0); - agent.set_domain(SUBSURFACE_FLAG, 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; - - char output_fname[1024]; - - // 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 = " << agent.get_timestepsize() << " s" << std::endl; - } - agent.advance(&mesh, step); - } - - // write out results - sprintf(output_fname, "results/agent-%06lu.csv", i); - FILE *result_fp = fopen(output_fname, "w"); - - 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", - 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 MPI_SUCCESS; - -} diff --git a/new/makefile b/new/makefile deleted file mode 100644 index 4a91e28..0000000 --- a/new/makefile +++ /dev/null @@ -1,33 +0,0 @@ -PRG=agent-provocateur - -#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 -LIB=-lm -lhdf5 -lhdf5_cpp - -OBJ=$(SRC:.cpp=.o) - -all:$(PRG) - -$(PRG): $(OBJ) - echo "[OK] linking objects" - $(CC) $(LFLAGS) $(CXXFLAGS) $(OBJ) -o $@ $(LIB) - -%.o: %.cpp - echo "[OK] building objects" - $(CC) $(LFLAGS) $(CXXFLAGS) -c $< -o $@ - -run: $(PRG) - ./$(PRG) - -clean: - rm -f $(PRG) $(OBJ) - -cleanall: - make clean - rm -f *.dat; rm -f *.png - -.PHONY: all clean cleanall - diff --git a/new/mesh.cpp b/new/mesh.cpp deleted file mode 100644 index 8112ede..0000000 --- a/new/mesh.cpp +++ /dev/null @@ -1,546 +0,0 @@ -// -*- mode: c++ -*- - -#include "mesh.hpp" - - -// ------------------------------------------------------------------------ -// -- PUBLIC METHODS -// ------------------------------------------------------------------------ - -void Mesh::build_surface_mesh(H5::H5File fp, H5::H5File ep) -{ - - get_mesh_metadata(fp, SURFACE_NCOL); - - handle_surface_node_data(fp); - - build_surface_cells(ep); - -} - -void Mesh::build_subsurface_mesh(H5::H5File fp) -{ - - get_mesh_metadata(fp, SUBSURFACE_NCOL); - - handle_subsurface_node_data(fp); - - build_subsurface_cells(); - -} - -void Mesh::build_environment(H5::H5File ep, int flag) -{ - - std::string groupname; - std::string datasetname; - - // Read in velocities in x direction - - if (flag == SUBSURFACE_FLAG) { - groupname = "darcy_velocity.cell.0"; - datasetname = SUBSURFACE_DATA_GROUP_NAME_H5; - } else if (flag == SURFACE_FLAG) { - groupname = "surface-velocity.cell.0"; - datasetname = SURFACE_DATA_GROUP_NAME_H5; - } - - read_velocity(ep, groupname, datasetname, flag, 0); - - // Read in velocities in y direction - - if (flag == SUBSURFACE_FLAG) { - groupname = "darcy_velocity.cell.1"; - datasetname = SUBSURFACE_DATA_GROUP_NAME_H5; - } else if (flag == SURFACE_FLAG) { - groupname = "surface-velocity.cell.1"; - datasetname = SURFACE_DATA_GROUP_NAME_H5; - } - - read_velocity(ep, groupname, datasetname, flag, 1); - - // Read velocities in z direction - - if (flag == SUBSURFACE_FLAG) { - groupname = "darcy_velocity.cell.2"; - datasetname = SUBSURFACE_DATA_GROUP_NAME_H5; - } else if (flag == SURFACE_FLAG) { - groupname = "surface-surface_subsurface_flux.cell.0"; - datasetname = SURFACE_DATA_GROUP_NAME_H5; - } - - read_velocity(ep, groupname, datasetname, flag, 2); - - - // Surface environment needs adjustments to the vertical velocity - // and needs to read in water depth - - if (flag == SURFACE_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); - - } -#if VERBOSE - else if (flag == SUBSURFACE_FLAG) { - - for (int i = 0; i < n_cells_subsurface; i ++) { - if (i % (n_cells_subsurface / MAXPRINT) == 0) - printf("[~~~~>] Darcy velocity in cell[%d] = { %e, %e, %e } m/s\n", i, - subsurface_cells[i].get_velocity(0), - subsurface_cells[i].get_velocity(1), - subsurface_cells[i].get_velocity(2)); - } - - } -#endif - -} - -// ------------------------------------------------------------------------ -// GETTERS AND SETTERS -// ------------------------------------------------------------------------ - -Mesh::Mesh() -{ - n_cells_surface = 0; - n_cells_subsurface = 0; - high_coordinates.resize(DIM); - low_coordinates.resize(DIM); -} - -int Mesh::get_n_cells(int flag) -{ - int n_cells; - - if (flag == SUBSURFACE_FLAG) { - n_cells = n_cells_subsurface; - } else { - n_cells = n_cells_surface; - } - - return n_cells; -} - -Cell* Mesh::get_cell(int index, int flag) -{ - - if (flag == SUBSURFACE_FLAG) { - return &(subsurface_cells[index]); - } else { - return &(surface_cells[index]); - } -} - -void Mesh::set_high_coordinate(int index, double value) -{ - high_coordinates[index] = value; -} - -double Mesh::get_high_coordinate(int index) -{ - return high_coordinates[index]; -} - -void Mesh::set_low_coordinate(int index, double value) -{ - low_coordinates[index] = value; -} - -double Mesh::get_low_coordinate(int index) -{ - return low_coordinates[index]; -} - -// ------------------------------------------------------------------------ -// -- PRIVATE METHODS -// ------------------------------------------------------------------------ - -void Mesh::handle_surface_node_data(H5::H5File fp) -{ - - H5::DataSet dataset_n = fp.openDataSet(MESH_NODES); - H5::DataType datatype_n = dataset_n.getDataType(); - H5::DataSpace dataspace_n = dataset_n.getSpace(); - - int rank = dataspace_n.getSimpleExtentNdims(); - - 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); - - nodes.resize(dims_n[0]); - - for (size_t i = 0; i < dims_n[0]; i ++) { - nodes[i].resize(dims_n[1]); - for (size_t j = 0; j < dims_n[1]; j ++) - nodes[i][j] = data_n[i * dims_n[1] + j]; - } - -} - -void Mesh::build_surface_cells(H5::H5File ep) -{ - - std::vector<double> coords1(DIM-1); - std::vector<double> coords2(DIM-1); - std::vector<double> coords3(DIM-1); - std::vector<double> coords_avg(DIM-1); - - read_elevation(ep); - - for (int i = 0; i < n_cells_surface; i ++) { - - int n1 = element_node_indices[i][1]; - int n2 = element_node_indices[i][2]; - int n3 = element_node_indices[i][3]; - - std::vector<double> edge_length(DIM); - - for (size_t k = 0; k < DIM; k ++) - edge_length[k] = DBL_MAX; - - - double smallest_distance; - double largest_distance; - - for (size_t j = 0; j < DIM-1; j ++) { - - coords1[j] = nodes[n1][j]; - coords2[j] = nodes[n2][j]; - coords3[j] = nodes[n3][j]; - - coords_avg[j] = (coords1[j] + coords2[j] + coords3[j]) / 3.0; - - edge_length[0] = fabs(coords1[j] - coords_avg[j]); - edge_length[1] = fabs(coords2[j] - coords_avg[j]); - edge_length[2] = fabs(coords3[j] - coords_avg[j]); - - smallest_distance = edge_length[0]; - largest_distance = edge_length[0]; - - for (size_t k = 1; k < DIM; k ++) { - smallest_distance = std::min(smallest_distance, edge_length[k]); - largest_distance = std::max(largest_distance, edge_length[k]); - } - - } - - surface_cells[i].set_characteristic_length(smallest_distance); - surface_cells[i].set_maximum_length(largest_distance); - - for (size_t k = 0; k < DIM-1; k ++) - surface_cells[i].set_coordinate(k, coords_avg[k]); - - surface_cells[i].set_cell_id((int) i); - - // find the coupled cell in the subsurface - std::vector<double> distance(DIM); - int coupled_neighbor = -1; // store current coupled neighbor here - double dz = DBL_MAX; - - for (int k = 0; k < n_cells_subsurface; k ++) { - - for (size_t j = 0; j < DIM; j ++) - distance[j] = subsurface_cells[k].get_coordinate(j) - surface_cells[i].get_coordinate(j); - - // Because the mesh is a stacked mesh, the only candidates for the subsurface - // neighbor are the cells located in the same column. These cells will have - // the same centroid in the xy-plane as the surface cell. The problem reduces - // to finding the nearest cell in the vertical direction in the column. - - double xy_distance = sqrt(distance[0]*distance[0] + distance[1]*distance[1]); - - if (xy_distance <= FUZZ) { - dz = std::min(dz, fabs(distance[2])); - coupled_neighbor = k; - if (dz <= 2.0 * subsurface_cells[k].get_characteristic_length()) - break; - } - - } - - 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 - - } - -} - -void Mesh::handle_subsurface_node_data(H5::H5File fp) -{ - - H5::DataSet dataset_n = fp.openDataSet(MESH_NODES); - H5::DataType datatype_n = dataset_n.getDataType(); - H5::DataSpace dataspace_n = dataset_n.getSpace(); - - int rank = dataspace_n.getSimpleExtentNdims(); - - 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); - - nodes.resize(dims_n[0]); - - for (size_t i = 0; i < dims_n[0]; i ++) { - nodes[i].resize(dims_n[1]); - for (size_t j = 0; j < dims_n[1]; j ++) - nodes[i][j] = data_n[i * dims_n[1] + j]; - } - -} - -void Mesh::build_subsurface_cells() -{ - - std::vector<double> coords1(DIM); - std::vector<double> coords2(DIM); - std::vector<double> coords3(DIM); - std::vector<double> coords4(DIM); - std::vector<double> coords5(DIM); - std::vector<double> coords6(DIM); - - std::vector<double> coords_avg(DIM); - - for (size_t i = 0; i < DIM; i ++) { - set_high_coordinate(i, DBL_MIN); - set_low_coordinate(i, DBL_MAX); - } - - double smallest_distance; - double largest_distance; - - for (int i = 0; i < n_cells_subsurface; i ++) { - - int n1 = element_node_indices[i][1]; - int n2 = element_node_indices[i][2]; - int n3 = element_node_indices[i][3]; - int n4 = element_node_indices[i][4]; - int n5 = element_node_indices[i][5]; - int n6 = element_node_indices[i][6]; - - std::vector<double> edge_length(N_VERTICES); - - for (size_t k = 0; k < N_VERTICES; k ++) - edge_length[k] = DBL_MAX; - - for (size_t j = 0; j < DIM; j++) { - - coords1[j] = nodes[n1][j]; - coords2[j] = nodes[n2][j]; - coords3[j] = nodes[n3][j]; - coords4[j] = nodes[n4][j]; - coords5[j] = nodes[n5][j]; - coords6[j] = nodes[n6][j]; - - coords_avg[j] = (coords1[j] + coords2[j] + coords3[j] + coords4[j] + coords5[j] + coords6[j]) / 6.0; - - // the characteristic length is the radius of the insphere of - // the wedge - edge_length[0] = fabs(coords1[j] - coords_avg[j]); - edge_length[1] = fabs(coords2[j] - coords_avg[j]); - edge_length[2] = fabs(coords3[j] - coords_avg[j]); - edge_length[3] = fabs(coords4[j] - coords_avg[j]); - edge_length[4] = fabs(coords5[j] - coords_avg[j]); - edge_length[5] = fabs(coords6[j] - coords_avg[j]); - - smallest_distance = edge_length[0]; - largest_distance = edge_length[0]; - - for (size_t k = 1; k < N_VERTICES; k ++) { - smallest_distance = std::min(smallest_distance, edge_length[k]); - largest_distance = std::max(largest_distance, edge_length[k]); - } - - } - - subsurface_cells[i].set_characteristic_length(smallest_distance); - subsurface_cells[i].set_maximum_length(largest_distance); - - for (size_t k = 0; k < DIM; k ++) { - - // update lowest and highest corner of mesh - high_coordinates[k] = std::max(coords_avg[k], high_coordinates[k]); - low_coordinates[k] = std::min(coords_avg[k], low_coordinates[k]); - - subsurface_cells[i].set_coordinate(k, coords_avg[k]); - - } - - subsurface_cells[i].set_cell_id((int) i); - - } - -} - -void Mesh::get_mesh_metadata(H5::H5File fp, int flag) -{ - - H5::DataSet dataset_e = fp.openDataSet(MESH_MIXED_ELEMENTS); - H5::DataType datatype_e = dataset_e.getDataType(); - H5::DataSpace dataspace_e = dataset_e.getSpace(); - - int rank = dataspace_e.getSimpleExtentNdims(); - - 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); - - size_t nx_e = flag; - size_t ny_e = dims_e[0] * dims_e[1] / nx_e; - - element_node_indices.resize(ny_e); - - for (size_t i = 0; i < ny_e; i ++) { - element_node_indices[i].resize(nx_e); - for (size_t j = 0; j < nx_e; j ++) - element_node_indices[i][j] = data_e[i * nx_e + j]; - } - - if (flag == SUBSURFACE_NCOL) { - n_cells_subsurface = ny_e; - subsurface_cells.resize(ny_e); - } else if (flag == SURFACE_NCOL) { - n_cells_surface = ny_e; - surface_cells.resize(ny_e); - } - -} - -void Mesh::read_velocity(H5::H5File fp, std::string groupname, std::string datasetname, int flag, int index) -{ - - H5::Group group = fp.openGroup(groupname.c_str()); - H5::DataSet dataset = group.openDataSet(datasetname.c_str()); - H5::DataType datatype = dataset.getDataType(); - H5::DataSpace dataspace = dataset.getSpace(); - - int rank = dataspace.getSimpleExtentNdims(); - - 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); - - if (flag == SUBSURFACE_FLAG) { - for (int i = 0; i < n_cells_subsurface; i ++) - subsurface_cells[i].set_velocity(index, data[i]); - } else if (flag == SURFACE_FLAG) { - for (int i = 0; i < n_cells_surface; i ++) - surface_cells[i].set_velocity(index, data[i]); - } - -} - -void Mesh::read_water_depth(H5::H5File fp) -{ - - H5::Group group = fp.openGroup("surface-ponded_depth.cell.0"); - H5::DataSet dataset = group.openDataSet(SURFACE_DATA_GROUP_NAME_H5); - H5::DataType datatype = dataset.getDataType(); - H5::DataSpace dataspace = dataset.getSpace(); - - int rank = dataspace.getSimpleExtentNdims(); - - 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 - } - -} - -void Mesh::read_elevation(H5::H5File fp) -{ - - H5::Group group = fp.openGroup("surface-elevation.cell.0"); - H5::DataSet dataset = group.openDataSet(SURFACE_DATA_GROUP_NAME_H5); - H5::DataType datatype = dataset.getDataType(); - H5::DataSpace dataspace = dataset.getSpace(); - - int rank = dataspace.getSimpleExtentNdims(); - - 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 - } - -} diff --git a/new/mesh.hpp b/new/mesh.hpp deleted file mode 100644 index 70f7c34..0000000 --- a/new/mesh.hpp +++ /dev/null @@ -1,59 +0,0 @@ -// -*- mode: c++ -*- - -#ifndef MESH_HPP -#define MESH_HPP - -#include "header.hpp" -#include "cell.hpp" - -class Mesh -{ - -public: - - Mesh(); - - void build_surface_mesh(H5::H5File fp, H5::H5File ep); - void build_subsurface_mesh(H5::H5File fp); - - void build_environment(H5::H5File ep, int flag); - void build_subsurface_environment(H5::H5File ep); - - Cell* get_cell(int index, int flag); - - void set_high_coordinate(int index, double value); - double get_high_coordinate(int index); - void set_low_coordinate(int index, double value); - double get_low_coordinate(int index); - - int get_n_cells(int flag); - -private: - - int n_cells_subsurface; - int n_cells_surface; - - std::vector<Cell> subsurface_cells; - std::vector<Cell> surface_cells; - - std::vector<double> high_coordinates; - std::vector<double> low_coordinates; - - std::vector<std::vector<double>> nodes; - std::vector<std::vector<int>> element_node_indices; - //std::vector<std::vector<std::vector<hsize_t>>> edge_ids; - - void handle_subsurface_node_data(H5::H5File fp); - void build_subsurface_cells(); - - void handle_surface_node_data(H5::H5File fp); - void build_surface_cells(H5::H5File fp); - - void get_mesh_metadata(H5::H5File fp, int ncol); - void read_velocity(H5::H5File fp, std::string groupname, std::string datasetname, int flag, int index); - void read_water_depth(H5::H5File fp); - void read_elevation(H5::H5File fp); - -}; - -#endif diff --git a/test.cpp b/test.cpp deleted file mode 100755 index e31010a..0000000 --- a/test.cpp +++ /dev/null @@ -1,140 +0,0 @@ -#include <H5Cpp.h> -#include <H5File.h> - -#include <iostream> - -#define MESH_FILE_NAME "test1/visdump_mesh.h5" -#define DATA_FILE_NAME "test1/visdump_data.h5" - -int main(void) -{ - std::cout << "Hello, world!" << std::endl; - - H5::H5File fp = H5::H5File(MESH_FILE_NAME, H5F_ACC_RDONLY); - H5::DataSet dataset_e = fp.openDataSet("/0/Mesh/MixedElements"); - H5::DataType datatype_e = dataset_e.getDataType(); - H5::DataSpace dataspace_e = dataset_e.getSpace(); - - int rank = dataspace_e.getSimpleExtentNdims(); - - std::cout << "Rank of element node map: " << rank << std::endl; - - hsize_t *dims_e = new hsize_t[rank]; - dataspace_e.getSimpleExtentDims(dims_e); - - for (int i = 0; i < (int) rank; i ++) - std::cout << "dims[" << i << "] of element node map: " << dims_e[i] << std::endl; - - int *data_e = new int[dims_e[0]*dims_e[1]]; - - size_t nx_e = 7; - size_t ny_e = dims_e[0] * dims_e[1] / 7; - - int **md_e = new int*[ny_e]; - - for (size_t i = 0; i < ny_e; i ++) - md_e[i] = data_e + i * nx_e; - - dataset_e.read(data_e, H5::PredType::STD_I32LE); - - size_t maxprint = 10; - - if (ny_e < maxprint) - maxprint = ny_e; - - for (size_t i = 0; i < maxprint; i ++) - { - printf("[ok]"); - for (size_t j = 0; j < nx_e; j ++) - printf(" %d ", md_e[i][j]); - printf("\n"); - } - - H5::DataSet dataset_n = fp.openDataSet("/0/Mesh/Nodes"); - H5::DataType datatype_n = dataset_n.getDataType(); - H5::DataSpace dataspace_n = dataset_n.getSpace(); - - rank = dataspace_n.getSimpleExtentNdims(); - - hsize_t *dims_n = new hsize_t[rank]; - dataspace_n.getSimpleExtentDims(dims_n); - - for (int i = 0; i < (int) rank; i ++) - std::cout << "dims[" << i << "] of nodes: " << dims_n[i] << std::endl; - - double *data_n = new double[dims_n[0]*dims_n[1]]; - double **md_n = new double*[dims_n[0]]; - - for (size_t i = 0; i < dims_n[0]; i ++) - md_n[i] = data_n + i * dims_n[1]; - - dataset_n.read(data_n, H5::PredType::IEEE_F64LE); - - maxprint = 10; - if (dims_n[0] < 10) - maxprint = dims_n[0]; - - for (size_t i = 0; i < maxprint; i ++) - { - printf("[ok]"); - for (size_t j = 0; j < dims_n[1]; j ++) - printf(" %f ", md_n[i][j]); - printf("\n"); - } - - double *coords1 = new double[dims_n[1]]; - double *coords2 = new double[dims_n[1]]; - double *coords3 = new double[dims_n[1]]; - double *coords4 = new double[dims_n[1]]; - double *coords5 = new double[dims_n[1]]; - double *coords6 = new double[dims_n[1]]; - - double *coords_avg = new double[3 * ny_e]; - - for (size_t i = 0; i < ny_e; i ++) { - - int n1 = md_e[i][1]; - int n2 = md_e[i][2]; - int n3 = md_e[i][3]; - int n4 = md_e[i][4]; - int n5 = md_e[i][5]; - int n6 = md_e[i][6]; - - for (size_t j = 0; j < dims_n[1]; j++) - { - - coords1[j] = md_n[n1][j]; - coords2[j] = md_n[n2][j]; - coords3[j] = md_n[n3][j]; - coords4[j] = md_n[n4][j]; - coords5[j] = md_n[n5][j]; - coords6[j] = md_n[n6][j]; - - coords_avg[dims_n[1] * i + j] = (coords1[j] + coords2[j] + coords3[j] + coords4[j] + coords5[j] + coords6[j]) / 6.0; - - } - - } - -#define GNUPLOT "gnuplot -persist" - - FILE *gp; - gp = popen(GNUPLOT, "w"); - fprintf(gp, "splot '-' w p lt -1 pt 6 title 'e and z'\n"); - - for (size_t i = 0; i < ny_e; i ++) - fprintf(gp, "%f %f %f\n", coords_avg[dims_n[1] * i], coords_avg[dims_n[1] * i + 1], coords_avg[dims_n[1] * i + 2]); - fprintf(gp, "\033"); - - fflush(gp); - fclose(gp); - - fp.close(); - - int ns[9]; - hsize_t **szudzik_ids = new hsize_t*[ny_e]; - - - - return 0; -} diff --git a/tools.cpp b/tools.cpp deleted file mode 100644 index 4a62950..0000000 --- a/tools.cpp +++ /dev/null @@ -1,11 +0,0 @@ -// -*- mode: c++; c-default-style: "linux" -*- - -#include "header.hpp" - -hsize_t szudzik(int n1, int n2) -{ - hsize_t k1 = (hsize_t) ((n1 < n2) ? n1 : n2); - hsize_t k2 = (hsize_t) ((n1 < n2) ? n2 : n1); - - return k2 * k2 + k1; -} diff --git a/tools.hpp b/tools.hpp deleted file mode 100644 index 9bea0f0..0000000 --- a/tools.hpp +++ /dev/null @@ -1,8 +0,0 @@ -// -*- mode: c++; c-default-style: "linux" -*- - -#ifndef TOOLS_HPP -#define TOOLS_HPP - -hsize_t szudzik(int n1, int n2); - -#endif |