GRASS GIS 8 Programmer's Manual 8.3.2(2024)-exported
Loading...
Searching...
No Matches
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
394 data = (N_solute_transport_data3d *)G_calloc(
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
446 data = (N_solute_transport_data2d *)G_calloc(
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);
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
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);
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
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
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:150
#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
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
int count
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
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
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
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
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
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_field_2d * N_alloc_gradient_field_2d(int cols, int rows)
Allocate a N_gradient_field_2d.
Definition n_gradient.c:896
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_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_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_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_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_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.
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_2d(void *solutedata, N_geom_data *geom, int col, int row)
This callback function creates the mass balance of a 5 point star.
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.
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