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