[GRASS-SVN] r40301 - in grass/trunk: lib/gpde raster/r.gwflow
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jan 7 09:30:28 EST 2010
Author: huhabla
Date: 2010-01-07 09:30:25 -0500 (Thu, 07 Jan 2010)
New Revision: 40301
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.sh
grass/trunk/raster/r.gwflow/valid_calc_excavation.sh
Log:
Water budged calculatin activated in r.gwflow
Update of scipts and documentation
Modified: grass/trunk/lib/gpde/N_gwflow.c
===================================================================
--- grass/trunk/lib/gpde/N_gwflow.c 2010-01-07 14:17:06 UTC (rev 40300)
+++ grass/trunk/lib/gpde/N_gwflow.c 2010-01-07 14:30:25 UTC (rev 40301)
@@ -612,9 +612,9 @@
}
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;
}
Modified: grass/trunk/raster/r.gwflow/main.c
===================================================================
--- grass/trunk/raster/r.gwflow/main.c 2010-01-07 14:17:06 UTC (rev 40300)
+++ grass/trunk/raster/r.gwflow/main.c 2010-01-07 14:30:25 UTC (rev 40301)
@@ -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 @@
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,15 +128,31 @@
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";
param.type->type = TYPE_STRING;
@@ -426,16 +443,25 @@
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 @@
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 @@
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 @@
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 @@
return;
}
+
/* *************************************************************** */
/* ***** create and solve the linear equation system ************* */
/* *************************************************************** */
@@ -529,24 +588,22 @@
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_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_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"));
Modified: grass/trunk/raster/r.gwflow/r.gwflow.html
===================================================================
--- grass/trunk/raster/r.gwflow/r.gwflow.html 2010-01-07 14:17:06 UTC (rev 40300)
+++ grass/trunk/raster/r.gwflow/r.gwflow.html 2010-01-07 14:30:25 UTC (rev 40301)
@@ -16,8 +16,9 @@
</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,8 +27,12 @@
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>
The groundwater flow calculation is based on Darcy's law and a
@@ -78,27 +83,27 @@
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
Modified: grass/trunk/raster/r.gwflow/valid_calc_7x7.sh
===================================================================
--- grass/trunk/raster/r.gwflow/valid_calc_7x7.sh 2010-01-07 14:17:06 UTC (rev 40300)
+++ grass/trunk/raster/r.gwflow/valid_calc_7x7.sh 2010-01-07 14:30:25 UTC (rev 40301)
@@ -19,29 +19,16 @@
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
Modified: grass/trunk/raster/r.gwflow/valid_calc_excavation.sh
===================================================================
--- grass/trunk/raster/r.gwflow/valid_calc_excavation.sh 2010-01-07 14:17:06 UTC (rev 40300)
+++ grass/trunk/raster/r.gwflow/valid_calc_excavation.sh 2010-01-07 14:30:25 UTC (rev 40301)
@@ -27,24 +27,4 @@
#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
More information about the grass-commit
mailing list