1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
// -*- mode: c++ -*-
#include "agent.hpp"
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;
}
// 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)
{
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];
domain[step] = OUTSIDE_FLAG;
break;
}
}
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
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_size = 1e6;
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;
}
|
