GRASS Programmer's Manual  6.4.3(2013)-r
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Macros Pages
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 /* *************************************************************** */
38 N_gwflow_data3d *N_alloc_gwflow_data3d(int cols, int rows, int depths,
39  int river, int drain)
40 {
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 /* *************************************************************** */
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 /* *************************************************************** */
148 N_gwflow_data2d *N_alloc_gwflow_data2d(int cols, int rows, int river,
149  int drain)
150 {
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 /* *************************************************************** */
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 /* *************************************************************** */
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;
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  /*specific yield */
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  /*specific yield */
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 /* ****************** N_callback_gwflow_2d *********************** */
353 /* *************************************************************** */
370 N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
371  int row)
372 {
373  double T_e = 0, T_w = 0, T_n = 0, T_s = 0;
374  double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
375  double dx, dy, Az;
376  double hc_x, hc_y;
377  double z, top;
378  double hc_xw, hc_yn;
379  double z_xw, z_yn;
380  double hc_xe, hc_ys;
381  double z_xe, z_ys;
382  double hc, hc_start;
383  double Ss, r, nf, q;
384  double C, W, E, N, S, V;
386  N_data_star *mat_pos;
387  double river_vect = 0; /*entry in vector */
388  double river_mat = 0; /*entry in matrix */
389  double drain_vect = 0; /*entry in vector */
390  double drain_mat = 0; /*entry in matrix */
391 
392  /*cast the void pointer to the right data structure */
393  data = (N_gwflow_data2d *) gwdata;
394 
395  dx = geom->dx;
396  dy = geom->dy;
397  Az = N_get_geom_data_area_of_cell(geom, row);
398 
399  /*read the data from the arrays */
400  hc_start = N_get_array_2d_d_value(data->phead_start, col, row);
401  hc = N_get_array_2d_d_value(data->phead, col, row);
402  top = N_get_array_2d_d_value(data->top, col, row);
403 
404 
405  if (hc > top) { /*If the aquifer is confined */
406  z = N_get_array_2d_d_value(data->top, col,
407  row) -
408  N_get_array_2d_d_value(data->bottom, col, row);
409  z_xw =
410  N_get_array_2d_d_value(data->top, col - 1,
411  row) -
412  N_get_array_2d_d_value(data->bottom, col - 1, row);
413  z_xe =
414  N_get_array_2d_d_value(data->top, col + 1,
415  row) -
416  N_get_array_2d_d_value(data->bottom, col + 1, row);
417  z_yn =
418  N_get_array_2d_d_value(data->top, col,
419  row - 1) -
420  N_get_array_2d_d_value(data->bottom, col, row - 1);
421  z_ys =
422  N_get_array_2d_d_value(data->top, col,
423  row + 1) -
424  N_get_array_2d_d_value(data->bottom, col, row + 1);
425  }
426  else { /* the aquifer is unconfined */
427 
428  /* If the aquifer is unconfied use an explicite scheme to solve
429  * the nonlinear equation. We use the phead from the first iteration */
430  z = N_get_array_2d_d_value(data->phead, col, row) -
431  N_get_array_2d_d_value(data->bottom, col, row);
432  z_xw = N_get_array_2d_d_value(data->phead, col - 1, row) -
433  N_get_array_2d_d_value(data->bottom, col - 1, row);
434  z_xe = N_get_array_2d_d_value(data->phead, col + 1, row) -
435  N_get_array_2d_d_value(data->bottom, col + 1, row);
436  z_yn = N_get_array_2d_d_value(data->phead, col, row - 1) -
437  N_get_array_2d_d_value(data->bottom, col, row - 1);
438  z_ys = N_get_array_2d_d_value(data->phead, col, row + 1) -
439  N_get_array_2d_d_value(data->bottom, col, row + 1);
440  }
441 
442  /*geometrical mean of cell height */
443  if (z_w > 0 || z_w < 0 || z_w == 0)
444  z_w = N_calc_arith_mean(z_xw, z);
445  else
446  z_w = z;
447  if (z_e > 0 || z_e < 0 || z_e == 0)
448  z_e = N_calc_arith_mean(z_xe, z);
449  else
450  z_e = z;
451  if (z_n > 0 || z_n < 0 || z_n == 0)
452  z_n = N_calc_arith_mean(z_yn, z);
453  else
454  z_n = z;
455  if (z_s > 0 || z_s < 0 || z_s == 0)
456  z_s = N_calc_arith_mean(z_ys, z);
457  else
458  z_s = z;
459 
460  /* Inner sources */
461  q = N_get_array_2d_d_value(data->q, col, row);
462  nf = N_get_array_2d_d_value(data->nf, col, row);
463 
464  /* specific yield of current cell face */
465  Ss = N_get_array_2d_d_value(data->s, col, row) * Az;
466  /* recharge */
467  r = N_get_array_2d_d_value(data->r, col, row) * Az;
468 
469  /*get the surrounding permeabilities */
470  hc_x = N_get_array_2d_d_value(data->hc_x, col, row);
471  hc_y = N_get_array_2d_d_value(data->hc_y, col, row);
472  hc_xw = N_get_array_2d_d_value(data->hc_x, col - 1, row);
473  hc_xe = N_get_array_2d_d_value(data->hc_x, col + 1, row);
474  hc_yn = N_get_array_2d_d_value(data->hc_y, col, row - 1);
475  hc_ys = N_get_array_2d_d_value(data->hc_y, col, row + 1);
476 
477  /* calculate the transmissivities */
478  T_w = N_calc_harmonic_mean(hc_xw, hc_x) * z_w;
479  T_e = N_calc_harmonic_mean(hc_xe, hc_x) * z_e;
480  T_n = N_calc_harmonic_mean(hc_yn, hc_y) * z_n;
481  T_s = N_calc_harmonic_mean(hc_ys, hc_y) * z_s;
482 
483  /*compute the river leakage, this is an explicit method */
484  if (data->river_leak &&
485  (N_get_array_2d_d_value(data->river_leak, col, row) != 0)) {
486  if (hc > N_get_array_2d_d_value(data->river_bed, col, row)) {
487  river_vect = N_get_array_2d_d_value(data->river_head, col, row) *
488  N_get_array_2d_d_value(data->river_leak, col, row);
489  river_mat = N_get_array_2d_d_value(data->river_leak, col, row);
490  }
491  else if (hc < N_get_array_2d_d_value(data->river_bed, col, row)) {
492  river_vect = (N_get_array_2d_d_value(data->river_head, col, row) -
493  N_get_array_2d_d_value(data->river_bed, col, row))
494  * N_get_array_2d_d_value(data->river_leak, col, row);
495  river_mat = 0;
496  }
497  }
498 
499  /*compute the drainage, this is an explicit method */
500  if (data->drain_leak &&
501  (N_get_array_2d_d_value(data->drain_leak, col, row) != 0)) {
502  if (hc > N_get_array_2d_d_value(data->drain_bed, col, row)) {
503  drain_vect = N_get_array_2d_d_value(data->drain_bed, col, row) *
504  N_get_array_2d_d_value(data->drain_leak, col, row);
505  drain_mat = N_get_array_2d_d_value(data->drain_leak, col, row);
506  }
507  else if (hc <= N_get_array_2d_d_value(data->drain_bed, col, row)) {
508  drain_vect = 0;
509  drain_mat = 0;
510  }
511  }
512 
513  /*mass balance center cell to western cell */
514  W = -1 * T_w * dy / dx;
515  /*mass balance center cell to eastern cell */
516  E = -1 * T_e * dy / dx;
517  /*mass balance center cell to northern cell */
518  N = -1 * T_n * dx / dy;
519  /*mass balance center cell to southern cell */
520  S = -1 * T_s * dx / dy;
521 
522  /*the diagonal entry of the matrix */
523  C = -1 * (W + E + N + S - Ss / data->dt - river_mat * Az -
524  drain_mat * Az);
525 
526  /*the entry in the right side b of Ax = b */
527  V = (q + hc_start * Ss / data->dt) + r + river_vect * Az +
528  drain_vect * Az;
529 
530  G_debug(5, "N_callback_gwflow_2d: called [%i][%i]", row, col);
531 
532  /*create the 5 point star entries */
533  mat_pos = N_create_5star(C, W, E, N, S, V);
534 
535  return mat_pos;
536 }