GRASS Programmer's Manual  6.4.2(2012)
test_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:      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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines