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: gwflow integration tests 00009 * 00010 * COPYRIGHT: (C) 2000 by the GRASS Development Team 00011 * 00012 * This program is free software under the GNU General Public 00013 * License (>=v2). Read the file COPYING that comes with GRASS 00014 * for details. 00015 * 00016 *****************************************************************************/ 00017 00018 #include <grass/gis.h> 00019 #include <grass/glocale.h> 00020 #include <grass/N_pde.h> 00021 #include <grass/N_gwflow.h> 00022 #include "test_gpde_lib.h" 00023 00024 /*redefine */ 00025 #define TEST_N_NUM_DEPTHS_LOCAL 2 00026 #define TEST_N_NUM_ROWS_LOCAL 3 00027 #define TEST_N_NUM_COLS_LOCAL 3 00028 00029 /* prototypes */ 00030 static N_gwflow_data2d *create_gwflow_data_2d(void); 00031 static N_gwflow_data3d *create_gwflow_data_3d(void); 00032 static int test_gwflow_2d(void); 00033 static int test_gwflow_3d(void); 00034 00035 00036 /* *************************************************************** */ 00037 /* Performe the gwflow integration tests ************************* */ 00038 /* *************************************************************** */ 00039 int integration_test_gwflow(void) 00040 { 00041 int sum = 0; 00042 00043 G_message(_("\n++ Running gwflow integration tests ++")); 00044 00045 G_message(_("\t 1. testing 2d gwflow")); 00046 sum += test_gwflow_2d(); 00047 00048 G_message(_("\t 2. testing 3d gwflow")); 00049 sum += test_gwflow_3d(); 00050 00051 if (sum > 0) 00052 G_warning(_("\n-- gwflow integration tests failure --")); 00053 else 00054 G_message(_("\n-- gwflow integration tests finished successfully --")); 00055 00056 return sum; 00057 } 00058 00059 00060 /* *************************************************************** */ 00061 /* Create valid groundwater flow data **************************** */ 00062 /* *************************************************************** */ 00063 N_gwflow_data3d *create_gwflow_data_3d(void) 00064 { 00065 N_gwflow_data3d *data; 00066 int i, j, k; 00067 00068 data = 00069 N_alloc_gwflow_data3d(TEST_N_NUM_COLS_LOCAL, TEST_N_NUM_ROWS_LOCAL, 00070 TEST_N_NUM_DEPTHS_LOCAL, 1, 1); 00071 00072 #pragma omp parallel for private (i, j, k) shared (data) 00073 for (k = 0; k < TEST_N_NUM_DEPTHS_LOCAL; k++) 00074 for (j = 0; j < TEST_N_NUM_ROWS_LOCAL; j++) { 00075 for (i = 0; i < TEST_N_NUM_COLS_LOCAL; i++) { 00076 00077 00078 if (j == 0) { 00079 N_put_array_3d_d_value(data->phead, i, j, k, 50); 00080 N_put_array_3d_d_value(data->phead_start, i, j, k, 50); 00081 N_put_array_3d_d_value(data->status, i, j, k, 2); 00082 } 00083 else { 00084 00085 N_put_array_3d_d_value(data->phead, i, j, k, 40); 00086 N_put_array_3d_d_value(data->phead_start, i, j, k, 40); 00087 N_put_array_3d_d_value(data->status, i, j, k, 1); 00088 } 00089 N_put_array_3d_d_value(data->hc_x, i, j, k, 0.0001); 00090 N_put_array_3d_d_value(data->hc_y, i, j, k, 0.0001); 00091 N_put_array_3d_d_value(data->hc_z, i, j, k, 0.0001); 00092 N_put_array_3d_d_value(data->q, i, j, k, 0.0); 00093 N_put_array_3d_d_value(data->s, i, j, k, 0.001); 00094 N_put_array_2d_d_value(data->r, i, j, 0.0); 00095 N_put_array_3d_d_value(data->nf, i, j, k, 0.1); 00096 } 00097 } 00098 00099 return data; 00100 } 00101 00102 00103 /* *************************************************************** */ 00104 /* Create valid groundwater flow data **************************** */ 00105 /* *************************************************************** */ 00106 N_gwflow_data2d *create_gwflow_data_2d(void) 00107 { 00108 int i, j; 00109 N_gwflow_data2d *data; 00110 00111 data = 00112 N_alloc_gwflow_data2d(TEST_N_NUM_COLS_LOCAL, TEST_N_NUM_ROWS_LOCAL, 1, 00113 1); 00114 00115 #pragma omp parallel for private (i, j) shared (data) 00116 for (j = 0; j < TEST_N_NUM_ROWS_LOCAL; j++) { 00117 for (i = 0; i < TEST_N_NUM_COLS_LOCAL; i++) { 00118 00119 if (j == 0) { 00120 N_put_array_2d_d_value(data->phead, i, j, 50); 00121 N_put_array_2d_d_value(data->phead_start, i, j, 50); 00122 N_put_array_2d_d_value(data->status, i, j, 2); 00123 } 00124 else { 00125 00126 N_put_array_2d_d_value(data->phead, i, j, 40); 00127 N_put_array_2d_d_value(data->phead_start, i, j, 40); 00128 N_put_array_2d_d_value(data->status, i, j, 1); 00129 } 00130 N_put_array_2d_d_value(data->hc_x, i, j, 30.0001); 00131 N_put_array_2d_d_value(data->hc_y, i, j, 30.0001); 00132 N_put_array_2d_d_value(data->q, i, j, 0.0); 00133 N_put_array_2d_d_value(data->s, i, j, 0.001); 00134 N_put_array_2d_d_value(data->r, i, j, 0.0); 00135 N_put_array_2d_d_value(data->nf, i, j, 0.1); 00136 N_put_array_2d_d_value(data->top, i, j, 20.0); 00137 N_put_array_2d_d_value(data->bottom, i, j, 0.0); 00138 } 00139 } 00140 00141 return data; 00142 } 00143 00144 /* *************************************************************** */ 00145 /* Test the groundwater flow in 3d with different solvers ******** */ 00146 /* *************************************************************** */ 00147 int test_gwflow_3d(void) 00148 { 00149 00150 00151 N_gwflow_data3d *data; 00152 N_geom_data *geom; 00153 N_les *les; 00154 N_les_callback_3d *call; 00155 00156 call = N_alloc_les_callback_3d(); 00157 N_set_les_callback_3d_func(call, (*N_callback_gwflow_3d)); /*gwflow 3d */ 00158 00159 data = create_gwflow_data_3d(); 00160 00161 data->dt = 86400; 00162 00163 geom = N_alloc_geom_data(); 00164 00165 geom->dx = 10; 00166 geom->dy = 15; 00167 geom->dz = 3; 00168 00169 geom->Az = 150; 00170 00171 geom->depths = TEST_N_NUM_DEPTHS_LOCAL; 00172 geom->rows = TEST_N_NUM_ROWS_LOCAL; 00173 geom->cols = TEST_N_NUM_COLS_LOCAL; 00174 00175 00176 /*Assemble the matrix */ 00177 /* 00178 */ 00179 /*CG*/ les = 00180 N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start, 00181 (void *)data, call); 00182 N_solver_cg(les, 100, 0.1e-8); 00183 N_print_les(les); 00184 N_free_les(les); 00185 00186 /*PCG N_DIAGONAL_PRECONDITION */ les = 00187 N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start, 00188 (void *)data, call); 00189 N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION); 00190 N_print_les(les); 00191 N_free_les(les); 00192 00193 /*PCG N_ROWSCALE_EUKLIDNORM_PRECONDITION */ les = 00194 N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start, 00195 (void *)data, call); 00196 N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION); 00197 N_print_les(les); 00198 N_free_les(les); 00199 00200 /*PCG N_ROWSCALE_ABSSUMNORM_PRECONDITION */ les = 00201 N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start, 00202 (void *)data, call); 00203 N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION); 00204 N_print_les(les); 00205 N_free_les(les); 00206 00207 00208 /*CG*/ les = 00209 N_assemble_les_3d_dirichlet(N_SPARSE_LES, geom, data->status, 00210 data->phead_start, (void *)data, call); 00211 N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start); 00212 N_solver_cg(les, 100, 0.1e-8); 00213 N_print_les(les); 00214 N_free_les(les); 00215 00216 00217 /*CG*/ les = 00218 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start, 00219 (void *)data, call); 00220 00221 N_solver_cg(les, 100, 0.1e-8); 00222 N_print_les(les); 00223 N_free_les(les); 00224 00225 /*PCG N_DIAGONAL_PRECONDITION */ les = 00226 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start, 00227 (void *)data, call); 00228 00229 N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION); 00230 N_print_les(les); 00231 N_free_les(les); 00232 00233 /*PCG N_ROWSCALE_EUKLIDNORM_PRECONDITION */ les = 00234 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start, 00235 (void *)data, call); 00236 00237 N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION); 00238 N_print_les(les); 00239 N_free_les(les); 00240 00241 /*PCG N_ROWSCALE_ABSSUMNORM_PRECONDITION */ les = 00242 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start, 00243 (void *)data, call); 00244 00245 N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION); 00246 N_print_les(les); 00247 N_free_les(les); 00248 00249 00250 /*CG*/ les = 00251 N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status, 00252 data->phead_start, (void *)data, call); 00253 N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start); 00254 N_solver_cg(les, 100, 0.1e-8); 00255 N_print_les(les); 00256 N_free_les(les); 00257 00258 00259 /*BICG*/ les = 00260 N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start, 00261 (void *)data, call); 00262 N_solver_bicgstab(les, 100, 0.1e-8); 00263 N_print_les(les); 00264 N_free_les(les); 00265 00266 /*BICG*/ les = 00267 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start, 00268 (void *)data, call); 00269 N_solver_bicgstab(les, 100, 0.1e-8); 00270 N_print_les(les); 00271 N_free_les(les); 00272 00273 /*BICG*/ les = 00274 N_assemble_les_3d_dirichlet(N_SPARSE_LES, geom, data->status, 00275 data->phead_start, (void *)data, call); 00276 N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start); 00277 N_solver_bicgstab(les, 100, 0.1e-8); 00278 N_print_les(les); 00279 N_free_les(les); 00280 00281 /*BICG*/ les = 00282 N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status, 00283 data->phead_start, (void *)data, call); 00284 N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start); 00285 N_solver_bicgstab(les, 100, 0.1e-8); 00286 N_print_les(les); 00287 N_free_les(les); 00288 00289 00290 /*GUASS*/ les = 00291 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start, 00292 (void *)data, call); 00293 N_solver_gauss(les); 00294 N_print_les(les); 00295 N_free_les(les); 00296 00297 /*LU*/ les = 00298 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start, 00299 (void *)data, call); 00300 N_solver_lu(les); 00301 N_print_les(les); 00302 N_free_les(les); 00303 00304 /*GUASS*/ les = 00305 N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status, 00306 data->phead_start, (void *)data, call); 00307 N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start); 00308 N_solver_gauss(les); 00309 N_print_les(les); 00310 N_free_les(les); 00311 00312 /*LU*/ les = 00313 N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status, 00314 data->phead_start, (void *)data, call); 00315 N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start); 00316 N_solver_lu(les); 00317 N_print_les(les); 00318 N_free_les(les); 00319 00320 /*Cholesky */ les = 00321 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start, 00322 (void *)data, call); 00323 N_solver_cholesky(les); 00324 N_print_les(les); 00325 N_free_les(les); 00326 00327 /*Cholesky */ les = 00328 N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status, 00329 data->phead_start, (void *)data, call); 00330 N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start); 00331 N_solver_cholesky(les); 00332 N_print_les(les); 00333 N_free_les(les); 00334 00335 N_free_gwflow_data3d(data); 00336 G_free(geom); 00337 G_free(call); 00338 00339 return 0; 00340 } 00341 00342 /* *************************************************************** */ 00343 int test_gwflow_2d(void) 00344 { 00345 N_gwflow_data2d *data; 00346 N_geom_data *geom; 00347 N_les *les; 00348 N_les_callback_2d *call; 00349 00350 /*set the callback */ 00351 call = N_alloc_les_callback_2d(); 00352 N_set_les_callback_2d_func(call, (*N_callback_gwflow_2d)); 00353 00354 data = create_gwflow_data_2d(); 00355 data->dt = 600; 00356 00357 geom = N_alloc_geom_data(); 00358 00359 geom->dx = 10; 00360 geom->dy = 15; 00361 00362 geom->Az = 150; 00363 00364 geom->rows = TEST_N_NUM_ROWS_LOCAL; 00365 geom->cols = TEST_N_NUM_COLS_LOCAL; 00366 00367 00368 /*Assemble the matrix */ 00369 /* 00370 */ 00371 /*CG*/ les = 00372 N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->phead_start, 00373 (void *)data, call); 00374 N_solver_cg(les, 100, 0.1e-8); 00375 N_print_les(les); 00376 N_free_les(les); 00377 00378 /*CG*/ les = 00379 N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, data->status, 00380 data->phead_start, (void *)data, call); 00381 N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start); 00382 N_solver_cg(les, 100, 0.1e-8); 00383 N_print_les(les); 00384 N_free_les(les); 00385 00386 /*CG*/ les = 00387 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start, 00388 (void *)data, call); 00389 N_solver_cg(les, 100, 0.1e-8); 00390 N_print_les(les); 00391 N_free_les(les); 00392 00393 00394 /*CG*/ les = 00395 N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, 00396 data->phead_start, (void *)data, call); 00397 N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start); 00398 N_solver_cg(les, 100, 0.1e-8); 00399 N_print_les(les); 00400 N_free_les(les); 00401 00402 /*PCG*/ les = 00403 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start, 00404 (void *)data, call); 00405 N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION); 00406 N_print_les(les); 00407 N_free_les(les); 00408 00409 00410 /*PCG*/ les = 00411 N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, 00412 data->phead_start, (void *)data, call); 00413 N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start); 00414 N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION); 00415 N_print_les(les); 00416 N_free_les(les); 00417 00418 00419 /*BICG*/ les = 00420 N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->phead_start, 00421 (void *)data, call); 00422 N_solver_bicgstab(les, 100, 0.1e-8); 00423 N_print_les(les); 00424 N_free_les(les); 00425 00426 /*BICG*/ les = 00427 N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, data->status, 00428 data->phead_start, (void *)data, call); 00429 N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start); 00430 N_solver_bicgstab(les, 100, 0.1e-8); 00431 N_print_les(les); 00432 N_free_les(les); 00433 00434 00435 /*BICG*/ les = 00436 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start, 00437 (void *)data, call); 00438 N_solver_bicgstab(les, 100, 0.1e-8); 00439 N_print_les(les); 00440 N_free_les(les); 00441 00442 /*BICG*/ les = 00443 N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, 00444 data->phead_start, (void *)data, call); 00445 N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start); 00446 N_solver_bicgstab(les, 100, 0.1e-8); 00447 N_print_les(les); 00448 N_free_les(les); 00449 00450 /*GAUSS*/ les = 00451 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start, 00452 (void *)data, call); 00453 N_solver_gauss(les); 00454 N_print_les(les); 00455 N_free_les(les); 00456 00457 /*GAUSS*/ les = 00458 N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, 00459 data->phead_start, (void *)data, call); 00460 N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start); 00461 N_solver_gauss(les); 00462 N_print_les(les); 00463 N_free_les(les); 00464 00465 00466 /*LU*/ les = 00467 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start, 00468 (void *)data, call); 00469 N_solver_lu(les); 00470 N_print_les(les); 00471 N_free_les(les); 00472 00473 /*LU*/ les = 00474 N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, 00475 data->phead_start, (void *)data, call); 00476 N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start); 00477 N_solver_lu(les); 00478 N_print_les(les); 00479 N_free_les(les); 00480 00481 /*Cholesky */ les = 00482 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start, 00483 (void *)data, call); 00484 N_solver_cholesky(les); 00485 N_print_les(les); 00486 N_free_les(les); 00487 00488 /*Cholesky */ les = 00489 N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, 00490 data->phead_start, (void *)data, call); 00491 N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start); 00492 N_solver_cholesky(les); 00493 N_print_les(les); 00494 N_free_les(les); 00495 00496 00497 00498 00499 N_free_gwflow_data2d(data); 00500 G_free(geom); 00501 G_free(call); 00502 00503 return 0; 00504 }