GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
test_gwflow.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: gwflow 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_gwflow.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_gwflow_data2d *create_gwflow_data_2d(void);
31 static N_gwflow_data3d *create_gwflow_data_3d(void);
32 static int test_gwflow_2d(void);
33 static int test_gwflow_3d(void);
34 
35 
36 /* *************************************************************** */
37 /* Performe the gwflow integration tests ************************* */
38 /* *************************************************************** */
40 {
41  int sum = 0;
42 
43  G_message(_("\n++ Running gwflow integration tests ++"));
44 
45  G_message(_("\t 1. testing 2d gwflow"));
46  sum += test_gwflow_2d();
47 
48  G_message(_("\t 2. testing 3d gwflow"));
49  sum += test_gwflow_3d();
50 
51  if (sum > 0)
52  G_warning(_("\n-- gwflow integration tests failure --"));
53  else
54  G_message(_("\n-- gwflow integration tests finished successfully --"));
55 
56  return sum;
57 }
58 
59 
60 /* *************************************************************** */
61 /* Create valid groundwater flow data **************************** */
62 /* *************************************************************** */
63 N_gwflow_data3d *create_gwflow_data_3d(void)
64 {
66  int i, j, k;
67 
68  data =
71 
72 #pragma omp parallel for private (i, j, k) shared (data)
73  for (k = 0; k < TEST_N_NUM_DEPTHS_LOCAL; k++)
74  for (j = 0; j < TEST_N_NUM_ROWS_LOCAL; j++) {
75  for (i = 0; i < TEST_N_NUM_COLS_LOCAL; i++) {
76 
77 
78  if (j == 0) {
79  N_put_array_3d_d_value(data->phead, i, j, k, 50);
80  N_put_array_3d_d_value(data->phead_start, i, j, k, 50);
81  N_put_array_3d_d_value(data->status, i, j, k, 2);
82  }
83  else {
84 
85  N_put_array_3d_d_value(data->phead, i, j, k, 40);
86  N_put_array_3d_d_value(data->phead_start, i, j, k, 40);
87  N_put_array_3d_d_value(data->status, i, j, k, 1);
88  }
89  N_put_array_3d_d_value(data->hc_x, i, j, k, 0.0001);
90  N_put_array_3d_d_value(data->hc_y, i, j, k, 0.0001);
91  N_put_array_3d_d_value(data->hc_z, i, j, k, 0.0001);
92  N_put_array_3d_d_value(data->q, i, j, k, 0.0);
93  N_put_array_3d_d_value(data->s, i, j, k, 0.001);
94  N_put_array_2d_d_value(data->r, i, j, 0.0);
95  N_put_array_3d_d_value(data->nf, i, j, k, 0.1);
96  }
97  }
98 
99  return data;
100 }
101 
102 
103 /* *************************************************************** */
104 /* Create valid groundwater flow data **************************** */
105 /* *************************************************************** */
106 N_gwflow_data2d *create_gwflow_data_2d(void)
107 {
108  int i, j;
110 
111  data =
112  N_alloc_gwflow_data2d(TEST_N_NUM_COLS_LOCAL, TEST_N_NUM_ROWS_LOCAL, 1,
113  1);
114 
115 #pragma omp parallel for private (i, j) shared (data)
116  for (j = 0; j < TEST_N_NUM_ROWS_LOCAL; j++) {
117  for (i = 0; i < TEST_N_NUM_COLS_LOCAL; i++) {
118 
119  if (j == 0) {
120  N_put_array_2d_d_value(data->phead, i, j, 50);
121  N_put_array_2d_d_value(data->phead_start, i, j, 50);
122  N_put_array_2d_d_value(data->status, i, j, 2);
123  }
124  else {
125 
126  N_put_array_2d_d_value(data->phead, i, j, 40);
127  N_put_array_2d_d_value(data->phead_start, i, j, 40);
128  N_put_array_2d_d_value(data->status, i, j, 1);
129  }
130  N_put_array_2d_d_value(data->hc_x, i, j, 30.0001);
131  N_put_array_2d_d_value(data->hc_y, i, j, 30.0001);
132  N_put_array_2d_d_value(data->q, i, j, 0.0);
133  N_put_array_2d_d_value(data->s, i, j, 0.001);
134  N_put_array_2d_d_value(data->r, i, j, 0.0);
135  N_put_array_2d_d_value(data->nf, i, j, 0.1);
136  N_put_array_2d_d_value(data->top, i, j, 20.0);
137  N_put_array_2d_d_value(data->bottom, i, j, 0.0);
138  }
139  }
140 
141  return data;
142 }
143 
144 /* *************************************************************** */
145 /* Test the groundwater flow in 3d with different solvers ******** */
146 /* *************************************************************** */
147 int test_gwflow_3d(void)
148 {
149 
150 
152  N_geom_data *geom;
153  N_les *les;
155 
156  call = N_alloc_les_callback_3d();
157  N_set_les_callback_3d_func(call, (*N_callback_gwflow_3d)); /*gwflow 3d */
158 
159  data = create_gwflow_data_3d();
160 
161  data->dt = 86400;
162 
163  geom = N_alloc_geom_data();
164 
165  geom->dx = 10;
166  geom->dy = 15;
167  geom->dz = 3;
168 
169  geom->Az = 150;
170 
172  geom->rows = TEST_N_NUM_ROWS_LOCAL;
173  geom->cols = TEST_N_NUM_COLS_LOCAL;
174 
175 
176  /*Assemble the matrix */
177  /*
178  */
179  /*CG*/ les =
180  N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start,
181  (void *)data, call);
182  N_solver_cg(les, 100, 0.1e-8);
183  N_print_les(les);
184  N_free_les(les);
185 
186  /*PCG N_DIAGONAL_PRECONDITION */ les =
187  N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start,
188  (void *)data, call);
189  N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
190  N_print_les(les);
191  N_free_les(les);
192 
193  /*PCG N_ROWSCALE_EUKLIDNORM_PRECONDITION */ les =
194  N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start,
195  (void *)data, call);
197  N_print_les(les);
198  N_free_les(les);
199 
200  /*PCG N_ROWSCALE_ABSSUMNORM_PRECONDITION */ les =
201  N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start,
202  (void *)data, call);
204  N_print_les(les);
205  N_free_les(les);
206 
207 
208  /*CG*/ les =
210  data->phead_start, (void *)data, call);
211  N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
212  N_solver_cg(les, 100, 0.1e-8);
213  N_print_les(les);
214  N_free_les(les);
215 
216 
217  /*CG*/ les =
218  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
219  (void *)data, call);
220 
221  N_solver_cg(les, 100, 0.1e-8);
222  N_print_les(les);
223  N_free_les(les);
224 
225  /*PCG N_DIAGONAL_PRECONDITION */ les =
226  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
227  (void *)data, call);
228 
229  N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
230  N_print_les(les);
231  N_free_les(les);
232 
233  /*PCG N_ROWSCALE_EUKLIDNORM_PRECONDITION */ les =
234  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
235  (void *)data, call);
236 
238  N_print_les(les);
239  N_free_les(les);
240 
241  /*PCG N_ROWSCALE_ABSSUMNORM_PRECONDITION */ les =
242  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
243  (void *)data, call);
244 
246  N_print_les(les);
247  N_free_les(les);
248 
249 
250  /*CG*/ les =
252  data->phead_start, (void *)data, call);
253  N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
254  N_solver_cg(les, 100, 0.1e-8);
255  N_print_les(les);
256  N_free_les(les);
257 
258 
259  /*BICG*/ les =
260  N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start,
261  (void *)data, call);
262  N_solver_bicgstab(les, 100, 0.1e-8);
263  N_print_les(les);
264  N_free_les(les);
265 
266  /*BICG*/ les =
267  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
268  (void *)data, call);
269  N_solver_bicgstab(les, 100, 0.1e-8);
270  N_print_les(les);
271  N_free_les(les);
272 
273  /*BICG*/ les =
275  data->phead_start, (void *)data, call);
276  N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
277  N_solver_bicgstab(les, 100, 0.1e-8);
278  N_print_les(les);
279  N_free_les(les);
280 
281  /*BICG*/ les =
283  data->phead_start, (void *)data, call);
284  N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
285  N_solver_bicgstab(les, 100, 0.1e-8);
286  N_print_les(les);
287  N_free_les(les);
288 
289 
290  /*GUASS*/ les =
291  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
292  (void *)data, call);
293  N_solver_gauss(les);
294  N_print_les(les);
295  N_free_les(les);
296 
297  /*LU*/ les =
298  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
299  (void *)data, call);
300  N_solver_lu(les);
301  N_print_les(les);
302  N_free_les(les);
303 
304  /*GUASS*/ les =
306  data->phead_start, (void *)data, call);
307  N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
308  N_solver_gauss(les);
309  N_print_les(les);
310  N_free_les(les);
311 
312  /*LU*/ les =
314  data->phead_start, (void *)data, call);
315  N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
316  N_solver_lu(les);
317  N_print_les(les);
318  N_free_les(les);
319 
320  /*Cholesky */ les =
321  N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
322  (void *)data, call);
323  N_solver_cholesky(les);
324  N_print_les(les);
325  N_free_les(les);
326 
327  /*Cholesky */ les =
329  data->phead_start, (void *)data, call);
330  N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
331  N_solver_cholesky(les);
332  N_print_les(les);
333  N_free_les(les);
334 
335  N_free_gwflow_data3d(data);
336  G_free(geom);
337  G_free(call);
338 
339  return 0;
340 }
341 
342 /* *************************************************************** */
343 int test_gwflow_2d(void)
344 {
346  N_geom_data *geom;
347  N_les *les;
349 
350  /*set the callback */
351  call = N_alloc_les_callback_2d();
353 
354  data = create_gwflow_data_2d();
355  data->dt = 600;
356 
357  geom = N_alloc_geom_data();
358 
359  geom->dx = 10;
360  geom->dy = 15;
361 
362  geom->Az = 150;
363 
364  geom->rows = TEST_N_NUM_ROWS_LOCAL;
365  geom->cols = TEST_N_NUM_COLS_LOCAL;
366 
367 
368  /*Assemble the matrix */
369  /*
370  */
371  /*CG*/ les =
372  N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->phead_start,
373  (void *)data, call);
374  N_solver_cg(les, 100, 0.1e-8);
375  N_print_les(les);
376  N_free_les(les);
377 
378  /*CG*/ les =
380  data->phead_start, (void *)data, call);
381  N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
382  N_solver_cg(les, 100, 0.1e-8);
383  N_print_les(les);
384  N_free_les(les);
385 
386  /*CG*/ les =
387  N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
388  (void *)data, call);
389  N_solver_cg(les, 100, 0.1e-8);
390  N_print_les(les);
391  N_free_les(les);
392 
393 
394  /*CG*/ les =
396  data->phead_start, (void *)data, call);
397  N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
398  N_solver_cg(les, 100, 0.1e-8);
399  N_print_les(les);
400  N_free_les(les);
401 
402  /*PCG*/ les =
403  N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
404  (void *)data, call);
405  N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
406  N_print_les(les);
407  N_free_les(les);
408 
409 
410  /*PCG*/ les =
412  data->phead_start, (void *)data, call);
413  N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
414  N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
415  N_print_les(les);
416  N_free_les(les);
417 
418 
419  /*BICG*/ les =
420  N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->phead_start,
421  (void *)data, call);
422  N_solver_bicgstab(les, 100, 0.1e-8);
423  N_print_les(les);
424  N_free_les(les);
425 
426  /*BICG*/ les =
428  data->phead_start, (void *)data, call);
429  N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
430  N_solver_bicgstab(les, 100, 0.1e-8);
431  N_print_les(les);
432  N_free_les(les);
433 
434 
435  /*BICG*/ les =
436  N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
437  (void *)data, call);
438  N_solver_bicgstab(les, 100, 0.1e-8);
439  N_print_les(les);
440  N_free_les(les);
441 
442  /*BICG*/ les =
444  data->phead_start, (void *)data, call);
445  N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
446  N_solver_bicgstab(les, 100, 0.1e-8);
447  N_print_les(les);
448  N_free_les(les);
449 
450  /*GAUSS*/ les =
451  N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
452  (void *)data, call);
453  N_solver_gauss(les);
454  N_print_les(les);
455  N_free_les(les);
456 
457  /*GAUSS*/ les =
459  data->phead_start, (void *)data, call);
460  N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
461  N_solver_gauss(les);
462  N_print_les(les);
463  N_free_les(les);
464 
465 
466  /*LU*/ les =
467  N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
468  (void *)data, call);
469  N_solver_lu(les);
470  N_print_les(les);
471  N_free_les(les);
472 
473  /*LU*/ les =
475  data->phead_start, (void *)data, call);
476  N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
477  N_solver_lu(les);
478  N_print_les(les);
479  N_free_les(les);
480 
481  /*Cholesky */ les =
482  N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
483  (void *)data, call);
484  N_solver_cholesky(les);
485  N_print_les(les);
486  N_free_les(les);
487 
488  /*Cholesky */ les =
490  data->phead_start, (void *)data, call);
491  N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
492  N_solver_cholesky(les);
493  N_print_les(les);
494  N_free_les(les);
495 
496 
497 
498 
499  N_free_gwflow_data2d(data);
500  G_free(geom);
501  G_free(call);
502 
503  return 0;
504 }