GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
n_gwflow.c
Go to the documentation of this file.
1
2/*****************************************************************************
3*
4* MODULE: Grass PDE Numerical Library
5* AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6* soerengebbert <at> gmx <dot> de
7*
8* PURPOSE: groundwater flow in porous media
9* part of the gpde library
10*
11* COPYRIGHT: (C) 2000 by the GRASS Development Team
12*
13* This program is free software under the GNU General Public
14* License (>=v2). Read the file COPYING that comes with GRASS
15* for details.
16*
17*****************************************************************************/
18
19#include <grass/N_gwflow.h>
20
21/* *************************************************************** */
22/* ***************** N_gwflow_data3d ***************************** */
23/* *************************************************************** */
24/*!
25 * \brief Alllocate memory for the groundwater calculation data structure in 3 dimensions
26 *
27 * The groundwater calculation data structure will be allocated including
28 * all appendant 3d and 2d arrays. The offset for the 3d arrays is one
29 * to establish homogeneous Neumann boundary conditions at the calculation area border.
30 * This data structure is used to create a linear equation system based on the computation of
31 * groundwater flow in porous media with the finite volume method.
32 *
33 * \param cols int
34 * \param rows int
35 * \param depths int
36 * \return N_gwflow_data3d *
37 * */
38N_gwflow_data3d *N_alloc_gwflow_data3d(int cols, int rows, int depths,
39 int river, int drain)
40{
41 N_gwflow_data3d *data;
42
43 data = (N_gwflow_data3d *) G_calloc(1, sizeof(N_gwflow_data3d));
44
45 data->phead = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
46 data->phead_start = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
47 data->status = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
48 data->hc_x = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
49 data->hc_y = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
50 data->hc_z = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
51 data->q = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
52 data->s = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
53 data->nf = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
54 data->r = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
55
56 if (river) {
57 data->river_head =
58 N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
59 data->river_leak =
60 N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
61 data->river_bed = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
62 }
63 else {
64 data->river_head = NULL;
65 data->river_leak = NULL;
66 data->river_bed = NULL;
67 }
68
69 if (drain) {
70 data->drain_leak =
71 N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
72 data->drain_bed = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
73 }
74 else {
75 data->drain_leak = NULL;
76 data->drain_bed = NULL;
77 }
78
79 return data;
80}
81
82/* *************************************************************** */
83/* ********************* N_free_gwflow_data3d ******************** */
84/* *************************************************************** */
85/*!
86 * \brief Release the memory of the groundwater flow data structure in three dimensions
87 *
88 * \param data N_gwflow_data3d *
89 * \return void *
90 * */
91
93{
94 if (data->phead)
96 if (data->phead_start)
98 if (data->status)
100 if (data->hc_x)
101 N_free_array_3d(data->hc_x);
102 if (data->hc_y)
103 N_free_array_3d(data->hc_y);
104 if (data->hc_z)
105 N_free_array_3d(data->hc_z);
106 if (data->q)
107 N_free_array_3d(data->q);
108 if (data->s)
109 N_free_array_3d(data->s);
110 if (data->nf)
111 N_free_array_3d(data->nf);
112 if (data->r)
113 N_free_array_2d(data->r);
114 if (data->river_head)
116 if (data->river_leak)
118 if (data->river_bed)
120 if (data->drain_leak)
122 if (data->drain_bed)
124
125 G_free(data);
126
127 data = NULL;
128
129 return;
130}
131
132/* *************************************************************** */
133/* ******************** N_alloc_gwflow_data2d ******************** */
134/* *************************************************************** */
135/*!
136 * \brief Alllocate memory for the groundwater calculation data structure in 2 dimensions
137 *
138 * The groundwater calculation data structure will be allocated including
139 * all appendant 2d arrays. The offset for the 3d arrays is one
140 * to establish homogeneous Neumann boundary conditions at the calculation area border.
141 * This data structure is used to create a linear equation system based on the computation of
142 * groundwater flow in porous media with the finite volume method.
143 *
144 * \param cols int
145 * \param rows int
146 * \return N_gwflow_data2d *
147 * */
148N_gwflow_data2d *N_alloc_gwflow_data2d(int cols, int rows, int river,
149 int drain)
150{
151 N_gwflow_data2d *data;
152
153 data = (N_gwflow_data2d *) G_calloc(1, sizeof(N_gwflow_data2d));
154
155 data->phead = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
156 data->phead_start = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
157 data->status = N_alloc_array_2d(cols, rows, 1, CELL_TYPE);
158 data->hc_x = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
159 data->hc_y = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
160 data->q = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
161 data->s = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
162 data->nf = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
163 data->r = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
164 data->top = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
165 data->bottom = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
166
167 if (river) {
168 data->river_head = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
169 data->river_leak = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
170 data->river_bed = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
171 }
172 else {
173 data->river_head = NULL;
174 data->river_leak = NULL;
175 data->river_bed = NULL;
176 }
177
178 if (drain) {
179 data->drain_leak = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
180 data->drain_bed = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
181 }
182 else {
183 data->drain_leak = NULL;
184 data->drain_bed = NULL;
185 }
186
187
188 return data;
189}
190
191/* *************************************************************** */
192/* ****************** N_free_gwflow_data2d *********************** */
193/* *************************************************************** */
194/*!
195 * \brief Release the memory of the groundwater flow data structure in two dimensions
196 *
197 * \param data N_gwflow_data2d *
198 * \return void
199 * */
201{
202 if (data->phead)
203 N_free_array_2d(data->phead);
204 if (data->phead_start)
206 if (data->status)
207 N_free_array_2d(data->status);
208 if (data->hc_x)
209 N_free_array_2d(data->hc_x);
210 if (data->hc_y)
211 N_free_array_2d(data->hc_y);
212 if (data->q)
213 N_free_array_2d(data->q);
214 if (data->s)
215 N_free_array_2d(data->s);
216 if (data->nf)
217 N_free_array_2d(data->nf);
218 if (data->r)
219 N_free_array_2d(data->r);
220 if (data->top)
221 N_free_array_2d(data->top);
222 if (data->bottom)
223 N_free_array_2d(data->bottom);
224 if (data->river_head)
226 if (data->river_leak)
228 if (data->river_bed)
230 if (data->drain_leak)
232 if (data->drain_bed)
234
235 G_free(data);
236
237 data = NULL;;
238
239 return;
240}
241
242/* *************************************************************** */
243/* ***************** N_callback_gwflow_3d ************************ */
244/* *************************************************************** */
245/*!
246 * \brief This callback function creates the mass balance of a 7 point star
247 *
248 * The mass balance is based on the common groundwater flow equation:
249 *
250 * \f[Ss \frac{\partial h}{\partial t} = \nabla {\bf K} \nabla h + q \f]
251 *
252 * This equation is discretizised with the finite volume method in three dimensions.
253 *
254 *
255 * \param gwdata N_gwflow_data3d *
256 * \param geom N_geom_data *
257 * \param col int
258 * \param row int
259 * \param depth int
260 * \return N_data_star *
261 *
262 * */
263N_data_star *N_callback_gwflow_3d(void *gwdata, N_geom_data * geom, int col,
264 int row, int depth)
265{
266 double hc_e = 0, hc_w = 0, hc_n = 0, hc_s = 0, hc_t = 0, hc_b = 0;
267 double dx, dy, dz, Ax, Ay, Az;
268 double hc_x, hc_y, hc_z;
269 double hc_xw, hc_yn, hc_zt;
270 double hc_xe, hc_ys, hc_zb;
271 double hc_start;
272 double Ss, r, nf, q;
273 double C, W, E, N, S, T, B, V;
274 N_data_star *mat_pos;
275 N_gwflow_data3d *data;
276
277 /*cast the void pointer to the right data structure */
278 data = (N_gwflow_data3d *) gwdata;
279
280 dx = geom->dx;
281 dy = geom->dy;
282 dz = geom->dz;
283 Az = N_get_geom_data_area_of_cell(geom, row);
284 Ay = geom->dx * geom->dz;
285 Ax = geom->dz * geom->dy;
286
287 /*read the data from the arrays */
288 hc_start = N_get_array_3d_d_value(data->phead_start, col, row, depth);
289
290 hc_x = N_get_array_3d_d_value(data->hc_x, col, row, depth);
291 hc_y = N_get_array_3d_d_value(data->hc_y, col, row, depth);
292 hc_z = N_get_array_3d_d_value(data->hc_z, col, row, depth);
293
294 hc_xw = N_get_array_3d_d_value(data->hc_x, col - 1, row, depth);
295 hc_xe = N_get_array_3d_d_value(data->hc_x, col + 1, row, depth);
296 hc_yn = N_get_array_3d_d_value(data->hc_y, col, row - 1, depth);
297 hc_ys = N_get_array_3d_d_value(data->hc_y, col, row + 1, depth);
298 hc_zt = N_get_array_3d_d_value(data->hc_z, col, row, depth + 1);
299 hc_zb = N_get_array_3d_d_value(data->hc_z, col, row, depth - 1);
300
301 hc_w = N_calc_harmonic_mean(hc_xw, hc_x);
302 hc_e = N_calc_harmonic_mean(hc_xe, hc_x);
303 hc_n = N_calc_harmonic_mean(hc_yn, hc_y);
304 hc_s = N_calc_harmonic_mean(hc_ys, hc_y);
305 hc_t = N_calc_harmonic_mean(hc_zt, hc_z);
306 hc_b = N_calc_harmonic_mean(hc_zb, hc_z);
307
308 /*inner sources */
309 q = N_get_array_3d_d_value(data->q, col, row, depth);
310 /*storativity */
311 Ss = N_get_array_3d_d_value(data->s, col, row, depth);
312 /*porosity */
313 nf = N_get_array_3d_d_value(data->nf, col, row, depth);
314
315 /*mass balance center cell to western cell */
316 W = -1 * Ax * hc_w / dx;
317 /*mass balance center cell to eastern cell */
318 E = -1 * Ax * hc_e / dx;
319 /*mass balance center cell to northern cell */
320 N = -1 * Ay * hc_n / dy;
321 /*mass balance center cell to southern cell */
322 S = -1 * Ay * hc_s / dy;
323 /*mass balance center cell to top cell */
324 T = -1 * Az * hc_t / dz;
325 /*mass balance center cell to bottom cell */
326 B = -1 * Az * hc_b / dz;
327
328 /*storativity */
329 Ss = Az * dz * Ss;
330
331 /*the diagonal entry of the matrix */
332 C = -1 * (W + E + N + S + T + B - Ss / data->dt * Az);
333
334 /*the entry in the right side b of Ax = b */
335 V = (q + hc_start * Ss / data->dt * Az);
336
337 /*only the top cells will have recharge */
338 if (depth == geom->depths - 2) {
339 r = N_get_array_2d_d_value(data->r, col, row);
340 V += r * Az;
341 }
342
343 G_debug(5, "N_callback_gwflow_3d: called [%i][%i][%i]", depth, col, row);
344
345 /*create the 7 point star entries */
346 mat_pos = N_create_7star(C, W, E, N, S, T, B, V);
347
348 return mat_pos;
349}
350
351
352/* *************************************************************** */
353/* ****************** N_gwflow_3d_calc_water_budget ************** */
354/* *************************************************************** */
355/*!
356 * \brief This function computes the water budget of the entire groundwater
357 *
358 * The water budget is calculated for each active and dirichlet cell from
359 * its surrounding neighbours. This is based on the 7 star mass balance computation
360 * of N_callback_gwflow_3d and the gradient of the water heights in the cells.
361 * The sum of the water budget of each active/dirichlet cell must be near zero
362 * due the effect of numerical inaccuracy of cpu's.
363 *
364 * \param gwdata N_gwflow_data3d *
365 * \param geom N_geom_data *
366 * \param budget N_array_3d
367 * \return void
368 *
369 * */
370void
372{
373 int z, y, x, stat;
374 double h, hc;
375 double val;
376 double sum;
377 N_data_star *dstar;
378
379 int rows = data->status->rows;
380 int cols = data->status->cols;
381 int depths = data->status->depths;
382 sum = 0;
383
384 for (z = 0; z < depths; z++) {
385 for (y = 0; y < rows; y++) {
386 G_percent(y, rows - 1, 10);
387 for (x = 0; x < cols; x++) {
388 stat = (int)N_get_array_3d_d_value(data->status, x, y, z);
389
390 val = 0.0;
391
392 if (stat != N_CELL_INACTIVE ) { /*all active/dirichlet cells */
393
394 /* Compute the flow parameter */
395 dstar = N_callback_gwflow_3d(data, geom, x, y, z);
396 /* Compute the gradient in each direction pointing from the center */
397 hc = N_get_array_3d_d_value(data->phead, x, y, z);
398
399 if((int)N_get_array_3d_d_value(data->status, x + 1, y , z) != N_CELL_INACTIVE) {
400 h = N_get_array_3d_d_value(data->phead, x + 1, y , z);
401 val += dstar->E * (hc - h);
402 }
403 if((int)N_get_array_3d_d_value(data->status, x - 1, y , z) != N_CELL_INACTIVE) {
404 h = N_get_array_3d_d_value(data->phead, x - 1, y , z);
405 val += dstar->W * (hc - h);
406 }
407 if((int)N_get_array_3d_d_value(data->status, x , y + 1, z) != N_CELL_INACTIVE) {
408 h = N_get_array_3d_d_value(data->phead, x , y + 1, z);
409 val += dstar->S * (hc - h);
410 }
411 if((int)N_get_array_3d_d_value(data->status, x , y - 1, z) != N_CELL_INACTIVE) {
412 h = N_get_array_3d_d_value(data->phead, x , y - 1, z);
413 val += dstar->N * (hc - h);
414 }
415 if((int)N_get_array_3d_d_value(data->status, x , y , z + 1) != N_CELL_INACTIVE) {
416 h = N_get_array_3d_d_value(data->phead, x , y , z + 1);
417 val += dstar->T * (hc - h);
418 }
419 if((int)N_get_array_3d_d_value(data->status, x , y , z - 1) != N_CELL_INACTIVE) {
420 h = N_get_array_3d_d_value(data->phead, x , y , z - 1);
421 val += dstar->B * (hc - h);
422 }
423 sum += val;
424
425 G_free(dstar);
426 }
427 else {
428 Rast_set_null_value(&val, 1, DCELL_TYPE);
429 }
430 N_put_array_3d_d_value(budget, x, y, z, val);
431 }
432 }
433 }
434
435 if(fabs(sum) < 0.0000000001)
436 G_message(_("The total sum of the water budget: %g\n"), sum);
437 else
438 G_warning(_("The total sum of the water budget is significantly larger then 0: %g\n"), sum);
439
440 return;
441}
442
443
444
445/* *************************************************************** */
446/* ****************** N_callback_gwflow_2d *********************** */
447/* *************************************************************** */
448/*!
449 * \brief This callback function creates the mass balance of a 5 point star
450 *
451 * The mass balance is based on the common groundwater flow equation:
452 *
453 * \f[Ss \frac{\partial h}{\partial t} = \nabla {\bf K} \nabla h + q \f]
454 *
455 * This equation is discretizised with the finite volume method in two dimensions.
456 *
457 * \param gwdata N_gwflow_data2d *
458 * \param geom N_geom_data *
459 * \param col int
460 * \param row int
461 * \return N_data_star *
462 *
463 * */
464N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
465 int row)
466{
467 double T_e = 0, T_w = 0, T_n = 0, T_s = 0;
468 double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
469 double dx, dy, Az;
470 double hc_x, hc_y;
471 double z, top;
472 double hc_xw, hc_yn;
473 double z_xw, z_yn;
474 double hc_xe, hc_ys;
475 double z_xe, z_ys;
476 double hc, hc_start;
477 double Ss, r, q;
478 double C, W, E, N, S, V;
479 N_gwflow_data2d *data;
480 N_data_star *mat_pos;
481 double river_vect = 0; /*entry in vector */
482 double river_mat = 0; /*entry in matrix */
483 double drain_vect = 0; /*entry in vector */
484 double drain_mat = 0; /*entry in matrix */
485
486 /*cast the void pointer to the right data structure */
487 data = (N_gwflow_data2d *) gwdata;
488
489 dx = geom->dx;
490 dy = geom->dy;
491 Az = N_get_geom_data_area_of_cell(geom, row);
492
493 /*read the data from the arrays */
494 hc_start = N_get_array_2d_d_value(data->phead_start, col, row);
495 hc = N_get_array_2d_d_value(data->phead, col, row);
496 top = N_get_array_2d_d_value(data->top, col, row);
497
498 /* Inner sources */
499 q = N_get_array_2d_d_value(data->q, col, row);
500
501 /* storativity or porosity of current cell face [-]*/
502 Ss = N_get_array_2d_d_value(data->s, col, row);
503 /* recharge */
504 r = N_get_array_2d_d_value(data->r, col, row) * Az;
505
506
507 if (hc > top) { /*If the aquifer is confined */
508 z = N_get_array_2d_d_value(data->top, col,
509 row) -
510 N_get_array_2d_d_value(data->bottom, col, row);
511 z_xw =
512 N_get_array_2d_d_value(data->top, col - 1,
513 row) -
514 N_get_array_2d_d_value(data->bottom, col - 1, row);
515 z_xe =
516 N_get_array_2d_d_value(data->top, col + 1,
517 row) -
518 N_get_array_2d_d_value(data->bottom, col + 1, row);
519 z_yn =
520 N_get_array_2d_d_value(data->top, col,
521 row - 1) -
522 N_get_array_2d_d_value(data->bottom, col, row - 1);
523 z_ys =
524 N_get_array_2d_d_value(data->top, col,
525 row + 1) -
526 N_get_array_2d_d_value(data->bottom, col, row + 1);
527 }
528 else { /* the aquifer is unconfined */
529
530 /* If the aquifer is unconfied use an explicite scheme to solve
531 * the nonlinear equation. We use the phead from the first iteration */
532 z = N_get_array_2d_d_value(data->phead, col, row) -
533 N_get_array_2d_d_value(data->bottom, col, row);
534 z_xw = N_get_array_2d_d_value(data->phead, col - 1, row) -
535 N_get_array_2d_d_value(data->bottom, col - 1, row);
536 z_xe = N_get_array_2d_d_value(data->phead, col + 1, row) -
537 N_get_array_2d_d_value(data->bottom, col + 1, row);
538 z_yn = N_get_array_2d_d_value(data->phead, col, row - 1) -
539 N_get_array_2d_d_value(data->bottom, col, row - 1);
540 z_ys = N_get_array_2d_d_value(data->phead, col, row + 1) -
541 N_get_array_2d_d_value(data->bottom, col, row + 1);
542 }
543
544 /*geometrical mean of cell height */
545 if (z_w > 0 || z_w < 0 || z_w == 0)
546 z_w = N_calc_arith_mean(z_xw, z);
547 else
548 z_w = z;
549 if (z_e > 0 || z_e < 0 || z_e == 0)
550 z_e = N_calc_arith_mean(z_xe, z);
551 else
552 z_e = z;
553 if (z_n > 0 || z_n < 0 || z_n == 0)
554 z_n = N_calc_arith_mean(z_yn, z);
555 else
556 z_n = z;
557 if (z_s > 0 || z_s < 0 || z_s == 0)
558 z_s = N_calc_arith_mean(z_ys, z);
559 else
560 z_s = z;
561
562 /*get the surrounding permeabilities */
563 hc_x = N_get_array_2d_d_value(data->hc_x, col, row);
564 hc_y = N_get_array_2d_d_value(data->hc_y, col, row);
565 hc_xw = N_get_array_2d_d_value(data->hc_x, col - 1, row);
566 hc_xe = N_get_array_2d_d_value(data->hc_x, col + 1, row);
567 hc_yn = N_get_array_2d_d_value(data->hc_y, col, row - 1);
568 hc_ys = N_get_array_2d_d_value(data->hc_y, col, row + 1);
569
570 /* calculate the transmissivities */
571 T_w = N_calc_harmonic_mean(hc_xw, hc_x) * z_w;
572 T_e = N_calc_harmonic_mean(hc_xe, hc_x) * z_e;
573 T_n = N_calc_harmonic_mean(hc_yn, hc_y) * z_n;
574 T_s = N_calc_harmonic_mean(hc_ys, hc_y) * z_s;
575
576 /* Compute the river leakage, this is an explicit method
577 * Influent and effluent flow is computed.
578 */
579 if (data->river_leak &&
580 (N_get_array_2d_d_value(data->river_leak, col, row) != 0) &&
581 N_get_array_2d_d_value(data->river_bed, col, row) <= top) {
582 /* Groundwater surface is above the river bed*/
583 if (hc > N_get_array_2d_d_value(data->river_bed, col, row)) {
584 river_vect = N_get_array_2d_d_value(data->river_head, col, row) *
585 N_get_array_2d_d_value(data->river_leak, col, row);
586 river_mat = N_get_array_2d_d_value(data->river_leak, col, row);
587 } /* Groundwater surface is below the river bed */
588 else if (hc < N_get_array_2d_d_value(data->river_bed, col, row)) {
589 river_vect = (N_get_array_2d_d_value(data->river_head, col, row) -
590 N_get_array_2d_d_value(data->river_bed, col, row))
591 * N_get_array_2d_d_value(data->river_leak, col, row);
592 river_mat = 0;
593 }
594 }
595
596 /* compute the drainage, this is an explicit method
597 * Drainage is only enabled, if the drain bed is lower the groundwater surface
598 */
599 if (data->drain_leak &&
600 (N_get_array_2d_d_value(data->drain_leak, col, row) != 0) &&
601 N_get_array_2d_d_value(data->drain_bed, col, row) <= top) {
602 if (hc > N_get_array_2d_d_value(data->drain_bed, col, row)) {
603 drain_vect = N_get_array_2d_d_value(data->drain_bed, col, row) *
604 N_get_array_2d_d_value(data->drain_leak, col, row);
605 drain_mat = N_get_array_2d_d_value(data->drain_leak, col, row);
606 }
607 else if (hc <= N_get_array_2d_d_value(data->drain_bed, col, row)) {
608 drain_vect = 0;
609 drain_mat = 0;
610 }
611 }
612
613 /*mass balance center cell to western cell */
614 W = -1 * T_w * dy / dx;
615 /*mass balance center cell to eastern cell */
616 E = -1 * T_e * dy / dx;
617 /*mass balance center cell to northern cell */
618 N = -1 * T_n * dx / dy;
619 /*mass balance center cell to southern cell */
620 S = -1 * T_s * dx / dy;
621
622 /*the diagonal entry of the matrix */
623 C = -1 * (W + E + N + S - Az *Ss / data->dt - river_mat * Az -
624 drain_mat * Az);
625
626 /*the entry in the right side b of Ax = b */
627 V = (q + hc_start * Az * Ss / data->dt) + r + river_vect * Az +
628 drain_vect * Az;
629
630 G_debug(5, "N_callback_gwflow_2d: called [%i][%i]", row, col);
631
632 /*create the 5 point star entries */
633 mat_pos = N_create_5star(C, W, E, N, S, V);
634
635 return mat_pos;
636}
637
638
639
640/* *************************************************************** */
641/* ****************** N_gwflow_2d_calc_water_budget ************** */
642/* *************************************************************** */
643/*!
644 * \brief This function computes the water budget of the entire groundwater
645 *
646 * The water budget is calculated for each active and dirichlet cell from
647 * its surrounding neighbours. This is based on the 5 star mass balance computation
648 * of N_callback_gwflow_2d and the gradient of the water heights in the cells.
649 * The sum of the water budget of each active/dirichlet cell must be near zero
650 * due the effect of numerical inaccuracy of cpu's.
651 *
652 * \param gwdata N_gwflow_data2d *
653 * \param geom N_geom_data *
654 * \param budget N_array_2d
655 * \return void
656 *
657 * */
658void
660{
661 int y, x, stat;
662 double h, hc;
663 double val;
664 double sum;
665 N_data_star *dstar;
666
667 int rows = data->status->rows;
668 int cols = data->status->cols;
669
670 sum = 0;
671
672 for (y = 0; y < rows; y++) {
673 G_percent(y, rows - 1, 10);
674 for (x = 0; x < cols; x++) {
675 stat = N_get_array_2d_c_value(data->status, x, y);
676
677 val = 0.0;
678
679 if (stat != N_CELL_INACTIVE ) { /*all active/dirichlet cells */
680
681 /* Compute the flow parameter */
682 dstar = N_callback_gwflow_2d(data, geom, x, y);
683 /* Compute the gradient in each direction pointing from the center */
684 hc = N_get_array_2d_d_value(data->phead, x, y);
685
686 if((int)N_get_array_2d_d_value(data->status, x + 1, y ) != N_CELL_INACTIVE) {
687 h = N_get_array_2d_d_value(data->phead, x + 1, y);
688 val += dstar->E * (hc - h);
689 }
690 if((int)N_get_array_2d_d_value(data->status, x - 1, y ) != N_CELL_INACTIVE) {
691 h = N_get_array_2d_d_value(data->phead, x - 1, y);
692 val += dstar->W * (hc - h);
693 }
694 if((int)N_get_array_2d_d_value(data->status, x , y + 1) != N_CELL_INACTIVE) {
695 h = N_get_array_2d_d_value(data->phead, x , y + 1);
696 val += dstar->S * (hc - h);
697 }
698 if((int)N_get_array_2d_d_value(data->status, x , y - 1) != N_CELL_INACTIVE) {
699 h = N_get_array_2d_d_value(data->phead, x , y - 1);
700 val += dstar->N * (hc - h);
701 }
702
703 sum += val;
704
705 G_free(dstar);
706 }
707 else {
708 Rast_set_null_value(&val, 1, DCELL_TYPE);
709 }
710 N_put_array_2d_d_value(budget, x, y, val);
711 }
712 }
713
714 if(fabs(sum) < 0.0000000001)
715 G_message(_("The total sum of the water budget: %g\n"), sum);
716 else
717 G_warning(_("The total sum of the water budget is significantly larger then 0: %g\n"), sum);
718
719 return;
720}
#define N_CELL_INACTIVE
Definition: N_pde.h:31
double N_calc_arith_mean(double a, double b)
Calculate the arithmetic mean of values a and b.
Definition: n_tools.c:33
double N_calc_harmonic_mean(double a, double b)
Calculate the harmonical mean of values a and b.
Definition: n_tools.c:119
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
#define NULL
Definition: ccmath.h:32
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
double r
void G_message(const char *msg,...)
Print a message to stderr.
Definition: gis/error.c:90
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
N_array_3d * N_alloc_array_3d(int cols, int rows, int depths, int offset, int type)
Allocate memory for a N_array_3d data structure.
Definition: n_arrays.c:726
CELL N_get_array_2d_c_value(N_array_2d *data, int col, int row)
Returns the value of type CELL at position col, row.
Definition: n_arrays.c:314
void N_free_array_3d(N_array_3d *data)
Release the memory of a N_array_3d.
Definition: n_arrays.c:780
DCELL N_get_array_2d_d_value(N_array_2d *data, int col, int row)
Returns the value of type DCELL at position col, row.
Definition: n_arrays.c:378
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d structure.
Definition: n_arrays.c:130
void N_put_array_3d_d_value(N_array_3d *data, int col, int row, int depth, double value)
Writes a double value to the N_array_3d struct at position col, row, depth.
Definition: n_arrays.c:1175
N_array_2d * N_alloc_array_2d(int cols, int rows, int offset, int type)
Allocate memory for a N_array_2d data structure.
Definition: n_arrays.c:72
double N_get_array_3d_d_value(N_array_3d *data, int col, int row, int depth)
This function returns the value of type float at position col, row, depth.
Definition: n_arrays.c:990
void N_put_array_2d_d_value(N_array_2d *data, int col, int row, DCELL value)
Writes a DCELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:584
double N_get_geom_data_area_of_cell(N_geom_data *geom, int row)
Get the areay size in square meter of one cell (x*y) at row.
Definition: n_geom.c:192
N_data_star * N_callback_gwflow_2d(void *gwdata, N_geom_data *geom, int col, int row)
This callback function creates the mass balance of a 5 point star.
Definition: n_gwflow.c:464
void N_gwflow_3d_calc_water_budget(N_gwflow_data3d *data, N_geom_data *geom, N_array_3d *budget)
This function computes the water budget of the entire groundwater.
Definition: n_gwflow.c:371
N_gwflow_data2d * N_alloc_gwflow_data2d(int cols, int rows, int river, int drain)
Alllocate memory for the groundwater calculation data structure in 2 dimensions.
Definition: n_gwflow.c:148
void N_gwflow_2d_calc_water_budget(N_gwflow_data2d *data, N_geom_data *geom, N_array_2d *budget)
This function computes the water budget of the entire groundwater.
Definition: n_gwflow.c:659
void N_free_gwflow_data3d(N_gwflow_data3d *data)
Release the memory of the groundwater flow data structure in three dimensions.
Definition: n_gwflow.c:92
N_gwflow_data3d * N_alloc_gwflow_data3d(int cols, int rows, int depths, int river, int drain)
Alllocate memory for the groundwater calculation data structure in 3 dimensions.
Definition: n_gwflow.c:38
void N_free_gwflow_data2d(N_gwflow_data2d *data)
Release the memory of the groundwater flow data structure in two dimensions.
Definition: n_gwflow.c:200
N_data_star * N_callback_gwflow_3d(void *gwdata, N_geom_data *geom, int col, int row, int depth)
This callback function creates the mass balance of a 7 point star.
Definition: n_gwflow.c:263
N_data_star * N_create_5star(double C, double W, double E, double N, double S, double V)
allocate and initialize a 5 point star data structure
N_data_star * N_create_7star(double C, double W, double E, double N, double S, double T, double B, double V)
allocate and initialize a 7 point star data structure
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition: percent.c:62
int cols
Definition: N_pde.h:134
int rows
Definition: N_pde.h:134
int rows
Definition: N_pde.h:168
int depths
Definition: N_pde.h:168
int cols
Definition: N_pde.h:168
Matrix entries for a mass balance 5/7/9 star system.
Definition: N_pde.h:273
double E
Definition: N_pde.h:276
double S
Definition: N_pde.h:276
double W
Definition: N_pde.h:276
double T
Definition: N_pde.h:278
double N
Definition: N_pde.h:276
double B
Definition: N_pde.h:280
Geometric information about the structured grid.
Definition: N_pde.h:104
double dx
Definition: N_pde.h:109
int depths
Definition: N_pde.h:115
double dy
Definition: N_pde.h:110
double dz
Definition: N_pde.h:111
This data structure contains all data needed to compute the groundwater mass balance in two dimension...
Definition: N_gwflow.h:68
N_array_2d * hc_y
Definition: N_gwflow.h:72
N_array_2d * nf
Definition: N_gwflow.h:76
N_array_2d * top
Definition: N_gwflow.h:88
N_array_2d * drain_leak
Definition: N_gwflow.h:84
N_array_2d * bottom
Definition: N_gwflow.h:89
N_array_2d * phead_start
Definition: N_gwflow.h:70
N_array_2d * river_leak
Definition: N_gwflow.h:79
N_array_2d * s
Definition: N_gwflow.h:75
N_array_2d * hc_x
Definition: N_gwflow.h:71
N_array_2d * r
Definition: N_gwflow.h:74
N_array_2d * drain_bed
Definition: N_gwflow.h:85
N_array_2d * q
Definition: N_gwflow.h:73
N_array_2d * river_bed
Definition: N_gwflow.h:81
N_array_2d * river_head
Definition: N_gwflow.h:80
N_array_2d * phead
Definition: N_gwflow.h:69
N_array_2d * status
Definition: N_gwflow.h:91
This data structure contains all data needed to compute the groundwater mass balance in three dimensi...
Definition: N_gwflow.h:36
N_array_3d * phead
Definition: N_gwflow.h:37
N_array_3d * hc_z
Definition: N_gwflow.h:41
N_array_3d * phead_start
Definition: N_gwflow.h:38
N_array_3d * hc_x
Definition: N_gwflow.h:39
N_array_3d * drain_leak
Definition: N_gwflow.h:53
N_array_3d * status
Definition: N_gwflow.h:56
N_array_3d * drain_bed
Definition: N_gwflow.h:54
N_array_3d * river_bed
Definition: N_gwflow.h:50
N_array_3d * river_head
Definition: N_gwflow.h:49
N_array_3d * s
Definition: N_gwflow.h:44
N_array_2d * r
Definition: N_gwflow.h:43
N_array_3d * q
Definition: N_gwflow.h:42
N_array_3d * hc_y
Definition: N_gwflow.h:40
N_array_3d * nf
Definition: N_gwflow.h:45
N_array_3d * river_leak
Definition: N_gwflow.h:48
#define x