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: solute_transport 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_solute_transport.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_solute_transport_data2d *create_solute_transport_data_2d(void); 00031 static N_solute_transport_data3d *create_solute_transport_data_3d(void); 00032 static int test_solute_transport_2d(void); 00033 static int test_solute_transport_3d(void); 00034 00035 00036 /* *************************************************************** */ 00037 /* Performe the solute_transport integration tests ************************* */ 00038 /* *************************************************************** */ 00039 int integration_test_solute_transport(void) 00040 { 00041 int sum = 0; 00042 00043 G_message(_("\n++ Running solute_transport integration tests ++")); 00044 00045 G_message(_("\t 1. testing 2d solute_transport")); 00046 sum += test_solute_transport_2d(); 00047 00048 G_message(_("\t 2. testing 3d solute_transport")); 00049 sum += test_solute_transport_3d(); 00050 00051 if (sum > 0) 00052 G_warning(_("\n-- solute_transport integration tests failure --")); 00053 else 00054 G_message(_("\n-- solute_transport integration tests finished successfully --")); 00055 00056 return sum; 00057 } 00058 00059 00060 /* *************************************************************** */ 00061 /* Create valid solute transport data **************************** */ 00062 /* *************************************************************** */ 00063 N_solute_transport_data3d *create_solute_transport_data_3d(void) 00064 { 00065 N_solute_transport_data3d *data; 00066 int i, j, k; 00067 00068 data = 00069 N_alloc_solute_transport_data3d(TEST_N_NUM_COLS_LOCAL, 00070 TEST_N_NUM_ROWS_LOCAL, 00071 TEST_N_NUM_DEPTHS_LOCAL); 00072 00073 #pragma omp parallel for private (i, j, k) shared (data) 00074 for (k = 0; k < TEST_N_NUM_DEPTHS_LOCAL; k++) 00075 for (j = 0; j < TEST_N_NUM_ROWS_LOCAL; j++) { 00076 for (i = 0; i < TEST_N_NUM_COLS_LOCAL; i++) { 00077 00078 00079 if (j == 0) { 00080 N_put_array_3d_d_value(data->c, i, j, k, 1); 00081 N_put_array_3d_d_value(data->c_start, i, j, k, 1); 00082 N_put_array_3d_d_value(data->status, i, j, k, 3); 00083 } 00084 else { 00085 00086 N_put_array_3d_d_value(data->c, i, j, k, 0); 00087 N_put_array_3d_d_value(data->c_start, i, j, k, 0); 00088 N_put_array_3d_d_value(data->status, i, j, k, 1); 00089 } 00090 N_put_array_3d_d_value(data->diff_x, i, j, k, 0.000001); 00091 N_put_array_3d_d_value(data->diff_y, i, j, k, 0.000001); 00092 N_put_array_3d_d_value(data->diff_z, i, j, k, 0.000001); 00093 N_put_array_3d_d_value(data->q, i, j, k, 0.0); 00094 N_put_array_3d_d_value(data->cs, i, j, k, 0.0); 00095 N_put_array_3d_d_value(data->R, i, j, k, 1.0); 00096 N_put_array_3d_d_value(data->nf, i, j, k, 0.1); 00097 if (j == 1 && i == 1 && k == 1) 00098 N_put_array_3d_d_value(data->cs, i, j, k, 5.0); 00099 00100 } 00101 } 00102 00103 return data; 00104 } 00105 00106 00107 /* *************************************************************** */ 00108 /* Create valid solute transport data **************************** */ 00109 /* *************************************************************** */ 00110 N_solute_transport_data2d *create_solute_transport_data_2d(void) 00111 { 00112 int i, j; 00113 N_solute_transport_data2d *data; 00114 00115 data = 00116 N_alloc_solute_transport_data2d(TEST_N_NUM_COLS_LOCAL, 00117 TEST_N_NUM_ROWS_LOCAL); 00118 00119 #pragma omp parallel for private (i, j) shared (data) 00120 for (j = 0; j < TEST_N_NUM_ROWS_LOCAL; j++) { 00121 for (i = 0; i < TEST_N_NUM_COLS_LOCAL; i++) { 00122 00123 if (j == 0) { 00124 N_put_array_2d_d_value(data->c, i, j, 0); 00125 N_put_array_2d_d_value(data->c_start, i, j, 0); 00126 N_put_array_2d_d_value(data->status, i, j, 2); 00127 } 00128 else { 00129 00130 N_put_array_2d_d_value(data->c, i, j, 0); 00131 N_put_array_2d_d_value(data->c_start, i, j, 0); 00132 N_put_array_2d_d_value(data->status, i, j, 1); 00133 } 00134 N_put_array_2d_d_value(data->diff_x, i, j, 0.000001); 00135 N_put_array_2d_d_value(data->diff_y, i, j, 0.000001); 00136 N_put_array_2d_d_value(data->cs, i, j, 0.0); 00137 N_put_array_2d_d_value(data->R, i, j, 1.0); 00138 N_put_array_2d_d_value(data->q, i, j, 0.0); 00139 N_put_array_2d_d_value(data->nf, i, j, 0.1); 00140 N_put_array_2d_d_value(data->top, i, j, 20.0); 00141 N_put_array_2d_d_value(data->bottom, i, j, 0.0); 00142 if (j == 1 && i == 1) 00143 N_put_array_2d_d_value(data->cs, i, j, 1.0); 00144 } 00145 } 00146 /*dispersivity length */ 00147 data->al = 0.2; 00148 data->at = 0.02; 00149 00150 00151 00152 00153 00154 return data; 00155 } 00156 00157 /* *************************************************************** */ 00158 /* Test the solute transport in 3d with different solvers ******** */ 00159 /* *************************************************************** */ 00160 int test_solute_transport_3d(void) 00161 { 00162 N_solute_transport_data3d *data; 00163 N_geom_data *geom; 00164 N_les *les; 00165 N_les_callback_3d *call; 00166 00167 call = N_alloc_les_callback_3d(); 00168 N_set_les_callback_3d_func(call, (*N_callback_solute_transport_3d)); /*solute_transport 3d */ 00169 00170 data = create_solute_transport_data_3d(); 00171 00172 N_calc_solute_transport_disptensor_3d(data); 00173 00174 data->dt = 86400; 00175 00176 geom = N_alloc_geom_data(); 00177 00178 geom->dx = 10; 00179 geom->dy = 15; 00180 geom->dz = 3; 00181 00182 geom->Az = 150; 00183 00184 geom->depths = TEST_N_NUM_DEPTHS_LOCAL; 00185 geom->rows = TEST_N_NUM_ROWS_LOCAL; 00186 geom->cols = TEST_N_NUM_COLS_LOCAL; 00187 /*Assemble the matrix */ 00188 /* 00189 */ 00190 /*Jacobi */ les = 00191 N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->c_start, 00192 (void *)data, call); 00193 N_solver_jacobi(les, 100, 1, 0.1e-8); 00194 N_print_les(les); 00195 N_free_les(les); 00196 00197 /*jacobi */ les = 00198 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start, 00199 (void *)data, call); 00200 N_solver_jacobi(les, 100, 1, 0.1e-8); 00201 N_print_les(les); 00202 N_free_les(les); 00203 00204 /*SOR*/ les = 00205 N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->c_start, 00206 (void *)data, call); 00207 N_solver_SOR(les, 100, 1, 0.1e-8); 00208 N_print_les(les); 00209 N_free_les(les); 00210 00211 /*SOR*/ les = 00212 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start, 00213 (void *)data, call); 00214 N_solver_SOR(les, 100, 1, 0.1e-8); 00215 N_print_les(les); 00216 N_free_les(les); 00217 00218 /*BICG*/ les = 00219 N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->c_start, 00220 (void *)data, call); 00221 N_solver_bicgstab(les, 100, 0.1e-8); 00222 N_print_les(les); 00223 N_free_les(les); 00224 00225 /*BICG*/ les = 00226 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start, 00227 (void *)data, call); 00228 N_solver_bicgstab(les, 100, 0.1e-8); 00229 N_print_les(les); 00230 N_free_les(les); 00231 00232 /*GUASS*/ les = 00233 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start, 00234 (void *)data, call); 00235 N_solver_gauss(les); 00236 N_print_les(les); 00237 N_free_les(les); 00238 00239 /*LU*/ les = 00240 N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start, 00241 (void *)data, call); 00242 N_solver_lu(les); 00243 N_print_les(les); 00244 N_free_les(les); 00245 00246 N_free_solute_transport_data3d(data); 00247 G_free(geom); 00248 G_free(call); 00249 00250 return 0; 00251 } 00252 00253 /* *************************************************************** */ 00254 /* Test the solute transport in 2d with different solvers ******** */ 00255 /* *************************************************************** */ 00256 int test_solute_transport_2d(void) 00257 { 00258 N_solute_transport_data2d *data = NULL; 00259 N_geom_data *geom = NULL; 00260 N_les *les = NULL; 00261 N_les_callback_2d *call = NULL; 00262 N_array_2d *pot, *relax = NULL; 00263 N_gradient_field_2d *field = NULL; 00264 int i, j; 00265 00266 /*set the callback */ 00267 call = N_alloc_les_callback_2d(); 00268 N_set_les_callback_2d_func(call, (*N_callback_solute_transport_2d)); 00269 00270 pot = 00271 N_alloc_array_2d(TEST_N_NUM_COLS_LOCAL, TEST_N_NUM_ROWS_LOCAL, 1, 00272 DCELL_TYPE); 00273 relax = 00274 N_alloc_array_2d(TEST_N_NUM_COLS_LOCAL, TEST_N_NUM_ROWS_LOCAL, 1, 00275 DCELL_TYPE); 00276 00277 data = create_solute_transport_data_2d(); 00278 00279 00280 data->dt = 600; 00281 00282 geom = N_alloc_geom_data(); 00283 00284 geom->dx = 10; 00285 geom->dy = 15; 00286 00287 geom->Az = 150; 00288 00289 geom->rows = TEST_N_NUM_ROWS_LOCAL; 00290 geom->cols = TEST_N_NUM_COLS_LOCAL; 00291 00292 00293 for (j = 0; j < TEST_N_NUM_ROWS_LOCAL; j++) { 00294 for (i = 0; i < TEST_N_NUM_COLS_LOCAL; i++) { 00295 N_put_array_2d_d_value(pot, i, j, j); 00296 N_put_array_2d_d_value(relax, i, j, 1); 00297 } 00298 } 00299 00300 field = N_compute_gradient_field_2d(pot, relax, relax, geom, NULL); 00301 N_copy_gradient_field_2d(field, data->grad); 00302 N_free_gradient_field_2d(field); 00303 00304 N_compute_gradient_field_2d(pot, relax, relax, geom, data->grad); 00305 /*The dispersivity tensor */ 00306 N_calc_solute_transport_disptensor_2d(data); 00307 00308 /*Assemble the matrix */ 00309 /* 00310 */ 00311 /*Jacobi */ les = 00312 N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c_start, 00313 (void *)data, call); 00314 N_solver_jacobi(les, 100, 1, 0.1e-8); 00315 N_print_les(les); 00316 N_free_les(les); 00317 00318 /*jacobi */ les = 00319 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start, 00320 (void *)data, call); 00321 N_solver_jacobi(les, 100, 1, 0.1e-8); 00322 N_print_les(les); 00323 N_free_les(les); 00324 00325 /*SOR*/ les = 00326 N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c_start, 00327 (void *)data, call); 00328 N_solver_SOR(les, 100, 1, 0.1e-8); 00329 N_print_les(les); 00330 N_free_les(les); 00331 00332 /*SOR*/ les = 00333 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start, 00334 (void *)data, call); 00335 N_solver_SOR(les, 100, 1, 0.1e-8); 00336 N_print_les(les); 00337 N_free_les(les); 00338 00339 /*BICG*/ les = 00340 N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c_start, 00341 (void *)data, call); 00342 N_solver_bicgstab(les, 100, 0.1e-8); 00343 N_print_les(les); 00344 N_free_les(les); 00345 00346 /*BICG*/ les = 00347 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start, 00348 (void *)data, call); 00349 N_solver_bicgstab(les, 100, 0.1e-8); 00350 N_print_les(les); 00351 N_free_les(les); 00352 00353 /*GAUSS*/ les = 00354 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start, 00355 (void *)data, call); 00356 N_solver_gauss(les); 00357 N_print_les(les); 00358 N_free_les(les); 00359 00360 /*LU*/ les = 00361 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start, 00362 (void *)data, call); 00363 N_solver_lu(les); 00364 N_print_les(les); 00365 N_free_les(les); 00366 00367 N_free_solute_transport_data2d(data); 00368 G_free(geom); 00369 G_free(call); 00370 00371 return 0; 00372 }