GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
test_solute_transport.c
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * MODULE: Grass PDE Numerical Library
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> gmx <dot> de
7 *
8 * PURPOSE: solute_transport integration tests
9 *
10 * COPYRIGHT: (C) 2000 by the GRASS Development Team
11 *
12 * This program is free software under the GNU General Public
13 * License (>=v2). Read the file COPYING that comes with GRASS
14 * for details.
15 *
16 *****************************************************************************/
17 
18 #include <grass/gis.h>
19 #include <grass/glocale.h>
20 #include <grass/N_pde.h>
21 #include <grass/N_solute_transport.h>
22 #include "test_gpde_lib.h"
23 
24 /*redefine */
25 #define TEST_N_NUM_DEPTHS_LOCAL 2
26 #define TEST_N_NUM_ROWS_LOCAL 3
27 #define TEST_N_NUM_COLS_LOCAL 3
28 
29 /* prototypes */
30 static N_solute_transport_data2d *create_solute_transport_data_2d(void);
31 static N_solute_transport_data3d *create_solute_transport_data_3d(void);
32 static int test_solute_transport_2d(void);
33 static int test_solute_transport_3d(void);
34 
35 
36 /* *************************************************************** */
37 /* Performe the solute_transport integration tests ************************* */
38 /* *************************************************************** */
40 {
41  int sum = 0;
42 
43  G_message(_("\n++ Running solute_transport integration tests ++"));
44 
45  G_message(_("\t 1. testing 2d solute_transport"));
46  sum += test_solute_transport_2d();
47 
48  G_message(_("\t 2. testing 3d solute_transport"));
49  sum += test_solute_transport_3d();
50 
51  if (sum > 0)
52  G_warning(_("\n-- solute_transport integration tests failure --"));
53  else
54  G_message(_("\n-- solute_transport integration tests finished successfully --"));
55 
56  return sum;
57 }
58 
59 
60 /* *************************************************************** */
61 /* Create valid solute transport data **************************** */
62 /* *************************************************************** */
63 N_solute_transport_data3d *create_solute_transport_data_3d(void)
64 {
66  int i, j, k;
67 
68  data =
72 
73 #pragma omp parallel for private (i, j, k) shared (data)
74  for (k = 0; k < TEST_N_NUM_DEPTHS_LOCAL; k++)
75  for (j = 0; j < TEST_N_NUM_ROWS_LOCAL; j++) {
76  for (i = 0; i < TEST_N_NUM_COLS_LOCAL; i++) {
77 
78 
79  if (j == 0) {
80  N_put_array_3d_d_value(data->c, i, j, k, 1);
81  N_put_array_3d_d_value(data->c_start, i, j, k, 1);
82  N_put_array_3d_d_value(data->status, i, j, k, 3);
83  }
84  else {
85 
86  N_put_array_3d_d_value(data->c, i, j, k, 0);
87  N_put_array_3d_d_value(data->c_start, i, j, k, 0);
88  N_put_array_3d_d_value(data->status, i, j, k, 1);
89  }
90  N_put_array_3d_d_value(data->diff_x, i, j, k, 0.000001);
91  N_put_array_3d_d_value(data->diff_y, i, j, k, 0.000001);
92  N_put_array_3d_d_value(data->diff_z, i, j, k, 0.000001);
93  N_put_array_3d_d_value(data->q, i, j, k, 0.0);
94  N_put_array_3d_d_value(data->cs, i, j, k, 0.0);
95  N_put_array_3d_d_value(data->R, i, j, k, 1.0);
96  N_put_array_3d_d_value(data->nf, i, j, k, 0.1);
97  if (j == 1 && i == 1 && k == 1)
98  N_put_array_3d_d_value(data->cs, i, j, k, 5.0);
99 
100  }
101  }
102 
103  return data;
104 }
105 
106 
107 /* *************************************************************** */
108 /* Create valid solute transport data **************************** */
109 /* *************************************************************** */
110 N_solute_transport_data2d *create_solute_transport_data_2d(void)
111 {
112  int i, j;
114 
115  data =
116  N_alloc_solute_transport_data2d(TEST_N_NUM_COLS_LOCAL,
117  TEST_N_NUM_ROWS_LOCAL);
118 
119 #pragma omp parallel for private (i, j) shared (data)
120  for (j = 0; j < TEST_N_NUM_ROWS_LOCAL; j++) {
121  for (i = 0; i < TEST_N_NUM_COLS_LOCAL; i++) {
122 
123  if (j == 0) {
124  N_put_array_2d_d_value(data->c, i, j, 0);
125  N_put_array_2d_d_value(data->c_start, i, j, 0);
126  N_put_array_2d_d_value(data->status, i, j, 2);
127  }
128  else {
129 
130  N_put_array_2d_d_value(data->c, i, j, 0);
131  N_put_array_2d_d_value(data->c_start, i, j, 0);
132  N_put_array_2d_d_value(data->status, i, j, 1);
133  }
134  N_put_array_2d_d_value(data->diff_x, i, j, 0.000001);
135  N_put_array_2d_d_value(data->diff_y, i, j, 0.000001);
136  N_put_array_2d_d_value(data->cs, i, j, 0.0);
137  N_put_array_2d_d_value(data->R, i, j, 1.0);
138  N_put_array_2d_d_value(data->q, i, j, 0.0);
139  N_put_array_2d_d_value(data->nf, i, j, 0.1);
140  N_put_array_2d_d_value(data->top, i, j, 20.0);
141  N_put_array_2d_d_value(data->bottom, i, j, 0.0);
142  if (j == 1 && i == 1)
143  N_put_array_2d_d_value(data->cs, i, j, 1.0);
144  }
145  }
146  /*dispersivity length */
147  data->al = 0.2;
148  data->at = 0.02;
149 
150 
151 
152 
153 
154  return data;
155 }
156 
157 /* *************************************************************** */
158 /* Test the solute transport in 3d with different solvers ******** */
159 /* *************************************************************** */
160 int test_solute_transport_3d(void)
161 {
163  N_geom_data *geom;
164  N_les *les;
166 
167  call = N_alloc_les_callback_3d();
168  N_set_les_callback_3d_func(call, (*N_callback_solute_transport_3d)); /*solute_transport 3d */
169 
170  data = create_solute_transport_data_3d();
171 
173 
174  data->dt = 86400;
175 
176  geom = N_alloc_geom_data();
177 
178  geom->dx = 10;
179  geom->dy = 15;
180  geom->dz = 3;
181 
182  geom->Az = 150;
183 
185  geom->rows = TEST_N_NUM_ROWS_LOCAL;
186  geom->cols = TEST_N_NUM_COLS_LOCAL;
187  /*Assemble the matrix */
188  /*
189  */
190  /*Jacobi */ les =
191  N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->c_start,
192  (void *)data, call);
193  N_solver_jacobi(les, 100, 1, 0.1e-8);
194  N_print_les(les);
195  N_free_les(les);
196 
197  /*jacobi */ les =
198  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start,
199  (void *)data, call);
200  N_solver_jacobi(les, 100, 1, 0.1e-8);
201  N_print_les(les);
202  N_free_les(les);
203 
204  /*SOR*/ les =
205  N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->c_start,
206  (void *)data, call);
207  N_solver_SOR(les, 100, 1, 0.1e-8);
208  N_print_les(les);
209  N_free_les(les);
210 
211  /*SOR*/ les =
212  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start,
213  (void *)data, call);
214  N_solver_SOR(les, 100, 1, 0.1e-8);
215  N_print_les(les);
216  N_free_les(les);
217 
218  /*BICG*/ les =
219  N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->c_start,
220  (void *)data, call);
221  N_solver_bicgstab(les, 100, 0.1e-8);
222  N_print_les(les);
223  N_free_les(les);
224 
225  /*BICG*/ les =
226  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start,
227  (void *)data, call);
228  N_solver_bicgstab(les, 100, 0.1e-8);
229  N_print_les(les);
230  N_free_les(les);
231 
232  /*GUASS*/ les =
233  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start,
234  (void *)data, call);
235  N_solver_gauss(les);
236  N_print_les(les);
237  N_free_les(les);
238 
239  /*LU*/ les =
240  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start,
241  (void *)data, call);
242  N_solver_lu(les);
243  N_print_les(les);
244  N_free_les(les);
245 
247  G_free(geom);
248  G_free(call);
249 
250  return 0;
251 }
252 
253 /* *************************************************************** */
254 /* Test the solute transport in 2d with different solvers ******** */
255 /* *************************************************************** */
256 int test_solute_transport_2d(void)
257 {
259  N_geom_data *geom = NULL;
260  N_les *les = NULL;
261  N_les_callback_2d *call = NULL;
262  N_array_2d *pot, *relax = NULL;
263  N_gradient_field_2d *field = NULL;
264  int i, j;
265 
266  /*set the callback */
267  call = N_alloc_les_callback_2d();
269 
270  pot =
271  N_alloc_array_2d(TEST_N_NUM_COLS_LOCAL, TEST_N_NUM_ROWS_LOCAL, 1,
272  DCELL_TYPE);
273  relax =
274  N_alloc_array_2d(TEST_N_NUM_COLS_LOCAL, TEST_N_NUM_ROWS_LOCAL, 1,
275  DCELL_TYPE);
276 
277  data = create_solute_transport_data_2d();
278 
279 
280  data->dt = 600;
281 
282  geom = N_alloc_geom_data();
283 
284  geom->dx = 10;
285  geom->dy = 15;
286 
287  geom->Az = 150;
288 
289  geom->rows = TEST_N_NUM_ROWS_LOCAL;
290  geom->cols = TEST_N_NUM_COLS_LOCAL;
291 
292 
293  for (j = 0; j < TEST_N_NUM_ROWS_LOCAL; j++) {
294  for (i = 0; i < TEST_N_NUM_COLS_LOCAL; i++) {
295  N_put_array_2d_d_value(pot, i, j, j);
296  N_put_array_2d_d_value(relax, i, j, 1);
297  }
298  }
299 
300  field = N_compute_gradient_field_2d(pot, relax, relax, geom, NULL);
301  N_copy_gradient_field_2d(field, data->grad);
303 
304  N_compute_gradient_field_2d(pot, relax, relax, geom, data->grad);
305  /*The dispersivity tensor */
307 
308  /*Assemble the matrix */
309  /*
310  */
311  /*Jacobi */ les =
312  N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c_start,
313  (void *)data, call);
314  N_solver_jacobi(les, 100, 1, 0.1e-8);
315  N_print_les(les);
316  N_free_les(les);
317 
318  /*jacobi */ les =
319  N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start,
320  (void *)data, call);
321  N_solver_jacobi(les, 100, 1, 0.1e-8);
322  N_print_les(les);
323  N_free_les(les);
324 
325  /*SOR*/ les =
326  N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c_start,
327  (void *)data, call);
328  N_solver_SOR(les, 100, 1, 0.1e-8);
329  N_print_les(les);
330  N_free_les(les);
331 
332  /*SOR*/ les =
333  N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start,
334  (void *)data, call);
335  N_solver_SOR(les, 100, 1, 0.1e-8);
336  N_print_les(les);
337  N_free_les(les);
338 
339  /*BICG*/ les =
340  N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c_start,
341  (void *)data, call);
342  N_solver_bicgstab(les, 100, 0.1e-8);
343  N_print_les(les);
344  N_free_les(les);
345 
346  /*BICG*/ les =
347  N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start,
348  (void *)data, call);
349  N_solver_bicgstab(les, 100, 0.1e-8);
350  N_print_les(les);
351  N_free_les(les);
352 
353  /*GAUSS*/ les =
354  N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start,
355  (void *)data, call);
356  N_solver_gauss(les);
357  N_print_les(les);
358  N_free_les(les);
359 
360  /*LU*/ les =
361  N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start,
362  (void *)data, call);
363  N_solver_lu(les);
364  N_print_les(les);
365  N_free_les(les);
366 
368  G_free(geom);
369  G_free(call);
370 
371  return 0;
372 }