GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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  * */
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)
95  N_free_array_3d(data->phead);
96  if (data->phead_start)
98  if (data->status)
99  N_free_array_3d(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)
119  N_free_array_3d(data->river_bed);
120  if (data->drain_leak)
122  if (data->drain_bed)
123  N_free_array_3d(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  * */
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)
229  N_free_array_2d(data->river_bed);
230  if (data->drain_leak)
232  if (data->drain_bed)
233  N_free_array_2d(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  * */
263 N_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  * */
370 void
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 {
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  * */
464 N_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  * */
658 void
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 {
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 CELL_TYPE
Definition: raster.h:11
N_array_2d * status
Definition: N_gwflow.h:91
void Rast_set_null_value(void *, int, RASTER_MAP_TYPE)
To set one or more raster values to null.
Definition: null_val.c:98
double N
Definition: N_pde.h:276
N_array_2d * top
Definition: N_gwflow.h:88
Matrix entries for a mass balance 5/7/9 star system.
Definition: N_pde.h:272
N_array_2d * s
Definition: N_gwflow.h:75
int rows
Definition: N_pde.h:134
double E
Definition: N_pde.h:276
N_array_2d * hc_y
Definition: N_gwflow.h: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
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
int depths
Definition: N_pde.h:168
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
N_array_2d * river_bed
Definition: N_gwflow.h:81
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 G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
N_array_3d * river_head
Definition: N_gwflow.h:49
This data structure contains all data needed to compute the groundwater mass balance in three dimensi...
Definition: N_gwflow.h:35
#define N
Definition: e_intersect.c:917
int cols
Definition: N_pde.h:134
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_array_2d * phead_start
Definition: N_gwflow.h:70
N_array_3d * river_leak
Definition: N_gwflow.h:48
N_array_2d * phead
Definition: N_gwflow.h:69
N_array_2d * q
Definition: N_gwflow.h:73
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_array_3d * s
Definition: N_gwflow.h:44
double T
Definition: N_pde.h:278
#define NULL
Definition: ccmath.h:32
N_array_3d * nf
Definition: N_gwflow.h:45
int depths
Definition: N_pde.h:115
#define x
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
#define G_calloc(m, n)
Definition: defs/gis.h:113
N_array_3d * river_bed
Definition: N_gwflow.h:50
double S
Definition: N_pde.h:276
double dy
Definition: N_pde.h:110
N_array_3d * hc_z
Definition: N_gwflow.h:41
N_array_2d * r
Definition: N_gwflow.h:74
#define W
Definition: ogsf.h:140
double top
Extent coordinates (top) - 3D data.
Definition: gis.h:477
N_array_2d * r
Definition: N_gwflow.h:43
void G_message(const char *,...) __attribute__((format(printf
Geometric information about the structured grid.
Definition: N_pde.h:103
double dz
Definition: N_pde.h:111
N_array_3d * status
Definition: N_gwflow.h:56
#define DCELL_TYPE
Definition: raster.h:13
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d structure.
Definition: n_arrays.c:130
N_array_3d * drain_leak
Definition: N_gwflow.h:53
int cols
Definition: N_pde.h:168
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_3d * phead_start
Definition: N_gwflow.h:38
N_array_3d * hc_y
Definition: N_gwflow.h:40
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
#define N_CELL_INACTIVE
Definition: N_pde.h:31
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_array_2d * hc_x
Definition: N_gwflow.h:71
int depths
number of depths for 3D data
Definition: gis.h:446
N_array_2d * river_head
Definition: N_gwflow.h:80
N_array_2d * nf
Definition: N_gwflow.h:76
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
double N_calc_arith_mean(double a, double b)
Calculate the arithmetic mean of values a and b.
Definition: n_tools.c:33
N_array_2d * drain_bed
Definition: N_gwflow.h:85
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_array_3d * hc_x
Definition: N_gwflow.h:39
void G_percent(long, long, int)
Print percent complete messages.
Definition: percent.c:62
N_array_2d * drain_leak
Definition: N_gwflow.h:84
int cols
Number of columns for 2D data.
Definition: gis.h:442
This data structure contains all data needed to compute the groundwater mass balance in two dimension...
Definition: N_gwflow.h:67
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
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
N_array_3d * drain_bed
Definition: N_gwflow.h:54
void G_warning(const char *,...) __attribute__((format(printf
N_array_3d * phead
Definition: N_gwflow.h:37
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
#define _(str)
Definition: glocale.h:10
N_array_2d * bottom
Definition: N_gwflow.h:89
N_array_3d * q
Definition: N_gwflow.h:42
double B
Definition: N_pde.h:280
double N_calc_harmonic_mean(double a, double b)
Calculate the harmonical mean of values a and b.
Definition: n_tools.c:119
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
double W
Definition: N_pde.h:276
double dx
Definition: N_pde.h:109
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_array_3d(N_array_3d *data)
Release the memory of a N_array_3d.
Definition: n_arrays.c:780
N_array_2d * river_leak
Definition: N_gwflow.h:79
int rows
Number of rows for 2D data.
Definition: gis.h:438
int G_debug(int, const char *,...) __attribute__((format(printf
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
double r
Definition: r_raster.c:39
int rows
Definition: N_pde.h:168