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