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: 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 }