Port of r.gwflow shell scripts to python

Update r.solute.transport examples + bugfixing

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

raster/r.gwflow/main.c

@@ -31,7 +31,7 @@
 typedef struct
 {
     struct Option *output, *phead, *status, *hc_x, *hc_y, *q, *s, *r, *top,
-	*bottom, *vector_x, *vector_y, *water_balance, *type, *dt, *maxit, *error, *solver, *sor,
+	*bottom, *vector_x, *vector_y, *budged, *type, *dt, *maxit, *error, *solver, *sor,
 	*river_head, *river_bed, *river_leak, *drain_bed, *drain_leak;
     struct Flag *sparse;
 } paramType;
@@ -43,12 +43,9 @@ static void set_params(void);		/*Fill the paramType structure */
 static void copy_result(N_array_2d * status, N_array_2d * phead_start,
 		 double *result, struct Cell_head *region,
 		 N_array_2d * target);
-static void calc_water_balance(N_gwflow_data2d * data, N_geom_data * geom, N_array_2d * balance);
 static N_les *create_solve_les(N_geom_data * geom, N_gwflow_data2d * data,
                         N_les_callback_2d * call, const char *solver, int maxit,
                         double error);
-static void copy_water_balance(N_gwflow_data2d * data, double *result,
-	    struct Cell_head *region, N_array_2d * target);
 /* ************************************************************************* */
 /* Set up the arguments we are expecting ********************************** */
 /* ************************************************************************* */
@@ -145,12 +142,12 @@ void set_params(void)
 	_("Calculate and store the groundwater filter velocity vector part in y direction [m/s]\n");
 
 
-    param.water_balance = G_define_option();
-    param.water_balance->key = "budged";
-    param.water_balance->type = TYPE_STRING;
-    param.water_balance->required = NO;
-    param.water_balance->gisprompt = "new,raster,raster";
-    param.water_balance->description =
+    param.budged = G_define_option();
+    param.budged->key = "budged";
+    param.budged->type = TYPE_STRING;
+    param.budged->required = NO;
+    param.budged->gisprompt = "new,raster,raster";
+    param.budged->description =
 	_("Store the groundwater budged for each cell\n");
 
     param.type = G_define_option();
@@ -455,9 +452,9 @@ int main(int argc, char *argv[])
     N_write_array_2d_to_rast(data->phead, param.output->answer);
 
     /*Write the water balance */
-    if(param.water_balance->answer)
+    if(param.budged->answer)
     {
-	N_write_array_2d_to_rast(budged, param.water_balance->answer);
+	N_write_array_2d_to_rast(budged, param.budged->answer);
     }
 
     /*Compute the the velocity field if required and write the result into three rast maps */
@@ -497,40 +494,6 @@ int main(int argc, char *argv[])
 }
 
 /* ************************************************************************* */
-/* this function copies the water balance into a N_array_2d struct           */
-/* ************************************************************************* */
-void
-copy_water_balance(N_gwflow_data2d * data, double *result,
-	    struct Cell_head *region, N_array_2d * target)
-{
-    int y, x, rows, cols, count, stat;
-    double d1 = 0;
-    DCELL val;
-
-    rows = region->rows;
-    cols = region->cols;
-
-    count = 0;
-    for (y = 0; y < rows; y++) {
-	G_percent(y, rows - 1, 10);
-	for (x = 0; x < cols; x++) {
-	    stat = (int)N_get_array_2d_d_value(data->status, x, y);
-	    if (stat == N_CELL_ACTIVE || stat == N_CELL_DIRICHLET) {
-		d1 = result[count];
-		val = (DCELL) d1;
-	    }
-	    else {
-		Rast_set_null_value(&val, 1, DCELL_TYPE);
-	    }
-	    N_put_array_2d_d_value(target, x, y, val);
-            count++;
-	}
-    }
-
-    return;
-}
-
-/* ************************************************************************* */
 /* this function copies the result into a N_array_2d struct                  */
 /* ************************************************************************* */
 void

raster/r.gwflow/valid_calc_7x7.py

