GRASS 8 Programmer's Manual  8.5.0dev(2025)-c070206eb1
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  * */
39 N_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)
94  N_free_array_3d(data->phead);
95  if (data->phead_start)
97  if (data->status)
98  N_free_array_3d(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)
118  N_free_array_3d(data->river_bed);
119  if (data->drain_leak)
121  if (data->drain_bed)
122  N_free_array_3d(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  * \param river
148  * \param drain
149  * \return N_gwflow_data2d *
150  * */
151 N_gwflow_data2d *N_alloc_gwflow_data2d(int cols, int rows, int river, int drain)
152 {
153  N_gwflow_data2d *data;
154 
155  data = (N_gwflow_data2d *)G_calloc(1, sizeof(N_gwflow_data2d));
156 
157  data->phead = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
158  data->phead_start = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
159  data->status = N_alloc_array_2d(cols, rows, 1, CELL_TYPE);
160  data->hc_x = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
161  data->hc_y = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
162  data->q = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
163  data->s = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
164  data->nf = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
165  data->r = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
166  data->top = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
167  data->bottom = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
168 
169  if (river) {
170  data->river_head = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
171  data->river_leak = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
172  data->river_bed = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
173  }
174  else {
175  data->river_head = NULL;
176  data->river_leak = NULL;
177  data->river_bed = NULL;
178  }
179 
180  if (drain) {
181  data->drain_leak = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
182  data->drain_bed = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
183  }
184  else {
185  data->drain_leak = NULL;
186  data->drain_bed = NULL;
187  }
188 
189  return data;
190 }
191 
192 /* *************************************************************** */
193 /* ****************** N_free_gwflow_data2d *********************** */
194 /* *************************************************************** */
195 /*!
196  * \brief Release the memory of the groundwater flow data structure in two
197  * dimensions
198  *
199  * \param data N_gwflow_data2d *
200  * \return void
201  * */
203 {
204  if (data->phead)
205  N_free_array_2d(data->phead);
206  if (data->phead_start)
208  if (data->status)
209  N_free_array_2d(data->status);
210  if (data->hc_x)
211  N_free_array_2d(data->hc_x);
212  if (data->hc_y)
213  N_free_array_2d(data->hc_y);
214  if (data->q)
215  N_free_array_2d(data->q);
216  if (data->s)
217  N_free_array_2d(data->s);
218  if (data->nf)
219  N_free_array_2d(data->nf);
220  if (data->r)
221  N_free_array_2d(data->r);
222  if (data->top)
223  N_free_array_2d(data->top);
224  if (data->bottom)
225  N_free_array_2d(data->bottom);
226  if (data->river_head)
228  if (data->river_leak)
230  if (data->river_bed)
231  N_free_array_2d(data->river_bed);
232  if (data->drain_leak)
234  if (data->drain_bed)
235  N_free_array_2d(data->drain_bed);
236 
237  G_free(data);
238 
239  data = NULL;
240  ;
241 
242  return;
243 }
244 
245 /* *************************************************************** */
246 /* ***************** N_callback_gwflow_3d ************************ */
247 /* *************************************************************** */
248 /*!
249  * \brief This callback function creates the mass balance of a 7 point star
250  *
251  * The mass balance is based on the common groundwater flow equation:
252  *
253  * \f[Ss \frac{\partial h}{\partial t} = \nabla {\bf K} \nabla h + q \f]
254  *
255  * This equation is discretizised with the finite volume method in three
256  * dimensions.
257  *
258  *
259  * \param gwdata N_gwflow_data3d *
260  * \param geom N_geom_data *
261  * \param col int
262  * \param row int
263  * \param depth int
264  * \return N_data_star *
265  *
266  * */
267 N_data_star *N_callback_gwflow_3d(void *gwdata, N_geom_data *geom, int col,
268  int row, int depth)
269 {
270  double hc_e = 0, hc_w = 0, hc_n = 0, hc_s = 0, hc_t = 0, hc_b = 0;
271  double dx, dy, dz, Ax, Ay, Az;
272  double hc_x, hc_y, hc_z;
273  double hc_xw, hc_yn, hc_zt;
274  double hc_xe, hc_ys, hc_zb;
275  double hc_start;
276  double Ss, r, /* nf, */ q;
277  double C, W, E, N, S, T, B, V;
278  N_data_star *mat_pos;
279  N_gwflow_data3d *data;
280 
281  /*cast the void pointer to the right data structure */
282  data = (N_gwflow_data3d *)gwdata;
283 
284  dx = geom->dx;
285  dy = geom->dy;
286  dz = geom->dz;
287  Az = N_get_geom_data_area_of_cell(geom, row);
288  Ay = geom->dx * geom->dz;
289  Ax = geom->dz * geom->dy;
290 
291  /*read the data from the arrays */
292  hc_start = N_get_array_3d_d_value(data->phead_start, col, row, depth);
293 
294  hc_x = N_get_array_3d_d_value(data->hc_x, col, row, depth);
295  hc_y = N_get_array_3d_d_value(data->hc_y, col, row, depth);
296  hc_z = N_get_array_3d_d_value(data->hc_z, col, row, depth);
297 
298  hc_xw = N_get_array_3d_d_value(data->hc_x, col - 1, row, depth);
299  hc_xe = N_get_array_3d_d_value(data->hc_x, col + 1, row, depth);
300  hc_yn = N_get_array_3d_d_value(data->hc_y, col, row - 1, depth);
301  hc_ys = N_get_array_3d_d_value(data->hc_y, col, row + 1, depth);
302  hc_zt = N_get_array_3d_d_value(data->hc_z, col, row, depth + 1);
303  hc_zb = N_get_array_3d_d_value(data->hc_z, col, row, depth - 1);
304 
305  hc_w = N_calc_harmonic_mean(hc_xw, hc_x);
306  hc_e = N_calc_harmonic_mean(hc_xe, hc_x);
307  hc_n = N_calc_harmonic_mean(hc_yn, hc_y);
308  hc_s = N_calc_harmonic_mean(hc_ys, hc_y);
309  hc_t = N_calc_harmonic_mean(hc_zt, hc_z);
310  hc_b = N_calc_harmonic_mean(hc_zb, hc_z);
311 
312  /*inner sources */
313  q = N_get_array_3d_d_value(data->q, col, row, depth);
314  /*storativity */
315  Ss = N_get_array_3d_d_value(data->s, col, row, depth);
316  /*porosity */
317  /* nf = N_get_array_3d_d_value(data->nf, col, row, depth); */
318 
319  /*mass balance center cell to western cell */
320  W = -1 * Ax * hc_w / dx;
321  /*mass balance center cell to eastern cell */
322  E = -1 * Ax * hc_e / dx;
323  /*mass balance center cell to northern cell */
324  N = -1 * Ay * hc_n / dy;
325  /*mass balance center cell to southern cell */
326  S = -1 * Ay * hc_s / dy;
327  /*mass balance center cell to top cell */
328  T = -1 * Az * hc_t / dz;
329  /*mass balance center cell to bottom cell */
330  B = -1 * Az * hc_b / dz;
331 
332  /*storativity */
333  Ss = Az * dz * Ss;
334 
335  /*the diagonal entry of the matrix */
336  C = -1 * (W + E + N + S + T + B - Ss / data->dt * Az);
337 
338  /*the entry in the right side b of Ax = b */
339  V = (q + hc_start * Ss / data->dt * Az);
340 
341  /*only the top cells will have recharge */
342  if (depth == geom->depths - 2) {
343  r = N_get_array_2d_d_value(data->r, col, row);
344  V += r * Az;
345  }
346 
347  G_debug(5, "N_callback_gwflow_3d: called [%i][%i][%i]", depth, col, row);
348 
349  /*create the 7 point star entries */
350  mat_pos = N_create_7star(C, W, E, N, S, T, B, V);
351 
352  return mat_pos;
353 }
354 
355 /* *************************************************************** */
356 /* ****************** N_gwflow_3d_calc_water_budget ************** */
357 /* *************************************************************** */
358 /*!
359  * \brief This function computes the water budget of the entire groundwater
360  *
361  * The water budget is calculated for each active and dirichlet cell from
362  * its surrounding neighbours. This is based on the 7 star mass balance
363  * computation of N_callback_gwflow_3d and the gradient of the water heights in
364  * the cells. The sum of the water budget of each active/dirichlet cell must be
365  * near zero due the effect of numerical inaccuracy of cpu's.
366  *
367  * \param gwdata N_gwflow_data3d *
368  * \param geom N_geom_data *
369  * \param budget N_array_3d
370  * \return void
371  *
372  * */
374  N_array_3d *budget)
375 {
376  int z, y, x, stat;
377  double h, hc;
378  double val;
379  double sum;
380  N_data_star *dstar;
381 
382  int rows = data->status->rows;
383  int cols = data->status->cols;
384  int depths = data->status->depths;
385 
386  sum = 0;
387 
388  for (z = 0; z < depths; z++) {
389  for (y = 0; y < rows; y++) {
390  G_percent(y, rows - 1, 10);
391  for (x = 0; x < cols; x++) {
392  stat = (int)N_get_array_3d_d_value(data->status, x, y, z);
393 
394  val = 0.0;
395 
396  if (stat != N_CELL_INACTIVE) { /*all active/dirichlet cells */
397 
398  /* Compute the flow parameter */
399  dstar = N_callback_gwflow_3d(data, geom, x, y, z);
400  /* Compute the gradient in each direction pointing from the
401  * center */
402  hc = N_get_array_3d_d_value(data->phead, x, y, z);
403 
404  if ((int)N_get_array_3d_d_value(data->status, x + 1, y,
405  z) != N_CELL_INACTIVE) {
406  h = N_get_array_3d_d_value(data->phead, x + 1, y, z);
407  val += dstar->E * (hc - h);
408  }
409  if ((int)N_get_array_3d_d_value(data->status, x - 1, y,
410  z) != N_CELL_INACTIVE) {
411  h = N_get_array_3d_d_value(data->phead, x - 1, y, z);
412  val += dstar->W * (hc - h);
413  }
414  if ((int)N_get_array_3d_d_value(data->status, x, y + 1,
415  z) != N_CELL_INACTIVE) {
416  h = N_get_array_3d_d_value(data->phead, x, y + 1, z);
417  val += dstar->S * (hc - h);
418  }
419  if ((int)N_get_array_3d_d_value(data->status, x, y - 1,
420  z) != N_CELL_INACTIVE) {
421  h = N_get_array_3d_d_value(data->phead, x, y - 1, z);
422  val += dstar->N * (hc - h);
423  }
424  if ((int)N_get_array_3d_d_value(data->status, x, y,
425  z + 1) != N_CELL_INACTIVE) {
426  h = N_get_array_3d_d_value(data->phead, x, y, z + 1);
427  val += dstar->T * (hc - h);
428  }
429  if ((int)N_get_array_3d_d_value(data->status, x, y,
430  z - 1) != N_CELL_INACTIVE) {
431  h = N_get_array_3d_d_value(data->phead, x, y, z - 1);
432  val += dstar->B * (hc - h);
433  }
434  sum += val;
435 
436  G_free(dstar);
437  }
438  else {
440  }
441  N_put_array_3d_d_value(budget, x, y, z, val);
442  }
443  }
444  }
445 
446  if (fabs(sum) < 0.0000000001)
447  G_message(_("The total sum of the water budget: %g\n"), sum);
448  else
449  G_warning(_("The total sum of the water budget is significantly larger "
450  "then 0: %g\n"),
451  sum);
452 
453  return;
454 }
455 
456 /* *************************************************************** */
457 /* ****************** N_callback_gwflow_2d *********************** */
458 /* *************************************************************** */
459 /*!
460  * \brief This callback function creates the mass balance of a 5 point star
461  *
462  * The mass balance is based on the common groundwater flow equation:
463  *
464  * \f[Ss \frac{\partial h}{\partial t} = \nabla {\bf K} \nabla h + q \f]
465  *
466  * This equation is discretizised with the finite volume method in two
467  * dimensions.
468  *
469  * \param gwdata N_gwflow_data2d *
470  * \param geom N_geom_data *
471  * \param col int
472  * \param row int
473  * \return N_data_star *
474  *
475  * */
476 N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data *geom, int col,
477  int row)
478 {
479  double T_e = 0, T_w = 0, T_n = 0, T_s = 0;
480  double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
481  double dx, dy, Az;
482  double hc_x, hc_y;
483  double z, top;
484  double hc_xw, hc_yn;
485  double z_xw, z_yn;
486  double hc_xe, hc_ys;
487  double z_xe, z_ys;
488  double hc, hc_start;
489  double Ss, r, q;
490  double C, W, E, N, S, V;
491  N_gwflow_data2d *data;
492  N_data_star *mat_pos;
493  double river_vect = 0; /*entry in vector */
494  double river_mat = 0; /*entry in matrix */
495  double drain_vect = 0; /*entry in vector */
496  double drain_mat = 0; /*entry in matrix */
497 
498  /*cast the void pointer to the right data structure */
499  data = (N_gwflow_data2d *)gwdata;
500 
501  dx = geom->dx;
502  dy = geom->dy;
503  Az = N_get_geom_data_area_of_cell(geom, row);
504 
505  /*read the data from the arrays */
506  hc_start = N_get_array_2d_d_value(data->phead_start, col, row);
507  hc = N_get_array_2d_d_value(data->phead, col, row);
508  top = N_get_array_2d_d_value(data->top, col, row);
509 
510  /* Inner sources */
511  q = N_get_array_2d_d_value(data->q, col, row);
512 
513  /* storativity or porosity of current cell face [-] */
514  Ss = N_get_array_2d_d_value(data->s, col, row);
515  /* recharge */
516  r = N_get_array_2d_d_value(data->r, col, row) * Az;
517 
518  if (hc > top) { /*If the aquifer is confined */
519  z = N_get_array_2d_d_value(data->top, col, row) -
520  N_get_array_2d_d_value(data->bottom, col, row);
521  z_xw = 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_xe = N_get_array_2d_d_value(data->top, col + 1, row) -
524  N_get_array_2d_d_value(data->bottom, col + 1, row);
525  z_yn = N_get_array_2d_d_value(data->top, col, row - 1) -
526  N_get_array_2d_d_value(data->bottom, col, row - 1);
527  z_ys = N_get_array_2d_d_value(data->top, col, row + 1) -
528  N_get_array_2d_d_value(data->bottom, col, row + 1);
529  }
530  else { /* the aquifer is unconfined */
531 
532  /* If the aquifer is unconfied use an explicit scheme to solve
533  * the nonlinear equation. We use the phead from the first iteration */
534  z = N_get_array_2d_d_value(data->phead, col, row) -
535  N_get_array_2d_d_value(data->bottom, col, row);
536  z_xw = 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_xe = N_get_array_2d_d_value(data->phead, col + 1, row) -
539  N_get_array_2d_d_value(data->bottom, col + 1, row);
540  z_yn = N_get_array_2d_d_value(data->phead, col, row - 1) -
541  N_get_array_2d_d_value(data->bottom, col, row - 1);
542  z_ys = N_get_array_2d_d_value(data->phead, col, row + 1) -
543  N_get_array_2d_d_value(data->bottom, col, row + 1);
544  }
545 
546  /*geometrical mean of cell height */
547  if (z_w > 0 || z_w < 0 || z_w == 0)
548  z_w = N_calc_arith_mean(z_xw, z);
549  else
550  z_w = z;
551  if (z_e > 0 || z_e < 0 || z_e == 0)
552  z_e = N_calc_arith_mean(z_xe, z);
553  else
554  z_e = z;
555  if (z_n > 0 || z_n < 0 || z_n == 0)
556  z_n = N_calc_arith_mean(z_yn, z);
557  else
558  z_n = z;
559  if (z_s > 0 || z_s < 0 || z_s == 0)
560  z_s = N_calc_arith_mean(z_ys, z);
561  else
562  z_s = z;
563 
564  /*get the surrounding permeabilities */
565  hc_x = N_get_array_2d_d_value(data->hc_x, col, row);
566  hc_y = N_get_array_2d_d_value(data->hc_y, col, row);
567  hc_xw = N_get_array_2d_d_value(data->hc_x, col - 1, row);
568  hc_xe = N_get_array_2d_d_value(data->hc_x, col + 1, row);
569  hc_yn = N_get_array_2d_d_value(data->hc_y, col, row - 1);
570  hc_ys = N_get_array_2d_d_value(data->hc_y, col, row + 1);
571 
572  /* calculate the transmissivities */
573  T_w = N_calc_harmonic_mean(hc_xw, hc_x) * z_w;
574  T_e = N_calc_harmonic_mean(hc_xe, hc_x) * z_e;
575  T_n = N_calc_harmonic_mean(hc_yn, hc_y) * z_n;
576  T_s = N_calc_harmonic_mean(hc_ys, hc_y) * z_s;
577 
578  /* Compute the river leakage, this is an explicit method
579  * Influent and effluent flow is computed.
580  */
581  if (data->river_leak &&
582  (N_get_array_2d_d_value(data->river_leak, col, row) != 0) &&
583  N_get_array_2d_d_value(data->river_bed, col, row) <= top) {
584  /* Groundwater surface is above the river bed */
585  if (hc > N_get_array_2d_d_value(data->river_bed, col, row)) {
586  river_vect = N_get_array_2d_d_value(data->river_head, col, row) *
587  N_get_array_2d_d_value(data->river_leak, col, row);
588  river_mat = N_get_array_2d_d_value(data->river_leak, col, row);
589  } /* Groundwater surface is below the river bed */
590  else if (hc < N_get_array_2d_d_value(data->river_bed, col, row)) {
591  river_vect = (N_get_array_2d_d_value(data->river_head, col, row) -
592  N_get_array_2d_d_value(data->river_bed, col, row)) *
593  N_get_array_2d_d_value(data->river_leak, col, row);
594  river_mat = 0;
595  }
596  }
597 
598  /* compute the drainage, this is an explicit method
599  * Drainage is only enabled, if the drain bed is lower the groundwater
600  * surface
601  */
602  if (data->drain_leak &&
603  (N_get_array_2d_d_value(data->drain_leak, col, row) != 0) &&
604  N_get_array_2d_d_value(data->drain_bed, col, row) <= top) {
605  if (hc > N_get_array_2d_d_value(data->drain_bed, col, row)) {
606  drain_vect = N_get_array_2d_d_value(data->drain_bed, col, row) *
607  N_get_array_2d_d_value(data->drain_leak, col, row);
608  drain_mat = N_get_array_2d_d_value(data->drain_leak, col, row);
609  }
610  else if (hc <= N_get_array_2d_d_value(data->drain_bed, col, row)) {
611  drain_vect = 0;
612  drain_mat = 0;
613  }
614  }
615 
616  /*mass balance center cell to western cell */
617  W = -1 * T_w * dy / dx;
618  /*mass balance center cell to eastern cell */
619  E = -1 * T_e * dy / dx;
620  /*mass balance center cell to northern cell */
621  N = -1 * T_n * dx / dy;
622  /*mass balance center cell to southern cell */
623  S = -1 * T_s * dx / dy;
624 
625  /*the diagonal entry of the matrix */
626  C = -1 *
627  (W + E + N + S - Az * Ss / data->dt - river_mat * Az - drain_mat * Az);
628 
629  /*the entry in the right side b of Ax = b */
630  V = (q + hc_start * Az * Ss / data->dt) + r + river_vect * Az +
631  drain_vect * Az;
632 
633  G_debug(5, "N_callback_gwflow_2d: called [%i][%i]", row, col);
634 
635  /*create the 5 point star entries */
636  mat_pos = N_create_5star(C, W, E, N, S, V);
637 
638  return mat_pos;
639 }
640 
641 /* *************************************************************** */
642 /* ****************** N_gwflow_2d_calc_water_budget ************** */
643 /* *************************************************************** */
644 /*!
645  * \brief This function computes the water budget of the entire groundwater
646  *
647  * The water budget is calculated for each active and dirichlet cell from
648  * its surrounding neighbours. This is based on the 5 star mass balance
649  * computation of N_callback_gwflow_2d and the gradient of the water heights in
650  * the cells. The sum of the water budget of each active/dirichlet cell must be
651  * near zero due the effect of numerical inaccuracy of cpu's.
652  *
653  * \param data N_gwflow_data2d *
654  * \param geom N_geom_data *
655  * \param budget N_array_2d
656  * \return void
657  *
658  * */
660  N_array_2d *budget)
661 {
662  int y, x, stat;
663  double h, hc;
664  double val;
665  double sum;
666  N_data_star *dstar;
667 
668  int rows = data->status->rows;
669  int cols = data->status->cols;
670 
671  sum = 0;
672 
673  for (y = 0; y < rows; y++) {
674  G_percent(y, rows - 1, 10);
675  for (x = 0; x < cols; x++) {
676  stat = N_get_array_2d_c_value(data->status, x, y);
677 
678  val = 0.0;
679 
680  if (stat != N_CELL_INACTIVE) { /*all active/dirichlet cells */
681 
682  /* Compute the flow parameter */
683  dstar = N_callback_gwflow_2d(data, geom, x, y);
684  /* Compute the gradient in each direction pointing from the
685  * center */
686  hc = N_get_array_2d_d_value(data->phead, x, y);
687 
688  if ((int)N_get_array_2d_d_value(data->status, x + 1, y) !=
689  N_CELL_INACTIVE) {
690  h = N_get_array_2d_d_value(data->phead, x + 1, y);
691  val += dstar->E * (hc - h);
692  }
693  if ((int)N_get_array_2d_d_value(data->status, x - 1, y) !=
694  N_CELL_INACTIVE) {
695  h = N_get_array_2d_d_value(data->phead, x - 1, y);
696  val += dstar->W * (hc - h);
697  }
698  if ((int)N_get_array_2d_d_value(data->status, x, y + 1) !=
699  N_CELL_INACTIVE) {
700  h = N_get_array_2d_d_value(data->phead, x, y + 1);
701  val += dstar->S * (hc - h);
702  }
703  if ((int)N_get_array_2d_d_value(data->status, x, y - 1) !=
704  N_CELL_INACTIVE) {
705  h = N_get_array_2d_d_value(data->phead, x, y - 1);
706  val += dstar->N * (hc - h);
707  }
708 
709  sum += val;
710 
711  G_free(dstar);
712  }
713  else {
715  }
716  N_put_array_2d_d_value(budget, x, y, val);
717  }
718  }
719 
720  if (fabs(sum) < 0.0000000001)
721  G_message(_("The total sum of the water budget: %g\n"), sum);
722  else
723  G_warning(_("The total sum of the water budget is significantly larger "
724  "then 0: %g\n"),
725  sum);
726 
727  return;
728 }
#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
#define NULL
Definition: ccmath.h:32
void G_percent(long, long, int)
Print percent complete messages.
Definition: percent.c:61
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:147
#define G_calloc(m, n)
Definition: defs/gis.h:95
void G_warning(const char *,...) __attribute__((format(printf
void G_message(const char *,...) __attribute__((format(printf
int G_debug(int, const char *,...) __attribute__((format(printf
void Rast_set_null_value(void *, int, RASTER_MAP_TYPE)
To set one or more raster values to null.
Definition: null_val.c:98
#define N
Definition: e_intersect.c:926
#define _(str)
Definition: glocale.h:10
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
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
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
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
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
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:373
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:476
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
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:151
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
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:267
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:202
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
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
#define W
Definition: ogsf.h:143
double r
Definition: r_raster.c:39
#define DCELL_TYPE
Definition: raster.h:13
#define CELL_TYPE
Definition: raster.h:11
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