GRASS GIS 7 Programmer's Manual  7.5.svn(2017)-r71817
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
n_solute_transport.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: solute transport in porous media
9 * part of the gpde library
10 *
11 * COPYRIGHT: (C) 2007 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 <math.h>
20 #include <grass/N_solute_transport.h>
21 
22 /* ************************************************************************* *
23  * ************************************************************************* *
24  * ************************************************************************* */
25 /*! \brief This is just a placeholder
26  *
27  * */
29  N_geom_data * geom, int col,
30  int row, int depth)
31 {
32  double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0, Df_t = 0, Df_b = 0;
33  double dx, dy, dz, Az;
34  double diff_x, diff_y, diff_z;
35  double diff_xw, diff_yn;
36  double diff_xe, diff_ys;
37  double diff_zt, diff_zb;
38  double cin = 0, cg, cg_start;
39  double R, nf, cs, q;
40  double C, W, E, N, S, T, B, V;
41  double vw = 0, ve = 0, vn = 0, vs = 0, vt = 0, vb = 0;
42  double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0, Ds_t = 0, Ds_b = 0;
43  double Dw = 0, De = 0, Dn = 0, Ds = 0, Dt = 0, Db = 0;
44  double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5, rt = 0.5, rb = 0.5;
45 
47  N_data_star *mat_pos;
48  N_gradient_3d grad;
49 
50  /*cast the void pointer to the right data structure */
51  data = (N_solute_transport_data3d *) solutedata;
52 
53  N_get_gradient_3d(data->grad, &grad, col, row, depth);
54 
55  dx = geom->dx;
56  dy = geom->dy;
57  dz = geom->dz;
58  Az = N_get_geom_data_area_of_cell(geom, row);
59 
60  /*read the data from the arrays */
61  cg_start = N_get_array_3d_d_value(data->c_start, col, row, depth);
62  cg = N_get_array_3d_d_value(data->c, col, row, depth);
63 
64  /*get the surrounding diffusion tensor entries */
65  diff_x = N_get_array_3d_d_value(data->diff_x, col, row, depth);
66  diff_y = N_get_array_3d_d_value(data->diff_y, col, row, depth);
67  diff_z = N_get_array_3d_d_value(data->diff_z, col, row, depth);
68  diff_xw = N_get_array_3d_d_value(data->diff_x, col - 1, row, depth);
69  diff_xe = N_get_array_3d_d_value(data->diff_x, col + 1, row, depth);
70  diff_yn = N_get_array_3d_d_value(data->diff_y, col, row - 1, depth);
71  diff_ys = N_get_array_3d_d_value(data->diff_y, col, row + 1, depth);
72  diff_zt = N_get_array_3d_d_value(data->diff_z, col, row, depth + 1);
73  diff_zb = N_get_array_3d_d_value(data->diff_z, col, row, depth - 1);
74 
75  /* calculate the diffusion on the cell borders using the harmonical mean */
76  Df_w = N_calc_harmonic_mean(diff_xw, diff_x);
77  Df_e = N_calc_harmonic_mean(diff_xe, diff_x);
78  Df_n = N_calc_harmonic_mean(diff_yn, diff_y);
79  Df_s = N_calc_harmonic_mean(diff_ys, diff_y);
80  Df_t = N_calc_harmonic_mean(diff_zt, diff_z);
81  Df_b = N_calc_harmonic_mean(diff_zb, diff_z);
82 
83  /* calculate the dispersion */
84  /*todo */
85 
86  /* calculate the velocity parts with full upwinding scheme */
87  vw = grad.WC;
88  ve = grad.EC;
89  vn = grad.NC;
90  vs = grad.SC;
91  vt = grad.TC;
92  vb = grad.BC;
93 
94  /* put the diffusion and dispersion together */
95  Dw = ((Df_w + Ds_w)) / dx;
96  De = ((Df_e + Ds_e)) / dx;
97  Dn = ((Df_n + Ds_n)) / dy;
98  Ds = ((Df_s + Ds_s)) / dy;
99  Dt = ((Df_t + Ds_t)) / dz;
100  Db = ((Df_b + Ds_b)) / dz;
101 
102  rw = N_exp_upwinding(-1 * vw, dx, Dw);
103  re = N_exp_upwinding(ve, dx, De);
104  rs = N_exp_upwinding(-1 * vs, dy, Ds);
105  rn = N_exp_upwinding(vn, dy, Dn);
106  rb = N_exp_upwinding(-1 * vb, dz, Dn);
107  rt = N_exp_upwinding(vt, dz, Dn);
108 
109  /*mass balance center cell to western cell */
110  W = -1 * (Dw) * dy * dz - vw * (1 - rw) * dy * dz;
111  /*mass balance center cell to eastern cell */
112  E = -1 * (De) * dy * dz + ve * (1 - re) * dy * dz;
113  /*mass balance center cell to southern cell */
114  S = -1 * (Ds) * dx * dz - vs * (1 - rs) * dx * dz;
115  /*mass balance center cell to northern cell */
116  N = -1 * (Dn) * dx * dz + vn * (1 - rn) * dx * dz;
117  /*mass balance center cell to bottom cell */
118  B = -1 * (Db) * Az - vb * (1 - rb) * Az;
119  /*mass balance center cell to top cell */
120  T = -1 * (Dt) * Az + vt * (1 - rt) * Az;
121 
122  /* Retardation */
123  R = N_get_array_3d_d_value(data->R, col, row, depth);
124  /* Inner sources */
125  cs = N_get_array_3d_d_value(data->cs, col, row, depth);
126  /* effective porosity */
127  nf = N_get_array_3d_d_value(data->nf, col, row, depth);
128  /* groundwater sources and sinks */
129  q = N_get_array_3d_d_value(data->q, col, row, depth);
130  /* concentration of influent water */
131  cin = N_get_array_3d_d_value(data->cin, col, row, depth);
132 
133  /*the diagonal entry of the matrix */
134  C = ((Dw - vw) * dy * dz +
135  (De + ve) * dy * dz +
136  (Ds - vs) * dx * dz +
137  (Dn + vn) * dx * dz +
138  (Db - vb) * Az + (Dt + vt) * Az + Az * dz * R / data->dt - q / nf);
139 
140  /*the entry in the right side b of Ax = b */
141  V = (cs + cg_start * Az * dz * R / data->dt - q / nf * cin);
142 
143  /*
144  * printf("nf %g\n", nf);
145  * printf("q %g\n", q);
146  * printf("cs %g\n", cs);
147  * printf("cin %g\n", cin);
148  * printf("cg %g\n", cg);
149  * printf("cg_start %g\n", cg_start);
150  * printf("Az %g\n", Az);
151  * printf("z %g\n", z);
152  * printf("R %g\n", R);
153  * printf("dt %g\n", data->dt);
154  */
155  G_debug(6, "N_callback_solute_transport_3d: called [%i][%i][%i]", row,
156  col, depth);
157 
158  /*create the 7 point star entries */
159  mat_pos = N_create_7star(C, W, E, N, S, T, B, V);
160 
161  return mat_pos;
162 }
163 
164 /* ************************************************************************* *
165  * ************************************************************************* *
166  * ************************************************************************* */
167 /*!
168  * \brief This callback function creates the mass balance of a 5 point star
169  *
170  * The mass balance is based on the common solute transport equation:
171  *
172  * \f[\frac{\partial c_g}{\partial t} R = \nabla \cdot ({\bf D} \nabla c_g - {\bf u} c_g) + \sigma + \frac{q}{n_f}(c_g - c_in) \f]
173  *
174  * This equation is discretizised with the finite volume method in two dimensions.
175  *
176  *
177  * \param solutedata * N_solute_transport_data2d - a void pointer to the data structure
178  * \param geom N_geom_data *
179  * \param col int
180  * \param row int
181  * \return N_data_star * - a five point data star
182  *
183  * */
185  N_geom_data * geom, int col,
186  int row)
187 {
188  double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0;
189  double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
190  double dx, dy, Az;
191  double diff_x, diff_y;
192  double disp_x, disp_y;
193  double z;
194  double diff_xw, diff_yn;
195  double disp_xw, disp_yn;
196  double z_xw, z_yn;
197  double diff_xe, diff_ys;
198  double disp_xe, disp_ys;
199  double z_xe, z_ys;
200  double cin = 0, cg, cg_start;
201  double R, nf, cs, q;
202  double C, W, E, N, S, V, NE, NW, SW, SE;
203  double vw = 0, ve = 0, vn = 0, vs = 0;
204  double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0;
205  double Dw = 0, De = 0, Dn = 0, Ds = 0;
206  double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5;
207 
209  N_data_star *mat_pos;
210  N_gradient_2d grad;
211 
212  /*cast the void pointer to the right data structure */
213  data = (N_solute_transport_data2d *) solutedata;
214 
215  N_get_gradient_2d(data->grad, &grad, col, row);
216 
217  dx = geom->dx;
218  dy = geom->dy;
219  Az = N_get_geom_data_area_of_cell(geom, row);
220 
221  /*read the data from the arrays */
222  cg_start = N_get_array_2d_d_value(data->c_start, col, row);
223  cg = N_get_array_2d_d_value(data->c, col, row);
224 
225  /* calculate the cell height */
226  z = N_get_array_2d_d_value(data->top, col,
227  row) -
228  N_get_array_2d_d_value(data->bottom, col, row);
229  z_xw =
230  N_get_array_2d_d_value(data->top, col - 1,
231  row) -
232  N_get_array_2d_d_value(data->bottom, col - 1, row);
233  z_xe =
234  N_get_array_2d_d_value(data->top, col + 1,
235  row) -
236  N_get_array_2d_d_value(data->bottom, col + 1, row);
237  z_yn =
238  N_get_array_2d_d_value(data->top, col,
239  row - 1) -
240  N_get_array_2d_d_value(data->bottom, col, row - 1);
241  z_ys =
242  N_get_array_2d_d_value(data->top, col,
243  row + 1) -
244  N_get_array_2d_d_value(data->bottom, col, row + 1);
245 
246  /*geometrical mean of cell height */
247  z_w = N_calc_geom_mean(z_xw, z);
248  z_e = N_calc_geom_mean(z_xe, z);
249  z_n = N_calc_geom_mean(z_yn, z);
250  z_s = N_calc_geom_mean(z_ys, z);
251 
252  /*get the surrounding diffusion tensor entries */
253  diff_x = N_get_array_2d_d_value(data->diff_x, col, row);
254  diff_y = N_get_array_2d_d_value(data->diff_y, col, row);
255  diff_xw = N_get_array_2d_d_value(data->diff_x, col - 1, row);
256  diff_xe = N_get_array_2d_d_value(data->diff_x, col + 1, row);
257  diff_yn = N_get_array_2d_d_value(data->diff_y, col, row - 1);
258  diff_ys = N_get_array_2d_d_value(data->diff_y, col, row + 1);
259 
260  /* calculate the diffusion at the cell borders using the harmonical mean */
261  Df_w = N_calc_harmonic_mean(diff_xw, diff_x);
262  Df_e = N_calc_harmonic_mean(diff_xe, diff_x);
263  Df_n = N_calc_harmonic_mean(diff_yn, diff_y);
264  Df_s = N_calc_harmonic_mean(diff_ys, diff_y);
265 
266  /* calculate the dispersion */
267  /*get the surrounding dispersion tensor entries */
268  disp_x = N_get_array_2d_d_value(data->disp_xx, col, row);
269  disp_y = N_get_array_2d_d_value(data->disp_yy, col, row);
270  if (N_get_array_2d_d_value(data->status, col - 1, row) ==
272  disp_xw = disp_x;
273  }
274  else {
275  disp_xw = N_get_array_2d_d_value(data->disp_xx, col - 1, row);
276  }
277  if (N_get_array_2d_d_value(data->status, col + 1, row) ==
279  disp_xe = disp_x;
280  }
281  else {
282  disp_xe = N_get_array_2d_d_value(data->disp_xx, col + 1, row);
283  }
284  if (N_get_array_2d_d_value(data->status, col, row - 1) ==
286  disp_yn = disp_y;
287  }
288  else {
289  disp_yn = N_get_array_2d_d_value(data->disp_yy, col, row - 1);
290  }
291  if (N_get_array_2d_d_value(data->status, col, row + 1) ==
293  disp_ys = disp_y;
294  }
295  else {
296  disp_ys = N_get_array_2d_d_value(data->disp_yy, col, row + 1);
297  }
298 
299  /* calculate the dispersion at the cell borders using the harmonical mean */
300  Ds_w = N_calc_harmonic_mean(disp_xw, disp_x);
301  Ds_e = N_calc_harmonic_mean(disp_xe, disp_x);
302  Ds_n = N_calc_harmonic_mean(disp_yn, disp_y);
303  Ds_s = N_calc_harmonic_mean(disp_ys, disp_y);
304 
305  /* put the diffusion and dispersion together */
306  Dw = ((Df_w + Ds_w)) / dx;
307  De = ((Df_e + Ds_e)) / dx;
308  Ds = ((Df_s + Ds_s)) / dy;
309  Dn = ((Df_n + Ds_n)) / dy;
310 
311  vw = -1.0 * grad.WC;
312  ve = grad.EC;
313  vs = -1.0 * grad.SC;
314  vn = grad.NC;
315 
316  if (data->stab == N_UPWIND_FULL) {
317  rw = N_full_upwinding(vw, dx, Dw);
318  re = N_full_upwinding(ve, dx, De);
319  rs = N_full_upwinding(vs, dy, Ds);
320  rn = N_full_upwinding(vn, dy, Dn);
321  }
322  else if (data->stab == N_UPWIND_EXP) {
323  rw = N_exp_upwinding(vw, dx, Dw);
324  re = N_exp_upwinding(ve, dx, De);
325  rs = N_exp_upwinding(vs, dy, Ds);
326  rn = N_exp_upwinding(vn, dy, Dn);
327  }
328 
329  /*mass balance center cell to western cell */
330  W = -1 * (Dw) * dy * z_w + vw * (1 - rw) * dy * z_w;
331  /*mass balance center cell to eastern cell */
332  E = -1 * (De) * dy * z_e + ve * (1 - re) * dy * z_e;
333  /*mass balance center cell to southern cell */
334  S = -1 * (Ds) * dx * z_s + vs * (1 - rs) * dx * z_s;
335  /*mass balance center cell to northern cell */
336  N = -1 * (Dn) * dx * z_n + vn * (1 - rn) * dx * z_n;
337 
338  NW = 0.0;
339  SW = 0.0;
340  NE = 0.0;
341  SE = 0.0;
342 
343  /* Retardation */
344  R = N_get_array_2d_d_value(data->R, col, row);
345  /* Inner sources */
346  cs = N_get_array_2d_d_value(data->cs, col, row);
347  /* effective porosity */
348  nf = N_get_array_2d_d_value(data->nf, col, row);
349  /* groundwater sources and sinks */
350  q = N_get_array_2d_d_value(data->q, col, row);
351  /* concentration of influent water */
352  cin = N_get_array_2d_d_value(data->cin, col, row);
353 
354  /*the diagonal entry of the matrix */
355  C = (Dw + vw * rw) * dy * z_w +
356  (De + ve * re) * dy * z_e +
357  (Ds + vs * rs) * dx * z_s +
358  (Dn + vn * rn) * dx * z_n + Az * z * R / data->dt - q / nf;
359 
360  /*the entry in the right side b of Ax = b */
361  V = (cs + cg_start * Az * z * R / data->dt + q / nf * cin);
362 
363  /*
364  fprintf(stderr, "nf %g\n", nf);
365  fprintf(stderr, "q %g\n", q);
366  fprintf(stderr, "cs %g\n", cs);
367  fprintf(stderr, "cin %g\n", cin);
368  fprintf(stderr, "cg %g\n", cg);
369  fprintf(stderr, "cg_start %g\n", cg_start);
370  fprintf(stderr, "Az %g\n", Az);
371  fprintf(stderr, "z %g\n", z);
372  fprintf(stderr, "R %g\n", R);
373  fprintf(stderr, "dt %g\n", data->dt);
374  */
375 
376  G_debug(6, "N_callback_solute_transport_2d: called [%i][%i]", row, col);
377 
378  /*create the 9 point star entries */
379  mat_pos = N_create_9star(C, W, E, N, S, NW, SW, NE, SE, V);
380 
381  return mat_pos;
382 }
383 
384 /* ************************************************************************* *
385  * ************************************************************************* *
386  * ************************************************************************* */
387 /*!
388  * \brief Alllocate memory for the solute transport data structure in three dimensions
389  *
390  * The solute transport data structure will be allocated including
391  * all appendant 3d arrays. The offset for the 3d arrays is one
392  * to establish homogeneous Neumann boundary conditions at the calculation area border.
393  * This data structure is used to create a linear equation system based on the computation of
394  * solute transport in porous media with the finite volume method.
395  *
396  * \param cols int
397  * \param rows int
398  * \param depths int
399  * \return N_solute_transport_data3d *
400  * */
401 
403  int depths)
404 {
406 
407  data =
408  (N_solute_transport_data3d *) G_calloc(1,
409  sizeof
411 
412  data->c = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
413  data->c_start = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
414  data->status = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
415  data->diff_x = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
416  data->diff_y = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
417  data->diff_z = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
418  data->q = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
419  data->cs = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
420  data->R = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
421  data->nf = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
422  data->cin = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
423 
424  /*Allocate the dispersivity tensor */
425  data->disp_xx = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
426  data->disp_yy = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
427  data->disp_zz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
428  data->disp_xy = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
429  data->disp_xz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
430  data->disp_yz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
431 
432 
433  data->grad = N_alloc_gradient_field_3d(cols, rows, depths);
434  data->stab = N_UPWIND_EXP;
435 
436  return data;
437 }
438 
439 /* ************************************************************************* *
440  * ************************************************************************* *
441  * ************************************************************************* */
442 /*!
443  * \brief Alllocate memory for the solute transport data structure in two dimensions
444  *
445  * The solute transport data structure will be allocated including
446  * all appendant 2d arrays. The offset for the 2d arrays is one
447  * to establish homogeneous Neumann boundary conditions at the calculation area border.
448  * This data structure is used to create a linear equation system based on the computation of
449  * solute transport in porous media with the finite volume method.
450  *
451  * \param cols int
452  * \param rows int
453  * \return N_solute_transport_data2d *
454  * */
455 
456 
458 {
460 
461  data =
462  (N_solute_transport_data2d *) G_calloc(1,
463  sizeof
465 
466  data->c = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
467  data->c_start = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
468  data->status = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
469  data->diff_x = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
470  data->diff_y = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
471  data->q = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
472  data->cs = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
473  data->R = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
474  data->nf = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
475  data->cin = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
476  data->top = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
477  data->bottom = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
478 
479  /*Allocate the dispersivity tensor */
480  data->disp_xx = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
481  data->disp_yy = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
482  data->disp_xy = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
483 
484  data->grad = N_alloc_gradient_field_2d(cols, rows);
485  data->stab = N_UPWIND_EXP;
486 
487  return data;
488 }
489 
490 /* ************************************************************************* *
491  * ************************************************************************* *
492  * ************************************************************************* */
493 /*!
494  * \brief Release the memory of the solute transport data structure in three dimensions
495  *
496  * \param data N_solute_transport_data2d *
497  * \return void *
498  * */
500 {
501  N_free_array_3d(data->c);
502  N_free_array_3d(data->c_start);
503  N_free_array_3d(data->status);
504  N_free_array_3d(data->diff_x);
505  N_free_array_3d(data->diff_y);
506  N_free_array_3d(data->diff_z);
507  N_free_array_3d(data->q);
508  N_free_array_3d(data->cs);
509  N_free_array_3d(data->R);
510  N_free_array_3d(data->nf);
511  N_free_array_3d(data->cin);
512 
513  N_free_array_3d(data->disp_xx);
514  N_free_array_3d(data->disp_yy);
515  N_free_array_3d(data->disp_zz);
516  N_free_array_3d(data->disp_xy);
517  N_free_array_3d(data->disp_xz);
518  N_free_array_3d(data->disp_yz);
519 
520  G_free(data);
521 
522  data = NULL;
523 
524  return;
525 }
526 
527 /* ************************************************************************* *
528  * ************************************************************************* *
529  * ************************************************************************* */
530 /*!
531  * \brief Release the memory of the solute transport data structure in two dimensions
532  *
533  * \param data N_solute_transport_data2d *
534  * \return void *
535  * */
537 {
538  N_free_array_2d(data->c);
539  N_free_array_2d(data->c_start);
540  N_free_array_2d(data->status);
541  N_free_array_2d(data->diff_x);
542  N_free_array_2d(data->diff_y);
543  N_free_array_2d(data->q);
544  N_free_array_2d(data->cs);
545  N_free_array_2d(data->R);
546  N_free_array_2d(data->nf);
547  N_free_array_2d(data->cin);
548  N_free_array_2d(data->top);
549  N_free_array_2d(data->bottom);
550 
551  N_free_array_2d(data->disp_xx);
552  N_free_array_2d(data->disp_yy);
553  N_free_array_2d(data->disp_xy);
554 
555  G_free(data);
556 
557  data = NULL;
558 
559  return;
560 }
561 
562 /*!
563  * \brief Compute the transmission boundary condition in 2d
564  *
565  * This function calculates the transmission boundary condition
566  * for each cell with status N_CELL_TRANSMISSION. The surrounding
567  * gradient field is used to verfiy the flow direction. If a flow
568  * goes into a cell, the concentration (data->c) from the neighbour cell is
569  * added to the transmission cell. If the flow from several neighbour
570  * cells goes into the cell, the concentration mean is calculated.
571  *
572  * The new concentrations are written into the data->c_start array,
573  * so they can be handled by the matrix assembling function.
574  *
575  * \param data N_solute_transport_data2d *
576  * \return void *
577  * */
579 {
580  int i, j, count = 1;
581  int cols, rows;
582  double c;
583  N_gradient_2d grad;
584 
585  cols = data->grad->cols;
586  rows = data->grad->rows;
587 
588  G_debug(2,
589  "N_calc_solute_transport_transmission_2d: calculating transmission boundary");
590 
591  for (j = 0; j < rows; j++) {
592  for (i = 0; i < cols; i++) {
593  if (N_get_array_2d_d_value(data->status, i, j) ==
595  count = 0;
596  /*get the gradient neighbours */
597  N_get_gradient_2d(data->grad, &grad, i, j);
598  c = 0;
599  /*
600  c = N_get_array_2d_d_value(data->c_start, i, j);
601  if(c > 0)
602  count++;
603  */
604 
605  if (grad.WC > 0 &&
606  !N_is_array_2d_value_null(data->c, i - 1, j)) {
607  c += N_get_array_2d_d_value(data->c, i - 1, j);
608  count++;
609  }
610  if (grad.EC < 0 &&
611  !N_is_array_2d_value_null(data->c, i + 1, j)) {
612  c += N_get_array_2d_d_value(data->c, i + 1, j);
613  count++;
614  }
615  if (grad.NC < 0 &&
616  !N_is_array_2d_value_null(data->c, i, j - 1)) {
617  c += N_get_array_2d_d_value(data->c, i, j - 1);
618  count++;
619  }
620  if (grad.SC > 0 &&
621  !N_is_array_2d_value_null(data->c, i, j + 1)) {
622  c += N_get_array_2d_d_value(data->c, i, j + 1);
623  count++;
624  }
625  if (count != 0)
626  c = c / (double)count;
627  /*make sure it is not NAN */
628  if (c > 0 || c == 0 || c < 0)
629  N_put_array_2d_d_value(data->c_start, i, j, c);
630  }
631  }
632  }
633 
634  return;
635 }
636 
637 /*!
638  * \brief Compute the dispersivity tensor based on the solute transport data in 2d
639  *
640  * The dispersivity tensor is stored in the data structure.
641  * To compute the dispersivity tensor, the dispersivity lentghs and the gradient field
642  * must be present.
643  *
644  * This is just a simple tensor computation which should be extended.
645  *
646  * \todo Change the tensor calculation to a mor realistic algorithm
647  *
648  * \param data N_solute_transport_data2d *
649  * \return void *
650  * */
652 {
653  int i, j;
654  int cols, rows;
655  double vx, vy, vv;
656  double disp_xx, disp_yy, disp_xy;
657  N_gradient_2d grad;
658 
659  cols = data->grad->cols;
660  rows = data->grad->rows;
661 
662  G_debug(2,
663  "N_calc_solute_transport_disptensor_2d: calculating the dispersivity tensor");
664 
665  for (j = 0; j < rows; j++) {
666  for (i = 0; i < cols; i++) {
667 
668  disp_xx = 0;
669  disp_yy = 0;
670  disp_xy = 0;
671 
672  /*get the gradient neighbours */
673  N_get_gradient_2d(data->grad, &grad, i, j);
674  vx = (grad.WC + grad.EC) / 2;
675  vy = (grad.NC + grad.SC) / 2;
676  vv = sqrt(vx * vx + vy * vy);
677 
678  if (vv != 0) {
679  disp_xx = data->al * vx * vx / vv + data->at * vy * vy / vv;
680  disp_yy = data->at * vx * vx / vv + data->al * vy * vy / vv;
681  disp_xy = (data->al - data->at) * vx * vy / vv;
682  }
683 
684  G_debug(5,
685  "N_calc_solute_transport_disptensor_2d: [%i][%i] disp_xx %g disp_yy %g disp_xy %g",
686  i, j, disp_xx, disp_yy, disp_xy);
687  N_put_array_2d_d_value(data->disp_xx, i, j, disp_xx);
688  N_put_array_2d_d_value(data->disp_yy, i, j, disp_yy);
689  N_put_array_2d_d_value(data->disp_xy, i, j, disp_xy);
690  }
691  }
692 
693  return;
694 }
695 
696 /*!
697  * \brief Compute the dispersivity tensor based on the solute transport data in 3d
698  *
699  * The dispersivity tensor is stored in the data structure.
700  * To compute the dispersivity tensor, the dispersivity lentghs and the gradient field
701  * must be present.
702  *
703  * This is just a simple tensor computation which should be extended.
704  *
705  * \todo Change the tensor calculation to a mor realistic algorithm
706  *
707  * \param data N_solute_transport_data3d *
708  * \return void *
709  * */
711 {
712  int i, j, k;
713  int cols, rows, depths;
714  double vx, vy, vz, vv;
715  double disp_xx, disp_yy, disp_zz, disp_xy, disp_xz, disp_yz;
716  N_gradient_3d grad;
717 
718  cols = data->grad->cols;
719  rows = data->grad->rows;
720  depths = data->grad->depths;
721 
722  G_debug(2,
723  "N_calc_solute_transport_disptensor_3d: calculating the dispersivity tensor");
724 
725  for (k = 0; k < depths; k++) {
726  for (j = 0; j < rows; j++) {
727  for (i = 0; i < cols; i++) {
728  disp_xx = 0;
729  disp_yy = 0;
730  disp_zz = 0;
731  disp_xy = 0;
732  disp_xz = 0;
733  disp_yz = 0;
734 
735  /*get the gradient neighbours */
736  N_get_gradient_3d(data->grad, &grad, i, j, k);
737  vx = (grad.WC + grad.EC) / 2;
738  vy = (grad.NC + grad.SC) / 2;
739  vz = (grad.BC + grad.TC) / 2;
740  vv = sqrt(vx * vx + vy * vy + vz * vz);
741 
742  if (vv != 0) {
743  disp_xx =
744  data->al * vx * vx / vv + data->at * vy * vy / vv +
745  data->at * vz * vz / vv;
746  disp_yy =
747  data->at * vx * vx / vv + data->al * vy * vy / vv +
748  data->at * vz * vz / vv;
749  disp_zz =
750  data->at * vx * vx / vv + data->at * vy * vy / vv +
751  data->al * vz * vz / vv;
752  disp_xy = (data->al - data->at) * vx * vy / vv;
753  disp_xz = (data->al - data->at) * vx * vz / vv;
754  disp_yz = (data->al - data->at) * vy * vz / vv;
755  }
756 
757  G_debug(5,
758  "N_calc_solute_transport_disptensor_3d: [%i][%i][%i] disp_xx %g disp_yy %g disp_zz %g disp_xy %g disp_xz %g disp_yz %g ",
759  i, j, k, disp_xx, disp_yy, disp_zz, disp_xy, disp_xz,
760  disp_yz);
761  N_put_array_3d_d_value(data->disp_xx, i, j, k, disp_xx);
762  N_put_array_3d_d_value(data->disp_yy, i, j, k, disp_yy);
763  N_put_array_3d_d_value(data->disp_zz, i, j, k, disp_zz);
764  N_put_array_3d_d_value(data->disp_xy, i, j, k, disp_xy);
765  N_put_array_3d_d_value(data->disp_xz, i, j, k, disp_xz);
766  N_put_array_3d_d_value(data->disp_yz, i, j, k, disp_yz);
767  }
768  }
769  }
770 
771  return;
772 }
N_gradient_2d * N_get_gradient_2d(N_gradient_field_2d *field, N_gradient_2d *gradient, int col, int row)
Return a N_gradient_2d structure calculated from the input gradient field at position [row][col]...
Definition: n_gradient.c:115
double EC
Definition: N_pde.h:419
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:149
N_data_star * N_create_9star(double C, double W, double E, double N, double S, double NW, double SW, double NE, double SE, double V)
allocate and initialize a 9 point star data structure
Matrix entries for a mass balance 5/7/9 star system.
Definition: N_pde.h:272
#define N_UPWIND_FULL
Definition: N_pde.h:53
double N_calc_geom_mean(double a, double b)
Calculate the geometrical mean of values a and b.
Definition: n_tools.c:76
#define N_CELL_TRANSMISSION
Definition: N_pde.h:34
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
Gradient between the cells in X and Y direction.
Definition: N_pde.h:408
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_gradient_field_3d * N_alloc_gradient_field_3d(int cols, int rows, int depths)
Allocate a N_gradient_field_3d.
Definition: n_gradient.c:1018
int count
#define NE
Definition: dataquad.h:30
#define N
Definition: e_intersect.c:923
void N_calc_solute_transport_transmission_2d(N_solute_transport_data2d *data)
Compute the transmission boundary condition in 2d.
double WC
Definition: N_pde.h:419
#define NW
Definition: dataquad.h:29
#define NULL
Definition: ccmath.h:32
double SC
Definition: N_pde.h:419
double dy
Definition: N_pde.h:110
#define W
Definition: ogsf.h:140
#define SW
Definition: dataquad.h:31
double NC
Definition: N_pde.h:419
N_gradient_field_3d * grad
void N_calc_solute_transport_disptensor_3d(N_solute_transport_data3d *data)
Compute the dispersivity tensor based on the solute transport data in 3d.
Geometric information about the structured grid.
Definition: N_pde.h:103
double dz
Definition: N_pde.h:111
#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
double TC
Definition: N_pde.h:419
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
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
double N_exp_upwinding(double sprod, double distance, double D)
exponential upwinding stabilization algorithm
Definition: n_upwind.c:63
double BC
Definition: N_pde.h:419
double EC
Definition: N_pde.h:411
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
int N_is_array_2d_value_null(N_array_2d *data, int col, int row)
Returns 1 if the value of N_array_2d struct at position col, row is of type null, otherwise 0...
Definition: n_arrays.c:231
Gradient between the cells in X, Y and Z direction.
Definition: N_pde.h:416
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
#define SE
Definition: dataquad.h:32
void N_calc_solute_transport_disptensor_2d(N_solute_transport_data2d *data)
Compute the dispersivity tensor based on the solute transport data in 2d.
void N_free_solute_transport_data2d(N_solute_transport_data2d *data)
Release the memory of the solute transport data structure in two dimensions.
void N_free_solute_transport_data3d(N_solute_transport_data3d *data)
Release the memory of the solute transport data structure in three dimensions.
N_gradient_field_2d * grad
#define N_UPWIND_EXP
Definition: N_pde.h:54
N_gradient_field_2d * N_alloc_gradient_field_2d(int cols, int rows)
Allocate a N_gradient_field_2d.
Definition: n_gradient.c:920
N_data_star * N_callback_solute_transport_2d(void *solutedata, N_geom_data *geom, int col, int row)
This callback function creates the mass balance of a 5 point star.
void N_put_array_2d_d_value(N_array_2d *data, int col, int row, DCELL value)
Writes a DCELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:584
double WC
Definition: N_pde.h:411
double N_calc_harmonic_mean(double a, double b)
Calculate the harmonical mean of values a and b.
Definition: n_tools.c:119
N_data_star * N_callback_solute_transport_3d(void *solutedata, N_geom_data *geom, int col, int row, int depth)
This is just a placeholder.
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 dx
Definition: N_pde.h:109
void N_free_array_3d(N_array_3d *data)
Release the memory of a N_array_3d.
Definition: n_arrays.c:780
double SC
Definition: N_pde.h:411
double NC
Definition: N_pde.h:411
double N_full_upwinding(double sprod, double distance, double D)
full upwinding stabilization algorithm
Definition: n_upwind.c:33
N_gradient_3d * N_get_gradient_3d(N_gradient_field_3d *field, N_gradient_3d *gradient, int col, int row, int depth)
Return a N_gradient_3d structure calculated from the input gradient field at position [depth][row][co...
Definition: n_gradient.c:248
N_solute_transport_data3d * N_alloc_solute_transport_data3d(int cols, int rows, int depths)
Alllocate memory for the solute transport data structure in three dimensions.
N_solute_transport_data2d * N_alloc_solute_transport_data2d(int cols, int rows)
Alllocate memory for the solute transport data structure in two dimensions.