1
2 /****************************************************************************
3 *
4 * MODULE: simwe library
5 * AUTHOR(S): Helena Mitasova, Jaro Hofierka, Lubos Mitas:
6 * PURPOSE: Hydrologic and sediment transport simulation (SIMWE)
7 *
8 * COPYRIGHT: (C) 2002 by the GRASS Development Team
9 *
10 * This program is free software under the GNU General Public
11 * License (>=v2). Read the file COPYING that comes with GRASS
12 * for details.
13 *
14 *****************************************************************************/
15
16 /* hydro.c (simlib), 20.nov.2002, JH */
17
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <math.h>
21 #include <grass/gis.h>
22 #include <grass/bitmap.h>
23 #include <grass/linkm.h>
24 #include <grass/glocale.h>
25
26 #include <grass/waterglobs.h>
27 #include <grass/simlib.h>
28 /*
29 * Soeren 8. Mar 2011 TODO:
30 * Put all these global variables into several meaningful structures and
31 * document use and purpose.
32 *
33 */
34
35 struct _points points;
36 struct point2D;
37 struct point3D;
38
39 char *elevin;
40 char *dxin;
41 char *dyin;
42 char *rain;
43 char *infil;
44 char *traps;
45 char *manin;
46 char *depth;
47 char *disch;
48 char *err;
49 char *outwalk;
50 char *observation;
51 char *logfile;
52 char *mapset;
53 char *mscale;
54 char *tserie;
55
56 char *wdepth;
57 char *detin;
58 char *tranin;
59 char *tauin;
60 char *tc;
61 char *et;
62 char *conc;
63 char *flux;
64 char *erdep;
65
66 char *rainval;
67 char *maninval;
68 char *infilval;
69
70 struct seed seed;
71
72 double xmin, ymin, xmax, ymax;
73 double mayy, miyy, maxx, mixx;
74 int mx, my;
75 int mx2, my2;
76
77 double bxmi, bymi, bxma, byma, bresx, bresy;
78 int maxwab;
79 double step, conv;
80
81 double frac;
82 double bxmi, bymi;
83
84 float **zz, **cchez;
85 double **v1, **v2, **slope;
86 double **gama, **gammas, **si, **inf, **sigma;
87 float **dc, **tau, **er, **ct, **trap;
88 float **dif;
89
90
91 /* int iflag[MAXW];*/
92 struct point3D *w;
93 struct point3D *stack;
94 struct point2D *vavg;
95
96 double hbeta;
97 double hhmax, sisum, vmean;
98 double infsum, infmean;
99 int maxw, maxwa, nwalk;
100 double rwalk, bresx, bresy, xrand, yrand;
101 double stepx, stepy, xp0, yp0;
102 double chmean, si0, deltap, deldif, cch, hhc, halpha;
103 double eps;
104 int maxwab, nstack;
105 int iterout, mx2o, my2o;
106 int miter, nwalka;
107 double timec;
108 int ts, timesec;
109
110 double rain_val;
111 double manin_val;
112 double infil_val;
113
114 struct History history; /* holds meta-data (title, comments,..) */
115
116 /* **************************************************** */
117 /* create walker representation of si */
118 /* ******************************************************** */
119 /* .......... iblock loop */
120
main_loop(void)121 void main_loop(void)
122 {
123
124 int i, ii, l, k;
125 int icoub, nmult;
126 int iw, iblock, lw;
127 int itime, iter1;
128 int nfiterh, nfiterw;
129 int mgen, mgen2, mgen3;
130 int nblock;
131 int icfl;
132 int mitfac;
133 /* int mitfac, p; */
134 double x, y;
135 double velx, vely, stxm, stym;
136 double factor, conn, gaux, gauy;
137 double d1, addac, decr;
138 double barea, sarea, walkwe;
139 double gen, gen2, wei2, wei3, wei, weifac;
140 float eff;
141
142 nblock = 1;
143 icoub = 0;
144 icfl = 0;
145 nstack = 0;
146
147 if (maxwa > (MAXW - mx * my)) {
148 mitfac = maxwa / (MAXW - mx * my);
149 nblock = mitfac + 1;
150 maxwa = maxwa / nblock;
151 }
152
153 /* Create the observation points */
154 create_observation_points();
155
156 G_debug(2, " maxwa, nblock %d %d", maxwa, nblock);
157
158 for (iblock = 1; iblock <= nblock; iblock++) {
159 ++icoub;
160
161 lw = 0;
162 walkwe = 0.;
163 barea = stepx * stepy;
164 sarea = bresx * bresy;
165 G_debug(2, " barea,sarea,rwalk,sisum: %f %f %f %f", barea, sarea,
166 rwalk, sisum);
167 /* write hh.walkers0 */
168
169 for (k = 0; k < my; k++) {
170 for (l = 0; l < mx; l++) { /* run thru the whole area */
171 if (zz[k][l] != UNDEF) {
172
173 x = xp0 + stepx * (double)(l);
174 y = yp0 + stepy * (double)(k);
175
176 gen = rwalk * si[k][l] / sisum;
177 mgen = (int)gen;
178 wei = gen / (double)(mgen + 1);
179
180 /*if (si[k][l] != 0.) { */
181 /* this stuff later for multiscale */
182
183 gen2 =
184 (double)maxwab *si[k][l] / (si0 *
185 (double)(mx2o * my2o));
186 gen2 = gen2 * (barea / sarea);
187 mgen2 = (int)gen2;
188 wei2 = gen2 / (double)(mgen2 + 1);
189 mgen3 =
190 (int)((double)mgen2 * wei2 / ((double)mgen * wei));
191 nmult = mgen3 + 1;
192 wei3 = gen2 / (double)((mgen + 1) * (mgen2 + 1));
193 weifac = wei3 / wei;
194 /* } else {
195 nmult = 1;
196 weifac = 1.;
197 fprintf(stderr, "\n zero rainfall excess in cell");
198 } */
199
200 /*G_debug(2, " gen,gen2,wei,wei2,mgen3,nmult: %f %f %f %f %d %d",gen,gen2,wei,wei2,mgen3,nmult);
201 */
202 for (iw = 1; iw <= mgen + 1; iw++) { /* assign walkers */
203 w[lw].x = x + stepx * (simwe_rand() - 0.5);
204 w[lw].y = y + stepy * (simwe_rand() - 0.5);
205 w[lw].m = wei;
206
207 walkwe += w[lw].m;
208 vavg[lw].x = v1[k][l];
209 vavg[lw].y = v2[k][l];
210 /* deactivated as unused, what was iflag for?
211 if (w[lw].x >= xmin && w[lw].y >= ymin &&
212 w[lw].x <= xmax && w[lw].y <= ymax) {
213 iflag[lw] = 0;
214 }
215 else {
216 iflag[lw] = 1;
217 }
218 */
219 lw++;
220 }
221 } /*DEFined area */
222 }
223 }
224 nwalk = lw;
225 G_debug(2, " nwalk, maxw %d %d", nwalk, MAXW);
226 G_debug(2, " walkwe (walk weight),frac %f %f", walkwe, frac);
227
228 stxm = stepx * (double)(mx + 1) - xmin;
229 stym = stepy * (double)(my + 1) - ymin;
230 nwalka = 0;
231 deldif = sqrt(deltap) * frac; /* diffuse factor */
232
233
234 factor = deltap * sisum / (rwalk * (double)nblock);
235
236 G_debug(2, " deldif,factor %f %e", deldif, factor);
237
238 /* ********************************************************** */
239 /* main loop over the projection time */
240 /* *********************************************************** */
241
242
243 G_debug(2, "main loop over the projection time... ");
244
245 for (i = 1; i <= miter; i++) { /* iteration loop depending on simulation time and deltap */
246 G_percent(i, miter, 1);
247 iter1 = i / iterout;
248 iter1 *= iterout;
249 if (iter1 == i) {
250 nfiterw = i / iterout + 10;
251 nfiterh = i / iterout + 40;
252 G_debug(2, "iblock=%d i=%d miter=%d nwalk=%d nwalka=%d",
253 iblock, i, miter, nwalk, nwalka);
254 }
255
256 if (nwalka == 0 && i > 1)
257 goto L_800;
258
259 /* ************************************************************ */
260 /* .... propagate one step */
261 /* ************************************************************ */
262
263 addac = factor;
264 conn = (double)nblock / (double)iblock;
265 if (i == 1) {
266 addac = factor * .5;
267 }
268 nwalka = 0;
269 nstack = 0;
270
271 #pragma omp parallel firstprivate(l,lw,k,decr,d1,hhc,velx,vely,eff,gaux,gauy)//nwalka
272 {
273 #if defined(_OPENMP)
274 int steps = (int)((((double)nwalk) / ((double) omp_get_num_threads())) + 0.5);
275 int tid = omp_get_thread_num();
276 int min_loop = tid * steps;
277 int max_loop = ((tid + 1) * steps) > nwalk ? nwalk : (tid + 1) * steps;
278
279 for (lw = min_loop; lw < max_loop; lw++) {
280 #else
281 for (lw = 0; lw < nwalk; lw++) {
282 #endif
283 if (w[lw].m > EPS) { /* check the walker weight */
284 ++nwalka;
285 l = (int)((w[lw].x + stxm) / stepx) - mx - 1;
286 k = (int)((w[lw].y + stym) / stepy) - my - 1;
287
288 if (l > mx - 1 || k > my - 1 || k < 0 || l < 0) {
289
290 G_debug(2, " k,l=%d,%d", k, l);
291 printf(" lw,w=%d %f %f", lw, w[lw].y, w[lw].m);
292 G_debug(2, " stxym=%f %f", stxm, stym);
293 printf(" step=%f %f", stepx, stepy);
294 G_debug(2, " m=%d %d", my, mx);
295 printf(" nwalka,nwalk=%d %d", nwalka, nwalk);
296 G_debug(2, " ");
297 }
298
299 if (zz[k][l] != UNDEF) {
300 if (infil != NULL) { /* infiltration part */
301 if (inf[k][l] - si[k][l] > 0.) {
302
303 decr = pow(addac * w[lw].m, 3. / 5.); /* decreasing factor in m */
304 if (inf[k][l] > decr) {
305 inf[k][l] -= decr; /* decrease infilt. in cell and eliminate the walker */
306 w[lw].m = 0.;
307 }
308 else {
309 w[lw].m -= pow(inf[k][l], 5. / 3.) / addac; /* use just proportional part of the walker weight */
310 inf[k][l] = 0.;
311
312 }
313 }
314 }
315
316 gama[k][l] += (addac * w[lw].m); /* add walker weigh to water depth or conc. */
317
318 d1 = gama[k][l] * conn;
319 #if defined(_OPENMP)
320 gasdev_for_paralel(&gaux, &gauy);
321 #else
322 gaux = gasdev();
323 gauy = gasdev();
324 #endif
325 hhc = pow(d1, 3. / 5.);
326
327 if (hhc > hhmax && wdepth == NULL) { /* increased diffusion if w.depth > hhmax */
328 dif[k][l] = (halpha + 1) * deldif;
329 velx = vavg[lw].x;
330 vely = vavg[lw].y;
331 }
332 else {
333 dif[k][l] = deldif;
334 velx = v1[k][l];
335 vely = v2[k][l];
336 }
337
338
339 if (traps != NULL && trap[k][l] != 0.) { /* traps */
340
341 eff = simwe_rand(); /* random generator */
342
343 if (eff <= trap[k][l]) {
344 velx = -0.1 * v1[k][l]; /* move it slightly back */
345 vely = -0.1 * v2[k][l];
346 }
347 }
348
349 w[lw].x += (velx + dif[k][l] * gaux); /* move the walker */
350 w[lw].y += (vely + dif[k][l] * gauy);
351
352 if (hhc > hhmax && wdepth == NULL) {
353 vavg[lw].x = hbeta * (vavg[lw].x + v1[k][l]);
354 vavg[lw].y = hbeta * (vavg[lw].y + v2[k][l]);
355 }
356
357 if (w[lw].x <= xmin || w[lw].y <= ymin || w[lw].x
358 >= xmax || w[lw].y >= ymax) {
359 w[lw].m = 1e-10; /* eliminate walker if it is out of area */
360 }
361 else {
362 if (wdepth != NULL) {
363 l = (int)((w[lw].x + stxm) / stepx) - mx - 1;
364 k = (int)((w[lw].y + stym) / stepy) - my - 1;
365 w[lw].m *= sigma[k][l];
366 }
367
368 } /* else */
369 } /*DEFined area */
370 else {
371 w[lw].m = 1e-10; /* eliminate walker if it is out of area */
372 }
373 }
374 } /* lw loop */
375 }
376 /* Changes made by Soeren 8. Mar 2011 to replace the site walker output implementation */
377 /* Save all walkers located within the computational region and with valid
378 z coordinates */
379 if (outwalk != NULL && (i == miter || i == iter1)) {
380 nstack = 0;
381
382 for (lw = 0; lw < nwalk; lw++) {
383 /* Compute the elevation raster map index */
384 l = (int)((w[lw].x + stxm) / stepx) - mx - 1;
385 k = (int)((w[lw].y + stym) / stepy) - my - 1;
386
387 /* Check for correct elevation raster map index */
388 if(l < 0 || l >= mx || k < 0 || k >= my)
389 continue;
390
391 if (w[lw].m > EPS && zz[k][l] != UNDEF) {
392
393 /* Save the 3d position of the walker */
394 stack[nstack].x = mixx / conv + w[lw].x / conv;
395 stack[nstack].y = miyy / conv + w[lw].y / conv;
396 stack[nstack].m = zz[k][l];
397
398 nstack++;
399 }
400 } /* lw loop */
401 }
402
403 if (i == iter1 && ts == 1) {
404 /* call output for iteration output */
405 if (erdep != NULL)
406 erod(gama); /* divergence of gama field */
407
408 conn = (double)nblock / (double)iblock;
409 itime = (int)(i * deltap * timec);
410 ii = output_data(itime, conn);
411 if (ii != 1)
412 G_fatal_error(_("Unable to write raster maps"));
413 }
414
415 /* Write the water depth each time step at an observation point */
416 if(points.is_open)
417 {
418 double value = 0.0;
419 int p;
420 fprintf(points.output, "%.6d ", i);
421 /* Write for each point */
422 for(p = 0; p < points.npoints; p++)
423 {
424 l = (int)((points.x[p] - mixx + stxm) / stepx) - mx - 1;
425 k = (int)((points.y[p] - miyy + stym) / stepy) - my - 1;
426
427 if (zz[k][l] != UNDEF) {
428
429 if (wdepth == NULL)
430 value = step * gama[k][l] * cchez[k][l];
431 else
432 value = gama[k][l] * slope[k][l];
433
434 fprintf(points.output, "%2.4f ", value);
435 } else {
436 /* Point is invalid, so a negative value is written */
437 fprintf(points.output, "%2.4f ", -1.0);
438 }
439 }
440 fprintf(points.output, "\n");
441 }
442 } /* miter */
443
444 L_800:
445 /* Soeren 8. Mar 2011: Why is this commented out?*/
446 /* if (iwrib != nblock) {
447 icount = icoub / iwrib;
448
449 if (icoub == (icount * iwrib)) {
450 ++icfl;
451 nflw = icfl + 50;
452 conn = (double) nblock / (double) iblock;
453
454 }
455 } */
456
457
458 if (err != NULL) {
459 for (k = 0; k < my; k++) {
460 for (l = 0; l < mx; l++) {
461 if (zz[k][l] != UNDEF) {
462 d1 = gama[k][l] * (double)conn;
463 gammas[k][l] += pow(d1, 3. / 5.);
464 } /* DEFined area */
465 }
466 }
467 }
468 if (erdep != NULL)
469 erod(gama);
470 }
471 /* ........ end of iblock loop */
472
473 /* Write final maps here because we know the last time stamp here */
474 if (ts == 0) {
475 conn = (double)nblock / (double)iblock;
476 itime = (int)(i * deltap * timec);
477 ii = output_data(itime, conn);
478 if (ii != 1)
479 G_fatal_error(_("Cannot write raster maps"));
480 }
481 /* Close the observation logfile */
482 if(points.is_open)
483 fclose(points.output);
484
485 points.is_open = 0;
486
487 }
488