From 1d1c5b27100ce499d1061397a2201eb7ca3bcedf Mon Sep 17 00:00:00 2001 From: Ilhan Özgen Xian Date: Fri, 21 May 2021 14:14:56 -0700 Subject: fix bugs in in- and exfiltration --- new/agent.cpp | 155 ++++++++++++++++++++++++++++++++++++++++++++------- new/agent.hpp | 3 +- new/cell.hpp | 2 +- new/header.hpp | 9 ++- new/main.cpp | 47 ++++++++++------ new/merge-results.sh | 8 ++- new/mesh.cpp | 109 +++++++++++++++++++++++++----------- new/mesh.hpp | 5 +- 8 files changed, 261 insertions(+), 77 deletions(-) mode change 100644 => 100755 new/merge-results.sh diff --git a/new/agent.cpp b/new/agent.cpp index a6bb9e9..fc8c1f1 100644 --- a/new/agent.cpp +++ b/new/agent.cpp @@ -2,6 +2,33 @@ #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); @@ -38,32 +65,108 @@ double Agent::get_timestepsize() return time_step_size; } -// std::vector> Agent::get_path() -// { -// return path; -// } - double Agent::get_path(size_t step, size_t index) { return path[step][index]; } -void Agent::advance(Mesh* mesh, size_t step) +int Agent::get_domain(size_t step) { + return domain[step]; +} - switch(domain[step-1]) { - 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]; +// ------------------------------------------------------------------------ +// -- PRIVATE METHODS +// ------------------------------------------------------------------------ + +void Agent::advance_surface(Mesh* mesh, size_t step) +{ + Cell *location = mesh->get_cell(cell_index, SURFACE_FLAG); + + std::vector velocity(DIM); + std::vector coordinates(DIM-1); + std::vector coordinates_neighbor(DIM-1); + std::vector 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; - break; + + + } 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); @@ -73,13 +176,13 @@ void Agent::advance_subsurface(Mesh* mesh, size_t step) std::vector coordinates_neighbor(DIM); std::vector distance(DIM); - domain[step] = SUBSURFACE_FLAG; // set domain for this step + 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 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) { @@ -88,7 +191,6 @@ void Agent::advance_subsurface(Mesh* mesh, size_t step) time_step_size = MAXTIMESTEP; } - // time_step_size = 1e6; time[step] = time[step-1] + time_step_size; // advect the agent @@ -119,11 +221,26 @@ void Agent::advance_subsurface(Mesh* mesh, size_t step) } } - + // 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()) + 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 index 19d3fb4..f99c2ab 100644 --- a/new/agent.hpp +++ b/new/agent.hpp @@ -18,9 +18,9 @@ public: 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); - //std::vector> get_path(); void advance(Mesh* mesh, size_t step); @@ -33,6 +33,7 @@ private: std::vector time; std::vector> path; + void advance_surface(Mesh *mesh, size_t step); void advance_subsurface(Mesh *mesh, size_t step); }; diff --git a/new/cell.hpp b/new/cell.hpp index 3981013..7c2cd83 100644 --- a/new/cell.hpp +++ b/new/cell.hpp @@ -41,7 +41,7 @@ public: private: int cell_id; - hsize_t coupled_neighbor; // id of the neighbor cell that's in the + int coupled_neighbor; // id of the neighbor cell that's in the // other domain, for example the surface // neighbor of a subsurface cell diff --git a/new/header.hpp b/new/header.hpp index b8ef032..9acc489 100644 --- a/new/header.hpp +++ b/new/header.hpp @@ -6,7 +6,7 @@ #include #include #include // string class -#include +//#include #include // min #include // DBL_MAX #include // fabs, sqrt @@ -14,8 +14,8 @@ #include #include -#define MAXPRINT 10 -#define MAXTIMESTEP 10.0 +#define MAXPRINT 5 +#define MAXTIMESTEP 1.0e5 #define EPS 1.0e-10 #define N_AGENT 100 @@ -31,6 +31,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" diff --git a/new/main.cpp b/new/main.cpp index 83aaf13..5363f89 100644 --- a/new/main.cpp +++ b/new/main.cpp @@ -9,29 +9,41 @@ int main(int argc, char* argv[]) 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); + mesh.build_subsurface_mesh(fp_subsurface); fp_subsurface.close(); - H5::H5File fp_surface = H5::H5File(SURFACE_MESH_FILE, H5F_ACC_RDONLY); - mesh.build_surface_mesh(fp_surface); + mesh.build_surface_mesh(fp_surface, dp_surface); fp_surface.close(); - - // FILE *mesh_fp = fopen("mesh.csv", "w"); - // fprintf(mesh_fp, "x (m), y (m), z (m), index\n"); -// Cell *c; -// for (int n = 0; n < mesh.get_n_cells(); n ++) { -// c = mesh.get_cell(n); -// fprintf(mesh_fp, "%f, %f, %f, %d\n", c->get_coordinate(0), -// c->get_coordinate(1), c->get_coordinate(2), n); -// } -// fclose(mesh_fp); + FILE *mesh_fp = fopen("mesh.csv", "w"); + fprintf(mesh_fp, "x (m), y (m), z (m), index\n"); + + Cell *c; + int len_surf = mesh.get_n_cells(SURFACE_FLAG); + for (int n = 0; n < len_surf; n ++) { + std::cout << n << "/" << len_surf << " "; + c = mesh.get_cell(n, SURFACE_FLAG); + std::cout << "read in c "; + fprintf(mesh_fp, "%f, %f, %f, %d\n", + c->get_coordinate(0), + c->get_coordinate(1), + c->get_coordinate(2), n); + std::cout << "done reading." << std::endl; + } - H5::H5File dp_subsurface = H5::H5File(SUBSURFACE_DATA_FILE, H5F_ACC_RDONLY); + std::cout << "closing file..." << std::endl; + + fclose(mesh_fp); + + std::cout << "reading the environment" << std::endl; + mesh.build_environment(dp_subsurface, SUBSURFACE_FLAG); dp_subsurface.close(); - H5::H5File dp_surface = H5::H5File(SURFACE_DATA_FILE, H5F_ACC_RDONLY); mesh.build_environment(dp_surface, SURFACE_FLAG); dp_subsurface.close(); @@ -76,13 +88,14 @@ int main(int argc, char* argv[]) 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)\n"); + 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\n", + fprintf(result_fp, "%f, %f, %f, %f, %d\n", agents[i].get_time(step), agents[i].get_path(step, 0), agents[i].get_path(step, 1), - agents[i].get_path(step, 2)); + agents[i].get_path(step, 2), + agents[i].get_domain(step)); } fclose(result_fp); diff --git a/new/merge-results.sh b/new/merge-results.sh old mode 100644 new mode 100755 index ed89062..5459897 --- a/new/merge-results.sh +++ b/new/merge-results.sh @@ -1,5 +1,11 @@ #!/usr/bin/env bash +echo "[ok] remove existing aggregation file" +rm -f results/agent-total.csv +echo "[ok] cat everything into agent-total.csv" cat results/*.csv > results/agent-total.csv +echo "[ok] cleaning up the file" gsed -i '/^time/d' results/agent-total.csv -gsed -i '1i time (s), x (m), y (m), z (m)' results/agent-total.csv +echo "[ok] inserting header" +gsed -i '1i time (s), x (m), y (m), z (m), domain' results/agent-total.csv +echo "[ok] done" diff --git a/new/mesh.cpp b/new/mesh.cpp index 6180ea8..6ae6c22 100644 --- a/new/mesh.cpp +++ b/new/mesh.cpp @@ -7,7 +7,7 @@ // -- PUBLIC METHODS // ------------------------------------------------------------------------ -void Mesh::build_surface_mesh(H5::H5File fp) +void Mesh::build_surface_mesh(H5::H5File fp, H5::H5File ep) { std::cout << "[~>] Building the surface mesh" << std::endl; @@ -15,7 +15,7 @@ void Mesh::build_surface_mesh(H5::H5File fp) handle_surface_node_data(fp); - build_surface_cells(); + build_surface_cells(ep); } @@ -92,7 +92,7 @@ void Mesh::build_environment(H5::H5File ep, int flag) surface_cells[i].get_velocity(2)); } - // Read in water depth (TODO) + // Read in water depth read_water_depth(ep); } else if (flag == SUBSURFACE_FLAG) { @@ -128,7 +128,7 @@ int Mesh::get_n_cells(int flag) if (flag == SUBSURFACE_FLAG) { n_cells = n_cells_subsurface; } else { - n_cells = n_cells_subsurface; + n_cells = n_cells_surface; } return n_cells; @@ -164,8 +164,6 @@ double Mesh::get_low_coordinate(int index) return low_coordinates[index]; } - - // ------------------------------------------------------------------------ // -- PRIVATE METHODS // ------------------------------------------------------------------------ @@ -198,16 +196,20 @@ void Mesh::handle_surface_node_data(H5::H5File fp) } -void Mesh::build_surface_cells() +void Mesh::build_surface_cells(H5::H5File ep) { - std::vector coords1(DIM); - std::vector coords2(DIM); - std::vector coords3(DIM); - std::vector coords_avg(DIM); + std::vector coords1(DIM-1); + std::vector coords2(DIM-1); + std::vector coords3(DIM-1); + std::vector coords_avg(DIM-1); + + read_elevation(ep); for (int i = 0; i < n_cells_surface; i ++) { + //std::cout << "build cell " << i << std::endl; + int n1 = element_node_indices[i][1]; int n2 = element_node_indices[i][2]; int n3 = element_node_indices[i][3]; @@ -217,7 +219,11 @@ void Mesh::build_surface_cells() for (size_t k = 0; k < DIM; k ++) edge_length[k] = DBL_MAX; - for (size_t j = 0; j < DIM; j ++) { + + 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]; @@ -228,40 +234,45 @@ void Mesh::build_surface_cells() 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]); - - } - double smallest_distance = edge_length[0]; - double 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]); + 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; k ++) + 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 distance(DIM); - hsize_t coupled_neighbor = -1; // store current coupled neighbor here - double dz = DBL_MAX; // find the nearest cell in the vertical direction + 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) - coords_avg[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 = (hsize_t) k; + coupled_neighbor = k; if (dz <= 2.0 * subsurface_cells[k].get_characteristic_length()) break; } @@ -326,6 +337,9 @@ void Mesh::build_subsurface_cells() 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]; @@ -359,15 +373,15 @@ void Mesh::build_subsurface_cells() 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]); - - } - double smallest_distance = edge_length[0]; - double largest_distance = edge_length[0]; - - for (size_t k = 1; k < N_VERTICES; k ++) { + 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); @@ -487,3 +501,32 @@ void Mesh::read_water_depth(H5::H5File fp) } } + +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); + + for (int i = 0; i < rank; i ++) + std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl; + + 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 (i%(n_cells_surface/5)==0) { + printf("[~~~~>] Surface elevation in cell[%d] = %e m\n", i, + surface_cells[i].get_coordinate(2)); + } + } + +} diff --git a/new/mesh.hpp b/new/mesh.hpp index dc26ac1..70f7c34 100644 --- a/new/mesh.hpp +++ b/new/mesh.hpp @@ -13,7 +13,7 @@ public: Mesh(); - void build_surface_mesh(H5::H5File fp); + 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); @@ -47,11 +47,12 @@ private: void build_subsurface_cells(); void handle_surface_node_data(H5::H5File fp); - void build_surface_cells(); + 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); }; -- cgit v1.2.3-13-gbd6f