[GRASS-SVN] r40538 - in grass/trunk: lib/gpde raster/r.gwflow raster3d/r3.gwflow

svn_grass at osgeo.org svn_grass at osgeo.org
Mon Jan 18 16:11:16 EST 2010


Author: huhabla
Date: 2010-01-18 16:11:12 -0500 (Mon, 18 Jan 2010)
New Revision: 40538

Modified:
   grass/trunk/lib/gpde/N_gwflow.c
   grass/trunk/raster/r.gwflow/main.c
   grass/trunk/raster/r.gwflow/r.gwflow.html
   grass/trunk/raster/r.gwflow/valid_calc_7x7.py
   grass/trunk/raster/r.gwflow/valid_calc_excavation.py
   grass/trunk/raster3d/r3.gwflow/main.c
   grass/trunk/raster3d/r3.gwflow/r3.gwflow.html
Log:
General groundwater flow module update.
Example bugfix and update.
Documentation update.

Modified: grass/trunk/lib/gpde/N_gwflow.c
===================================================================
--- grass/trunk/lib/gpde/N_gwflow.c	2010-01-18 19:39:53 UTC (rev 40537)
+++ grass/trunk/lib/gpde/N_gwflow.c	2010-01-18 21:11:12 UTC (rev 40538)
@@ -364,7 +364,7 @@
  * \param gwdata N_gwflow_data3d *
  * \param geom N_geom_data *
  * \param budget N_array_3d
- * \returnvoid
+ * \return void
  *
  * */
 void
@@ -396,12 +396,12 @@
                     /* 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 @@
                         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 @@
  * \param gwdata N_gwflow_data2d *
  * \param geom N_geom_data *
  * \param budget N_array_2d
- * \returnvoid
+ * \return void
  *
  * */
 void

Modified: grass/trunk/raster/r.gwflow/main.c
===================================================================
--- grass/trunk/raster/r.gwflow/main.c	2010-01-18 19:39:53 UTC (rev 40537)
+++ grass/trunk/raster/r.gwflow/main.c	2010-01-18 21:11:12 UTC (rev 40538)
@@ -33,7 +33,7 @@
     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 @@
     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 @@
 
     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 @@
     /*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 *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);

Modified: grass/trunk/raster/r.gwflow/r.gwflow.html
===================================================================
--- grass/trunk/raster/r.gwflow/r.gwflow.html	2010-01-18 19:39:53 UTC (rev 40537)
+++ grass/trunk/raster/r.gwflow/r.gwflow.html	2010-01-18 21:11:12 UTC (rev 40538)
@@ -67,9 +67,9 @@
 <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 @@
 <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 @@
 #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="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

Modified: grass/trunk/raster/r.gwflow/valid_calc_7x7.py
===================================================================
--- grass/trunk/raster/r.gwflow/valid_calc_7x7.py	2010-01-18 19:39:53 UTC (rev 40537)
+++ grass/trunk/raster/r.gwflow/valid_calc_7x7.py	2010-01-18 21:11:12 UTC (rev 40538)
@@ -29,15 +29,15 @@
 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
 

Modified: grass/trunk/raster/r.gwflow/valid_calc_excavation.py
===================================================================
--- grass/trunk/raster/r.gwflow/valid_calc_excavation.py	2010-01-18 19:39:53 UTC (rev 40537)
+++ grass/trunk/raster/r.gwflow/valid_calc_excavation.py	2010-01-18 21:11:12 UTC (rev 40538)
@@ -34,6 +34,6 @@
 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")

Modified: grass/trunk/raster3d/r3.gwflow/main.c
===================================================================
--- grass/trunk/raster3d/r3.gwflow/main.c	2010-01-18 19:39:53 UTC (rev 40537)
+++ grass/trunk/raster3d/r3.gwflow/main.c	2010-01-18 21:11:12 UTC (rev 40538)
@@ -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,16 +117,38 @@
     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);
     param.error = N_define_standard_option(N_OPT_ITERATION_ERROR);
@@ -137,9 +159,10 @@
     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 @@
 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 @@
     /*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 @@
     }
 
     /*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 @@
 		 &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,15 +346,14 @@
 
 	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);
 	if (ycomp)

Modified: grass/trunk/raster3d/r3.gwflow/r3.gwflow.html
===================================================================
--- grass/trunk/raster3d/r3.gwflow/r3.gwflow.html	2010-01-18 19:39:53 UTC (rev 40537)
+++ grass/trunk/raster3d/r3.gwflow/r3.gwflow.html	2010-01-18 21:11:12 UTC (rev 40538)
@@ -53,11 +53,13 @@
 <br>
 <br>
 The groundwater flow equation can be solved with several 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.
 
-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).
-
 <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.
@@ -67,18 +69,18 @@
 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>



More information about the grass-commit mailing list