1
2 /****************************************************************************
3 *
4 * MODULE: r.solute.transport
5 *
6 * AUTHOR(S): Original author
7 * Soeren Gebbert soerengebbert <at> gmx <dot> de
8 * 27 11 2006 Berlin
9 * PURPOSE: Calculates transient two dimensional solute transport
10 * in porous media
11 *
12 * COPYRIGHT: (C) 2006-2009 by Soeren Gebbert, and the GRASS Development Team
13 *
14 * This program is free software under the GNU General Public
15 * License (>=v2). Read the file COPYING that comes with GRASS
16 * for details.
17 *
18 *****************************************************************************/
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <string.h>
22 #include <math.h>
23 #include <grass/gis.h>
24 #include <grass/raster.h>
25 #include <grass/glocale.h>
26 #include <grass/gmath.h>
27 #include <grass/N_pde.h>
28 #include <grass/N_solute_transport.h>
29
30
31 /*- Parameters and global variables -----------------------------------------*/
32 typedef struct
33 {
34 struct Option *output, *phead, *hc_x, *hc_y,
35 *c, *status, *diff_x, *diff_y, *q, *cs, *r, *top, *nf, *cin,
36 *bottom, *vector_x, *vector_y, *type, *dt, *maxit, *error, *solver, *sor,
37 *al, *at, *loops, *stab;
38 struct Flag *full_les;
39 struct Flag *cfl;
40 } paramType;
41
42 paramType param; /*Parameters */
43
44 /*- prototypes --------------------------------------------------------------*/
45 void set_params(); /*Fill the paramType structure */
46 void copy_result(N_array_2d * status, N_array_2d * c_start, double *result,
47 struct Cell_head *region, N_array_2d * target, int tflag);
48 N_les *create_solve_les(N_geom_data * geom, N_solute_transport_data2d * data,
49 N_les_callback_2d * call, const char *solver, int maxit,
50 double error, double sor);
51
52 /* ************************************************************************* */
53 /* Set up the arguments we are expecting ********************************** */
54 /* ************************************************************************* */
set_params()55 void set_params()
56 {
57 param.c = G_define_standard_option(G_OPT_R_INPUT);
58 param.c->key = "c";
59 param.c->description = _("The initial concentration in [kg/m^3]");
60
61 param.phead = G_define_standard_option(G_OPT_R_INPUT);
62 param.phead->key = "phead";
63 param.phead->description = _("The piezometric head in [m]");
64
65 param.hc_x = G_define_standard_option(G_OPT_R_INPUT);
66 param.hc_x->key = "hc_x";
67 param.hc_x->description =
68 _("The x-part of the hydraulic conductivity tensor in [m/s]");
69
70 param.hc_y = G_define_standard_option(G_OPT_R_INPUT);
71 param.hc_y->key = "hc_y";
72 param.hc_y->description =
73 _("The y-part of the hydraulic conductivity tensor in [m/s]");
74
75
76 param.status = G_define_standard_option(G_OPT_R_INPUT);
77 param.status->key = "status";
78 param.status->description =
79 _("The status for each cell, = 0 - inactive cell, 1 - active cell, "
80 "2 - dirichlet- and 3 - transfer boundary condition");
81
82 param.diff_x = G_define_standard_option(G_OPT_R_INPUT);
83 param.diff_x->key = "diff_x";
84 param.diff_x->description =
85 _("The x-part of the diffusion tensor in [m^2/s]");
86
87 param.diff_y = G_define_standard_option(G_OPT_R_INPUT);
88 param.diff_y->key = "diff_y";
89 param.diff_y->description =
90 _("The y-part of the diffusion tensor in [m^2/s]");
91
92 param.q = G_define_standard_option(G_OPT_R_INPUT);
93 param.q->key = "q";
94 param.q->guisection = _("Water flow");
95 param.q->required = NO;
96 param.q->description = _("Groundwater sources and sinks in [m^3/s]");
97
98 param.cin = G_define_standard_option(G_OPT_R_INPUT);
99 param.cin->key = "cin";
100 param.cin->required = NO;
101 param.cin->gisprompt = "old,raster,raster";
102 param.cin->guisection = _("Water flow");
103 param.cin->description = _("Concentration sources and sinks bounded to a "
104 "water source or sink in [kg/s]");
105
106
107 param.cs = G_define_standard_option(G_OPT_R_INPUT);
108 param.cs->key = "cs";
109 param.cs->type = TYPE_STRING;
110 param.cs->required = YES;
111 param.cs->gisprompt = "old,raster,raster";
112 param.cs->description = _("Concentration of inner sources and inner sinks in [kg/s] "
113 "(i.e. a chemical reaction)");
114
115 param.r = G_define_standard_option(G_OPT_R_INPUT);
116 param.r->key = "rd";
117 param.r->description = _("Retardation factor [-]");
118
119 param.nf = G_define_standard_option(G_OPT_R_INPUT);
120 param.nf->key = "nf";
121 param.nf->description = _("Effective porosity [-]");
122
123 param.top = G_define_standard_option(G_OPT_R_INPUT);
124 param.top->key = "top";
125 param.top->description = _("Top surface of the aquifer in [m]");
126
127 param.bottom = G_define_standard_option(G_OPT_R_INPUT);
128 param.bottom->key = "bottom";
129 param.bottom->description = _("Bottom surface of the aquifer in [m]");
130
131 param.output = G_define_standard_option(G_OPT_R_OUTPUT);
132 param.output->description = _("The resulting concentration of the numerical solute "
133 "transport calculation will be written to this map. [kg/m^3]");
134
135 param.vector_x = G_define_standard_option(G_OPT_R_OUTPUT);
136 param.vector_x->key = "vx";
137 param.vector_x->required = NO;
138 param.vector_x->guisection = _("Water flow");
139 param.vector_x->description =
140 _("Calculate and store the groundwater filter velocity vector part in x direction [m/s]\n");
141
142 param.vector_y = G_define_standard_option(G_OPT_R_OUTPUT);
143 param.vector_y->key = "vy";
144 param.vector_y->required = NO;
145 param.vector_y->guisection = _("Water flow");
146 param.vector_y->description =
147 _("Calculate and store the groundwater filter velocity vector part in y direction [m/s]\n");
148
149 param.dt = N_define_standard_option(N_OPT_CALC_TIME);
150 param.maxit = N_define_standard_option(N_OPT_MAX_ITERATIONS);
151 param.error = N_define_standard_option(N_OPT_ITERATION_ERROR);
152 param.solver = N_define_standard_option(N_OPT_SOLVER_UNSYMM);
153 param.sor = N_define_standard_option(N_OPT_SOR_VALUE);
154
155 param.al = G_define_option();
156 param.al->key = "al";
157 param.al->type = TYPE_DOUBLE;
158 param.al->required = NO;
159 param.al->answer = "0.0";
160 param.al->description =
161 _("The longditudinal dispersivity length. [m]");
162
163 param.at = G_define_option();
164 param.at->key = "at";
165 param.at->type = TYPE_DOUBLE;
166 param.at->required = NO;
167 param.at->answer = "0.0";
168 param.at->description =
169 _("The transversal dispersivity length. [m]");
170
171 param.loops = G_define_option();
172 param.loops->key = "loops";
173 param.loops->type = TYPE_DOUBLE;
174 param.loops->required = NO;
175 param.loops->answer = "1";
176 param.loops->description =
177 _("Use this number of time loops if the CFL flag is off. The timestep will become dt/loops.");
178
179 param.stab = G_define_option();
180 param.stab->key = "stab";
181 param.stab->type = TYPE_STRING;
182 param.stab->required = NO;
183 param.stab->answer = "full";
184 param.stab->options = "full,exp";
185 param.stab->guisection = "Stabelization";
186 param.stab->description =
187 _("Set the flow stabilizing scheme (full or exponential upwinding).");
188
189 param.full_les = G_define_flag();
190 param.full_les->key = 'f';
191 param.full_les->guisection = "Solver";
192 param.full_les->description = _("Use a full filled quadratic linear equation system,"
193 " default is a sparse linear equation system.");
194
195 param.cfl = G_define_flag();
196 param.cfl->key = 'c';
197 param.cfl->guisection = "Stabelization";
198 param.cfl->description =
199 _("Use the Courant-Friedrichs-Lewy criteria for time step calculation");
200 }
201
202 /* ************************************************************************* */
203 /* Main function *********************************************************** */
204 /* ************************************************************************* */
main(int argc,char * argv[])205 int main(int argc, char *argv[])
206 {
207 struct GModule *module = NULL;
208 N_solute_transport_data2d *data = NULL;
209 N_geom_data *geom = NULL;
210 N_les *les = NULL;
211 N_les_callback_2d *call = NULL;
212 struct Cell_head region;
213 double error, sor;
214 char *solver;
215 int x, y, stat, i, maxit = 1;
216 double loops = 1;
217 N_array_2d *xcomp = NULL;
218 N_array_2d *ycomp = NULL;
219 N_array_2d *hc_x = NULL;
220 N_array_2d *hc_y = NULL;
221 N_array_2d *phead = NULL;
222
223 double time_step, cfl, length, time_loops, time_sum;
224
225 /* Initialize GRASS */
226 G_gisinit(argv[0]);
227
228 module = G_define_module();
229 G_add_keyword(_("raster"));
230 G_add_keyword(_("hydrology"));
231 G_add_keyword(_("solute transport"));
232 module->description =
233 _("Numerical calculation program for transient, confined and unconfined "
234 "solute transport in two dimensions");
235
236 /* Get parameters from user */
237 set_params();
238
239 if (G_parser(argc, argv))
240 exit(EXIT_FAILURE);
241
242 /* Make sure that the current projection is not lat/long */
243 if ((G_projection() == PROJECTION_LL))
244 G_fatal_error(_("Lat/Long location is not supported by %s. Please reproject map first."),
245 G_program_name());
246
247 /*Set the maximum iterations */
248 sscanf(param.maxit->answer, "%i", &(maxit));
249 /*Set the calculation error break criteria */
250 sscanf(param.error->answer, "%lf", &(error));
251 sscanf(param.sor->answer, "%lf", &(sor));
252 /*number of loops*/
253 sscanf(param.loops->answer, "%lf", &(loops));
254 /*Set the solver */
255 solver = param.solver->answer;
256
257 if (strcmp(solver, G_MATH_SOLVER_DIRECT_LU) == 0 && !param.full_les->answer)
258 G_fatal_error(_("The direct LU solver do not work with sparse matrices"));
259 if (strcmp(solver, G_MATH_SOLVER_DIRECT_GAUSS) == 0 && !param.full_les->answer)
260 G_fatal_error(_("The direct Gauss solver do not work with sparse matrices"));
261
262
263 /*get the current region */
264 G_get_set_window(®ion);
265
266 /*allocate the geometry structure for geometry and area calculation */
267 geom = N_init_geom_data_2d(®ion, geom);
268
269 /*Set the function callback to the groundwater flow function */
270 call = N_alloc_les_callback_2d();
271 N_set_les_callback_2d_func(call, (*N_callback_solute_transport_2d)); /*solute_transport 2d */
272
273 /*Allocate the groundwater flow data structure */
274 data = N_alloc_solute_transport_data2d(geom->cols, geom->rows);
275
276 /*Set the stabilizing scheme*/
277 if (strncmp("full", param.stab->answer, 4) == 0) {
278 data->stab = N_UPWIND_FULL;
279 }
280 if (strncmp("exp", param.stab->answer, 3) == 0) {
281 data->stab = N_UPWIND_EXP;
282 }
283
284 /*the dispersivity lengths*/
285 sscanf(param.al->answer, "%lf", &(data->al));
286 sscanf(param.at->answer, "%lf", &(data->at));
287
288 /*Set the calculation time */
289 sscanf(param.dt->answer, "%lf", &(data->dt));
290
291 /*read all input maps into the memory and take care of the
292 * null values.*/
293 N_read_rast_to_array_2d(param.c->answer, data->c);
294 N_convert_array_2d_null_to_zero(data->c);
295 N_read_rast_to_array_2d(param.c->answer, data->c_start);
296 N_convert_array_2d_null_to_zero(data->c_start);
297 N_read_rast_to_array_2d(param.status->answer, data->status);
298 N_convert_array_2d_null_to_zero(data->status);
299 N_read_rast_to_array_2d(param.diff_x->answer, data->diff_x);
300 N_convert_array_2d_null_to_zero(data->diff_x);
301 N_read_rast_to_array_2d(param.diff_y->answer, data->diff_y);
302 N_convert_array_2d_null_to_zero(data->diff_y);
303 N_read_rast_to_array_2d(param.q->answer, data->q);
304 N_convert_array_2d_null_to_zero(data->q);
305 N_read_rast_to_array_2d(param.nf->answer, data->nf);
306 N_convert_array_2d_null_to_zero(data->nf);
307 N_read_rast_to_array_2d(param.cs->answer, data->cs);
308 N_convert_array_2d_null_to_zero(data->cs);
309 N_read_rast_to_array_2d(param.top->answer, data->top);
310 N_convert_array_2d_null_to_zero(data->top);
311 N_read_rast_to_array_2d(param.bottom->answer, data->bottom);
312 N_convert_array_2d_null_to_zero(data->bottom);
313 N_read_rast_to_array_2d(param.r->answer, data->R);
314 N_convert_array_2d_null_to_zero(data->R);
315
316 if(param.cin->answer) {
317 N_read_rast_to_array_2d(param.cin->answer, data->cin);
318 N_convert_array_2d_null_to_zero(data->cin);
319 }
320
321 /*initiate the values for velocity calculation*/
322 hc_x = N_alloc_array_2d(geom->cols, geom->rows, 1, DCELL_TYPE);
323 hc_x = N_read_rast_to_array_2d(param.hc_x->answer, hc_x);
324 N_convert_array_2d_null_to_zero(hc_x);
325 hc_y = N_alloc_array_2d(geom->cols, geom->rows, 1, DCELL_TYPE);
326 hc_y = N_read_rast_to_array_2d(param.hc_y->answer, hc_y);
327 N_convert_array_2d_null_to_zero(hc_y);
328 phead = N_alloc_array_2d(geom->cols, geom->rows, 1, DCELL_TYPE);
329 phead = N_read_rast_to_array_2d(param.phead->answer, phead);
330 N_convert_array_2d_null_to_zero(phead);
331
332 /* Set the inactive values to zero, to assure a no flow boundary */
333 for (y = 0; y < geom->rows; y++) {
334 for (x = 0; x < geom->cols; x++) {
335 stat = (int)N_get_array_2d_d_value(data->status, x, y);
336 if (stat == N_CELL_INACTIVE) { /*only inactive cells */
337 N_put_array_2d_d_value(data->diff_x, x, y, 0);
338 N_put_array_2d_d_value(data->diff_y, x, y, 0);
339 N_put_array_2d_d_value(data->cs, x, y, 0);
340 N_put_array_2d_d_value(data->q, x, y, 0);
341 }
342 }
343 }
344
345 /*compute the velocities */
346 N_math_array_2d(hc_x, data->nf, hc_x, N_ARRAY_DIV);
347 N_math_array_2d(hc_y, data->nf, hc_y, N_ARRAY_DIV);
348 N_compute_gradient_field_2d(phead, hc_x, hc_y, geom, data->grad);
349
350 /*Now compute the dispersivity tensor*/
351 N_calc_solute_transport_disptensor_2d(data);
352
353 /***************************************/
354 /*the Courant-Friedrichs-Lewy criteria */
355 /*Compute the correct time step */
356 if (geom->dx > geom->dy)
357 length = geom->dx;
358 else
359 length = geom->dy;
360
361 if (fabs(data->grad->max) > fabs(data->grad->min)) {
362 cfl = (double)data->dt * fabs(data->grad->max) / length;
363 time_step = 1*length / fabs(data->grad->max);
364 }
365 else {
366 cfl = (double)data->dt * fabs(data->grad->min) / length;
367 time_step = 1*length / fabs(data->grad->min);
368 }
369
370 G_message(_("The Courant-Friedrichs-Lewy criteria is %g it should be within [0:1]"), cfl);
371 G_message(_("The largest stable time step is %g"), time_step);
372
373 /*Set the number of inner loops and the time step*/
374 if (data->dt > time_step && param.cfl->answer) {
375 /*safe the user time step */
376 time_sum = data->dt;
377 time_loops = data->dt / time_step;
378 time_loops = floor(time_loops) + 1;
379 data->dt = data->dt / time_loops;
380 G_message(_("Number of inner loops is %g"), time_loops);
381 G_message(_("Time step for each loop %g"), data->dt);
382 }
383 else {
384 if(data->dt > time_step)
385 G_warning(_("The time step is to large: %gs. The largest time step should be of size %gs."), data->dt, time_step);
386
387 time_loops = loops;
388 data->dt = data->dt / loops;
389 }
390
391 N_free_array_2d(phead);
392 N_free_array_2d(hc_x);
393 N_free_array_2d(hc_y);
394
395 /*Compute for each time step*/
396 for (i = 0; i < time_loops; i++) {
397 G_message(_("Time step %i with time sum %g"), i + 1, (i + 1)*data->dt);
398
399 /*assemble the linear equation system and solve it */
400 les = create_solve_les(geom, data, call, solver, maxit, error, sor);
401
402 /* copy the result into the c array for output */
403 copy_result(data->status, data->c_start, les->x, ®ion, data->c, 1);
404 N_convert_array_2d_null_to_zero(data->c_start);
405
406 if (les)
407 N_free_les(les);
408
409 /*Set the start array*/
410 N_copy_array_2d(data->c, data->c_start);
411 /*Set the transmission boundary*/
412 N_calc_solute_transport_transmission_2d(data);
413
414 }
415
416 /*write the result to the output file */
417 N_write_array_2d_to_rast(data->c, param.output->answer);
418
419 /*Compute the the velocity field if required and write the result into three rast maps */
420 if (param.vector_x->answer || param.vector_y->answer) {
421 xcomp = N_alloc_array_2d(geom->cols, geom->rows, 1, DCELL_TYPE);
422 ycomp = N_alloc_array_2d(geom->cols, geom->rows, 1, DCELL_TYPE);
423
424 N_compute_gradient_field_components_2d(data->grad, xcomp, ycomp);
425
426 if (param.vector_x->answer)
427 N_write_array_2d_to_rast(xcomp, param.vector_x->answer);
428 if (param.vector_y->answer)
429 N_write_array_2d_to_rast(ycomp, param.vector_y->answer);
430
431 if (xcomp)
432 N_free_array_2d(xcomp);
433 if (ycomp)
434 N_free_array_2d(ycomp);
435 }
436
437
438 if (data)
439 N_free_solute_transport_data2d(data);
440 if (geom)
441 N_free_geom_data(geom);
442 if (call)
443 G_free(call);
444
445 return (EXIT_SUCCESS);
446 }
447
448
449 /* ************************************************************************* */
450 /* this function copies the result from the x vector to a N_array_2d array * */
451 /* ************************************************************************* */
452 void
copy_result(N_array_2d * status,N_array_2d * c_start,double * result,struct Cell_head * region,N_array_2d * target,int tflag)453 copy_result(N_array_2d * status, N_array_2d * c_start, double *result,
454 struct Cell_head *region, N_array_2d * target, int tflag)
455 {
456 int y, x, rows, cols, count, stat;
457 double d1 = 0;
458 DCELL val;
459
460 rows = region->rows;
461 cols = region->cols;
462
463 count = 0;
464 for (y = 0; y < rows; y++) {
465 G_percent(y, rows - 1, 10);
466 for (x = 0; x < cols; x++) {
467 stat = (int)N_get_array_2d_d_value(status, x, y);
468 if (stat == N_CELL_ACTIVE) { /*only active cells */
469 d1 = result[count];
470 val = (DCELL) d1;
471 count++;
472 }
473 else if (stat == N_CELL_DIRICHLET) { /*dirichlet cells */
474 d1 = N_get_array_2d_d_value(c_start, x, y);
475 val = (DCELL) d1;
476 }
477 else if (tflag == 1 && stat == N_CELL_TRANSMISSION) {/*transmission cells*/
478 d1 = N_get_array_2d_d_value(c_start, x, y);
479 val = (DCELL) d1;
480 }
481 else {
482 Rast_set_null_value(&val, 1, DCELL_TYPE);
483 }
484 N_put_array_2d_d_value(target, x, y, val);
485 }
486 }
487
488 return;
489 }
490
491 /* *************************************************************** */
492 /* ***** create and solve the linear equation system ************* */
493 /* *************************************************************** */
create_solve_les(N_geom_data * geom,N_solute_transport_data2d * data,N_les_callback_2d * call,const char * solver,int maxit,double error,double sor)494 N_les *create_solve_les(N_geom_data * geom, N_solute_transport_data2d * data,
495 N_les_callback_2d * call, const char *solver, int maxit,
496 double error, double sor)
497 {
498
499 N_les *les;
500
501 /*assemble the linear equation system */
502 if (param.full_les->answer)
503 les =
504 N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c,
505 (void *)data, call);
506 else
507 les =
508 N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c,
509 (void *)data, call);
510
511 /*solve the equation system */
512 if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_JACOBI) == 0)
513 {
514 if (!param.full_les->answer)
515 G_math_solver_sparse_jacobi(les->Asp, les->x, les->b, les->rows, maxit, sor, error);
516 else
517 G_math_solver_jacobi(les->A, les->x, les->b, les->rows, maxit, sor, error);
518 }
519
520 if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_SOR) == 0)
521 {
522 if (!param.full_les->answer)
523 G_math_solver_sparse_gs(les->Asp, les->x, les->b, les->rows, maxit, sor, error);
524 else
525 G_math_solver_gs(les->A, les->x, les->b, les->rows, maxit, sor, error);
526 }
527
528 if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_BICGSTAB) == 0)
529 {
530 if (!param.full_les->answer)
531 G_math_solver_sparse_bicgstab(les->Asp, les->x, les->b, les->rows, maxit, error);
532 else
533 G_math_solver_bicgstab(les->A, les->x, les->b, les->rows, maxit, error);
534 }
535
536 if (strcmp(solver, G_MATH_SOLVER_DIRECT_LU) == 0)
537 G_math_solver_lu(les->A, les->x, les->b, les->rows);
538
539 if (strcmp(solver, G_MATH_SOLVER_DIRECT_GAUSS) == 0)
540 G_math_solver_gauss(les->A, les->x, les->b, les->rows);
541
542 if (les == NULL)
543 G_fatal_error(_("Could not create and solve the linear equation system"));
544
545 return les;
546 }
547
548