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