← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to raster/r.solute.transport/main.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.solute.transport/main.c
 
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
/* ************************************************************************* */
 
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
/* ************************************************************************* */
 
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(&region);
 
265
 
 
266
    /*allocate the geometry structure for geometry and area calculation */
 
267
    geom = N_init_geom_data_2d(&region, 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, &region, 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
 
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
/* *************************************************************** */
 
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