Water budged calculatin activated in r.gwflow

Update of scipts and documentation


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

lib/gpde/N_gwflow.c

@@ -612,9 +612,9 @@ N_gwflow_2d_calc_water_budged(N_gwflow_data2d * data, N_geom_data * geom, N_arra
     }
 
     if(fabs(sum) < 0.0000000001)
-        G_message("The total sum of the water balance: %g\n", sum);
+        G_message("The total sum of the water budged: %g\n", sum);
     else
-        G_warning("The total sum of the water balance is significant larger then 0: %g\n", sum);
+        G_warning("The total sum of the water budged is significant larger then 0: %g\n", sum);
 
     return;
 }

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, *type, *dt, *maxit, *error, *solver, *sor,
+	*bottom, *vector_x, *vector_y, *water_balance, *type, *dt, *maxit, *error, *solver, *sor,
 	*river_head, *river_bed, *river_leak, *drain_bed, *drain_leak;
     struct Flag *sparse;
 } paramType;
@@ -39,15 +39,16 @@ typedef struct
 paramType param;		/*Parameters */
 
 /*- prototypes --------------------------------------------------------------*/
-void set_params(void);		/*Fill the paramType structure */
-void copy_result(N_array_2d * status, N_array_2d * phead_start,
+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);
-
-N_les *create_solve_les(N_geom_data * geom, N_gwflow_data2d * data,
+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 ********************************** */
 /* ************************************************************************* */
@@ -127,14 +128,30 @@ void set_params(void)
     param.output->gisprompt = "new,raster,raster";
     param.output->description = _("The map storing the numerical result [m]");
 
-    param.vector = G_define_option();
-    param.vector->key = "velocity";
-    param.vector->type = TYPE_STRING;
-    param.vector->required = NO;
-    param.vector->gisprompt = "new,raster,raster";
-    param.vector->description =
-	_("Calculate the groundwater filter velocity vector field [m/s]\n"
-	  "and write the x, and y components to maps named name_[xy]");
+    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.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 =
+	_("Store the groundwater budged for each cell\n");
 
     param.type = G_define_option();
     param.type->key = "type";
@@ -426,16 +443,25 @@ int main(int argc, char *argv[])
 	    free(tmp_vect);
     }
 
-    /*write the result to the output file */
-    N_write_array_2d_to_rast(data->phead, param.output->answer);
-
     /*release the memory */
     if (les)
 	N_free_les(les);
 
+    /* Compute the water budged for each cell */
+    N_array_2d *budged = N_alloc_array_2d(geom->cols, geom->rows, 1, DCELL_TYPE);
+    N_gwflow_2d_calc_water_budged(data, geom, budged);
+
+    /*write the result to the output file */
+    N_write_array_2d_to_rast(data->phead, param.output->answer);
+
+    /*Write the water balance */
+    if(param.water_balance->answer)
+    {
+	N_write_array_2d_to_rast(budged, param.water_balance->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) {
 	field =
 	    N_compute_gradient_field_2d(data->phead, data->hc_x, data->hc_y,
 					geom, NULL);
@@ -445,10 +471,8 @@ int main(int argc, char *argv[])
 
 	N_compute_gradient_field_components_2d(field, 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);
+	N_write_array_2d_to_rast(xcomp, param.vector_x->answer);
+	N_write_array_2d_to_rast(ycomp, param.vector_y->answer);
 	if (buff)
 	    G_free(buff);
 
@@ -460,7 +484,8 @@ int main(int argc, char *argv[])
 	    N_free_gradient_field_2d(field);
     }
 
-
+    if(budged)
+        N_free_array_2d(budged);
     if (data)
 	N_free_gwflow_data2d(data);
     if (geom)
@@ -471,9 +496,42 @@ int main(int argc, char *argv[])
     return (EXIT_SUCCESS);
 }
 
+/* ************************************************************************* */
+/* 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 from the x vector to a N_array_2d struct  */
+/* this function copies the result into a N_array_2d struct                  */
 /* ************************************************************************* */
 void
 copy_result(N_array_2d * status, N_array_2d * phead_start, double *result,
@@ -510,6 +568,7 @@ copy_result(N_array_2d * status, N_array_2d * phead_start, double *result,
 
     return;
 }
+
 /* *************************************************************** */
 /* ***** create and solve the linear equation system ************* */
 /* *************************************************************** */
@@ -529,24 +588,22 @@ N_les *create_solve_les(N_geom_data * geom, N_gwflow_data2d * data,
     N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead);
 
     /*solve the linear equation system */
-    if(les && les->type == N_NORMAL_LES)
-    {
-    if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_CG) == 0)
-        G_math_solver_cg(les->A, les->x, les->b, les->rows, maxit, error);
+    if(les && les->type == N_NORMAL_LES) {
+        if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_CG) == 0)
+            G_math_solver_cg(les->A, les->x, les->b, les->rows, maxit, error);
 
-    if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_PCG) == 0)
-        G_math_solver_pcg(les->A, les->x, les->b, les->rows, maxit, error, G_MATH_DIAGONAL_PRECONDITION);
+        if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_PCG) == 0)
+            G_math_solver_pcg(les->A, les->x, les->b, les->rows, maxit, error, G_MATH_DIAGONAL_PRECONDITION);
 
