GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /***************************************************************************** 00003 * 00004 * MODULE: Grass PDE Numerical Library 00005 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006 00006 * soerengebbert <at> gmx <dot> de 00007 * 00008 * PURPOSE: solute transport in porous media 00009 * part of the gpde library 00010 * 00011 * COPYRIGHT: (C) 2007 by the GRASS Development Team 00012 * 00013 * This program is free software under the GNU General Public 00014 * License (>=v2). Read the file COPYING that comes with GRASS 00015 * for details. 00016 * 00017 *****************************************************************************/ 00018 00019 #include "grass/N_solute_transport.h" 00020 #include <math.h> 00021 00022 /* ************************************************************************* * 00023 * ************************************************************************* * 00024 * ************************************************************************* */ 00028 N_data_star *N_callback_solute_transport_3d(void *solutedata, 00029 N_geom_data * geom, int col, 00030 int row, int depth) 00031 { 00032 double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0, Df_t = 0, Df_b = 0; 00033 double dx, dy, dz, Az; 00034 double diff_x, diff_y, diff_z; 00035 double diff_xw, diff_yn; 00036 double diff_xe, diff_ys; 00037 double diff_zt, diff_zb; 00038 double cin = 0, cg, cg_start; 00039 double R, nf, cs, q; 00040 double C, W, E, N, S, T, B, V; 00041 double vw = 0, ve = 0, vn = 0, vs = 0, vt = 0, vb = 0; 00042 double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0, Ds_t = 0, Ds_b = 0; 00043 double Dw = 0, De = 0, Dn = 0, Ds = 0, Dt = 0, Db = 0; 00044 double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5, rt = 0.5, rb = 0.5; 00045 00046 N_solute_transport_data3d *data = NULL; 00047 N_data_star *mat_pos; 00048 N_gradient_3d grad; 00049 00050 /*cast the void pointer to the right data structure */ 00051 data = (N_solute_transport_data3d *) solutedata; 00052 00053 N_get_gradient_3d(data->grad, &grad, col, row, depth); 00054 00055 dx = geom->dx; 00056 dy = geom->dy; 00057 dz = geom->dz; 00058 Az = N_get_geom_data_area_of_cell(geom, row); 00059 00060 /*read the data from the arrays */ 00061 cg_start = N_get_array_3d_d_value(data->c_start, col, row, depth); 00062 cg = N_get_array_3d_d_value(data->c, col, row, depth); 00063 00064 /*get the surrounding diffusion tensor entries */ 00065 diff_x = N_get_array_3d_d_value(data->diff_x, col, row, depth); 00066 diff_y = N_get_array_3d_d_value(data->diff_y, col, row, depth); 00067 diff_z = N_get_array_3d_d_value(data->diff_z, col, row, depth); 00068 diff_xw = N_get_array_3d_d_value(data->diff_x, col - 1, row, depth); 00069 diff_xe = N_get_array_3d_d_value(data->diff_x, col + 1, row, depth); 00070 diff_yn = N_get_array_3d_d_value(data->diff_y, col, row - 1, depth); 00071 diff_ys = N_get_array_3d_d_value(data->diff_y, col, row + 1, depth); 00072 diff_zt = N_get_array_3d_d_value(data->diff_z, col, row, depth + 1); 00073 diff_zb = N_get_array_3d_d_value(data->diff_z, col, row, depth - 1); 00074 00075 /* calculate the diffusion on the cell borders using the harmonical mean */ 00076 Df_w = N_calc_harmonic_mean(diff_xw, diff_x); 00077 Df_e = N_calc_harmonic_mean(diff_xe, diff_x); 00078 Df_n = N_calc_harmonic_mean(diff_yn, diff_y); 00079 Df_s = N_calc_harmonic_mean(diff_ys, diff_y); 00080 Df_t = N_calc_harmonic_mean(diff_zt, diff_z); 00081 Df_b = N_calc_harmonic_mean(diff_zb, diff_z); 00082 00083 /* calculate the dispersion */ 00084 /*todo */ 00085 00086 /* calculate the velocity parts with full upwinding scheme */ 00087 vw = grad.WC; 00088 ve = grad.EC; 00089 vn = grad.NC; 00090 vs = grad.SC; 00091 vt = grad.TC; 00092 vb = grad.BC; 00093 00094 /* put the diffusion and dispersion together */ 00095 Dw = ((Df_w + Ds_w)) / dx; 00096 De = ((Df_e + Ds_e)) / dx; 00097 Dn = ((Df_n + Ds_n)) / dy; 00098 Ds = ((Df_s + Ds_s)) / dy; 00099 Dt = ((Df_t + Ds_t)) / dz; 00100 Db = ((Df_b + Ds_b)) / dz; 00101 00102 rw = N_exp_upwinding(-1 * vw, dx, Dw); 00103 re = N_exp_upwinding(ve, dx, De); 00104 rs = N_exp_upwinding(-1 * vs, dy, Ds); 00105 rn = N_exp_upwinding(vn, dy, Dn); 00106 rb = N_exp_upwinding(-1 * vb, dz, Dn); 00107 rt = N_exp_upwinding(vt, dz, Dn); 00108 00109 /*mass balance center cell to western cell */ 00110 W = -1 * (Dw) * dy * dz - vw * (1 - rw) * dy * dz; 00111 /*mass balance center cell to eastern cell */ 00112 E = -1 * (De) * dy * dz + ve * (1 - re) * dy * dz; 00113 /*mass balance center cell to southern cell */ 00114 S = -1 * (Ds) * dx * dz - vs * (1 - rs) * dx * dz; 00115 /*mass balance center cell to northern cell */ 00116 N = -1 * (Dn) * dx * dz + vn * (1 - rn) * dx * dz; 00117 /*mass balance center cell to bottom cell */ 00118 B = -1 * (Db) * Az - vb * (1 - rb) * Az; 00119 /*mass balance center cell to top cell */ 00120 T = -1 * (Dt) * Az + vt * (1 - rt) * Az; 00121 00122 /* Retardation */ 00123 R = N_get_array_3d_d_value(data->R, col, row, depth); 00124 /* Inner sources */ 00125 cs = N_get_array_3d_d_value(data->cs, col, row, depth); 00126 /* effective porosity */ 00127 nf = N_get_array_3d_d_value(data->nf, col, row, depth); 00128 /* groundwater sources and sinks */ 00129 q = N_get_array_3d_d_value(data->q, col, row, depth); 00130 /* concentration of influent water */ 00131 cin = N_get_array_3d_d_value(data->cin, col, row, depth); 00132 00133 /*the diagonal entry of the matrix */ 00134 C = ((Dw - vw) * dy * dz + 00135 (De + ve) * dy * dz + 00136 (Ds - vs) * dx * dz + 00137 (Dn + vn) * dx * dz + 00138 (Db - vb) * Az + (Dt + vt) * Az + Az * dz * R / data->dt - q / nf); 00139 00140 /*the entry in the right side b of Ax = b */ 00141 V = (cs + cg_start * Az * dz * R / data->dt - q / nf * cin); 00142 00143 /* 00144 * printf("nf %g\n", nf); 00145 * printf("q %g\n", q); 00146 * printf("cs %g\n", cs); 00147 * printf("cin %g\n", cin); 00148 * printf("cg %g\n", cg); 00149 * printf("cg_start %g\n", cg_start); 00150 * printf("Az %g\n", Az); 00151 * printf("z %g\n", z); 00152 * printf("R %g\n", R); 00153 * printf("dt %g\n", data->dt); 00154 */ 00155 G_debug(6, "N_callback_solute_transport_3d: called [%i][%i][%i]", row, 00156 col, depth); 00157 00158 /*create the 7 point star entries */ 00159 mat_pos = N_create_7star(C, W, E, N, S, T, B, V); 00160 00161 return mat_pos; 00162 } 00163 00164 /* ************************************************************************* * 00165 * ************************************************************************* * 00166 * ************************************************************************* */ 00184 N_data_star *N_callback_solute_transport_2d(void *solutedata, 00185 N_geom_data * geom, int col, 00186 int row) 00187 { 00188 double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0; 00189 double z_e = 0, z_w = 0, z_n = 0, z_s = 0; 00190 double dx, dy, Az; 00191 double diff_x, diff_y; 00192 double disp_x, disp_y; 00193 double z; 00194 double diff_xw, diff_yn; 00195 double disp_xw, disp_yn; 00196 double z_xw, z_yn; 00197 double diff_xe, diff_ys; 00198 double disp_xe, disp_ys; 00199 double z_xe, z_ys; 00200 double cin = 0, cg, cg_start; 00201 double R, nf, cs, q; 00202 double C, W, E, N, S, V, NE, NW, SW, SE; 00203 double vw = 0, ve = 0, vn = 0, vs = 0; 00204 double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0; 00205 double Dw = 0, De = 0, Dn = 0, Ds = 0; 00206 double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5; 00207 00208 N_solute_transport_data2d *data = NULL; 00209 N_data_star *mat_pos; 00210 N_gradient_2d grad; 00211 00212 /*cast the void pointer to the right data structure */ 00213 data = (N_solute_transport_data2d *) solutedata; 00214 00215 N_get_gradient_2d(data->grad, &grad, col, row); 00216 00217 dx = geom->dx; 00218 dy = geom->dy; 00219 Az = N_get_geom_data_area_of_cell(geom, row); 00220 00221 /*read the data from the arrays */ 00222 cg_start = N_get_array_2d_d_value(data->c_start, col, row); 00223 cg = N_get_array_2d_d_value(data->c, col, row); 00224 00225 /* calculate the cell height */ 00226 z = N_get_array_2d_d_value(data->top, col, 00227 row) - 00228 N_get_array_2d_d_value(data->bottom, col, row); 00229 z_xw = 00230 N_get_array_2d_d_value(data->top, col - 1, 00231 row) - 00232 N_get_array_2d_d_value(data->bottom, col - 1, row); 00233 z_xe = 00234 N_get_array_2d_d_value(data->top, col + 1, 00235 row) - 00236 N_get_array_2d_d_value(data->bottom, col + 1, row); 00237 z_yn = 00238 N_get_array_2d_d_value(data->top, col, 00239 row - 1) - 00240 N_get_array_2d_d_value(data->bottom, col, row - 1); 00241 z_ys = 00242 N_get_array_2d_d_value(data->top, col, 00243 row + 1) - 00244 N_get_array_2d_d_value(data->bottom, col, row + 1); 00245 00246 /*geometrical mean of cell height */ 00247 z_w = N_calc_geom_mean(z_xw, z); 00248 z_e = N_calc_geom_mean(z_xe, z); 00249 z_n = N_calc_geom_mean(z_yn, z); 00250 z_s = N_calc_geom_mean(z_ys, z); 00251 00252 /*get the surrounding diffusion tensor entries */ 00253 diff_x = N_get_array_2d_d_value(data->diff_x, col, row); 00254 diff_y = N_get_array_2d_d_value(data->diff_y, col, row); 00255 diff_xw = N_get_array_2d_d_value(data->diff_x, col - 1, row); 00256 diff_xe = N_get_array_2d_d_value(data->diff_x, col + 1, row); 00257 diff_yn = N_get_array_2d_d_value(data->diff_y, col, row - 1); 00258 diff_ys = N_get_array_2d_d_value(data->diff_y, col, row + 1); 00259 00260 /* calculate the diffusion at the cell borders using the harmonical mean */ 00261 Df_w = N_calc_harmonic_mean(diff_xw, diff_x); 00262 Df_e = N_calc_harmonic_mean(diff_xe, diff_x); 00263 Df_n = N_calc_harmonic_mean(diff_yn, diff_y); 00264 Df_s = N_calc_harmonic_mean(diff_ys, diff_y); 00265 00266 /* calculate the dispersion */ 00267 /*get the surrounding dispersion tensor entries */ 00268 disp_x = N_get_array_2d_d_value(data->disp_xx, col, row); 00269 disp_y = N_get_array_2d_d_value(data->disp_yy, col, row); 00270 if (N_get_array_2d_d_value(data->status, col - 1, row) == 00271 N_CELL_TRANSMISSION) { 00272 disp_xw = disp_x; 00273 } 00274 else { 00275 disp_xw = N_get_array_2d_d_value(data->disp_xx, col - 1, row); 00276 } 00277 if (N_get_array_2d_d_value(data->status, col + 1, row) == 00278 N_CELL_TRANSMISSION) { 00279 disp_xe = disp_x; 00280 } 00281 else { 00282 disp_xe = N_get_array_2d_d_value(data->disp_xx, col + 1, row); 00283 } 00284 if (N_get_array_2d_d_value(data->status, col, row - 1) == 00285 N_CELL_TRANSMISSION) { 00286 disp_yn = disp_y; 00287 } 00288 else { 00289 disp_yn = N_get_array_2d_d_value(data->disp_yy, col, row - 1); 00290 } 00291 if (N_get_array_2d_d_value(data->status, col, row + 1) == 00292 N_CELL_TRANSMISSION) { 00293 disp_ys = disp_y; 00294 } 00295 else { 00296 disp_ys = N_get_array_2d_d_value(data->disp_yy, col, row + 1); 00297 } 00298 00299 /* calculate the dispersion at the cell borders using the harmonical mean */ 00300 Ds_w = N_calc_harmonic_mean(disp_xw, disp_x); 00301 Ds_e = N_calc_harmonic_mean(disp_xe, disp_x); 00302 Ds_n = N_calc_harmonic_mean(disp_yn, disp_y); 00303 Ds_s = N_calc_harmonic_mean(disp_ys, disp_y); 00304 00305 /* put the diffusion and dispersion together */ 00306 Dw = ((Df_w + Ds_w)) / dx; 00307 De = ((Df_e + Ds_e)) / dx; 00308 Ds = ((Df_s + Ds_s)) / dy; 00309 Dn = ((Df_n + Ds_n)) / dy; 00310 00311 vw = -1.0 * grad.WC; 00312 ve = grad.EC; 00313 vs = -1.0 * grad.SC; 00314 vn = grad.NC; 00315 00316 if (data->stab == N_UPWIND_FULL) { 00317 rw = N_full_upwinding(vw, dx, Dw); 00318 re = N_full_upwinding(ve, dx, De); 00319 rs = N_full_upwinding(vs, dy, Ds); 00320 rn = N_full_upwinding(vn, dy, Dn); 00321 } 00322 else if (data->stab == N_UPWIND_EXP) { 00323 rw = N_exp_upwinding(vw, dx, Dw); 00324 re = N_exp_upwinding(ve, dx, De); 00325 rs = N_exp_upwinding(vs, dy, Ds); 00326 rn = N_exp_upwinding(vn, dy, Dn); 00327 } 00328 00329 /*mass balance center cell to western cell */ 00330 W = -1 * (Dw) * dy * z_w + vw * (1 - rw) * dy * z_w; 00331 /*mass balance center cell to eastern cell */ 00332 E = -1 * (De) * dy * z_e + ve * (1 - re) * dy * z_e; 00333 /*mass balance center cell to southern cell */ 00334 S = -1 * (Ds) * dx * z_s + vs * (1 - rs) * dx * z_s; 00335 /*mass balance center cell to northern cell */ 00336 N = -1 * (Dn) * dx * z_n + vn * (1 - rn) * dx * z_n; 00337 00338 NW = 0.0; 00339 SW = 0.0; 00340 NE = 0.0; 00341 SE = 0.0; 00342 00343 /* Retardation */ 00344 R = N_get_array_2d_d_value(data->R, col, row); 00345 /* Inner sources */ 00346 cs = N_get_array_2d_d_value(data->cs, col, row); 00347 /* effective porosity */ 00348 nf = N_get_array_2d_d_value(data->nf, col, row); 00349 /* groundwater sources and sinks */ 00350 q = N_get_array_2d_d_value(data->q, col, row); 00351 /* concentration of influent water */ 00352 cin = N_get_array_2d_d_value(data->cin, col, row); 00353 00354 /*the diagonal entry of the matrix */ 00355 C = (Dw + vw * rw) * dy * z_w + 00356 (De + ve * re) * dy * z_e + 00357 (Ds + vs * rs) * dx * z_s + 00358 (Dn + vn * rn) * dx * z_n + Az * z * R / data->dt - q / nf; 00359 00360 /*the entry in the right side b of Ax = b */ 00361 V = (cs + cg_start * Az * z * R / data->dt + q / nf * cin); 00362 00363 /* 00364 fprintf(stderr, "nf %g\n", nf); 00365 fprintf(stderr, "q %g\n", q); 00366 fprintf(stderr, "cs %g\n", cs); 00367 fprintf(stderr, "cin %g\n", cin); 00368 fprintf(stderr, "cg %g\n", cg); 00369 fprintf(stderr, "cg_start %g\n", cg_start); 00370 fprintf(stderr, "Az %g\n", Az); 00371 fprintf(stderr, "z %g\n", z); 00372 fprintf(stderr, "R %g\n", R); 00373 fprintf(stderr, "dt %g\n", data->dt); 00374 */ 00375 00376 G_debug(6, "N_callback_solute_transport_2d: called [%i][%i]", row, col); 00377 00378 /*create the 9 point star entries */ 00379 mat_pos = N_create_9star(C, W, E, N, S, NW, SW, NE, SE, V); 00380 00381 return mat_pos; 00382 } 00383 00384 /* ************************************************************************* * 00385 * ************************************************************************* * 00386 * ************************************************************************* */ 00402 N_solute_transport_data3d *N_alloc_solute_transport_data3d(int cols, int rows, 00403 int depths) 00404 { 00405 N_solute_transport_data3d *data = NULL; 00406 00407 data = 00408 (N_solute_transport_data3d *) G_calloc(1, 00409 sizeof 00410 (N_solute_transport_data3d)); 00411 00412 data->c = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00413 data->c_start = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00414 data->status = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00415 data->diff_x = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00416 data->diff_y = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00417 data->diff_z = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00418 data->q = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00419 data->cs = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00420 data->R = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00421 data->nf = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00422 data->cin = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00423 00424 /*Allocate the dispersivity tensor */ 00425 data->disp_xx = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00426 data->disp_yy = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00427 data->disp_zz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00428 data->disp_xy = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00429 data->disp_xz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00430 data->disp_yz = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE); 00431 00432 00433 data->grad = N_alloc_gradient_field_3d(cols, rows, depths); 00434 data->stab = N_UPWIND_EXP; 00435 00436 return data; 00437 } 00438 00439 /* ************************************************************************* * 00440 * ************************************************************************* * 00441 * ************************************************************************* */ 00457 N_solute_transport_data2d *N_alloc_solute_transport_data2d(int cols, int rows) 00458 { 00459 N_solute_transport_data2d *data = NULL; 00460 00461 data = 00462 (N_solute_transport_data2d *) G_calloc(1, 00463 sizeof 00464 (N_solute_transport_data2d)); 00465 00466 data->c = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00467 data->c_start = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00468 data->status = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00469 data->diff_x = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00470 data->diff_y = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00471 data->q = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00472 data->cs = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00473 data->R = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00474 data->nf = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00475 data->cin = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00476 data->top = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00477 data->bottom = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00478 00479 /*Allocate the dispersivity tensor */ 00480 data->disp_xx = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00481 data->disp_yy = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00482 data->disp_xy = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE); 00483 00484 data->grad = N_alloc_gradient_field_2d(cols, rows); 00485 data->stab = N_UPWIND_EXP; 00486 00487 return data; 00488 } 00489 00490 /* ************************************************************************* * 00491 * ************************************************************************* * 00492 * ************************************************************************* */ 00499 void N_free_solute_transport_data3d(N_solute_transport_data3d * data) 00500 { 00501 N_free_array_3d(data->c); 00502 N_free_array_3d(data->c_start); 00503 N_free_array_3d(data->status); 00504 N_free_array_3d(data->diff_x); 00505 N_free_array_3d(data->diff_y); 00506 N_free_array_3d(data->diff_z); 00507 N_free_array_3d(data->q); 00508 N_free_array_3d(data->cs); 00509 N_free_array_3d(data->R); 00510 N_free_array_3d(data->nf); 00511 N_free_array_3d(data->cin); 00512 00513 N_free_array_3d(data->disp_xx); 00514 N_free_array_3d(data->disp_yy); 00515 N_free_array_3d(data->disp_zz); 00516 N_free_array_3d(data->disp_xy); 00517 N_free_array_3d(data->disp_xz); 00518 N_free_array_3d(data->disp_yz); 00519 00520 G_free(data); 00521 00522 data = NULL; 00523 00524 return; 00525 } 00526 00527 /* ************************************************************************* * 00528 * ************************************************************************* * 00529 * ************************************************************************* */ 00536 void N_free_solute_transport_data2d(N_solute_transport_data2d * data) 00537 { 00538 N_free_array_2d(data->c); 00539 N_free_array_2d(data->c_start); 00540 N_free_array_2d(data->status); 00541 N_free_array_2d(data->diff_x); 00542 N_free_array_2d(data->diff_y); 00543 N_free_array_2d(data->q); 00544 N_free_array_2d(data->cs); 00545 N_free_array_2d(data->R); 00546 N_free_array_2d(data->nf); 00547 N_free_array_2d(data->cin); 00548 N_free_array_2d(data->top); 00549 N_free_array_2d(data->bottom); 00550 00551 N_free_array_2d(data->disp_xx); 00552 N_free_array_2d(data->disp_yy); 00553 N_free_array_2d(data->disp_xy); 00554 00555 G_free(data); 00556 00557 data = NULL; 00558 00559 return; 00560 } 00561 00578 void N_calc_solute_transport_transmission_2d(N_solute_transport_data2d * data) 00579 { 00580 int i, j, count = 1; 00581 int cols, rows; 00582 double c; 00583 N_gradient_2d grad; 00584 00585 cols = data->grad->cols; 00586 rows = data->grad->rows; 00587 00588 G_debug(2, 00589 "N_calc_solute_transport_transmission_2d: calculating transmission boundary"); 00590 00591 for (j = 0; j < rows; j++) { 00592 for (i = 0; i < cols; i++) { 00593 if (N_get_array_2d_d_value(data->status, i, j) == 00594 N_CELL_TRANSMISSION) { 00595 count = 0; 00596 /*get the gradient neighbours */ 00597 N_get_gradient_2d(data->grad, &grad, i, j); 00598 c = 0; 00599 /* 00600 c = N_get_array_2d_d_value(data->c_start, i, j); 00601 if(c > 0) 00602 count++; 00603 */ 00604 00605 if (grad.WC > 0 && 00606 !N_is_array_2d_value_null(data->c, i - 1, j)) { 00607 c += N_get_array_2d_d_value(data->c, i - 1, j); 00608 count++; 00609 } 00610 if (grad.EC < 0 && 00611 !N_is_array_2d_value_null(data->c, i + 1, j)) { 00612 c += N_get_array_2d_d_value(data->c, i + 1, j); 00613 count++; 00614 } 00615 if (grad.NC < 0 && 00616 !N_is_array_2d_value_null(data->c, i, j - 1)) { 00617 c += N_get_array_2d_d_value(data->c, i, j - 1); 00618 count++; 00619 } 00620 if (grad.SC > 0 && 00621 !N_is_array_2d_value_null(data->c, i, j + 1)) { 00622 c += N_get_array_2d_d_value(data->c, i, j + 1); 00623 count++; 00624 } 00625 if (count != 0) 00626 c = c / (double)count; 00627 /*make sure it is not NAN */ 00628 if (c > 0 || c == 0 || c < 0) 00629 N_put_array_2d_d_value(data->c_start, i, j, c); 00630 } 00631 } 00632 } 00633 00634 return; 00635 } 00636 00651 void N_calc_solute_transport_disptensor_2d(N_solute_transport_data2d * data) 00652 { 00653 int i, j; 00654 int cols, rows; 00655 double vx, vy, vv; 00656 double disp_xx, disp_yy, disp_xy; 00657 N_gradient_2d grad; 00658 00659 cols = data->grad->cols; 00660 rows = data->grad->rows; 00661 00662 G_debug(2, 00663 "N_calc_solute_transport_disptensor_2d: calculating the dispersivity tensor"); 00664 00665 for (j = 0; j < rows; j++) { 00666 for (i = 0; i < cols; i++) { 00667 00668 disp_xx = 0; 00669 disp_yy = 0; 00670 disp_xy = 0; 00671 00672 /*get the gradient neighbours */ 00673 N_get_gradient_2d(data->grad, &grad, i, j); 00674 vx = (grad.WC + grad.EC) / 2; 00675 vy = (grad.NC + grad.SC) / 2; 00676 vv = sqrt(vx * vx + vy * vy); 00677 00678 if (vv != 0) { 00679 disp_xx = data->al * vx * vx / vv + data->at * vy * vy / vv; 00680 disp_yy = data->at * vx * vx / vv + data->al * vy * vy / vv; 00681 disp_xy = (data->al - data->at) * vx * vy / vv; 00682 } 00683 00684 G_debug(5, 00685 "N_calc_solute_transport_disptensor_2d: [%i][%i] disp_xx %g disp_yy %g disp_xy %g", 00686 i, j, disp_xx, disp_yy, disp_xy); 00687 N_put_array_2d_d_value(data->disp_xx, i, j, disp_xx); 00688 N_put_array_2d_d_value(data->disp_yy, i, j, disp_yy); 00689 N_put_array_2d_d_value(data->disp_xy, i, j, disp_xy); 00690 } 00691 } 00692 00693 return; 00694 } 00695 00710 void N_calc_solute_transport_disptensor_3d(N_solute_transport_data3d * data) 00711 { 00712 int i, j, k; 00713 int cols, rows, depths; 00714 double vx, vy, vz, vv; 00715 double disp_xx, disp_yy, disp_zz, disp_xy, disp_xz, disp_yz; 00716 N_gradient_3d grad; 00717 00718 cols = data->grad->cols; 00719 rows = data->grad->rows; 00720 depths = data->grad->depths; 00721 00722 G_debug(2, 00723 "N_calc_solute_transport_disptensor_3d: calculating the dispersivity tensor"); 00724 00725 for (k = 0; k < depths; k++) { 00726 for (j = 0; j < rows; j++) { 00727 for (i = 0; i < cols; i++) { 00728 disp_xx = 0; 00729 disp_yy = 0; 00730 disp_zz = 0; 00731 disp_xy = 0; 00732 disp_xz = 0; 00733 disp_yz = 0; 00734 00735 /*get the gradient neighbours */ 00736 N_get_gradient_3d(data->grad, &grad, i, j, k); 00737 vx = (grad.WC + grad.EC) / 2; 00738 vy = (grad.NC + grad.SC) / 2; 00739 vz = (grad.BC + grad.TC) / 2; 00740 vv = sqrt(vx * vx + vy * vy + vz * vz); 00741 00742 if (vv != 0) { 00743 disp_xx = 00744 data->al * vx * vx / vv + data->at * vy * vy / vv + 00745 data->at * vz * vz / vv; 00746 disp_yy = 00747 data->at * vx * vx / vv + data->al * vy * vy / vv + 00748 data->at * vz * vz / vv; 00749 disp_zz = 00750 data->at * vx * vx / vv + data->at * vy * vy / vv + 00751 data->al * vz * vz / vv; 00752 disp_xy = (data->al - data->at) * vx * vy / vv; 00753 disp_xz = (data->al - data->at) * vx * vz / vv; 00754 disp_yz = (data->al - data->at) * vy * vz / vv; 00755 } 00756 00757 G_debug(5, 00758 "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 ", 00759 i, j, k, disp_xx, disp_yy, disp_zz, disp_xy, disp_xz, 00760 disp_yz); 00761 N_put_array_3d_d_value(data->disp_xx, i, j, k, disp_xx); 00762 N_put_array_3d_d_value(data->disp_yy, i, j, k, disp_yy); 00763 N_put_array_3d_d_value(data->disp_zz, i, j, k, disp_zz); 00764 N_put_array_3d_d_value(data->disp_xy, i, j, k, disp_xy); 00765 N_put_array_3d_d_value(data->disp_xz, i, j, k, disp_xz); 00766 N_put_array_3d_d_value(data->disp_yz, i, j, k, disp_yz); 00767 } 00768 } 00769 } 00770 00771 return; 00772 }