@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+# Python script to verify r.gwflow calculation, this calculation is based on
+# the example at page 133 of the following book:
+#	author = "Kinzelbach, W. and Rausch, R.",
+#	title = "Grundwassermodellierung",
+#	publisher = "Gebr{\"u}der Borntraeger (Berlin, Stuttgart)",
+#	year = "1995"
+#
+import sys
+import os
+import grass.script as grass
+
+# Overwrite existing maps
+grass.run_command("g.gisenv", set="OVERWRITE=1")
+
+grass.message(_("Set the region"))
+
+# The area is 2000m x 1000m with a cell size of 25m x 25m
+grass.run_command("g.region", res=100, n=700, s=0, w=0, e=700)
+
+grass.run_command("r.mapcalc", expression="phead=50")
+grass.run_command("r.mapcalc", expression="status=if(col() == 1 || col() == 7 , 2, 1)")
+grass.run_command("r.mapcalc", expression="well=if((row() == 4 && col() == 4), -0.1, 0)")
+grass.run_command("r.mapcalc", expression="hydcond=0.0005")
+grass.run_command("r.mapcalc", expression="recharge=0")
+grass.run_command("r.mapcalc", expression="top_conf=20")
+grass.run_command("r.mapcalc", expression="bottom=0")
+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",\
+ status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="syield",\
+ r="recharge", output="gwresult_conf", dt=5, type="confined", budged="water_budged")
+
+count=0
+# loop over the timesteps
+for i in range(20):
+    grass.run_command("r.gwflow", 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_" + str(count), dt=5, type="confined", budged="water_budged")
+    count += 500
+

raster/r.gwflow/valid_calc_7x7.sh

@@ -1,34 +0,0 @@
-#!/bin/sh
-# Shellscript to verify r.gwflow calculation, this calculation is based on
-# the example at page 133 of the following book:
-#	author = "Kinzelbach, W. and Rausch, R.",
-#	title = "Grundwassermodellierung",
-#	publisher = "Gebr{\"u}der Borntraeger (Berlin, Stuttgart)",
-#	year = "1995"
-#
-# set the region
-g.region res=100 n=700 s=0 w=0 e=700
-
-r.mapcalc --o expression="phead=50"
-r.mapcalc --o expression="status=if(col() == 1 || col() == 7 , 2, 1)"
-r.mapcalc --o expression="well=if((row() == 4 && col() == 4), -0.1, 0)"
-r.mapcalc --o expression="hydcond=0.0005"
-r.mapcalc --o expression="recharge=0"
-r.mapcalc --o expression="top_conf=20"
-r.mapcalc --o expression="bottom=0"
-r.mapcalc --o expression="syield=0.0001"
-r.mapcalc --o expression="null=0.0"
-
-#First compute the initial groundwater flow
-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=500 type=confined budged=water_budged
-
-count=500
-# loop over the timesteps
-while [ `expr $count \< 10000` -eq 1 ] ; do
-  r.gwflow --o solver=cg 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 dt=500 type=confined budged=water_budged
-  count=`expr $count + 500`
-done

raster/r.gwflow/valid_calc_excavation.py

@@ -0,0 +1,39 @@
+#!/usr/bin/env python
+# Shellscript to verify r.gwflow calculation, this calculation is based on
+# the example at page 167 of the following book:
+#	author = "Kinzelbach, W. and Rausch, R.",
+#	title = "Grundwassermodellierung",
+#	publisher = "Gebr{\"u}der Borntraeger (Berlin, Stuttgart)",
+#	year = "1995"
+#
+import sys
+import os
+import grass.script as grass
+
+# Overwrite existing maps
+grass.run_command("g.gisenv", set="OVERWRITE=1")
+
+grass.message(_("Set the region"))
+
+# The area is 2000m x 1000m with a cell size of 25m x 25m
+grass.run_command("g.region", res=50, n=950, s=0, w=0, e=2000)
+
+grass.run_command("r.mapcalc", expression="phead= if(row() == 19, 5, 3)")
+grass.run_command("r.mapcalc", expression="status=if((col() == 1 && row() == 13) ||\
+                     (col() == 1 && row() == 14) ||\
+		     (col() == 2 && row() == 13) ||\
+		     (col() == 2 && row() == 14) ||\
+		     (row() == 19), 2, 1)")
+
+grass.run_command("r.mapcalc", expression="well=0.0")
+grass.run_command("r.mapcalc", expression="hydcond=0.001")
+grass.run_command("r.mapcalc", expression="recharge=0.000000006")
+grass.run_command("r.mapcalc", expression="top=20")
+grass.run_command("r.mapcalc", expression="bottom=0")
+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", \
+ status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="syield", \
+ r="recharge", output="gwresult", dt=864000000000, type="unconfined", budged="water_budged")

raster/r.gwflow/valid_calc_excavation.sh

@@ -1,30 +0,0 @@
-#!/bin/sh
-# Shellscript to verify r.gwflow calculation, this calculation is based on
-# the example at page 167 of the following book:
-#	author = "Kinzelbach, W. and Rausch, R.",
-#	title = "Grundwassermodellierung",
-#	publisher = "Gebr{\"u}der Borntraeger (Berlin, Stuttgart)",
-#	year = "1995"
-#
-# set the region
-g.region res=50 n=950 s=0 w=0 e=2000
-
-r.mapcalc --o expression="phead= if(row() == 19, 5, 3)"
-r.mapcalc --o expression="status=if((col() == 1 && row() == 13) ||\
-                     (col() == 1 && row() == 14) ||\
-		     (col() == 2 && row() == 13) ||\
-		     (col() == 2 && row() == 14) ||\
-		     (row() == 19), 2, 1)"
-
-r.mapcalc --o expression="well=0.0"
-r.mapcalc --o expression="hydcond=0.001"
-r.mapcalc --o expression="recharge=0.000000006"
-r.mapcalc --o expression="top=20"
-r.mapcalc --o expression="bottom=0"
-r.mapcalc --o expression="syield=0.001"
-r.mapcalc --o expression="null=0.0"
-
-#compute a steady state groundwater flow
-r.gwflow --o 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 budged=water_budged

raster/r.solute.transport/example.py

@@ -8,7 +8,7 @@ import grass.script as grass
 # Overwrite existing maps
 grass.run_command("g.gisenv", set="OVERWRITE=1")
 
-grass.message(_("Set the region"))
+grass.message("Set the region")
 
 # The area is 200m x 100m with a cell size of 1m x 1m
 grass.run_command("g.region", res=1, res3=1, t=10, b=0, n=100, s=0, w=0, e=200)
@@ -49,5 +49,5 @@ for dt in range(30):
     grass.run_command("r.solute.transport", "s", solver="bicgstab", top="top_conf",\
     bottom="bottom", phead="gwresult_conf", status="tstatus", hc_x="hydcond", hc_y="hydcond",\
     r="R", cs="cs", q="well", nf="poros", output="stresult_conf_" + str(dt + 1), dt=864000, diff_x="diff",\
-    diff_y="diff", c="stresult_conf_" + str(dt), al=0.1, at=0.01)
+    diff_y="diff", c="stresult_conf_" + str(dt), al=0.1, at=0.01, vx="vx", vy="vy")
 

