General groundwater flow module update.

Example bugfix and update.
Documentation update.

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@40538 15284696-431f-4ddb-bdfa-cd5b030d7da7

lib/gpde/N_gwflow.c

@@ -364,7 +364,7 @@ N_data_star *N_callback_gwflow_3d(void *gwdata, N_geom_data * geom, int col,
  * \param gwdata N_gwflow_data3d *
  * \param geom N_geom_data *
  * \param budget N_array_3d
- * \returnvoid
+ * \return void
  *
  * */
 void
@@ -396,12 +396,12 @@ N_gwflow_3d_calc_water_budget(N_gwflow_data3d * data, N_geom_data * geom, N_arra
                     /* Compute the gradient in each direction pointing from the center */
                     hc = N_get_array_3d_d_value(data->phead, x, y, z);
 
-                    if((int)N_get_array_3d_d_value(data->status, x + 1, y, z ) != N_CELL_INACTIVE) {
-                        h = N_get_array_3d_d_value(data->phead,  x + 1, y, z );
+                    if((int)N_get_array_3d_d_value(data->status, x + 1, y    , z) != N_CELL_INACTIVE) {
+                        h = N_get_array_3d_d_value(data->phead,  x + 1, y    , z);
                         val += dstar->E * (hc - h);
                     }
-                    if((int)N_get_array_3d_d_value(data->status, x - 1, y, z ) != N_CELL_INACTIVE) {
-                        h = N_get_array_3d_d_value(data->phead,  x - 1, y, z);
+                    if((int)N_get_array_3d_d_value(data->status, x - 1, y    , z) != N_CELL_INACTIVE) {
+                        h = N_get_array_3d_d_value(data->phead,  x - 1, y    , z);
                         val += dstar->W * (hc - h);
                     }
                     if((int)N_get_array_3d_d_value(data->status, x    , y + 1, z) != N_CELL_INACTIVE) {
@@ -412,12 +412,12 @@ N_gwflow_3d_calc_water_budget(N_gwflow_data3d * data, N_geom_data * geom, N_arra
                         h = N_get_array_3d_d_value(data->phead,  x    , y - 1, z);
                         val += dstar->N * (hc - h);
                     }
-                    if((int)N_get_array_3d_d_value(data->status, x    , y, z + 1) != N_CELL_INACTIVE) {
-                        h = N_get_array_3d_d_value(data->phead,  x    , y, z + 1);
+                    if((int)N_get_array_3d_d_value(data->status, x    , y    , z + 1) != N_CELL_INACTIVE) {
+                        h = N_get_array_3d_d_value(data->phead,  x    , y    , z + 1);
                         val += dstar->T * (hc - h);
                     }
-                    if((int)N_get_array_3d_d_value(data->status, x    , y, z - 1) != N_CELL_INACTIVE) {
-                        h = N_get_array_3d_d_value(data->phead,  x    , y, z - 1);
+                    if((int)N_get_array_3d_d_value(data->status, x    , y    , z - 1) != N_CELL_INACTIVE) {
+                        h = N_get_array_3d_d_value(data->phead,  x    , y    , z - 1);
                         val += dstar->B * (hc - h);
                     }
                     sum += val;
@@ -652,7 +652,7 @@ N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
  * \param gwdata N_gwflow_data2d *
  * \param geom N_geom_data *
  * \param budget N_array_2d
- * \returnvoid
+ * \return void
  *
  * */
 void

raster/r.gwflow/main.c

@@ -33,7 +33,7 @@ typedef struct
     struct Option *output, *phead, *status, *hc_x, *hc_y, *q, *s, *r, *top,
 	*bottom, *vector_x, *vector_y, *budget, *type, *dt, *maxit, *error, *solver, *sor,
 	*river_head, *river_bed, *river_leak, *drain_bed, *drain_leak;
-    struct Flag *sparse;
+    struct Flag *full_les;
 } paramType;
 
 paramType param;		/*Parameters */
@@ -201,11 +201,12 @@ void set_params(void)
     param.maxit = N_define_standard_option(N_OPT_MAX_ITERATIONS);
     param.error = N_define_standard_option(N_OPT_ITERATION_ERROR);
     param.solver = N_define_standard_option(N_OPT_SOLVER_SYMM);
+    param.solver->options = "cg,pcg,cholesky";
 
-    param.sparse = G_define_flag();
-    param.sparse->key = 's';
-    param.sparse->description =
-	_("Use a sparse matrix, only available with iterative solvers");
+    param.full_les = G_define_flag();
+    param.full_les->key = 'f';
+    param.full_les->description = _("Use a full filled quadratic linear equation system,"
+            " default is a sparse linear equation system.");
 
 }
 
@@ -237,6 +238,7 @@ int main(int argc, char *argv[])
 
     module = G_define_module();
     G_add_keyword(_("raster"));
+    G_add_keyword(_("groundwater flow"));
     module->description =
 	_("Numerical calculation program for transient, confined and unconfined groundwater flow in two dimensions.");
 
@@ -282,12 +284,9 @@ int main(int argc, char *argv[])
     /*set the solver */
     solver = param.solver->answer;
 
-    if (strcmp(solver, G_MATH_SOLVER_DIRECT_LU) == 0 && param.sparse->answer)
-	G_fatal_error(_("The direct LU solver do not work with sparse matrices"));
-    if (strcmp(solver, G_MATH_SOLVER_DIRECT_GAUSS) == 0 && param.sparse->answer)
-	G_fatal_error(_("The direct Gauss solver do not work with sparse matrices"));
-    if (strcmp(solver, G_MATH_SOLVER_DIRECT_CHOLESKY) == 0 && param.sparse->answer)
-	G_fatal_error(_("The direct cholesky solver do not work with sparse matrices"));
+    if (strcmp(solver, G_MATH_SOLVER_DIRECT_CHOLESKY) == 0 && !param.full_les->answer)
+	G_fatal_error(_("The cholesky solver dos not work with sparse matrices.\n"
+                "You may choose a full filled quadratic matrix, flag -f. "));
 
 
     /*get the current region */
@@ -543,7 +542,7 @@ N_les *create_solve_les(N_geom_data * geom, N_gwflow_data2d * data,
     N_les *les;
 
     /*assemble the linear equation system */
-    if (param.sparse->answer)
+    if (!param.full_les->answer)
         les = N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, data->status, data->phead, (void *)data, call);
     else
         les = N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead, (void *)data, call);

raster/r.gwflow/r.gwflow.html

@@ -67,9 +67,9 @@ the following states are supportet:
 <br>
 <br>
 The groundwater flow equation can be solved with several solvers.
-Two iterative solvers with sparse and quadratic matrices support are implemented.
-The conjugate gradients (cg) method and the biconjugate gradients-stabilized (bicgstab) method. 
-Aditionally a direct Gauss solver and LU solver are available. Those direct solvers
+An iterative solvers with sparse and quadratic matrices support is implemented.
+The conjugate gradients method with (pcg) and without (cg) precondition.
+Aditionally a direct Cholesky solver is available. This direct solver
 only work with normal quadratic matrices, so be careful using them with large maps 
 (maps of size 10.000 cells will need more than one gigabyte of RAM).
 Always prefer a sparse matrix solver.
@@ -77,6 +77,7 @@ Always prefer a sparse matrix solver.
 <h2>EXAMPLE</h2>
 Use this small script to create a working
 groundwater flow area and data. Make sure you are not in a lat/lon projection.
+It includes drainage and river input as well.
 
 <div class="code"><pre>
 # set the region accordingly
@@ -85,7 +86,7 @@ g.region res=25 res3=25 t=100 b=0 n=1000 s=0 w=0 e=1000
 #now create the input raster maps for confined and unconfined aquifers
 r.mapcalc --o expression="phead=if(row() == 1 , 50, 40)"
 r.mapcalc --o expression="status=if(row() == 1 , 2, 1)"
-r.mapcalc --o expression="well=if(row() == 20 && col() == 20 , -0.001, 0)"
+r.mapcalc --o expression="well=if(row() == 20 && col() == 20 , -0.01, 0)"
 r.mapcalc --o expression="hydcond=0.00025"
 r.mapcalc --o expression="recharge=0"
 r.mapcalc --o expression="top_conf=20.0"
@@ -95,15 +96,28 @@ r.mapcalc --o expression="null=0.0"
 r.mapcalc --o expression="poros=0.15"
 r.mapcalc --o expression="syield=0.0001"
 
-#confined groundwater flow with cg solver and sparse matrix
-r.gwflow --o -s solver=cg top=top_conf bottom=bottom phead=phead status=status \
+# The maps of the river
+r.mapcalc --o expression="river_bed=if(col() == 35 , 48, null())"
+r.mapcalc --o expression="river_head=if(col() == 35 , 49, null())"
+r.mapcalc --o expression="river_leak=if(col() == 35 , 0.0001, null())"
+
+# The maps of the drainage
+r.mapcalc --o expression="drain_bed=if(col() == 5 , 48, null())"
+r.mapcalc --o expression="drain_leak=if(col() == 5 , 0.01, null())"
+
+#confined groundwater flow with cg solver and sparse matrix, river and drain
+#do not work with this confined aquifer (top == 20m)
+r.gwflow --o solver=cg top=top_conf bottom=bottom phead=phead status=status \
 hc_x=hydcond hc_y=hydcond q=well s=syield r=recharge output=gwresult_conf \
-dt=8640000 type=confined vx=gwresult_conf_velocity_x vy=gwresult_conf_velocity_y
+dt=8640000 type=confined vx=gwresult_conf_velocity_x vy=gwresult_conf_velocity_y budget=budget_conf
 
-#unconfined groundwater flow with cg solver and sparse matrix
-r.gwflow --o -s solver=cg top=top_unconf bottom=bottom phead=phead \
+#unconfined groundwater flow with cg solver and sparse matrix, river and drain are enabled
+r.gwflow --o solver=cg top=top_unconf bottom=bottom phead=phead \
 status=status hc_x=hydcond hc_y=hydcond q=well s=poros r=recharge \
-output=gwresult_unconf dt=8640000 type=unconfined vx=gwresult_unconf_velocity_x vy=gwresult_unconf_velocity_y
+river_bed=river_bed river_head=river_head river_leak=river_leak \
+drain_bed=drain_bed drain_leak=drain_leak \
+output=gwresult_unconf dt=8640000 type=unconfined vx=gwresult_unconf_velocity_x \
+budget=budget_unconf vy=gwresult_unconf_velocity_y
 
 # The data can be visulaized with paraview when exported with r.out.vtk
 r.out.vtk -p in=gwresult_conf,status vector=gwresult_conf_velocity_x,gwresult_conf_velocity_y,null out=/tmp/gwdata_conf2d.vtk

raster/r.gwflow/valid_calc_7x7.py

@@ -29,15 +29,15 @@ grass.run_command("r.mapcalc", expression="syield=0.0001")
 grass.run_command("r.mapcalc", expression="null=0.0")
 
 #First compute the initial groundwater flow
-grass.run_command("r.gwflow", solver="cg", top="top_conf", bottom="bottom", phead="phead",\
+grass.run_command("r.gwflow", "f", solver="cholesky", top="top_conf", bottom="bottom", phead="phead",\
  status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="syield",\
- r="recharge", output="gwresult_conf", dt=5, type="confined", budget="water_budget")
+ r="recharge", output="gwresult_conf", dt=500, type="confined", budget="water_budget")
 
-count=0
+count=500
 # loop over the timesteps
 for i in range(20):
-    grass.run_command("r.gwflow", solver="cg", top="top_conf", bottom="bottom", phead="phead",\
+    grass.run_command("r.gwflow", "f", solver="cholesky", top="top_conf", bottom="bottom", phead="gwresult_conf",\
      status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="syield",\
-     r="recharge", output="gwresult_conf_" + str(count), dt=5, type="confined", budget="water_budget")
+     r="recharge", output="gwresult_conf", dt=500, type="confined", budget="water_budget")
     count += 500
 

raster/r.gwflow/valid_calc_excavation.py

@@ -34,6 +34,6 @@ grass.run_command("r.mapcalc", expression="syield=0.001")
 grass.run_command("r.mapcalc", expression="null=0.0")
 
 #compute a steady state groundwater flow
-grass.run_command("r.gwflow", solver="cholesky", top="top", bottom="bottom", phead="phead", \
+grass.run_command("r.gwflow", "f", solver="cholesky", top="top", bottom="bottom", phead="phead", \
  status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="syield", \
  r="recharge", output="gwresult", dt=864000000000, type="unconfined", budget="water_budget")

raster3d/r3.gwflow/main.c

@@ -30,9 +30,9 @@
 typedef struct
 {
     struct Option *output, *phead, *status, *hc_x, *hc_y, *hc_z, *q, *s, *r,
-	*vector, *dt, *maxit, *error, *solver;
+	*vector_x, *vector_y, *vector_z, *budget, *dt, *maxit, *error, *solver;
     struct Flag *mask;
-    struct Flag *sparse;
+    struct Flag *full_les;
 } paramType;
 
 paramType param;		/*Parameters */
@@ -117,15 +117,37 @@ void set_params(void)
     param.output->gisprompt = "new,grid3,3d-raster";
     param.output->description = _("The piezometric head result of the numerical calculation will be written to this map");
 
-    param.vector = G_define_option();
-    param.vector->key = "velocity";
-    param.vector->type = TYPE_STRING;
-    param.vector->required = NO;
-    param.vector->gisprompt = "new,grid3,3d-raster";
-    param.vector->description = _("Calculate the groundwater distance velocity vector field \n"
-	                          "and write the x, y, and z components to maps named name_[xyz].\n"
-	                          "Name is basename for the new raster3d maps");
-
+    param.vector_x = G_define_option();
+    param.vector_x->key = "vx";
+    param.vector_x->type = TYPE_STRING;
+    param.vector_x->required = NO;
+    param.vector_x->gisprompt = "new,grid3,3d-raster";
+    param.vector_x->description =
+	_("Calculate and store the groundwater filter velocity vector part in x direction [m/s]\n");
+
+    param.vector_y = G_define_option();
+    param.vector_y->key = "vy";
+    param.vector_y->type = TYPE_STRING;
+    param.vector_y->required = NO;
+    param.vector_y->gisprompt = "new,grid3,3d-raster";
+    param.vector_y->description =
+	_("Calculate and store the groundwater filter velocity vector part in y direction [m/s]\n");
+
+    param.vector_z = G_define_option();
+    param.vector_z->key = "vz";
+    param.vector_z->type = TYPE_STRING;
+    param.vector_z->required = NO;
+    param.vector_z->gisprompt = "new,grid3,3d-raster";
+    param.vector_z->description =
+	_("Calculate and store the groundwater filter velocity vector part in y direction [m/s]\n");
+
+    param.budget = G_define_option();
+    param.budget->key = "budget";
+    param.budget->type = TYPE_STRING;
+    param.budget->required = NO;
+    param.budget->gisprompt = "new,grid3,3d-raster";
+    param.budget->description =
+	_("Store the groundwater budget for each cell\n");
 
     param.dt = N_define_standard_option(N_OPT_CALC_TIME);
     param.maxit = N_define_standard_option(N_OPT_MAX_ITERATIONS);
@@ -137,9 +159,10 @@ void set_params(void)
     param.mask->key = 'm';
     param.mask->description = _("Use G3D mask (if exists)");
 
-    param.sparse = G_define_flag();
-    param.sparse->key = 's';
-    param.sparse->description = _("Use a sparse linear equation system, only available with iterative solvers");
+    param.full_les = G_define_flag();
+    param.full_les->key = 'f';
+    param.full_les->description = _("Use a full filled quadratic linear equation system,"
+            " default is a sparse linear equation system.");
 }
 
 /* ************************************************************************* */
@@ -148,41 +171,29 @@ void set_params(void)
 int main(int argc, char *argv[])
 {
     struct GModule *module = NULL;
-
     N_gwflow_data3d *data = NULL;
-
     N_geom_data *geom = NULL;
-
     N_les *les = NULL;
-
     N_les_callback_3d *call = NULL;
-
     G3D_Region region;
-
     N_gradient_field_3d *field = NULL;
-
     N_array_3d *xcomp = NULL;
-
     N_array_3d *ycomp = NULL;
-
     N_array_3d *zcomp = NULL;
-
     double error;
-
     int maxit;
-
     const char *solver;
-
     int x, y, z, stat;
 
-    char *buff = NULL;
-
     /* Initialize GRASS */
     G_gisinit(argv[0]);
 
     module = G_define_module();
     G_add_keyword(_("raster3d"));
     G_add_keyword(_("voxel"));
+    G_add_keyword(_("groundwater"));
+    G_add_keyword(_("numeric"));
+    G_add_keyword(_("simulation"));
     module->description = _("Numerical calculation program for transient, confined groundwater flow in three dimensions");
 
     /* Get parameters from user */
@@ -199,8 +210,9 @@ int main(int argc, char *argv[])
     /*Set the solver */
     solver = param.solver->answer;
 
-    if (strcmp(solver, G_MATH_SOLVER_DIRECT_CHOLESKY) == 0 && param.sparse->answer)
-	G_fatal_error(_("The direct cholesky solver do not work with sparse matrices"));
+    if (strcmp(solver, G_MATH_SOLVER_DIRECT_CHOLESKY) == 0 && !param.full_les->answer)
+	G_fatal_error(_("The cholesky solver dos not work with sparse matrices.\n"
+                "You may choose a full filled quadratic matrix, flag -f. "));
 
 
 
@@ -264,7 +276,7 @@ int main(int argc, char *argv[])
     }
 
     /*assemble the linear equation system */
-    if (param.sparse->answer) {
+    if (!param.full_les->answer) {
 	les =
 	    N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead,
 			      (void *)data, call);
@@ -306,8 +318,18 @@ int main(int argc, char *argv[])
 		 &region, param.output->answer);
     N_free_les(les);
 
+    /* Compute the water budget for each cell */
+    N_array_3d *budget = N_alloc_array_3d(geom->cols, geom->rows, geom->depths, 1, DCELL_TYPE);
+    N_gwflow_3d_calc_water_budget(data, geom, budget);
+    
+    /*Write the water balance */
+    if(param.budget->answer)
+    {
+	N_write_array_3d_to_rast3d(budget, param.budget->answer, 1);
+    }
+
     /*Compute the the velocity field if required and write the result into three rast3d maps */
-    if (param.vector->answer) {
+    if (param.vector_x->answer || param.vector_y->answer || param.vector_z->answer) {
 	field =
 	    N_compute_gradient_field_3d(data->phead, data->hc_x, data->hc_y,
 					data->hc_z, geom, NULL);
@@ -324,14 +346,13 @@ int main(int argc, char *argv[])
 
 	N_compute_gradient_field_components_3d(field, xcomp, ycomp, zcomp);
 
-	G_asprintf(&buff, "%s_x", param.vector->answer);
-	N_write_array_3d_to_rast3d(xcomp, buff, 1);
-	G_asprintf(&buff, "%s_y", param.vector->answer);
-	N_write_array_3d_to_rast3d(ycomp, buff, 1);
-	G_asprintf(&buff, "%s_z", param.vector->answer);
-	N_write_array_3d_to_rast3d(zcomp, buff, 1);
-	if (buff)
-	    G_free(buff);
+
+        if(param.vector_x->answer)
+            N_write_array_3d_to_rast3d(xcomp, param.vector_x->answer, 1);
+        if(param.vector_y->answer)
+            N_write_array_3d_to_rast3d(ycomp, param.vector_y->answer, 1);
+        if(param.vector_z->answer)
+            N_write_array_3d_to_rast3d(zcomp, param.vector_z->answer, 1);
 
 	if (xcomp)
 	    N_free_array_3d(xcomp);

raster3d/r3.gwflow/r3.gwflow.html

@@ -53,10 +53,12 @@ the following cell states are supportet:
 <br>
 <br>
 The groundwater flow equation can be solved with several solvers.
-
-Aditionally a direct Gauss solver and LU solver are available. Those direct solvers
-only work with quadratic matrices, so be careful using them with large maps 
-(maps of size 10.000 cells will need more than one gigabyte of ram).
+An iterative solvers with sparse and quadratic matrices support is implemented.
+The conjugate gradients method with (pcg) and without (cg) precondition.
+Aditionally a direct Cholesky solver is available. This direct solver
+only work with normal quadratic matrices, so be careful using them with large maps
+(maps of size 10.000 cells will need more than one gigabyte of RAM).
+Always prefer a sparse matrix solver.
 
 <H2>EXAMPLE</H2>
 Use this small script to create a working
@@ -67,18 +69,18 @@ groundwater flow area and data. Make sure you are not in a lat/lon projection.
 g.region res=25 res3=25 t=100 b=0 n=1000 s=0 w=0 e=1000
 
 #now create the input raster maps for a confined aquifer
-r3.mapcalc "phead = if(row() == 1 && depth() == 4, 50, 40)"
-r3.mapcalc "status = if(row() == 1 && depth() == 4, 2, 1)"
-r3.mapcalc "well = if(row() == 20 && col() == 20 , -0.00025, 0)"
-r3.mapcalc "hydcond = 0.00025"
-r3.mapcalc "syield = 0.0001"
-r.mapcalc  "recharge = 0.0"
+r3.mapcalc --o expression="phead = if(row() == 1 && depth() == 4, 50, 40)"
+r3.mapcalc --o expression="status = if(row() == 1 && depth() == 4, 2, 1)"
+r3.mapcalc --o expression="well = if(row() == 20 && col() == 20 && depth() == 2, -0.25, 0)"
+r3.mapcalc --o expression="hydcond = 0.00025"
+r3.mapcalc --o expression="syield = 0.0001"
+r.mapcalc  --o expression="recharge = 0.0"
 
-r3.gwflow --o -s solver=cg phead=phead status=status hc_x=hydcond hc_y=hydcond  \
-hc_z=hydcond q=well s=syield r=recharge output=gwresult dt=8640000 velocity=gwresult_velocity
+r3.gwflow --o solver=cg phead=phead status=status hc_x=hydcond hc_y=hydcond  \
+hc_z=hydcond q=well s=syield r=recharge output=gwresult dt=8640000 vx=vx vy=vy vz=vz budget=budget
 
 # The data can be visulaized with paraview when exported with r3.out.vtk
-r3.out.vtk -p in=gwresult,status vector=gwresult_velocity_x,gwresult_velocity_y,gwresult_velocity_z out=/tmp/gwdata3d.vtk
+r3.out.vtk -p in=gwresult,status,budget vector=vx,vy,vz out=/tmp/gwdata3d.vtk
 
 #now load the data into paraview
 </pre></div>