-    if (strcmp(solver, G_MATH_SOLVER_DIRECT_CHOLESKY) == 0)
-        G_math_solver_cholesky(les->A, les->x, les->b, les->rows, les->rows);
-    } else if (les && les->type == N_SPARSE_LES)
-    {
-    if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_CG) == 0)
-        G_math_solver_sparse_cg(les->Asp, les->x, les->b, les->rows, maxit, error);
-
-    if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_PCG) == 0)
-        G_math_solver_sparse_pcg(les->Asp, les->x, les->b, les->rows, maxit, error, G_MATH_DIAGONAL_PRECONDITION);
+        if (strcmp(solver, G_MATH_SOLVER_DIRECT_CHOLESKY) == 0)
+            G_math_solver_cholesky(les->A, les->x, les->b, les->rows, les->rows);
+    }
+    else if (les && les->type == N_SPARSE_LES) {
+        if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_CG) == 0)
+            G_math_solver_sparse_cg(les->Asp, les->x, les->b, les->rows, maxit, error);
 
+        if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_PCG) == 0)
+            G_math_solver_sparse_pcg(les->Asp, les->x, les->b, les->rows, maxit, error, G_MATH_DIAGONAL_PRECONDITION);
     }
     if (les == NULL)
         G_fatal_error(_("Unable to create and solve the linear equation system"));

raster/r.gwflow/r.gwflow.html

@@ -16,8 +16,9 @@ raster maps.
 </center>
 <p>
 
-<em>r.gwflow</em> calculates the piezometric head and optionally the 
-filter velocity field, based on the hydraulic conductivity and the piezometric head. 
+<em>r.gwflow</em> calculates the piezometric head and optionally 
+the water budged and the filter velocity field,
+based on the hydraulic conductivity and the piezometric head. 
 The vector components can be visualized with paraview if they are exported
 with <em>r.out.vtk</em>.
 <br>
@@ -26,7 +27,11 @@ The groundwater flow will always be calculated transient.
 If you want to calculate stady state, set the timestep 
 to a large number (billions of seconds) or set the 
 specific yield/ effective porosity raster maps to zero.
-
+<br>
+<br>
+The water budged is calculated for each non inactive cell. The
+sum of the budget for each non inactive cell must be near zero.
+This is an indicator of the quality of the numerical result.
 
 <h2>NOTES</h2>
 
@@ -78,27 +83,27 @@ 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 confined and unconfined aquifers
-r.mapcalc "phead=if(row() == 1 , 50, 40)"
-r.mapcalc "status=if(row() == 1 , 2, 1)"
-r.mapcalc "well=if(row() == 20 && col() == 20 , -0.001, 0)"
-r.mapcalc "hydcond=0.00025"
-r.mapcalc "recharge=0"
-r.mapcalc "top_conf=20.0"
-r.mapcalc "top_unconf=70.0"
-r.mapcalc "bottom=0.0"
-r.mapcalc "null=0.0"
-r.mapcalc "poros=0.15"
-r.mapcalc "syield=0.0001"
+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="hydcond=0.00025"
+r.mapcalc --o expression="recharge=0"
+r.mapcalc --o expression="top_conf=20.0"
+r.mapcalc --o expression="top_unconf=70.0"
+r.mapcalc --o expression="bottom=0.0"
+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 \
 hc_x=hydcond hc_y=hydcond q=well s=syield r=recharge output=gwresult_conf \
-dt=8640000 type=confined velocity=gwresult_conf_velocity
+dt=8640000 type=confined vx=gwresult_conf_velocity_x vy=gwresult_conf_velocity_y
 
 #unconfined groundwater flow with cg solver and sparse matrix
 r.gwflow --o -s 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 velocity=gwresult_unconf_velocity
+output=gwresult_unconf dt=8640000 type=unconfined vx=gwresult_unconf_velocity_x 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.sh

@@ -19,29 +19,16 @@ r.mapcalc --o expression="bottom=0"
 r.mapcalc --o expression="syield=0.0001"
 r.mapcalc --o expression="null=0.0"
 
-#First compute the steady state groundwater flow
-r.gwflow --o solver=cholesky top=top_conf bottom=bottom phead=phead\
+#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 
+ 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=cholesky top=top_conf bottom=bottom phead=gwresult_conf\
+  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
+     r=recharge output=gwresult_conf dt=500 type=confined budged=water_budged
   count=`expr $count + 500`
 done
-
-export GRASS_WIDTH=640
-export GRASS_HEIGHT=480
-
-#export as png and convert into eps and pdf
-export GRASS_TRUECOLOR=TRUE
-export GRASS_PNGFILE=valid_calc_7x7.png
-d.rast gwresult_conf
-d.rast.num gwresult_conf dp=2
-d.barscale at=1,10 &
-echo "Groundwater flow 10.000s" | d.text size=6 color=black
-convert valid_calc_7x7.png valid_calc_7x7.eps
-convert valid_calc_7x7.png valid_calc_7x7.pdf

raster/r.gwflow/valid_calc_excavation.sh

@@ -27,24 +27,4 @@ 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 
-
-# create contour lines
-r.contour input=gwresult output=gwresult_contour step=0.2 --o
-#create flow lines
-r.flow elevin=gwresult flout=gwresult_flow skip=3 --o
-
-export GRASS_WIDTH=640
-export GRASS_HEIGHT=480
-#export as png and convert into eps and pdf
-export GRASS_TRUECOLOR=TRUE
-export GRASS_PNGFILE=Excavation_pit.png
-d.rast gwresult
-d.vect gwresult_flow color=grey
-d.vect gwresult_contour color=black display=attr,shape attrcol=level lsize=16 lcolor=black
-d.legend at=8,12,15,85 map=gwresult 
-d.barscale at=1,10 &
-echo "Groundwater flow steady state" | d.text size=6 color=black
-convert Excavation_pit.png Excavation_pit.eps
-convert Excavation_pit.png Excavation_pit.pdf
-
+ r=recharge output=gwresult dt=864000000000 type=unconfined budged=water_budged