GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
n_gwflow.c
Go to the documentation of this file.
1/*****************************************************************************
2 *
3 * MODULE: Grass PDE Numerical Library
4 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
5 * soerengebbert <at> gmx <dot> de
6 *
7 * PURPOSE: groundwater flow in porous media
8 * part of the gpde library
9 *
10 * COPYRIGHT: (C) 2000 by the GRASS Development Team
11 *
12 * This program is free software under the GNU General Public
13 * License (>=v2). Read the file COPYING that comes with GRASS
14 * for details.
15 *
16 *****************************************************************************/
17
18#include <grass/N_gwflow.h>
19
20/* *************************************************************** */
21/* ***************** N_gwflow_data3d ***************************** */
22/* *************************************************************** */
23/*!
24 * \brief Allocate memory for the groundwater calculation data structure in 3
25 * 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
30 * border. This data structure is used to create a linear equation system based
31 * on the computation of groundwater flow in porous media with the finite volume
32 * method.
33 *
34 * \param cols int
35 * \param rows int
36 * \param depths int
37 * \return N_gwflow_data3d *
38 * */
39N_gwflow_data3d *N_alloc_gwflow_data3d(int cols, int rows, int depths,
40 int river, int drain)
41{
42 N_gwflow_data3d *data;
43
44 data = (N_gwflow_data3d *)G_calloc(1, sizeof(N_gwflow_data3d));
45
46 data->phead = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
47 data->phead_start = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
48 data->status = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
49 data->hc_x = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
50 data->hc_y = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
51 data->hc_z = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
52 data->q = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
53 data->s = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
54 data->nf = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
55 data->r = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
56
57 if (river) {
58 data->river_head = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
59 data->river_leak = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
60 data->river_bed = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
61 }
62 else {
63 data->river_head = NULL;
64 data->river_leak = NULL;
65 data->river_bed = NULL;
66 }
67
68 if (drain) {
69 data->drain_leak = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
70 data->drain_bed = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
71 }
72 else {
73 data->drain_leak = NULL;
74 data->drain_bed = NULL;
75 }
76
77 return data;
78}
79
80/* *************************************************************** */
81/* ********************* N_free_gwflow_data3d ******************** */
82/* *************************************************************** */
83/*!
84 * \brief Release the memory of the groundwater flow data structure in three
85 * dimensions
86 *
87 * \param data N_gwflow_data3d *
88 * \return void *
89 * */
90
92{
93 if (data->phead)
95 if (data->phead_start)
97 if (data->status)
99 if (data->hc_x)
100 N_free_array_3d(data->hc_x);
101 if (data->hc_y)
102 N_free_array_3d(data->hc_y);
103 if (data->hc_z)
104 N_free_array_3d(data->hc_z);
105 if (data->q)
106 N_free_array_3d(data->q);
107 if (data->s)
108 N_free_array_3d(data->s);
109 if (data->nf)
110 N_free_array_3d(data->nf);
111 if (data->r)
112 N_free_array_2d(data->r);
113 if (data->river_head)
115 if (data->river_leak)
117 if (data->river_bed)
119 if (data->drain_leak)
121 if (data->drain_bed)
123
124 G_free(data);
125
126 data = NULL;
127
128 return;
129}
130
131/* *************************************************************** */
132/* ******************** N_alloc_gwflow_data2d ******************** */
133/* *************************************************************** */
134/*!
135 * \brief Allocate memory for the groundwater calculation data structure in 2
136 * 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
141 * border. This data structure is used to create a linear equation system based
142 * on the computation of groundwater flow in porous media with the finite volume
143 * method.
144 *
145 * \param cols int
146 * \param rows int
147 * \return N_gwflow_data2d *
148 * */
149N_gwflow_data2d *N_alloc_gwflow_data2d(int cols, int rows, int river, 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 return data;
188}
189
190/* *************************************************************** */
191/* ****************** N_free_gwflow_data2d *********************** */
192/* *************************************************************** */
193/*!
194 * \brief Release the memory of the groundwater flow data structure in two
195 * 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
240 return;
241}
242
243/* *************************************************************** */
244/* ***************** N_callback_gwflow_3d ************************ */
245/* *************************************************************** */
246/*!
247 * \brief This callback function creates the mass balance of a 7 point star
248 *
249 * The mass balance is based on the common groundwater flow equation:
250 *
251 * \f[Ss \frac{\partial h}{\partial t} = \nabla {\bf K} \nabla h + q \f]
252 *
253 * This equation is discretizised with the finite volume method in three
254 * dimensions.
255 *
256 *
257 * \param gwdata N_gwflow_data3d *
258 * \param geom N_geom_data *
259 * \param col int
260 * \param row int
261 * \param depth int
262 * \return N_data_star *
263 *
264 * */
265N_data_star *N_callback_gwflow_3d(void *gwdata, N_geom_data *geom, int col,
266 int row, int depth)
267{
268 double hc_e = 0, hc_w = 0, hc_n = 0, hc_s = 0, hc_t = 0, hc_b = 0;
269 double dx, dy, dz, Ax, Ay, Az;
270 double hc_x, hc_y, hc_z;
271 double hc_xw, hc_yn, hc_zt;
272 double hc_xe, hc_ys, hc_zb;
273 double hc_start;
274 double Ss, r, /* nf, */ q;
275 double C, W, E, N, S, T, B, V;
276 N_data_star *mat_pos;
277 N_gwflow_data3d *data;
278
279 /*cast the void pointer to the right data structure */
280 data = (N_gwflow_data3d *)gwdata;
281
282 dx = geom->dx;
283 dy = geom->dy;
284 dz = geom->dz;
285 Az = N_get_geom_data_area_of_cell(geom, row);
286 Ay = geom->dx * geom->dz;
287 Ax = geom->dz * geom->dy;
288
289 /*read the data from the arrays */
290 hc_start = N_get_array_3d_d_value(data->phead_start, col, row, depth);
291
292 hc_x = N_get_array_3d_d_value(data->hc_x, col, row, depth);
293 hc_y = N_get_array_3d_d_value(data->hc_y, col, row, depth);
294 hc_z = N_get_array_3d_d_value(data->hc_z, col, row, depth);
295
296 hc_xw = N_get_array_3d_d_value(data->hc_x, col - 1, row, depth);
297 hc_xe = N_get_array_3d_d_value(data->hc_x, col + 1, row, depth);
298 hc_yn = N_get_array_3d_d_value(data->hc_y, col, row - 1, depth);
299 hc_ys = N_get_array_3d_d_value(data->hc_y, col, row + 1, depth);
300 hc_zt = N_get_array_3d_d_value(data->hc_z, col, row, depth + 1);
301 hc_zb = N_get_array_3d_d_value(data->hc_z, col, row, depth - 1);
302
303 hc_w = N_calc_harmonic_mean(hc_xw, hc_x);
304 hc_e = N_calc_harmonic_mean(hc_xe, hc_x);
305 hc_n = N_calc_harmonic_mean(hc_yn, hc_y);
306 hc_s = N_calc_harmonic_mean(hc_ys, hc_y);
307 hc_t = N_calc_harmonic_mean(hc_zt, hc_z);
308 hc_b = N_calc_harmonic_mean(hc_zb, hc_z);
309
310 /*inner sources */
311 q = N_get_array_3d_d_value(data->q, col, row, depth);
312 /*storativity */
313 Ss = N_get_array_3d_d_value(data->s, col, row, depth);
314 /*porosity */
315 /* nf = N_get_array_3d_d_value(data->nf, col, row, depth); */
316
317 /*mass balance center cell to western cell */
318 W = -1 * Ax * hc_w / dx;
319 /*mass balance center cell to eastern cell */
320 E = -1 * Ax * hc_e / dx;
321 /*mass balance center cell to northern cell */
322 N = -1 * Ay * hc_n / dy;
323 /*mass balance center cell to southern cell */
324 S = -1 * Ay * hc_s / dy;
325 /*mass balance center cell to top cell */
326 T = -1 * Az * hc_t / dz;
327 /*mass balance center cell to bottom cell */
328 B = -1 * Az * hc_b / dz;
329
330 /*storativity */
331 Ss = Az * dz * Ss;
332
333 /*the diagonal entry of the matrix */
334 C = -1 * (W + E + N + S + T + B - Ss / data->dt * Az);
335
336 /*the entry in the right side b of Ax = b */
337 V = (q + hc_start * Ss / data->dt * Az);
338
339 /*only the top cells will have recharge */
340 if (depth == geom->depths - 2) {
341 r = N_get_array_2d_d_value(data->r, col, row);
342 V += r * Az;
343 }
344
345 G_debug(5, "N_callback_gwflow_3d: called [%i][%i][%i]", depth, col, row);
346
347 /*create the 7 point star entries */
348 mat_pos = N_create_7star(C, W, E, N, S, T, B, V);
349
350 return mat_pos;
351}
352
353/* *************************************************************** */
354/* ****************** N_gwflow_3d_calc_water_budget ************** */
355/* *************************************************************** */
356/*!
357 * \brief This function computes the water budget of the entire groundwater
358 *
359 * The water budget is calculated for each active and dirichlet cell from
360 * its surrounding neighbours. This is based on the 7 star mass balance
361 * computation of N_callback_gwflow_3d and the gradient of the water heights in
362 * the cells. The sum of the water budget of each active/dirichlet cell must be
363 * near zero due the effect of numerical inaccuracy of cpu's.
364 *
365 * \param gwdata N_gwflow_data3d *
366 * \param geom N_geom_data *
367 * \param budget N_array_3d
368 * \return void
369 *
370 * */
372 N_array_3d *budget)
373{
374 int z, y, x, stat;
375 double h, hc;
376 double val;
377 double sum;
378 N_data_star *dstar;
379
380 int rows = data->status->rows;
381 int cols = data->status->cols;
382 int depths = data->status->depths;
383
384 sum = 0;
385
386 for (z = 0; z < depths; z++) {
387 for (y = 0; y < rows; y++) {
388 G_percent(y, rows - 1, 10);
389 for (x = 0; x < cols; x++) {
390 stat = (int)N_get_array_3d_d_value(data->status, x, y, z);
391
392 val = 0.0;
393
394 if (stat != N_CELL_INACTIVE) { /*all active/dirichlet cells */
395
396 /* Compute the flow parameter */
397 dstar = N_callback_gwflow_3d(data, geom, x, y, z);
398 /* Compute the gradient in each direction pointing from the
399 * center */
400 hc = N_get_array_3d_d_value(data->phead, x, y, z);
401
402 if ((int)N_get_array_3d_d_value(data->status, x + 1, y,
403 z) != N_CELL_INACTIVE) {
404 h = N_get_array_3d_d_value(data->phead, x + 1, y, z);
405 val += dstar->E * (hc - h);
406 }
407 if ((int)N_get_array_3d_d_value(data->status, x - 1, y,
408 z) != N_CELL_INACTIVE) {
409 h = N_get_array_3d_d_value(data->phead, x - 1, y, z);
410 val += dstar->W * (hc - h);
411 }
412 if ((int)N_get_array_3d_d_value(data->status, x, y + 1,
413 z) != N_CELL_INACTIVE) {
414 h = N_get_array_3d_d_value(data->phead, x, y + 1, z);
415 val += dstar->S * (hc - h);
416 }
417 if ((int)N_get_array_3d_d_value(data->status, x, y - 1,
418 z) != N_CELL_INACTIVE) {
419 h = N_get_array_3d_d_value(data->phead, x, y - 1, z);
420 val += dstar->N * (hc - h);
421 }
422 if ((int)N_get_array_3d_d_value(data->status, x, y,
423 z + 1) != N_CELL_INACTIVE) {
424 h = N_get_array_3d_d_value(data->phead, x, y, z + 1);
425 val += dstar->T * (hc - h);
426 }
427 if ((int)N_get_array_3d_d_value(data->status, x, y,
428 z - 1) != N_CELL_INACTIVE) {
429 h = N_get_array_3d_d_value(data->phead, x, y, z - 1);
430 val += dstar->B * (hc - h);
431 }
432 sum += val;
433
434 G_free(dstar);
435 }
436 else {
437 Rast_set_null_value(&val, 1, DCELL_TYPE);
438 }
439 N_put_array_3d_d_value(budget, x, y, z, val);
440 }
441 }
442 }
443
444 if (fabs(sum) < 0.0000000001)
445 G_message(_("The total sum of the water budget: %g\n"), sum);
446 else
447 G_warning(_("The total sum of the water budget is significantly larger "
448 "then 0: %g\n"),
449 sum);
450
451 return;
452}
453
454/* *************************************************************** */
455/* ****************** N_callback_gwflow_2d *********************** */
456/* *************************************************************** */
457/*!
458 * \brief This callback function creates the mass balance of a 5 point star
459 *
460 * The mass balance is based on the common groundwater flow equation:
461 *
462 * \f[Ss \frac{\partial h}{\partial t} = \nabla {\bf K} \nabla h + q \f]
463 *
464 * This equation is discretizised with the finite volume method in two
465 * dimensions.
466 *
467 * \param gwdata N_gwflow_data2d *
468 * \param geom N_geom_data *
469 * \param col int
470 * \param row int
471 * \return N_data_star *
472 *
473 * */
474N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data *geom, int col,
475 int row)
476{
477 double T_e = 0, T_w = 0, T_n = 0, T_s = 0;
478 double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
479 double dx, dy, Az;
480 double hc_x, hc_y;
481 double z, top;
482 double hc_xw, hc_yn;
483 double z_xw, z_yn;
484 double hc_xe, hc_ys;
485 double z_xe, z_ys;
486 double hc, hc_start;
487 double Ss, r, q;
488 double C, W, E, N, S, V;
489 N_gwflow_data2d *data;
490 N_data_star *mat_pos;
491 double river_vect = 0; /*entry in vector */
492 double river_mat = 0; /*entry in matrix */
493 double drain_vect = 0; /*entry in vector */
494 double drain_mat = 0; /*entry in matrix */
495
496 /*cast the void pointer to the right data structure */
497 data = (N_gwflow_data2d *)gwdata;
498
499 dx = geom->dx;
500 dy = geom->dy;
501 Az = N_get_geom_data_area_of_cell(geom, row);
502
503 /*read the data from the arrays */
504 hc_start = N_get_array_2d_d_value(data->phead_start, col, row);
505 hc = N_get_array_2d_d_value(data->phead, col, row);
506 top = N_get_array_2d_d_value(data->top, col, row);
507
508 /* Inner sources */
509 q = N_get_array_2d_d_value(data->q, col, row);
510
511 /* storativity or porosity of current cell face [-] */
512 Ss = N_get_array_2d_d_value(data->s, col, row);
513 /* recharge */
514 r = N_get_array_2d_d_value(data->r, col, row) * Az;
515
516 if (hc > top) { /*If the aquifer is confined */
517 z = N_get_array_2d_d_value(data->top, col, row) -
518 N_get_array_2d_d_value(data->bottom, col, row);
519 z_xw = N_get_array_2d_d_value(data->top, col - 1, row) -
520 N_get_array_2d_d_value(data->bottom, col - 1, row);
521 z_xe = N_get_array_2d_d_value(data->top, col + 1, row) -
522 N_get_array_2d_d_value(data->bottom, col + 1, row);
523 z_yn = N_get_array_2d_d_value(data->top, col, row - 1) -
524 N_get_array_2d_d_value(data->bottom, col, row - 1);
525 z_ys = N_get_array_2d_d_value(data->top, col, 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 explicit 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
598 * surface
599 */
600 if (data->drain_leak &&
601 (N_get_array_2d_d_value(data->drain_leak, col, row) != 0) &&
602 N_get_array_2d_d_value(data->drain_bed, col, row) <= top) {
603 if (hc > N_get_array_2d_d_value(data->drain_bed, col, row)) {
604 drain_vect = N_get_array_2d_d_value(data->drain_bed, col, row) *
605 N_get_array_2d_d_value(data->drain_leak, col, row);
606 drain_mat = N_get_array_2d_d_value(data->drain_leak, col, row);
607 }
608 else if (hc <= N_get_array_2d_d_value(data->drain_bed, col, row)) {
609 drain_vect = 0;
610 drain_mat = 0;
611 }
612 }
613
614 /*mass balance center cell to western cell */
615 W = -1 * T_w * dy / dx;
616 /*mass balance center cell to eastern cell */
617 E = -1 * T_e * dy / dx;
618 /*mass balance center cell to northern cell */
619 N = -1 * T_n * dx / dy;
620 /*mass balance center cell to southern cell */
621 S = -1 * T_s * dx / dy;
622
623 /*the diagonal entry of the matrix */
624 C = -1 *
625 (W + E + N + S - Az * Ss / data->dt - river_mat * Az - drain_mat * Az);
626
627 /*the entry in the right side b of Ax = b */
628 V = (q + hc_start * Az * Ss / data->dt) + r + river_vect * Az +
629 drain_vect * Az;
630
631 G_debug(5, "N_callback_gwflow_2d: called [%i][%i]", row, col);
632
633 /*create the 5 point star entries */
634 mat_pos = N_create_5star(C, W, E, N, S, V);
635
636 return mat_pos;
637}
638
639/* *************************************************************** */
640/* ****************** N_gwflow_2d_calc_water_budget ************** */
641/* *************************************************************** */
642/*!
643 * \brief This function computes the water budget of the entire groundwater
644 *
645 * The water budget is calculated for each active and dirichlet cell from
646 * its surrounding neighbours. This is based on the 5 star mass balance
647 * computation of N_callback_gwflow_2d and the gradient of the water heights in
648 * the cells. The sum of the water budget of each active/dirichlet cell must be
649 * near zero due the effect of numerical inaccuracy of cpu's.
650 *
651 * \param gwdata N_gwflow_data2d *
652 * \param geom N_geom_data *
653 * \param budget N_array_2d
654 * \return void
655 *
656 * */
658 N_array_2d *budget)
659{
660 int y, x, stat;
661 double h, hc;
662 double val;
663 double sum;
664 N_data_star *dstar;
665
666 int rows = data->status->rows;
667 int cols = data->status->cols;
668
669 sum = 0;
670
671 for (y = 0; y < rows; y++) {
672 G_percent(y, rows - 1, 10);
673 for (x = 0; x < cols; x++) {
674 stat = N_get_array_2d_c_value(data->status, x, y);
675
676 val = 0.0;
677
678 if (stat != N_CELL_INACTIVE) { /*all active/dirichlet cells */
679
680 /* Compute the flow parameter */
681 dstar = N_callback_gwflow_2d(data, geom, x, y);
682 /* Compute the gradient in each direction pointing from the
683 * 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) !=
688 h = N_get_array_2d_d_value(data->phead, x + 1, y);
689 val += dstar->E * (hc - h);
690 }
691 if ((int)N_get_array_2d_d_value(data->status, x - 1, y) !=
693 h = N_get_array_2d_d_value(data->phead, x - 1, y);
694 val += dstar->W * (hc - h);
695 }
696 if ((int)N_get_array_2d_d_value(data->status, x, y + 1) !=
698 h = N_get_array_2d_d_value(data->phead, x, y + 1);
699 val += dstar->S * (hc - h);
700 }
701 if ((int)N_get_array_2d_d_value(data->status, x, y - 1) !=
703 h = N_get_array_2d_d_value(data->phead, x, y - 1);
704 val += dstar->N * (hc - h);
705 }
706
707 sum += val;
708
709 G_free(dstar);
710 }
711 else {
712 Rast_set_null_value(&val, 1, DCELL_TYPE);
713 }
714 N_put_array_2d_d_value(budget, x, y, val);
715 }
716 }
717
718 if (fabs(sum) < 0.0000000001)
719 G_message(_("The total sum of the water budget: %g\n"), sum);
720 else
721 G_warning(_("The total sum of the water budget is significantly larger "
722 "then 0: %g\n"),
723 sum);
724
725 return;
726}
#define N_CELL_INACTIVE
Definition N_pde.h:30
double N_calc_arith_mean(double a, double b)
Calculate the arithmetic mean of values a and b.
Definition n_tools.c:31
double N_calc_harmonic_mean(double a, double b)
Calculate the harmonical mean of values a and b.
Definition n_tools.c:115
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
#define NULL
Definition ccmath.h:32
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
double r
void G_message(const char *msg,...)
Print a message to stderr.
Definition gis/error.c:89
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
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:719
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:774
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:380
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d structure.
Definition n_arrays.c:132
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:1148
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:75
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:979
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:576
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:196
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:474
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)
Allocate memory for the groundwater calculation data structure in 2 dimensions.
Definition n_gwflow.c:149
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:657
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:91
N_gwflow_data3d * N_alloc_gwflow_data3d(int cols, int rows, int depths, int river, int drain)
Allocate memory for the groundwater calculation data structure in 3 dimensions.
Definition n_gwflow.c:39
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:265
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:61
int cols
Definition N_pde.h:134
int rows
Definition N_pde.h:134
int rows
Definition N_pde.h:177
int depths
Definition N_pde.h:177
int cols
Definition N_pde.h:177
Matrix entries for a mass balance 5/7/9 star system.
Definition N_pde.h:295
double E
Definition N_pde.h:298
double S
Definition N_pde.h:298
double W
Definition N_pde.h:298
double T
Definition N_pde.h:300
double N
Definition N_pde.h:298
double B
Definition N_pde.h:302
Geometric information about the structured grid.
Definition N_pde.h:101
double dx
Definition N_pde.h:108
int depths
Definition N_pde.h:114
double dy
Definition N_pde.h:109
double dz
Definition N_pde.h:110
This data structure contains all data needed to compute the groundwater mass balance in two dimension...
Definition N_gwflow.h:65
N_array_2d * hc_y
Definition N_gwflow.h:69
N_array_2d * nf
Definition N_gwflow.h:73
N_array_2d * top
Definition N_gwflow.h:84
N_array_2d * drain_leak
Definition N_gwflow.h:81
N_array_2d * bottom
Definition N_gwflow.h:85
N_array_2d * phead_start
Definition N_gwflow.h:67
N_array_2d * river_leak
Definition N_gwflow.h:76
N_array_2d * s
Definition N_gwflow.h:72
N_array_2d * hc_x
Definition N_gwflow.h:68
N_array_2d * r
Definition N_gwflow.h:71
N_array_2d * drain_bed
Definition N_gwflow.h:82
N_array_2d * q
Definition N_gwflow.h:70
N_array_2d * river_bed
Definition N_gwflow.h:78
N_array_2d * river_head
Definition N_gwflow.h:77
N_array_2d * phead
Definition N_gwflow.h:66
N_array_2d * status
Definition N_gwflow.h:87
This data structure contains all data needed to compute the groundwater mass balance in three dimensi...
Definition N_gwflow.h:34
N_array_3d * phead
Definition N_gwflow.h:35
N_array_3d * hc_z
Definition N_gwflow.h:39
N_array_3d * phead_start
Definition N_gwflow.h:36
N_array_3d * hc_x
Definition N_gwflow.h:37
N_array_3d * drain_leak
Definition N_gwflow.h:51
N_array_3d * status
Definition N_gwflow.h:54
N_array_3d * drain_bed
Definition N_gwflow.h:52
N_array_3d * river_bed
Definition N_gwflow.h:48
N_array_3d * river_head
Definition N_gwflow.h:47
N_array_3d * s
Definition N_gwflow.h:42
N_array_2d * r
Definition N_gwflow.h:41
N_array_3d * q
Definition N_gwflow.h:40
N_array_3d * hc_y
Definition N_gwflow.h:38
N_array_3d * nf
Definition N_gwflow.h:43
N_array_3d * river_leak
Definition N_gwflow.h:46
#define x