GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 "grass/N_solute_transport.h"
20 #include <math.h>
21 
22 /* ************************************************************************* *
23  * ************************************************************************* *
24  * ************************************************************************* */
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 
34  double dx, dy, dz, Az;
35 
36  double diff_x, diff_y, diff_z;
37 
38  double diff_xw, diff_yn;
39 
40  double diff_xe, diff_ys;
41 
42  double diff_zt, diff_zb;
43 
44  double cin = 0, cg, cg_start;
45 
46  double R, nf, cs, q;
47 
48  double C, W, E, N, S, T, B, V;
49 
50  double vw = 0, ve = 0, vn = 0, vs = 0, vt = 0, vb = 0;
51 
52  double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0, Ds_t = 0, Ds_b = 0;
53 
54  double Dw = 0, De = 0, Dn = 0, Ds = 0, Dt = 0, Db = 0;
55 
56  double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5, rt = 0.5, rb = 0.5;
57 
59 
60  N_data_star *mat_pos;
61 
62  N_gradient_3d grad;
63 
64  /*cast the void pointer to the right data structure */
65  data = (N_solute_transport_data3d *) solutedata;
66 
67  N_get_gradient_3d(data->grad, &grad, col, row, depth);
68 
69  dx = geom->dx;
70  dy = geom->dy;
71  dz = geom->dz;
72  Az = N_get_geom_data_area_of_cell(geom, row);
73 
74  /*read the data from the arrays */
75  cg_start = N_get_array_3d_d_value(data->c_start, col, row, depth);
76  cg = N_get_array_3d_d_value(data->c, col, row, depth);
77 
78  /*get the surrounding diffusion tensor entries */
79  diff_x = N_get_array_3d_d_value(data->diff_x, col, row, depth);
80  diff_y = N_get_array_3d_d_value(data->diff_y, col, row, depth);
81  diff_z = N_get_array_3d_d_value(data->diff_z, col, row, depth);
82  diff_xw = N_get_array_3d_d_value(data->diff_x, col - 1, row, depth);
83  diff_xe = N_get_array_3d_d_value(data->diff_x, col + 1, row, depth);
84  diff_yn = N_get_array_3d_d_value(data->diff_y, col, row - 1, depth);
85  diff_ys = N_get_array_3d_d_value(data->diff_y, col, row + 1, depth);
86  diff_zt = N_get_array_3d_d_value(data->diff_z, col, row, depth + 1);
87  diff_zb = N_get_array_3d_d_value(data->diff_z, col, row, depth - 1);
88 
89  /* calculate the diffusion on the cell borders using the harmonical mean */
90  Df_w = N_calc_harmonic_mean(diff_xw, diff_x);
91  Df_e = N_calc_harmonic_mean(diff_xe, diff_x);
92  Df_n = N_calc_harmonic_mean(diff_yn, diff_y);
93  Df_s = N_calc_harmonic_mean(diff_ys, diff_y);
94  Df_t = N_calc_harmonic_mean(diff_zt, diff_z);
95  Df_b = N_calc_harmonic_mean(diff_zb, diff_z);
96 
97  /* calculate the dispersion */
98  /*todo */
99 
100  /* calculate the velocity parts with full upwinding scheme */
101  vw = grad.WC;
102  ve = grad.EC;
103  vn = grad.NC;
104  vs = grad.SC;
105  vt = grad.TC;
106  vb = grad.BC;
107 
108  /* put the diffusion and dispersion together */
109  Dw = ((Df_w + Ds_w)) / dx;
110  De = ((Df_e + Ds_e)) / dx;
111  Dn = ((Df_n + Ds_n)) / dy;
112  Ds = ((Df_s + Ds_s)) / dy;
113  Dt = ((Df_t + Ds_t)) / dz;
114  Db = ((Df_b + Ds_b)) / dz;
115 
116  rw = N_exp_upwinding(-1 * vw, dx, Dw);
117  re = N_exp_upwinding(ve, dx, De);
118  rs = N_exp_upwinding(-1 * vs, dy, Ds);
119  rn = N_exp_upwinding(vn, dy, Dn);
120  rb = N_exp_upwinding(-1 * vb, dz, Dn);
121  rt = N_exp_upwinding(vt, dz, Dn);
122 
123  /*mass balance center cell to western cell */
124  W = -1 * (Dw) * dy * dz - vw * (1 - rw) * dy * dz;
125  /*mass balance center cell to eastern cell */
126  E = -1 * (De) * dy * dz + ve * (1 - re) * dy * dz;
127  /*mass balance center cell to southern cell */
128  S = -1 * (Ds) * dx * dz - vs * (1 - rs) * dx * dz;
129  /*mass balance center cell to northern cell */
130  N = -1 * (Dn) * dx * dz + vn * (1 - rn) * dx * dz;
131  /*mass balance center cell to bottom cell */
132  B = -1 * (Db) * Az - vb * (1 - rb) * Az;
133  /*mass balance center cell to top cell */
134  T = -1 * (Dt) * Az + vt * (1 - rt) * Az;
135 
136  /* Retardation */
137  R = N_get_array_3d_d_value(data->R, col, row, depth);
138  /* Inner sources */
139  cs = N_get_array_3d_d_value(data->cs, col, row, depth);
140  /* effective porosity */
141  nf = N_get_array_3d_d_value(data->nf, col, row, depth);
142  /* groundwater sources and sinks */
143  q = N_get_array_3d_d_value(data->q, col, row, depth);
144  /* concentration of influent water */
145  cin = N_get_array_3d_d_value(data->cin, col, row, depth);
146 
147  /*the diagonal entry of the matrix */
148  C = ((Dw - vw) * dy * dz +
149  (De + ve) * dy * dz +
150  (Ds - vs) * dx * dz +
151  (Dn + vn) * dx * dz +
152  (Db - vb) * Az + (Dt + vt) * Az + Az * dz * R / data->dt - q / nf);
153 
154  /*the entry in the right side b of Ax = b */
155  V = (cs + cg_start * Az * dz * R / data->dt - q / nf * cin);
156 
157  /*
158  * printf("nf %g\n", nf);
159  * printf("q %g\n", q);
160  * printf("cs %g\n", cs);
161  * printf("cin %g\n", cin);
162  * printf("cg %g\n", cg);
163  * printf("cg_start %g\n", cg_start);
164  * printf("Az %g\n", Az);
165  * printf("z %g\n", z);
166  * printf("R %g\n", R);
167  * printf("dt %g\n", data->dt);
168  */
169  G_debug(6, "N_callback_solute_transport_3d: called [%i][%i][%i]", row,
170  col, depth);
171 
172  /*create the 7 point star entries */
173  mat_pos = N_create_7star(C, W, E, N, S, T, B, V);
174 
175  return mat_pos;
176 }
177 
178 /* ************************************************************************* *
179  * ************************************************************************* *
180  * ************************************************************************* */
199  N_geom_data * geom, int col,
200  int row)
201 {
202  double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0;
203 
204  double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
205 
206  double dx, dy, Az;
207 
208  double diff_x, diff_y;
209 
210  double disp_x, disp_y;
211 
212  double z;
213 
214  double diff_xw, diff_yn;
215 
216  double disp_xw, disp_yn;
217 
218  double z_xw, z_yn;
219 
220  double diff_xe, diff_ys;
221 
222  double disp_xe, disp_ys;
223 
224  double z_xe, z_ys;
225 
226  double cin = 0, cg, cg_start;
227 
228  double R, nf, cs, q;
229 
230  double C, W, E, N, S, V, NE, NW, SW, SE;
231 
232  double vw = 0, ve = 0, vn = 0, vs = 0;
233 
234  double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0;
235 
236  double Dw = 0, De = 0, Dn = 0, Ds = 0;
237 
238  double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5;
239 
241 
242  N_data_star *mat_pos;
243 
244  N_gradient_2d grad;
245 
246  /*cast the void pointer to the right data structure */
247  data = (N_solute_transport_data2d *) solutedata;
248 
249  N_get_gradient_2d(data->grad, &grad, col, row);
250 
251  dx = geom->dx;
252  dy = geom->dy;
253  Az = N_get_geom_data_area_of_cell(geom, row);
254 
255  /*read the data from the arrays */
256  cg_start = N_get_array_2d_d_value(data->c_start, col, row);
257  cg = N_get_array_2d_d_value(data->c, col, row);
258 
259  /* calculate the cell height */
260  z = N_get_array_2d_d_value(data->top, col,
261  row) -
262  N_get_array_2d_d_value(data->bottom, col, row);
263  z_xw =
264  N_get_array_2d_d_value(data->top, col - 1,
265  row) -
266  N_get_array_2d_d_value(data->bottom, col - 1, row);
267  z_xe =
268  N_get_array_2d_d_value(data->top, col + 1,
269  row) -
270  N_get_array_2d_d_value(data->bottom, col + 1, row);
271  z_yn =
272  N_get_array_2d_d_value(data->top, col,
273  row - 1) -
274  N_get_array_2d_d_value(data->bottom, col, row - 1);
275  z_ys =
276  N_get_array_2d_d_value(data->top, col,
277  row + 1) -
278  N_get_array_2d_d_value(data->bottom, col, row + 1);
279 
280  /*geometrical mean of cell height */
281  z_w = N_calc_geom_mean(z_xw, z);
282  z_e = N_calc_geom_mean(z_xe, z);
283  z_n = N_calc_geom_mean(z_yn, z);
284  z_s = N_calc_geom_mean(z_ys, z);
285 
286  /*get the surrounding diffusion tensor entries */
287  diff_x = N_get_array_2d_d_value(data->diff_x, col, row);
288  diff_y = N_get_array_2d_d_value(data->diff_y, col, row);
289  diff_xw = N_get_array_2d_d_value(data->diff_x, col - 1, row);
290  diff_xe = N_get_array_2d_d_value(data->diff_x, col + 1, row);
291  diff_yn = N_get_array_2d_d_value(data->diff_y, col, row - 1);
292  diff_ys = N_get_array_2d_d_value(data->diff_y, col, row + 1);
293 
294  /* calculate the diffusion at the cell borders using the harmonical mean */
295  Df_w = N_calc_harmonic_mean(diff_xw, diff_x);
296  Df_e = N_calc_harmonic_mean(diff_xe, diff_x);
297  Df_n = N_calc_harmonic_mean(diff_yn, diff_y);
298  Df_s = N_calc_harmonic_mean(diff_ys, diff_y);
299 
300  /* calculate the dispersion */
301  /*get the surrounding dispersion tensor entries */
302  disp_x = N_get_array_2d_d_value(data->disp_xx, col, row);
303  disp_y = N_get_array_2d_d_value(data->disp_yy, col, row);
304  if (N_get_array_2d_d_value(data->status, col - 1, row) ==
306  disp_xw = disp_x;
307  }
308  else {
309  disp_xw = N_get_array_2d_d_value(data->disp_xx, col - 1, row);
310  }
311  if (N_get_array_2d_d_value(data->status, col + 1, row) ==
313  disp_xe = disp_x;
314  }
315  else {
316  disp_xe = N_get_array_2d_d_value(data->disp_xx, col + 1, row);
317  }
318  if (N_get_array_2d_d_value(data->status, col, row - 1) ==
320  disp_yn = disp_y;
321  }
322  else {
323  disp_yn = N_get_array_2d_d_value(data->disp_yy, col, row - 1);
324  }
325  if (N_get_array_2d_d_value(data->status, col, row + 1) ==
327  disp_ys = disp_y;
328  }
329  else {
330  disp_ys = N_get_array_2d_d_value(data->disp_yy, col, row + 1);
331  }
332 
333  /* calculate the dispersion at the cell borders using the harmonical mean */
334  Ds_w = N_calc_harmonic_mean(disp_xw, disp_x);
335  Ds_e = N_calc_harmonic_mean(disp_xe, disp_x);
336  Ds_n = N_calc_harmonic_mean(disp_yn, disp_y);
337  Ds_s = N_calc_harmonic_mean(disp_ys, disp_y);
338 
339  /* put the diffusion and dispersion together */
340  Dw = ((Df_w + Ds_w)) / dx;
341  De = ((Df_e + Ds_e)) / dx;
342  Ds = ((Df_s + Ds_s)) / dy;
343  Dn = ((Df_n + Ds_n)) / dy;
344 
345  vw = -1.0 * grad.WC;
346  ve = grad.EC;
347  vs = -1.0 * grad.SC;
348  vn = grad.NC;
349 
350  if (data->stab == N_UPWIND_FULL) {
351  rw = N_full_upwinding(vw, dx, Dw);
352  re = N_full_upwinding(ve, dx, De);
353  rs = N_full_upwinding(vs, dy, Ds);
354  rn = N_full_upwinding(vn, dy, Dn);
355  }
356  else if (data->stab == N_UPWIND_EXP) {
357  rw = N_exp_upwinding(vw, dx, Dw);
358  re = N_exp_upwinding(ve, dx, De);
359  rs = N_exp_upwinding(vs, dy, Ds);
360  rn = N_exp_upwinding(vn, dy, Dn);
361  }
362 
363  /*mass balance center cell to western cell */
364  W = -1 * (Dw) * dy * z_w + vw * (1 - rw) * dy * z_w;
365  /*mass balance center cell to eastern cell */
366  E = -1 * (De) * dy * z_e + ve * (1 - re) * dy * z_e;
367  /*mass balance center cell to southern cell */
368  S = -1 * (Ds) * dx * z_s + vs * (1 - rs) * dx * z_s;
369  /*mass balance center cell to northern cell */
370  N = -1 * (Dn) * dx * z_n + vn * (1 - rn) * dx * z_n;
371 
372  NW = 0.0;
373  SW = 0.0;
374  NE = 0.0;
375  SE = 0.0;
376 
377  /* Retardation */
378  R = N_get_array_2d_d_value(data->R, col, row);
379  /* Inner sources */
380  cs = N_get_array_2d_d_value(data->cs, col, row);
381  /* effective porosity */
382  nf = N_get_array_2d_d_value(data->nf, col, row);
383  /* groundwater sources and sinks */
384  q = N_get_array_2d_d_value(data->q, col, row);
385  /* concentration of influent water */
386  cin = N_get_array_2d_d_value(data->cin, col, row);
387 
388  /*the diagonal entry of the matrix */
389  C = (Dw + vw * rw) * dy * z_w +
390  (De + ve * re) * dy * z_e +
391  (Ds + vs * rs) * dx * z_s +
392  (Dn + vn * rn) * dx * z_n + Az * z * R / data->dt - q / nf;
393 
394  /*the entry in the right side b of Ax = b */
395  V = (cs + cg_start * Az * z * R / data->dt + q / nf * cin);
396 
397  /*
398  fprintf(stderr, "nf %g\n", nf);
399  fprintf(stderr, "q %g\n", q);
400  fprintf(stderr, "cs %g\n", cs);
401  fprintf(stderr, "cin %g\n", cin);
402  fprintf(stderr, "cg %g\n", cg);
403  fprintf(stderr, "cg_start %g\n", cg_start);
404  fprintf(stderr, "Az %g\n", Az);
405  fprintf(stderr, "z %g\n", z);
406  fprintf(stderr, "R %g\n", R);
407  fprintf(stderr, "dt %g\n", data->dt);
408  */
409 
410  G_debug(6, "N_callback_solute_transport_2d: called [%i][%i]", row, col);
411 
412  /*create the 9 point star entries */
413  mat_pos = N_create_9star(C, W, E, N, S, NW, SW, NE, SE, V);
414 
415  return mat_pos;
416 }
417 
418 /* ************************************************************************* *
419  * ************************************************************************* *
420  * ************************************************************************* */
437  int depths)
438 {
440 
441  data =
442  (N_solute_transport_data3d *) G_calloc(1,
443  sizeof
445 
446  data->c = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
447  data->c_start = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
448  data->status = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
449  data->diff_x = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
450  data->diff_y = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
451  data->diff_z = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
452  data->q = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
453  data->cs = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
454  data->R = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
455  data->nf = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
456  data->cin = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
457 
458  /*Allocate the dispersivity tensor */
459  data->disp_xx = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
460  data->disp_yy = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
461  data->disp_zz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
462  data->disp_xy = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
463  data->disp_xz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
464  data->disp_yz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
465 
466 
467  data->grad = N_alloc_gradient_field_3d(cols, rows, depths);
468  data->stab = N_UPWIND_EXP;
469 
470  return data;
471 }
472 
473 /* ************************************************************************* *
474  * ************************************************************************* *
475  * ************************************************************************* */
492 {
494 
495  data =
496  (N_solute_transport_data2d *) G_calloc(1,
497  sizeof
499 
500  data->c = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
501  data->c_start = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
502  data->status = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
503  data->diff_x = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
504  data->diff_y = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
505  data->q = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
506  data->cs = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
507  data->R = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
508  data->nf = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
509  data->cin = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
510  data->top = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
511  data->bottom = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
512 
513  /*Allocate the dispersivity tensor */
514  data->disp_xx = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
515  data->disp_yy = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
516  data->disp_xy = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
517 
518  data->grad = N_alloc_gradient_field_2d(cols, rows);
519  data->stab = N_UPWIND_EXP;
520 
521  return data;
522 }
523 
524 /* ************************************************************************* *
525  * ************************************************************************* *
526  * ************************************************************************* */
534 {
535  N_free_array_3d(data->c);
536  N_free_array_3d(data->c_start);
537  N_free_array_3d(data->status);
538  N_free_array_3d(data->diff_x);
539  N_free_array_3d(data->diff_y);
540  N_free_array_3d(data->diff_z);
541  N_free_array_3d(data->q);
542  N_free_array_3d(data->cs);
543  N_free_array_3d(data->R);
544  N_free_array_3d(data->nf);
545  N_free_array_3d(data->cin);
546 
547  N_free_array_3d(data->disp_xx);
548  N_free_array_3d(data->disp_yy);
549  N_free_array_3d(data->disp_zz);
550  N_free_array_3d(data->disp_xy);
551  N_free_array_3d(data->disp_xz);
552  N_free_array_3d(data->disp_yz);
553 
554  G_free(data);
555 
556  data = NULL;
557 
558  return;
559 }
560 
561 /* ************************************************************************* *
562  * ************************************************************************* *
563  * ************************************************************************* */
571 {
572  N_free_array_2d(data->c);
573  N_free_array_2d(data->c_start);
574  N_free_array_2d(data->status);
575  N_free_array_2d(data->diff_x);
576  N_free_array_2d(data->diff_y);
577  N_free_array_2d(data->q);
578  N_free_array_2d(data->cs);
579  N_free_array_2d(data->R);
580  N_free_array_2d(data->nf);
581  N_free_array_2d(data->cin);
582  N_free_array_2d(data->top);
583  N_free_array_2d(data->bottom);
584 
585  N_free_array_2d(data->disp_xx);
586  N_free_array_2d(data->disp_yy);
587  N_free_array_2d(data->disp_xy);
588 
589  G_free(data);
590 
591  data = NULL;
592 
593  return;
594 }
595 
613 {
614  int i, j, count = 1;
615 
616  int cols, rows;
617 
618  double c;
619 
620  N_gradient_2d grad;
621 
622  cols = data->grad->cols;
623  rows = data->grad->rows;
624 
625  G_debug(2,
626  "N_calc_solute_transport_transmission_2d: calculating transmission boundary");
627 
628  for (j = 0; j < rows; j++) {
629  for (i = 0; i < cols; i++) {
630  if (N_get_array_2d_d_value(data->status, i, j) ==
632  count = 0;
633  /*get the gradient neighbours */
634  N_get_gradient_2d(data->grad, &grad, i, j);
635  c = 0;
636  /*
637  c = N_get_array_2d_d_value(data->c_start, i, j);
638  if(c > 0)
639  count++;
640  */
641 
642  if (grad.WC > 0 &&
643  !N_is_array_2d_value_null(data->c, i - 1, j)) {
644  c += N_get_array_2d_d_value(data->c, i - 1, j);
645  count++;
646  }
647  if (grad.EC < 0 &&
648  !N_is_array_2d_value_null(data->c, i + 1, j)) {
649  c += N_get_array_2d_d_value(data->c, i + 1, j);
650  count++;
651  }
652  if (grad.NC < 0 &&
653  !N_is_array_2d_value_null(data->c, i, j - 1)) {
654  c += N_get_array_2d_d_value(data->c, i, j - 1);
655  count++;
656  }
657  if (grad.SC > 0 &&
658  !N_is_array_2d_value_null(data->c, i, j + 1)) {
659  c += N_get_array_2d_d_value(data->c, i, j + 1);
660  count++;
661  }
662  if (count != 0)
663  c = c / (double)count;
664  /*make sure it is not NAN */
665  if (c > 0 || c == 0 || c < 0)
666  N_put_array_2d_d_value(data->c_start, i, j, c);
667  }
668  }
669  }
670 
671  return;
672 }
673 
689 {
690  int i, j;
691 
692  int cols, rows;
693 
694  double vx, vy, vv;
695 
696  double disp_xx, disp_yy, disp_xy;
697 
698  N_gradient_2d grad;
699 
700  cols = data->grad->cols;
701  rows = data->grad->rows;
702 
703  G_debug(2,
704  "N_calc_solute_transport_disptensor_2d: calculating the dispersivity tensor");
705 
706  for (j = 0; j < rows; j++) {
707  for (i = 0; i < cols; i++) {
708 
709  disp_xx = 0;
710  disp_yy = 0;
711  disp_xy = 0;
712 
713  /*get the gradient neighbours */
714  N_get_gradient_2d(data->grad, &grad, i, j);
715  vx = (grad.WC + grad.EC) / 2;
716  vy = (grad.NC + grad.SC) / 2;
717  vv = sqrt(vx * vx + vy * vy);
718 
719  if (vv != 0) {
720  disp_xx = data->al * vx * vx / vv + data->at * vy * vy / vv;
721  disp_yy = data->at * vx * vx / vv + data->al * vy * vy / vv;
722  disp_xy = (data->al - data->at) * vx * vy / vv;
723  }
724 
725  G_debug(5,
726  "N_calc_solute_transport_disptensor_2d: [%i][%i] disp_xx %g disp_yy %g disp_xy %g",
727  i, j, disp_xx, disp_yy, disp_xy);
728  N_put_array_2d_d_value(data->disp_xx, i, j, disp_xx);
729  N_put_array_2d_d_value(data->disp_yy, i, j, disp_yy);
730  N_put_array_2d_d_value(data->disp_xy, i, j, disp_xy);
731  }
732  }
733 
734  return;
735 }
736 
752 {
753  int i, j, k;
754 
755  int cols, rows, depths;
756 
757  double vx, vy, vz, vv;
758 
759  double disp_xx, disp_yy, disp_zz, disp_xy, disp_xz, disp_yz;
760 
761  N_gradient_3d grad;
762 
763  cols = data->grad->cols;
764  rows = data->grad->rows;
765  depths = data->grad->depths;
766 
767  G_debug(2,
768  "N_calc_solute_transport_disptensor_3d: calculating the dispersivity tensor");
769 
770  for (k = 0; k < depths; k++) {
771  for (j = 0; j < rows; j++) {
772  for (i = 0; i < cols; i++) {
773  disp_xx = 0;
774  disp_yy = 0;
775  disp_zz = 0;
776  disp_xy = 0;
777  disp_xz = 0;
778  disp_yz = 0;
779 
780  /*get the gradient neighbours */
781  N_get_gradient_3d(data->grad, &grad, i, j, k);
782  vx = (grad.WC + grad.EC) / 2;
783  vy = (grad.NC + grad.SC) / 2;
784  vz = (grad.BC + grad.TC) / 2;
785  vv = sqrt(vx * vx + vy * vy + vz * vz);
786 
787  if (vv != 0) {
788  disp_xx =
789  data->al * vx * vx / vv + data->at * vy * vy / vv +
790  data->at * vz * vz / vv;
791  disp_yy =
792  data->at * vx * vx / vv + data->al * vy * vy / vv +
793  data->at * vz * vz / vv;
794  disp_zz =
795  data->at * vx * vx / vv + data->at * vy * vy / vv +
796  data->al * vz * vz / vv;
797  disp_xy = (data->al - data->at) * vx * vy / vv;
798  disp_xz = (data->al - data->at) * vx * vz / vv;
799  disp_yz = (data->al - data->at) * vy * vz / vv;
800  }
801 
802  G_debug(5,
803  "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 ",
804  i, j, k, disp_xx, disp_yy, disp_zz, disp_xy, disp_xz,
805  disp_yz);
806  N_put_array_3d_d_value(data->disp_xx, i, j, k, disp_xx);
807  N_put_array_3d_d_value(data->disp_yy, i, j, k, disp_yy);
808  N_put_array_3d_d_value(data->disp_zz, i, j, k, disp_zz);
809  N_put_array_3d_d_value(data->disp_xy, i, j, k, disp_xy);
810  N_put_array_3d_d_value(data->disp_xz, i, j, k, disp_xz);
811  N_put_array_3d_d_value(data->disp_yz, i, j, k, disp_yz);
812  }
813  }
814  }
815 
816  return;
817 }
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d structure.
Definition: N_arrays.c:127
double EC
Definition: N_pde.h:536
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.
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
Matrix entries for a mass balance 5/7/9 star system.
Definition: N_pde.h:342
#define N_UPWIND_FULL
Definition: N_pde.h:67
tuple q
Definition: forms.py:2019
double N_calc_geom_mean(double a, double b)
Calculate the geometrical mean of values a and b.
Definition: N_tools.c:77
#define N_CELL_TRANSMISSION
Definition: N_pde.h:48
void N_free_array_3d(N_array_3d *data)
Release the memory of a N_array_3d.
Definition: N_arrays.c:777
Gradient between the cells in X and Y direction.
Definition: N_pde.h:525
void N_calc_solute_transport_transmission_2d(N_solute_transport_data2d *data)
Compute the transmission boundary condition in 2d.
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:581
int count
#define NE
Definition: dataquad.h:18
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:69
double WC
Definition: N_pde.h:536
#define NW
Definition: dataquad.h:17
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:723
double SC
Definition: N_pde.h:536
double dy
Definition: N_pde.h:134
#define SW
Definition: dataquad.h:19
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:249
double NC
Definition: N_pde.h:536
N_gradient_field_3d * grad
#define C
Definition: intr_char.c:17
Geometric information about the structured grid.
Definition: N_pde.h:127
double dz
Definition: N_pde.h:135
N_gradient_field_2d * N_alloc_gradient_field_2d(int cols, int rows)
Allocate a N_gradient_field_2d.
Definition: N_gradient.c:928
tuple data
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
void N_calc_solute_transport_disptensor_3d(N_solute_transport_data3d *data)
Compute the dispersivity tensor based on the solute transport data in 3d.
double TC
Definition: N_pde.h:536
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:536
double EC
Definition: N_pde.h:528
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
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.
Gradient between the cells in X, Y and Z direction.
Definition: N_pde.h:533
void N_free_solute_transport_data3d(N_solute_transport_data3d *data)
Release the memory of the solute transport data structure in three dimensions.
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:199
#define SE
Definition: dataquad.h:20
N_gradient_field_2d * grad
#define N_UPWIND_EXP
Definition: N_pde.h:68
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 postion col, row is of type null, otherwise 0...
Definition: N_arrays.c:228
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.
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:987
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
return NULL
Definition: dbfopen.c:1394
tuple cols
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: gis/debug.c:51
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:1026
double WC
Definition: N_pde.h:528
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:375
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:1172
#define N
Definition: inverse.c:8
double N_calc_harmonic_mean(double a, double b)
Calculate the harmonical mean of values a and b.
Definition: N_tools.c:121
double dx
Definition: N_pde.h:133
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.
double SC
Definition: N_pde.h:528
double NC
Definition: N_pde.h:528
N_solute_transport_data2d * N_alloc_solute_transport_data2d(int cols, int rows)
Alllocate memory for the solute transport data structure in two dimensions.
double N_full_upwinding(double sprod, double distance, double D)
full upwinding stabilization algorithm
Definition: N_upwind.c:33