diff options
Diffstat (limited to 'new/agent.cpp')
| -rw-r--r-- | new/agent.cpp | 155 |
1 files changed, 136 insertions, 19 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); + } + } } |