path: root/agent.cpp

// -*- 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
// ------------------------------------------------------------------------

void Agent::init(int n_iteration)
{
  time.resize(n_iteration);
  domain.resize(n_iteration);
  path.resize(n_iteration);
  for (int 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 > 15.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 > 15.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);
    }
  }
  
}