raster/r.solute.transport/main.c

@@ -33,7 +33,7 @@ typedef struct
 {
     struct Option *output, *phead, *hc_x, *hc_y,
 	*c, *status, *diff_x, *diff_y, *q, *cs, *r, *top, *nf, *cin,
-	*bottom, *vector, *type, *dt, *maxit, *error, *solver, *sor,
+	*bottom, *vector_x, *vector_y, *bugded, *type, *dt, *maxit, *error, *solver, *sor,
 	*al, *at, *loops, *stab;
     struct Flag *sparse;
     struct Flag *cfl;
@@ -130,12 +130,22 @@ void set_params()
     param.output->description =	_("The resulting concentration of the numerical solute "
             "transport calculation will be written to this map. [kg/m^3]");
 
-    param.vector = G_define_standard_option(G_OPT_R_OUTPUT);
-    param.vector->key = "velocity";
-    param.vector->required = NO;
-    param.vector->description =
-            _("Calculate the groundwater distance velocity vector field and "
-            "write the x, and y components to maps named name_(xy), [m/s]");
+
+    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,raster,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,raster,raster";
+    param.vector_y->description =
+	_("Calculate and store the groundwater filter velocity vector part in y direction [m/s]\n");
 
     param.dt = N_define_standard_option(N_OPT_CALC_TIME);
     param.maxit = N_define_standard_option(N_OPT_MAX_ITERATIONS);
@@ -208,7 +218,6 @@ int main(int argc, char *argv[])
     N_array_2d *hc_x = NULL;
     N_array_2d *hc_y = NULL;
     N_array_2d *phead = NULL;
-    char *buff;
 
     double time_step, cfl, length, time_loops, time_sum;
 
@@ -331,8 +340,6 @@ int main(int argc, char *argv[])
     N_math_array_2d(hc_y, data->nf, hc_y, N_ARRAY_DIV);
     N_compute_gradient_field_2d(phead, hc_x, hc_y, geom, data->grad);
 
-    N_print_gradient_field_2d_info(data->grad);
-
     /*Now compute the dispersivity tensor*/
     N_calc_solute_transport_disptensor_2d(data);
 
@@ -403,18 +410,16 @@ int main(int argc, char *argv[])
     N_write_array_2d_to_rast(data->c, param.output->answer);
 
     /*Compute the the velocity field if required and write the result into three rast maps */
-    if (param.vector->answer) {
+    if (param.vector_x->answer || param.vector_y->answer) {
 	xcomp = N_alloc_array_2d(geom->cols, geom->rows, 1, DCELL_TYPE);
 	ycomp = N_alloc_array_2d(geom->cols, geom->rows, 1, DCELL_TYPE);
 
 	N_compute_gradient_field_components_2d(data->grad, xcomp, ycomp);
 
-	G_asprintf(&buff, "%s_x", param.vector->answer);
-	N_write_array_2d_to_rast(xcomp, buff);
-	G_asprintf(&buff, "%s_y", param.vector->answer);
-	N_write_array_2d_to_rast(ycomp, buff);
-	if (buff)
-	    G_free(buff);
+        if (param.vector_x->answer)
+            N_write_array_2d_to_rast(xcomp, param.vector_x->answer);
+        if (param.vector_y->answer)
+            N_write_array_2d_to_rast(ycomp, param.vector_y->answer);
 
 	if (xcomp)
 	    N_free_array_2d(xcomp);

raster/r.solute.transport/seguin_verify.py

@@ -58,7 +58,7 @@ grass.run_command("r.mapcalc", expression="c_1=if(col() == 20 && row() == 20 , 0
 # One inner sources (result of chemical reaction)
 grass.run_command("r.mapcalc", expression="cs_1=if(col() == 20 && row() == 20 , 0.0013888, 0.0)")
 # No pollution by well
-grass.run_command("r.mapcalc", expression="cin_1=0.0)")
+grass.run_command("r.mapcalc", expression="cin_1=0.0")
 # We have a transfer boundary condition
 grass.run_command("r.mapcalc", expression="tstatus_1=if(col() == 1 || col() == 80 , 3, 1)")
 # No diffusion coefficient known for the solution
@@ -74,7 +74,7 @@ AT=10
 grass.run_command("r.solute.transport", "c", "s", error=0.000000000000001, maxit=1000, solver="bicgstab",\
   top="top_conf_1", bottom="bottom_1", phead="gwresult_conf_1", status="tstatus_1", hc_x="hydcond_1",\
   hc_y="hydcond_1", r="R_1", cs="cs_1", q="well_1", nf="poros_1", output="stresult_conf_1", dt=86400000,\
-  diff_x="diff_1", diff_y="diff_1", cin="cin_1", c="c_1", al=AL, at=AT, velocity="stresult_conf_vel_1")
+  diff_x="diff_1", diff_y="diff_1", cin="cin_1", c="c_1", al=AL, at=AT, vx="stresult_conf_vel_1_x", vy="stresult_conf_vel_1_y")
 
 # The second computation uses different porosity for higher groundwater velocity
 
@@ -83,6 +83,6 @@ grass.run_command("r.mapcalc", expression="poros_2=1")
 
 # Compute the solute transport using the above defined dispersivity coefficients for a timestep of 1000d
 grass.run_command("r.solute.transport", "c", "s", error=0.000000000000001, maxit=1000, solver="bicgstab",\
-  top="top_conf_1", bottom="bottom_1", phead="gwresult_conf_2", status="tstatus_1", hc_x="hydcond_1",\
+  top="top_conf_1", bottom="bottom_1", phead="gwresult_conf_1", status="tstatus_1", hc_x="hydcond_1",\
   hc_y="hydcond_1", r="R_1", cs="cs_1", q="well_1", nf="poros_2", output="stresult_conf_2", dt=86400000,\
-  diff_x="diff_1", diff_y="diff_1", cin="cin_1", c="c_1", al=AL, at=AT, velocity="stresult_conf_vel_2")
+  diff_x="diff_1", diff_y="diff_1", cin="cin_1", c="c_1", al=AL, at=AT, vx="stresult_conf_vel_2_x", vy="stresult_conf_vel_2_y")

raster/r.solute.transport/seguin_verify_well.py

@@ -75,7 +75,7 @@ AT=5
 grass.run_command("r.solute.transport", "c", "s", error=0.000000000000001, maxit=1000, solver="bicgstab",\
   top="top_conf_1", bottom="bottom_1", phead="gwresult_conf_1", status="tstatus_1", hc_x="hydcond_1",\
   hc_y="hydcond_1", r="R_1", cs="cs_1", q="well_1", nf="poros_1", output="stresult_conf_1", dt=21600000,\
-  diff_x="diff_1", diff_y="diff_1", cin="cin_1", c="c_1", al=AL, at=AT, velocity="stresult_conf_vel_1")
+  diff_x="diff_1", diff_y="diff_1", cin="cin_1", c="c_1", al=AL, at=AT)
 
 # The second computation uses different dispersivities
 
@@ -86,5 +86,5 @@ AT=1
 grass.run_command("r.solute.transport", "c", "s", error=0.000000000000001, maxit=1000, solver="bicgstab",\
   top="top_conf_1", bottom="bottom_1", phead="gwresult_conf_1", status="tstatus_1", hc_x="hydcond_1",\
   hc_y="hydcond_1", r="R_1", cs="cs_1", q="well_1", nf="poros_1", output="stresult_conf_2", dt=21600000,\
-  diff_x="diff_1", diff_y="diff_1", cin="cin_1", c="c_1", al=AL, at=AT, velocity="stresult_conf_vel_2")
+  diff_x="diff_1", diff_y="diff_1", cin="cin_1", c="c_1", al=AL, at=AT)