GRASS Programmer's Manual  6.4.2(2012)
N_solute_transport.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines