GRASS Programmer's Manual
6.4.2(2012)
Main Page
Related Pages
Namespaces
Data Structures
Files
File List
Globals
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
/* *************************************************************** */
39
int
integration_test_gwflow
(
void
)
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
{
65
N_gwflow_data3d
*
data
;
66
int
i, j, k;
67
68
data =
69
N_alloc_gwflow_data3d
(
TEST_N_NUM_COLS_LOCAL
,
TEST_N_NUM_ROWS_LOCAL
,
70
TEST_N_NUM_DEPTHS_LOCAL
, 1, 1);
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;
109
N_gwflow_data2d
*
data
;
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
151
N_gwflow_data3d
*
data
;
152
N_geom_data
*geom;
153
N_les
*les;
154
N_les_callback_3d
*
call
;
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
171
geom->
depths
=
TEST_N_NUM_DEPTHS_LOCAL
;
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);
196
N_solver_pcg
(les, 100, 0.1e-8,
N_ROWSCALE_EUKLIDNORM_PRECONDITION
);
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);
203
N_solver_pcg
(les, 100, 0.1e-8,
N_ROWSCALE_ABSSUMNORM_PRECONDITION
);
204
N_print_les
(les);
205
N_free_les
(les);
206
207
208
/*CG*/
les =
209
N_assemble_les_3d_dirichlet
(
N_SPARSE_LES
, geom, data->
status
,
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
237
N_solver_pcg
(les, 100, 0.1e-8,
N_ROWSCALE_EUKLIDNORM_PRECONDITION
);
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
245
N_solver_pcg
(les, 100, 0.1e-8,
N_ROWSCALE_ABSSUMNORM_PRECONDITION
);
246
N_print_les
(les);
247
N_free_les
(les);
248
249
250
/*CG*/
les =
251
N_assemble_les_3d_dirichlet
(
N_NORMAL_LES
, geom, data->
status
,
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 =
274
N_assemble_les_3d_dirichlet
(
N_SPARSE_LES
, geom, data->
status
,
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 =
282
N_assemble_les_3d_dirichlet
(
N_NORMAL_LES
, geom, data->
status
,
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 =
305
N_assemble_les_3d_dirichlet
(
N_NORMAL_LES
, geom, data->
status
,
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 =
313
N_assemble_les_3d_dirichlet
(
N_NORMAL_LES
, geom, data->
status
,
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 =
328
N_assemble_les_3d_dirichlet
(
N_NORMAL_LES
, geom, data->
status
,
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
{
345
N_gwflow_data2d
*
data
;
346
N_geom_data
*geom;
347
N_les
*les;
348
N_les_callback_2d
*
call
;
349
350
/*set the callback */
351
call =
N_alloc_les_callback_2d
();
352
N_set_les_callback_2d_func
(call, (*
N_callback_gwflow_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 =
379
N_assemble_les_2d_dirichlet
(
N_SPARSE_LES
, geom, data->
status
,
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 =
395
N_assemble_les_2d_dirichlet
(
N_NORMAL_LES
, geom, data->
status
,
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 =
411
N_assemble_les_2d_dirichlet
(
N_NORMAL_LES
, geom, data->
status
,
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 =
427
N_assemble_les_2d_dirichlet
(
N_SPARSE_LES
, geom, data->
status
,
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 =
443
N_assemble_les_2d_dirichlet
(
N_NORMAL_LES
, geom, data->
status
,
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 =
458
N_assemble_les_2d_dirichlet
(
N_NORMAL_LES
, geom, data->
status
,
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 =
474
N_assemble_les_2d_dirichlet
(
N_NORMAL_LES
, geom, data->
status
,
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 =
489
N_assemble_les_2d_dirichlet
(
N_NORMAL_LES
, geom, data->
status
,
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
}
lib
gpde
test
test_gwflow.c
Generated on Wed Jun 6 2012 14:04:29 for GRASS Programmer's Manual by
1.8.1