[GRASS-SVN] r40304 - in grass/trunk/raster: r.gwflow
r.solute.transport
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jan 7 11:54:23 EST 2010
Author: huhabla
Date: 2010-01-07 11:54:22 -0500 (Thu, 07 Jan 2010)
New Revision: 40304
Added:
grass/trunk/raster/r.gwflow/valid_calc_7x7.py
grass/trunk/raster/r.gwflow/valid_calc_excavation.py
Removed:
grass/trunk/raster/r.gwflow/valid_calc_7x7.sh
grass/trunk/raster/r.gwflow/valid_calc_excavation.sh
Modified:
grass/trunk/raster/r.gwflow/main.c
grass/trunk/raster/r.solute.transport/example.py
grass/trunk/raster/r.solute.transport/main.c
grass/trunk/raster/r.solute.transport/seguin_verify.py
grass/trunk/raster/r.solute.transport/seguin_verify_well.py
Log:
Port of r.gwflow shell scripts to python
Update r.solute.transport examples + bugfixing
Modified: grass/trunk/raster/r.gwflow/main.c
===================================================================
--- grass/trunk/raster/r.gwflow/main.c 2010-01-07 14:46:08 UTC (rev 40303)
+++ grass/trunk/raster/r.gwflow/main.c 2010-01-07 16:54:22 UTC (rev 40304)
@@ -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 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 @@
_("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 @@
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 @@
}
/* ************************************************************************* */
-/* 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
Added: grass/trunk/raster/r.gwflow/valid_calc_7x7.py
===================================================================
--- grass/trunk/raster/r.gwflow/valid_calc_7x7.py (rev 0)
+++ grass/trunk/raster/r.gwflow/valid_calc_7x7.py 2010-01-07 16:54:22 UTC (rev 40304)
@@ -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
+
Deleted: grass/trunk/raster/r.gwflow/valid_calc_7x7.sh
===================================================================
--- grass/trunk/raster/r.gwflow/valid_calc_7x7.sh 2010-01-07 14:46:08 UTC (rev 40303)
+++ grass/trunk/raster/r.gwflow/valid_calc_7x7.sh 2010-01-07 16:54:22 UTC (rev 40304)
@@ -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
Added: grass/trunk/raster/r.gwflow/valid_calc_excavation.py
===================================================================
--- grass/trunk/raster/r.gwflow/valid_calc_excavation.py (rev 0)
+++ grass/trunk/raster/r.gwflow/valid_calc_excavation.py 2010-01-07 16:54:22 UTC (rev 40304)
@@ -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")
Deleted: grass/trunk/raster/r.gwflow/valid_calc_excavation.sh
===================================================================
--- grass/trunk/raster/r.gwflow/valid_calc_excavation.sh 2010-01-07 14:46:08 UTC (rev 40303)
+++ grass/trunk/raster/r.gwflow/valid_calc_excavation.sh 2010-01-07 16:54:22 UTC (rev 40304)
@@ -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
Modified: grass/trunk/raster/r.solute.transport/example.py
===================================================================
--- grass/trunk/raster/r.solute.transport/example.py 2010-01-07 14:46:08 UTC (rev 40303)
+++ grass/trunk/raster/r.solute.transport/example.py 2010-01-07 16:54:22 UTC (rev 40304)
@@ -8,7 +8,7 @@
# 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 @@
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")
Modified: grass/trunk/raster/r.solute.transport/main.c
===================================================================
--- grass/trunk/raster/r.solute.transport/main.c 2010-01-07 14:46:08 UTC (rev 40303)
+++ grass/trunk/raster/r.solute.transport/main.c 2010-01-07 16:54:22 UTC (rev 40304)
@@ -33,7 +33,7 @@
{
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,13 +130,23 @@
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);
param.error = N_define_standard_option(N_OPT_ITERATION_ERROR);
@@ -208,7 +218,6 @@
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 @@
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 @@
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);
Modified: grass/trunk/raster/r.solute.transport/seguin_verify.py
===================================================================
--- grass/trunk/raster/r.solute.transport/seguin_verify.py 2010-01-07 14:46:08 UTC (rev 40303)
+++ grass/trunk/raster/r.solute.transport/seguin_verify.py 2010-01-07 16:54:22 UTC (rev 40304)
@@ -58,7 +58,7 @@
# 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 @@
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 @@
# 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")
Modified: grass/trunk/raster/r.solute.transport/seguin_verify_well.py
===================================================================
--- grass/trunk/raster/r.solute.transport/seguin_verify_well.py 2010-01-07 14:46:08 UTC (rev 40303)
+++ grass/trunk/raster/r.solute.transport/seguin_verify_well.py 2010-01-07 16:54:22 UTC (rev 40304)
@@ -75,7 +75,7 @@
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 @@
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)
More information about the grass-commit
mailing list