GRASS Programmer's Manual  6.4.2(2012)
N_gwflow.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:      groundwater flow in porous media 
00009 *               part of the gpde library
00010 *
00011 * COPYRIGHT:    (C) 2000 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_gwflow.h"
00020 
00021 /* *************************************************************** */
00022 /* ***************** N_gwflow_data3d ***************************** */
00023 /* *************************************************************** */
00038 N_gwflow_data3d *N_alloc_gwflow_data3d(int cols, int rows, int depths,
00039                                        int river, int drain)
00040 {
00041     N_gwflow_data3d *data;
00042 
00043     data = (N_gwflow_data3d *) G_calloc(1, sizeof(N_gwflow_data3d));
00044 
00045     data->phead = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00046     data->phead_start = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00047     data->status = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00048     data->hc_x = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00049     data->hc_y = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00050     data->hc_z = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00051     data->q = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00052     data->s = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00053     data->nf = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00054     data->r = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00055 
00056     if (river) {
00057         data->river_head =
00058             N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00059         data->river_leak =
00060             N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00061         data->river_bed = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00062     }
00063     else {
00064         data->river_head = NULL;
00065         data->river_leak = NULL;
00066         data->river_bed = NULL;
00067     }
00068 
00069     if (drain) {
00070         data->drain_leak =
00071             N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00072         data->drain_bed = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
00073     }
00074     else {
00075         data->drain_leak = NULL;
00076         data->drain_bed = NULL;
00077     }
00078 
00079     return data;
00080 }
00081 
00082 /* *************************************************************** */
00083 /* ********************* N_free_gwflow_data3d ******************** */
00084 /* *************************************************************** */
00092 void N_free_gwflow_data3d(N_gwflow_data3d * data)
00093 {
00094     if (data->phead)
00095         N_free_array_3d(data->phead);
00096     if (data->phead_start)
00097         N_free_array_3d(data->phead_start);
00098     if (data->status)
00099         N_free_array_3d(data->status);
00100     if (data->hc_x)
00101         N_free_array_3d(data->hc_x);
00102     if (data->hc_y)
00103         N_free_array_3d(data->hc_y);
00104     if (data->hc_z)
00105         N_free_array_3d(data->hc_z);
00106     if (data->q)
00107         N_free_array_3d(data->q);
00108     if (data->s)
00109         N_free_array_3d(data->s);
00110     if (data->nf)
00111         N_free_array_3d(data->nf);
00112     if (data->r)
00113         N_free_array_2d(data->r);
00114     if (data->river_head)
00115         N_free_array_3d(data->river_head);
00116     if (data->river_leak)
00117         N_free_array_3d(data->river_leak);
00118     if (data->river_bed)
00119         N_free_array_3d(data->river_bed);
00120     if (data->drain_leak)
00121         N_free_array_3d(data->drain_leak);
00122     if (data->drain_bed)
00123         N_free_array_3d(data->drain_bed);
00124 
00125     G_free(data);
00126 
00127     data = NULL;
00128 
00129     return;
00130 }
00131 
00132 /* *************************************************************** */
00133 /* ******************** N_alloc_gwflow_data2d ******************** */
00134 /* *************************************************************** */
00148 N_gwflow_data2d *N_alloc_gwflow_data2d(int cols, int rows, int river,
00149                                        int drain)
00150 {
00151     N_gwflow_data2d *data;
00152 
00153     data = (N_gwflow_data2d *) G_calloc(1, sizeof(N_gwflow_data2d));
00154 
00155     data->phead = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00156     data->phead_start = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00157     data->status = N_alloc_array_2d(cols, rows, 1, CELL_TYPE);
00158     data->hc_x = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00159     data->hc_y = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00160     data->q = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00161     data->s = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00162     data->nf = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00163     data->r = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00164     data->top = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00165     data->bottom = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00166 
00167     if (river) {
00168         data->river_head = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00169         data->river_leak = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00170         data->river_bed = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00171     }
00172     else {
00173         data->river_head = NULL;
00174         data->river_leak = NULL;
00175         data->river_bed = NULL;
00176     }
00177 
00178     if (drain) {
00179         data->drain_leak = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00180         data->drain_bed = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
00181     }
00182     else {
00183         data->drain_leak = NULL;
00184         data->drain_bed = NULL;
00185     }
00186 
00187 
00188     return data;
00189 }
00190 
00191 /* *************************************************************** */
00192 /* ****************** N_free_gwflow_data2d *********************** */
00193 /* *************************************************************** */
00200 void N_free_gwflow_data2d(N_gwflow_data2d * data)
00201 {
00202     if (data->phead)
00203         N_free_array_2d(data->phead);
00204     if (data->phead_start)
00205         N_free_array_2d(data->phead_start);
00206     if (data->status)
00207         N_free_array_2d(data->status);
00208     if (data->hc_x)
00209         N_free_array_2d(data->hc_x);
00210     if (data->hc_y)
00211         N_free_array_2d(data->hc_y);
00212     if (data->q)
00213         N_free_array_2d(data->q);
00214     if (data->s)
00215         N_free_array_2d(data->s);
00216     if (data->nf)
00217         N_free_array_2d(data->nf);
00218     if (data->r)
00219         N_free_array_2d(data->r);
00220     if (data->top)
00221         N_free_array_2d(data->top);
00222     if (data->bottom)
00223         N_free_array_2d(data->bottom);
00224     if (data->river_head)
00225         N_free_array_2d(data->river_head);
00226     if (data->river_leak)
00227         N_free_array_2d(data->river_leak);
00228     if (data->river_bed)
00229         N_free_array_2d(data->river_bed);
00230     if (data->drain_leak)
00231         N_free_array_2d(data->drain_leak);
00232     if (data->drain_bed)
00233         N_free_array_2d(data->drain_bed);
00234 
00235     G_free(data);
00236 
00237     data = NULL;;
00238 
00239     return;
00240 }
00241 
00242 /* *************************************************************** */
00243 /* ***************** N_callback_gwflow_3d ************************ */
00244 /* *************************************************************** */
00263 N_data_star *N_callback_gwflow_3d(void *gwdata, N_geom_data * geom, int col,
00264                                   int row, int depth)
00265 {
00266     double hc_e = 0, hc_w = 0, hc_n = 0, hc_s = 0, hc_t = 0, hc_b = 0;
00267     double dx, dy, dz, Ax, Ay, Az;
00268     double hc_x, hc_y, hc_z;
00269     double hc_xw, hc_yn, hc_zt;
00270     double hc_xe, hc_ys, hc_zb;
00271     double hc_start;
00272     double Ss, r, nf, q;
00273     double C, W, E, N, S, T, B, V;
00274     N_data_star *mat_pos;
00275     N_gwflow_data3d *data;
00276 
00277     /*cast the void pointer to the right data structure */
00278     data = (N_gwflow_data3d *) gwdata;
00279 
00280     dx = geom->dx;
00281     dy = geom->dy;
00282     dz = geom->dz;
00283     Az = N_get_geom_data_area_of_cell(geom, row);
00284     Ay = geom->dx * geom->dz;
00285     Ax = geom->dz * geom->dy;
00286 
00287     /*read the data from the arrays */
00288     hc_start = N_get_array_3d_d_value(data->phead_start, col, row, depth);
00289 
00290     hc_x = N_get_array_3d_d_value(data->hc_x, col, row, depth);
00291     hc_y = N_get_array_3d_d_value(data->hc_y, col, row, depth);
00292     hc_z = N_get_array_3d_d_value(data->hc_z, col, row, depth);
00293 
00294     hc_xw = N_get_array_3d_d_value(data->hc_x, col - 1, row, depth);
00295     hc_xe = N_get_array_3d_d_value(data->hc_x, col + 1, row, depth);
00296     hc_yn = N_get_array_3d_d_value(data->hc_y, col, row - 1, depth);
00297     hc_ys = N_get_array_3d_d_value(data->hc_y, col, row + 1, depth);
00298     hc_zt = N_get_array_3d_d_value(data->hc_z, col, row, depth + 1);
00299     hc_zb = N_get_array_3d_d_value(data->hc_z, col, row, depth - 1);
00300 
00301     hc_w = N_calc_harmonic_mean(hc_xw, hc_x);
00302     hc_e = N_calc_harmonic_mean(hc_xe, hc_x);
00303     hc_n = N_calc_harmonic_mean(hc_yn, hc_y);
00304     hc_s = N_calc_harmonic_mean(hc_ys, hc_y);
00305     hc_t = N_calc_harmonic_mean(hc_zt, hc_z);
00306     hc_b = N_calc_harmonic_mean(hc_zb, hc_z);
00307 
00308     /*inner sources */
00309     q = N_get_array_3d_d_value(data->q, col, row, depth);
00310     /*specific yield */
00311     Ss = N_get_array_3d_d_value(data->s, col, row, depth);
00312     /*porosity */
00313     nf = N_get_array_3d_d_value(data->nf, col, row, depth);
00314 
00315     /*mass balance center cell to western cell */
00316     W = -1 * Ax * hc_w / dx;
00317     /*mass balance center cell to eastern cell */
00318     E = -1 * Ax * hc_e / dx;
00319     /*mass balance center cell to northern cell */
00320     N = -1 * Ay * hc_n / dy;
00321     /*mass balance center cell to southern cell */
00322     S = -1 * Ay * hc_s / dy;
00323     /*mass balance center cell to top cell */
00324     T = -1 * Az * hc_t / dz;
00325     /*mass balance center cell to bottom cell */
00326     B = -1 * Az * hc_b / dz;
00327 
00328     /*specific yield */
00329     Ss = Az * dz * Ss;
00330 
00331     /*the diagonal entry of the matrix */
00332     C = -1 * (W + E + N + S + T + B - Ss / data->dt * Az);
00333 
00334     /*the entry in the right side b of Ax = b */
00335     V = (q + hc_start * Ss / data->dt * Az);
00336 
00337     /*only the top cells will have recharge */
00338     if (depth == geom->depths - 2) {
00339         r = N_get_array_2d_d_value(data->r, col, row);
00340         V += r * Az;
00341     }
00342 
00343     G_debug(5, "N_callback_gwflow_3d: called [%i][%i][%i]", depth, col, row);
00344 
00345     /*create the 7 point star entries */
00346     mat_pos = N_create_7star(C, W, E, N, S, T, B, V);
00347 
00348     return mat_pos;
00349 }
00350 
00351 /* *************************************************************** */
00352 /* ****************** N_callback_gwflow_2d *********************** */
00353 /* *************************************************************** */
00370 N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
00371                                   int row)
00372 {
00373     double T_e = 0, T_w = 0, T_n = 0, T_s = 0;
00374     double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
00375     double dx, dy, Az;
00376     double hc_x, hc_y;
00377     double z, top;
00378     double hc_xw, hc_yn;
00379     double z_xw, z_yn;
00380     double hc_xe, hc_ys;
00381     double z_xe, z_ys;
00382     double hc, hc_start;
00383     double Ss, r, nf, q;
00384     double C, W, E, N, S, V;
00385     N_gwflow_data2d *data;
00386     N_data_star *mat_pos;
00387     double river_vect = 0;      /*entry in vector */
00388     double river_mat = 0;       /*entry in matrix */
00389     double drain_vect = 0;      /*entry in vector */
00390     double drain_mat = 0;       /*entry in matrix */
00391 
00392     /*cast the void pointer to the right data structure */
00393     data = (N_gwflow_data2d *) gwdata;
00394 
00395     dx = geom->dx;
00396     dy = geom->dy;
00397     Az = N_get_geom_data_area_of_cell(geom, row);
00398 
00399     /*read the data from the arrays */
00400     hc_start = N_get_array_2d_d_value(data->phead_start, col, row);
00401     hc = N_get_array_2d_d_value(data->phead, col, row);
00402     top = N_get_array_2d_d_value(data->top, col, row);
00403 
00404 
00405     if (hc > top) {             /*If the aquifer is confined */
00406         z = N_get_array_2d_d_value(data->top, col,
00407                                    row) -
00408             N_get_array_2d_d_value(data->bottom, col, row);
00409         z_xw =
00410             N_get_array_2d_d_value(data->top, col - 1,
00411                                    row) -
00412             N_get_array_2d_d_value(data->bottom, col - 1, row);
00413         z_xe =
00414             N_get_array_2d_d_value(data->top, col + 1,
00415                                    row) -
00416             N_get_array_2d_d_value(data->bottom, col + 1, row);
00417         z_yn =
00418             N_get_array_2d_d_value(data->top, col,
00419                                    row - 1) -
00420             N_get_array_2d_d_value(data->bottom, col, row - 1);
00421         z_ys =
00422             N_get_array_2d_d_value(data->top, col,
00423                                    row + 1) -
00424             N_get_array_2d_d_value(data->bottom, col, row + 1);
00425     }
00426     else {                      /* the aquifer is unconfined */
00427 
00428         /* If the aquifer is unconfied use an explicite scheme to solve
00429          * the nonlinear equation. We use the phead from the first iteration */
00430         z = N_get_array_2d_d_value(data->phead, col, row) -
00431             N_get_array_2d_d_value(data->bottom, col, row);
00432         z_xw = N_get_array_2d_d_value(data->phead, col - 1, row) -
00433             N_get_array_2d_d_value(data->bottom, col - 1, row);
00434         z_xe = N_get_array_2d_d_value(data->phead, col + 1, row) -
00435             N_get_array_2d_d_value(data->bottom, col + 1, row);
00436         z_yn = N_get_array_2d_d_value(data->phead, col, row - 1) -
00437             N_get_array_2d_d_value(data->bottom, col, row - 1);
00438         z_ys = N_get_array_2d_d_value(data->phead, col, row + 1) -
00439             N_get_array_2d_d_value(data->bottom, col, row + 1);
00440     }
00441 
00442     /*geometrical mean of cell height */
00443     if (z_w > 0 || z_w < 0 || z_w == 0)
00444         z_w = N_calc_arith_mean(z_xw, z);
00445     else
00446         z_w = z;
00447     if (z_e > 0 || z_e < 0 || z_e == 0)
00448         z_e = N_calc_arith_mean(z_xe, z);
00449     else
00450         z_e = z;
00451     if (z_n > 0 || z_n < 0 || z_n == 0)
00452         z_n = N_calc_arith_mean(z_yn, z);
00453     else
00454         z_n = z;
00455     if (z_s > 0 || z_s < 0 || z_s == 0)
00456         z_s = N_calc_arith_mean(z_ys, z);
00457     else
00458         z_s = z;
00459 
00460     /* Inner sources */
00461     q = N_get_array_2d_d_value(data->q, col, row);
00462     nf = N_get_array_2d_d_value(data->nf, col, row);
00463 
00464     /* specific yield  of current cell face */
00465     Ss = N_get_array_2d_d_value(data->s, col, row) * Az;
00466     /* recharge */
00467     r = N_get_array_2d_d_value(data->r, col, row) * Az;
00468 
00469     /*get the surrounding permeabilities */
00470     hc_x = N_get_array_2d_d_value(data->hc_x, col, row);
00471     hc_y = N_get_array_2d_d_value(data->hc_y, col, row);
00472     hc_xw = N_get_array_2d_d_value(data->hc_x, col - 1, row);
00473     hc_xe = N_get_array_2d_d_value(data->hc_x, col + 1, row);
00474     hc_yn = N_get_array_2d_d_value(data->hc_y, col, row - 1);
00475     hc_ys = N_get_array_2d_d_value(data->hc_y, col, row + 1);
00476 
00477     /* calculate the transmissivities */
00478     T_w = N_calc_harmonic_mean(hc_xw, hc_x) * z_w;
00479     T_e = N_calc_harmonic_mean(hc_xe, hc_x) * z_e;
00480     T_n = N_calc_harmonic_mean(hc_yn, hc_y) * z_n;
00481     T_s = N_calc_harmonic_mean(hc_ys, hc_y) * z_s;
00482 
00483     /*compute the river leakage, this is an explicit method */
00484     if (data->river_leak &&
00485         (N_get_array_2d_d_value(data->river_leak, col, row) != 0)) {
00486         if (hc > N_get_array_2d_d_value(data->river_bed, col, row)) {
00487             river_vect = N_get_array_2d_d_value(data->river_head, col, row) *
00488                 N_get_array_2d_d_value(data->river_leak, col, row);
00489             river_mat = N_get_array_2d_d_value(data->river_leak, col, row);
00490         }
00491         else if (hc < N_get_array_2d_d_value(data->river_bed, col, row)) {
00492             river_vect = (N_get_array_2d_d_value(data->river_head, col, row) -
00493                           N_get_array_2d_d_value(data->river_bed, col, row))
00494                 * N_get_array_2d_d_value(data->river_leak, col, row);
00495             river_mat = 0;
00496         }
00497     }
00498 
00499     /*compute the drainage, this is an explicit method */
00500     if (data->drain_leak &&
00501         (N_get_array_2d_d_value(data->drain_leak, col, row) != 0)) {
00502         if (hc > N_get_array_2d_d_value(data->drain_bed, col, row)) {
00503             drain_vect = N_get_array_2d_d_value(data->drain_bed, col, row) *
00504                 N_get_array_2d_d_value(data->drain_leak, col, row);
00505             drain_mat = N_get_array_2d_d_value(data->drain_leak, col, row);
00506         }
00507         else if (hc <= N_get_array_2d_d_value(data->drain_bed, col, row)) {
00508             drain_vect = 0;
00509             drain_mat = 0;
00510         }
00511     }
00512 
00513     /*mass balance center cell to western cell */
00514     W = -1 * T_w * dy / dx;
00515     /*mass balance center cell to eastern cell */
00516     E = -1 * T_e * dy / dx;
00517     /*mass balance center cell to northern cell */
00518     N = -1 * T_n * dx / dy;
00519     /*mass balance center cell to southern cell */
00520     S = -1 * T_s * dx / dy;
00521 
00522     /*the diagonal entry of the matrix */
00523     C = -1 * (W + E + N + S - Ss / data->dt - river_mat * Az -
00524               drain_mat * Az);
00525 
00526     /*the entry in the right side b of Ax = b */
00527     V = (q + hc_start * Ss / data->dt) + r + river_vect * Az +
00528         drain_vect * Az;
00529 
00530     G_debug(5, "N_callback_gwflow_2d: called [%i][%i]", row, col);
00531 
00532     /*create the 5 point star entries */
00533     mat_pos = N_create_5star(C, W, E, N, S, V);
00534 
00535     return mat_pos;
00536 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines