Skip to content

Commit

Permalink
r.sim: dynamic allocation of walkers
Browse files Browse the repository at this point in the history
  • Loading branch information
@petrasovaa @neteler
petrasovaa authored and neteler committed Nov 18, 2020
1 parent e4147fb commit 7a71955
Show file tree
Hide file tree
Showing 5 changed files with 88 additions and 50 deletions.
89 changes: 45 additions & 44 deletions raster/r.sim/simlib/hydro.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,8 @@
*/

struct _points points;
struct point2D;
struct point3D;

char *elevin;
char *dxin;
Expand Down Expand Up @@ -85,10 +87,11 @@ double **gama, **gammas, **si, **inf, **sigma;
float **dc, **tau, **er, **ct, **trap;
float **dif;

/* suspected BUG below: fist array subscripts go from 1 to MAXW
* correct: from 0 to MAXW - 1, e.g for (lw = 0; lw < MAXW; lw++) */
double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3];
int iflag[MAXW];

/* int iflag[MAXW];*/
struct point3D *w;
struct point3D *stack;
struct point2D *vavg;

double hbeta;
double hhmax, sisum, vmean;
Expand Down Expand Up @@ -197,24 +200,22 @@ void main_loop(void)
/*G_debug(2, " gen,gen2,wei,wei2,mgen3,nmult: %f %f %f %f %d %d",gen,gen2,wei,wei2,mgen3,nmult);
*/
for (iw = 1; iw <= mgen + 1; iw++) { /* assign walkers */

if (lw >= MAXW) /* max valid value is MAXW - 1, not MAXW */
G_fatal_error(_("nwalk (%d) > maxw (%d)!"), lw, MAXW);

w[lw][0] = x + stepx * (simwe_rand() - 0.5);
w[lw][1] = y + stepy * (simwe_rand() - 0.5);
w[lw][2] = wei;

walkwe += w[lw][2];
vavg[lw][0] = v1[k][l];
vavg[lw][1] = v2[k][l];
if (w[lw][0] >= xmin && w[lw][1] >= ymin &&
w[lw][0] <= xmax && w[lw][1] <= ymax) {
w[lw].x = x + stepx * (simwe_rand() - 0.5);
w[lw].y = y + stepy * (simwe_rand() - 0.5);
w[lw].m = wei;

walkwe += w[lw].m;
vavg[lw].x = v1[k][l];
vavg[lw].y = v2[k][l];
/* deactivated as unused, what was iflag for?
if (w[lw].x >= xmin && w[lw].y >= ymin &&
w[lw].x <= xmax && w[lw].y <= ymax) {
iflag[lw] = 0;
}
else {
iflag[lw] = 1;
}
*/
lw++;
}
} /*DEFined area */
Expand Down Expand Up @@ -279,15 +280,15 @@ void main_loop(void)
#else
for (lw = 0; lw < nwalk; lw++) {
#endif
if (w[lw][2] > EPS) { /* check the walker weight */
if (w[lw].m > EPS) { /* check the walker weight */
++nwalka;
l = (int)((w[lw][0] + stxm) / stepx) - mx - 1;
k = (int)((w[lw][1] + stym) / stepy) - my - 1;
l = (int)((w[lw].x + stxm) / stepx) - mx - 1;
k = (int)((w[lw].y + stym) / stepy) - my - 1;

if (l > mx - 1 || k > my - 1 || k < 0 || l < 0) {

G_debug(2, " k,l=%d,%d", k, l);
printf(" lw,w=%d %f %f", lw, w[lw][1], w[lw][2]);
printf(" lw,w=%d %f %f", lw, w[lw].y, w[lw].m);
G_debug(2, " stxym=%f %f", stxm, stym);
printf(" step=%f %f", stepx, stepy);
G_debug(2, " m=%d %d", my, mx);
Expand All @@ -299,20 +300,20 @@ void main_loop(void)
if (infil != NULL) { /* infiltration part */
if (inf[k][l] - si[k][l] > 0.) {

decr = pow(addac * w[lw][2], 3. / 5.); /* decreasing factor in m */
decr = pow(addac * w[lw].m, 3. / 5.); /* decreasing factor in m */
if (inf[k][l] > decr) {
inf[k][l] -= decr; /* decrease infilt. in cell and eliminate the walker */
w[lw][2] = 0.;
w[lw].m = 0.;
}
else {
w[lw][2] -= pow(inf[k][l], 5. / 3.) / addac; /* use just proportional part of the walker weight */
w[lw].m -= pow(inf[k][l], 5. / 3.) / addac; /* use just proportional part of the walker weight */
inf[k][l] = 0.;

}
}
}

gama[k][l] += (addac * w[lw][2]); /* add walker weigh to water depth or conc. */
gama[k][l] += (addac * w[lw].m); /* add walker weigh to water depth or conc. */

d1 = gama[k][l] * conn;
#if defined(_OPENMP)
Expand All @@ -325,8 +326,8 @@ void main_loop(void)

if (hhc > hhmax && wdepth == NULL) { /* increased diffusion if w.depth > hhmax */
dif[k][l] = (halpha + 1) * deldif;
velx = vavg[lw][0];
vely = vavg[lw][1];
velx = vavg[lw].x;
vely = vavg[lw].y;
}
else {
dif[k][l] = deldif;
Expand All @@ -345,29 +346,29 @@ void main_loop(void)
}
}

w[lw][0] += (velx + dif[k][l] * gaux); /* move the walker */
w[lw][1] += (vely + dif[k][l] * gauy);
w[lw].x += (velx + dif[k][l] * gaux); /* move the walker */
w[lw].y += (vely + dif[k][l] * gauy);

if (hhc > hhmax && wdepth == NULL) {
vavg[lw][0] = hbeta * (vavg[lw][0] + v1[k][l]);
vavg[lw][1] = hbeta * (vavg[lw][1] + v2[k][l]);
vavg[lw].x = hbeta * (vavg[lw].x + v1[k][l]);
vavg[lw].y = hbeta * (vavg[lw].y + v2[k][l]);
}

if (w[lw][0] <= xmin || w[lw][1] <= ymin || w[lw][0]
>= xmax || w[lw][1] >= ymax) {
w[lw][2] = 1e-10; /* eliminate walker if it is out of area */
if (w[lw].x <= xmin || w[lw].y <= ymin || w[lw].x
>= xmax || w[lw].y >= ymax) {
w[lw].m = 1e-10; /* eliminate walker if it is out of area */
}
else {
if (wdepth != NULL) {
l = (int)((w[lw][0] + stxm) / stepx) - mx - 1;
k = (int)((w[lw][1] + stym) / stepy) - my - 1;
w[lw][2] *= sigma[k][l];
l = (int)((w[lw].x + stxm) / stepx) - mx - 1;
k = (int)((w[lw].y + stym) / stepy) - my - 1;
w[lw].m *= sigma[k][l];
}

} /* else */
} /*DEFined area */
else {
w[lw][2] = 1e-10; /* eliminate walker if it is out of area */
w[lw].m = 1e-10; /* eliminate walker if it is out of area */
}
}
} /* lw loop */
Expand All @@ -380,19 +381,19 @@ void main_loop(void)

for (lw = 0; lw < nwalk; lw++) {
/* Compute the elevation raster map index */
l = (int)((w[lw][0] + stxm) / stepx) - mx - 1;
k = (int)((w[lw][1] + stym) / stepy) - my - 1;
l = (int)((w[lw].x + stxm) / stepx) - mx - 1;
k = (int)((w[lw].y + stym) / stepy) - my - 1;

/* Check for correct elevation raster map index */
if(l < 0 || l >= mx || k < 0 || k >= my)
continue;

if (w[lw][2] > EPS && zz[k][l] != UNDEF) {
if (w[lw].m > EPS && zz[k][l] != UNDEF) {

/* Save the 3d position of the walker */
stack[nstack][0] = mixx / conv + w[lw][0] / conv;
stack[nstack][1] = miyy / conv + w[lw][1] / conv;
stack[nstack][2] = zz[k][l];
stack[nstack].x = mixx / conv + w[lw].x / conv;
stack[nstack].y = miyy / conv + w[lw].y / conv;
stack[nstack].m = zz[k][l];

nstack++;
}
Expand Down
16 changes: 15 additions & 1 deletion raster/r.sim/simlib/input.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -268,6 +268,16 @@ void init_grids_sediment()
erod(si);
}

void alloc_walkers(int max_walkers)
{
G_debug(1, "beginning memory allocation for walkers");

w = (struct point3D *)G_calloc(max_walkers, sizeof(struct point3D));
vavg = (struct point2D *)G_calloc(max_walkers, sizeof(struct point2D));
if (outwalk != NULL)
stack = (struct point3D *)G_calloc(max_walkers, sizeof(struct point3D));
}

/* ************************************************************** */
/* GRASS input procedures, allocations */
/* *************************************************************** */
Expand All @@ -283,6 +293,7 @@ void init_grids_sediment()
int input_data(void)
{
int rows = my, cols = mx; /* my and mx are global variables */
int max_walkers;
double unitconv = 0.000000278; /* mm/hr to m/s */
int if_rain = 0;

Expand Down Expand Up @@ -359,7 +370,10 @@ int input_data(void)
gama = read_double_raster_map(rows, cols, wdepth, 1.0);
copy_matrix_undef_double_to_float_values(rows, cols, gama, zz);
}

/* allocate walkers */
max_walkers = maxwa + mx * my;
alloc_walkers(max_walkers);

/* Array for gradient checking */
slope = create_double_matrix(rows, cols, 0.0);

Expand Down
15 changes: 12 additions & 3 deletions raster/r.sim/simlib/output.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,14 @@

#include <grass/waterglobs.h>

void free_walkers()
{
G_free(w);
G_free(vavg);
if (outwalk != NULL)
G_free(stack);
}


static void output_walker_as_vector(int tt_minutes, int ndigit, struct TimeStamp *timestamp);

Expand Down Expand Up @@ -49,9 +57,9 @@ void output_walker_as_vector(int tt_minutes, int ndigit, struct TimeStamp *times
Cats = Vect_new_cats_struct();

for (i = 0; i < nstack; i++) {
x = stack[i][0];
y = stack[i][1];
z = stack[i][2];
x = stack[i].x;
y = stack[i].y;
z = stack[i].m;

Vect_reset_line(Points);
Vect_reset_cats(Cats);
Expand Down Expand Up @@ -759,6 +767,7 @@ int output_et()
Rast_free_colors(&colors);
/* } */
}
free_walkers();

return 1;
}
1 change: 1 addition & 0 deletions raster/r.sim/simlib/simlib.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,7 @@ int input_data(void);
int grad_check(void);
void main_loop(void);
int output_data(int, double);
void free_walkers();

struct options
{
Expand Down
17 changes: 15 additions & 2 deletions raster/r.sim/simlib/waterglobs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,18 @@ struct _points
int is_open; /* Set to 1 if open, 0 if closed */
};

struct point2D
{
double x;
double y;
};
struct point3D
{
double x;
double y;
double m;
};

extern struct _points points;
extern void erod(double **);
extern int output_et(void);
Expand Down Expand Up @@ -86,8 +98,9 @@ extern double **gama, **gammas, **si, **inf, **sigma;
extern float **dc, **tau, **er, **ct, **trap;
extern float **dif;

extern double vavg[MAXW][2], stack[MAXW][3], w[MAXW][3];
extern int iflag[MAXW];
extern struct point3D *w;
extern struct point3D *stack;
extern struct point2D *vavg;

extern double hbeta;
extern double hhmax, sisum, vmean;
Expand Down

0 comments on commit 7a71955

Please sign in to comment.