fix bugs in in- and exfiltration

8 files changed, 261 insertions, 77 deletions

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<std::vector<double>> 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<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;
-    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<double> coordinates_neighbor(DIM);
   std::vector<double> 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<std::vector<double>> get_path();
 
   void advance(Mesh* mesh, size_t step);
 
@@ -33,6 +33,7 @@ private:
   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);
   
 };
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 <iostream>
 #include <vector>
 #include <string>    // string class
-#include <omp.h>
+//#include <omp.h>
 #include <algorithm> // min
 #include <float.h>   // DBL_MAX
 #include <math.h>    // fabs, sqrt
@@ -14,8 +14,8 @@
 #include <H5Cpp.h>
 #include <H5File.h>
 
-#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
index ed89062..5459897 100644..100755
--- 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<double> coords1(DIM);
-  std::vector<double> coords2(DIM);
-  std::vector<double> coords3(DIM);
-  std::vector<double> coords_avg(DIM);
+  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 ++) {
 
+    //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<double> 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);
     
 };