[GRASS-SVN] r40628 - in grass/trunk: raster/r.gwflow
raster/r.solute.transport raster3d/r3.gwflow
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Jan 23 08:49:09 EST 2010
Author: huhabla
Date: 2010-01-23 08:49:09 -0500 (Sat, 23 Jan 2010)
New Revision: 40628
Modified:
grass/trunk/raster/r.gwflow/main.c
grass/trunk/raster/r.gwflow/r.gwflow.html
grass/trunk/raster/r.solute.transport/example.py
grass/trunk/raster/r.solute.transport/main.c
grass/trunk/raster/r.solute.transport/r.solute.transport.html
grass/trunk/raster/r.solute.transport/seguin_verify.py
grass/trunk/raster/r.solute.transport/seguin_verify_well.py
grass/trunk/raster3d/r3.gwflow/main.c
grass/trunk/raster3d/r3.gwflow/r3.gwflow.html
Log:
Major documentation update.
Added link to diploma thesis.
Small code cleanup.
Examples updated and extended.
Modified: grass/trunk/raster/r.gwflow/main.c
===================================================================
--- grass/trunk/raster/r.gwflow/main.c 2010-01-23 13:47:02 UTC (rev 40627)
+++ grass/trunk/raster/r.gwflow/main.c 2010-01-23 13:49:09 UTC (rev 40628)
@@ -148,7 +148,7 @@
param.budget->required = NO;
param.budget->gisprompt = "new,raster,raster";
param.budget->description =
- _("Store the groundwater budget for each cell\n");
+ _("Store the groundwater budget for each cell [m^3/s]\n");
param.type = G_define_option();
param.type->key = "type";
@@ -315,7 +315,7 @@
/*Set the calculation time */
sscanf(param.dt->answer, "%lf", &(data->dt));
- G_message("Calculation time: %g", data->dt);
+ G_message(_("Calculation time: %g"), data->dt);
/*read all input maps into the memory and take care of the
* null values.*/
Modified: grass/trunk/raster/r.gwflow/r.gwflow.html
===================================================================
--- grass/trunk/raster/r.gwflow/r.gwflow.html 2010-01-23 13:47:02 UTC (rev 40627)
+++ grass/trunk/raster/r.gwflow/r.gwflow.html 2010-01-23 13:49:09 UTC (rev 40628)
@@ -1,11 +1,13 @@
<h2>DESCRIPTION</h2>
-This numerical program calculates transient, confined and
+This numerical program calculates implicit transient, confined and
unconfined groundwater flow in two dimensions based on
-raster maps and the current region resolution.
+raster maps and the current region settings.
All initial and boundary conditions must be provided as
-raster maps.
-
+raster maps. The unit in the location must be meters.
<p>
+This module is sensitive to mask settings. All cells which are outside the mask
+are ignored and handled as no flow boundaries.
+<p>
<center>
<img src=r_gwflow_concept.png border=0><br>
<table border=0 width=700>
@@ -24,9 +26,9 @@
<br>
<br>
The groundwater flow will always be calculated transient.
-If you want to calculate stady state, set the timestep
+For stady state computation set the timestep
to a large number (billions of seconds) or set the
-specific yield/ effective porosity raster maps to zero.
+specific yield/ effective porosity raster map to zero.
<br>
<br>
The water budget is calculated for each non inactive cell. The
@@ -35,12 +37,17 @@
<h2>NOTES</h2>
-The groundwater flow calculation is based on Darcy's law and a
-finite volume discretization. The solved groundwater flow partial
+The groundwater flow calculation is based on Darcy's law and a numerical implicit
+finite volume discretization. The discretization results in a symmetric and positive definit
+linear equation system in form of <i>Ax = b</i>, which must be solved. The groundwater flow partial
differential equation is of the following form:
<p>
-(dh/dt)*Ss = Kxx * (d^2h/dx^2) + Kyy * (d^2h/dy^2) + q
+(dh/dt)*S = div (K grad h) + q
+<p>
+In detail for 2 dimensions:
+<p>
+(dh/dt)*S = Kxx * (d^2h/dx^2) + Kyy * (d^2h/dy^2) + q
<ul>
<li>h -- the piezometric head im [m]</li>
@@ -48,7 +55,7 @@
<li>S -- the specific yield [1/m]</li>
<li>Kxx -- the hydraulic conductivity tensor part in x direction in [m/s]</li>
<li>Kyy -- the hydraulic conductivity tensor part in y direction in [m/s]</li>
-<li>q - inner source in meter per second [1/s]</li>
+<li>q - inner source/sink in meter per second [1/s]</li>
</ul>
<br>
@@ -63,10 +70,14 @@
<li>1 == active - this cell is used for groundwater floaw calculation, inner sources and recharge can be defined for those cells</li>
<li>2 == Dirichlet - cells of this type will have a fixed piezometric head value which do not change over the time </li>
</ul>
-
<br>
<br>
-The groundwater flow equation can be solved with several solvers.
+Note that all required raster maps are read into main memory. Additionally the
+linear equation system will be allocated, so the memory consumption of this
+module rapidely grow with the size of the input maps.
+<br>
+<br>
+The resulting linear equation system <i>Ax = b</i> 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
@@ -130,10 +141,15 @@
<h2>SEE ALSO</h2>
+<em><a href="r.solute.transport.html">r.solute.transport</a></em><br>
<em><a href="r3.gwflow.html">r3.gwflow</a></em><br>
<em><a href="r.out.vtk.html">r.out.vtk</a></em><br>
<h2>AUTHOR</h2>
-Soeren Gebbert
+Sören Gebbert
+<p>
+This work is based on the Diploma Thesis of Sören Gebbert available
+<a href="http://www.hydrogeologie.tu-berlin.de/fileadmin/fg66/_hydro/Diplomarbeiten/2007_Diplomarbeit_Soeren_Gebbert.pdf">here</a>
+at Technical University Berlin in Germany.
<p><i>Last changed: $Date$</i>
Modified: grass/trunk/raster/r.solute.transport/example.py
===================================================================
--- grass/trunk/raster/r.solute.transport/example.py 2010-01-23 13:47:02 UTC (rev 40627)
+++ grass/trunk/raster/r.solute.transport/example.py 2010-01-23 13:49:09 UTC (rev 40628)
@@ -19,7 +19,6 @@
grass.run_command("r.mapcalc", expression="hydcond=0.00005")
grass.run_command("r.mapcalc", expression="recharge=0")
grass.run_command("r.mapcalc", expression="top_conf=20")
-grass.run_command("r.mapcalc", expression="top_unconf=60")
grass.run_command("r.mapcalc", expression="bottom=0")
grass.run_command("r.mapcalc", expression="poros=0.17")
grass.run_command("r.mapcalc", expression="syield=0.0001")
@@ -27,7 +26,7 @@
#
grass.message(_("Compute a steady state groundwater flow"))
-grass.run_command("r.gwflow", "s", solver="cg", top="top_conf", bottom="bottom", phead="phead",\
+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=8640000000000, type="confined")
@@ -39,14 +38,14 @@
grass.run_command("r.mapcalc", expression="R=1.0")
# Compute the initial state
-grass.run_command("r.solute.transport", "s", solver="bicgstab", top="top_conf",\
+grass.run_command("r.solute.transport", 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_0", dt=3600, diff_x="diff",\
diff_y="diff", c="c", al=0.1, at=0.01)
# Compute the solute transport for 300 days in 10 day steps
for dt in range(30):
- grass.run_command("r.solute.transport", "s", solver="bicgstab", top="top_conf",\
+ grass.run_command("r.solute.transport", 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, vx="vx", vy="vy")
Modified: grass/trunk/raster/r.solute.transport/main.c
===================================================================
--- grass/trunk/raster/r.solute.transport/main.c 2010-01-23 13:47:02 UTC (rev 40627)
+++ grass/trunk/raster/r.solute.transport/main.c 2010-01-23 13:49:09 UTC (rev 40628)
@@ -35,7 +35,7 @@
*c, *status, *diff_x, *diff_y, *q, *cs, *r, *top, *nf, *cin,
*bottom, *vector_x, *vector_y, *bugded, *type, *dt, *maxit, *error, *solver, *sor,
*al, *at, *loops, *stab;
- struct Flag *sparse;
+ struct Flag *full_les;
struct Flag *cfl;
} paramType;
@@ -186,12 +186,11 @@
param.stab->description =
_("Set the flow stabilizing scheme (full or exponential upwinding).");
+ 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.");
- param.sparse = G_define_flag();
- param.sparse->key = 's';
- param.sparse->description =
- _("Use a sparse linear equation system, only available with iterative solvers");
-
param.cfl = G_define_flag();
param.cfl->key = 'c';
param.cfl->description =
@@ -247,9 +246,9 @@
/*Set the solver */
solver = param.solver->answer;
- if (strcmp(solver, G_MATH_SOLVER_DIRECT_LU) == 0 && param.sparse->answer)
+ if (strcmp(solver, G_MATH_SOLVER_DIRECT_LU) == 0 && !param.full_les->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)
+ if (strcmp(solver, G_MATH_SOLVER_DIRECT_GAUSS) == 0 && !param.full_les->answer)
G_fatal_error(_("The direct Gauss solver do not work with sparse matrices"));
@@ -492,19 +491,19 @@
N_les *les;
/*assemble the linear equation system */
- if (param.sparse->answer)
+ if (param.full_les->answer)
les =
- N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c,
+ N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c,
(void *)data, call);
else
les =
- N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c,
+ N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c,
(void *)data, call);
/*solve the equation system */
if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_JACOBI) == 0)
{
- if (param.sparse->answer)
+ if (!param.full_les->answer)
G_math_solver_sparse_jacobi(les->Asp, les->x, les->b, les->rows, maxit, sor, error);
else
G_math_solver_jacobi(les->A, les->x, les->b, les->rows, maxit, sor, error);
@@ -512,7 +511,7 @@
if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_SOR) == 0)
{
- if (param.sparse->answer)
+ if (!param.full_les->answer)
G_math_solver_sparse_gs(les->Asp, les->x, les->b, les->rows, maxit, sor, error);
else
G_math_solver_gs(les->A, les->x, les->b, les->rows, maxit, sor, error);
@@ -520,7 +519,7 @@
if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_BICGSTAB) == 0)
{
- if (param.sparse->answer)
+ if (!param.full_les->answer)
G_math_solver_sparse_bicgstab(les->Asp, les->x, les->b, les->rows, maxit, error);
else
G_math_solver_bicgstab(les->A, les->x, les->b, les->rows, maxit, error);
Modified: grass/trunk/raster/r.solute.transport/r.solute.transport.html
===================================================================
--- grass/trunk/raster/r.solute.transport/r.solute.transport.html 2010-01-23 13:47:02 UTC (rev 40627)
+++ grass/trunk/raster/r.solute.transport/r.solute.transport.html 2010-01-23 13:49:09 UTC (rev 40628)
@@ -1,37 +1,42 @@
<H2>DESCRIPTION</H2>
-This numerical program calculates transient
-solute transport in porous media in the saturated zone of an aquifer in two dimensions based on
-raster maps and the current region resolution.
-All initial- and boundary-conditions must be provided as
-raster maps.
+This numerical program calculates numerical implicit transient and steady state
+solute transport in porous media in the saturated zone of an aquifer. The computation is based on
+raster maps and the current region settings. All initial- and boundary-conditions must be provided as
+raster maps. The unit in the location must be meters.
<br>
+<p>
+This module is sensitive to mask settings. All cells which are outside the mask
+are ignored and handled as no flow boundaries.
<br>
This module calculates the concentration of the solution and optional the
velocity field, based on the hydraulic conductivity,
-the effective porosity and the
-initial piecometric heads.
+the effective porosity and the initial piecometric heads.
The vector components can be visualized with paraview if they are exported
with r.out.vtk.
<br>
<br>
-Use <A HREF="r.gwflow.html">r.gwflow</A> to compute the piezometric hights
-of the aqufer. The piezometric heights and the hydraulic conductivity
+Use <A HREF="r.gwflow.html">r.gwflow</A> to compute the piezometric heights
+of the aquifer. The piezometric heights and the hydraulic conductivity
are used to compute the flow direction and the mean velocity of the groundwater.
-They are the base of the concentration computation.
+This is the base of the solute transport computation.
<br>
<br>
The solute transport will always be calculated transient.
-If you want to calculate stady state, set the timestep
+For stady state computation set the timestep
to a large number (billions of seconds).
<br>
<br>
-To reduce the numerical dispersion you have to use small time steps.
-You can choose additionally between full and exponential upwinding to reduce
-the numerical dispersion.
+To reduce the numerical dispersion, which is a consequence of the convection term and
+the finite volume discretization, you can use small time steps and choose between full
+and exponential upwinding.
<H2>NOTES</H2>
-The solute transport partial
+The solute transport calculation is based on a diffusion/convection partial differential equation and
+a numerical implicit finite volume discretization. Specific for this kind of differential
+equation is the combination of a diffusion/dispersion term and a convection term.
+The discretization results in an unsymmetric linear equation system in form of <i>Ax = b</i>,
+which must be solved. The solute transport partial
differential equation is of the following form:
<p>
@@ -39,6 +44,7 @@
<ul>
<li>c -- the concentration [kg/m^3]</li>
+<li>u -- vector of mean groundwater flow velocity</li>
<li>dt -- the time step for transient calculation in seconds [s]</li>
<li>R -- the linear retardation coefficient [-]</li>
<li>D -- the diffusion and dispersion tensor [m^2/s]</li>
@@ -50,10 +56,10 @@
<br>
<br>
-Two different boundary conditions are implemented,
-the Dirichlet and Neumann conditions.
-The calculation and boundary status of single cells can be set with the status map,
-the following states are supportet:
+Three different boundary conditions are implemented,
+the Dirichlet, Transmission and Neumann conditions.
+The calculation and boundary status of single cells can be set with the status map.
+The following states are supportet:
<ul>
<li>0 == inactive - the cell with status 0 will not be calulated, active cells will have a no flow boundary to an inactive cell</li>
@@ -70,15 +76,16 @@
<br>
<br>
-The solute transport equation can be solved with several solvers.
+The resulting linear equation system <i>Ax = b</i> can be solved with several solvers.
Several iterative solvers with unsymmetric sparse and quadratic matrices support are implemented.
The jacobi method, the Gauss-Seidel method and the biconjugate gradients-stabilized (bicgstab) method.
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).
+Always prefer a sparse matrix solver.
<H2>EXAMPLE</H2>
-Use this small script to create a working
+Use this small python script to create a working
groundwater flow / solute transport area and data.
Make sure you are not in a lat/lon projection.
@@ -104,7 +111,6 @@
grass.run_command("r.mapcalc", expression="hydcond=0.00005")
grass.run_command("r.mapcalc", expression="recharge=0")
grass.run_command("r.mapcalc", expression="top_conf=20")
-grass.run_command("r.mapcalc", expression="top_unconf=60")
grass.run_command("r.mapcalc", expression="bottom=0")
grass.run_command("r.mapcalc", expression="poros=0.17")
grass.run_command("r.mapcalc", expression="syield=0.0001")
@@ -112,7 +118,7 @@
grass.message(_("Compute a steady state groundwater flow"))
-grass.run_command("r.gwflow", "s", solver="cg", top="top_conf", bottom="bottom", phead="phead",\
+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=8640000000000, type="confined")
@@ -124,14 +130,14 @@
grass.run_command("r.mapcalc", expression="R=1.0")
# Compute the initial state
-grass.run_command("r.solute.transport", "s", solver="bicgstab", top="top_conf",\
+grass.run_command("r.solute.transport", 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_0", dt=3600, diff_x="diff",\
diff_y="diff", c="c", al=0.1, at=0.01)
# Compute the solute transport for 300 days in 10 day steps
for dt in range(30):
- grass.run_command("r.solute.transport", "s", solver="bicgstab", top="top_conf",\
+ grass.run_command("r.solute.transport", 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)
@@ -145,7 +151,12 @@
<EM><A HREF="r3.gwflow.html">r3.gwflow</A></EM><br>
<EM><A HREF="r.out.vtk.html">r.out.vtk</A></EM><br>
-<H2>AUTHOR</H2>
-Soeren Gebbert
+<h2>AUTHOR</h2>
+Sören Gebbert
+<p>
+This work is based on the Diploma Thesis of Sören Gebbert available
+<a href="http://www.hydrogeologie.tu-berlin.de/fileadmin/fg66/_hydro/Diplomarbeiten/2007_Diplomarbeit_Soeren_Gebbert.pdf">here</a>
+at Technical University Berlin in Germany.
+
<p><i>Last changed: $Date: 2007/02/15 23:27:58 $</i>
Modified: grass/trunk/raster/r.solute.transport/seguin_verify.py
===================================================================
--- grass/trunk/raster/r.solute.transport/seguin_verify.py 2010-01-23 13:47:02 UTC (rev 40627)
+++ grass/trunk/raster/r.solute.transport/seguin_verify.py 2010-01-23 13:49:09 UTC (rev 40628)
@@ -46,7 +46,7 @@
grass.message("First compute a steady state groundwater flow")
# Compute the steady state groundwater flow
-grass.run_command("r.gwflow", "s", solver="cg", top="top_conf_1", bottom="bottom_1", phead="phead_1",\
+grass.run_command("r.gwflow", solver="cg", top="top_conf_1", bottom="bottom_1", phead="phead_1",\
status="status_1", hc_x="hydcond_1", hc_y="hydcond_1", \
q="well_1", s="syield_1", r="recharge_1", output="gwresult_conf_1",\
dt=8640000000000, type="confined")
@@ -71,7 +71,7 @@
AT=10
# 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",\
+grass.run_command("r.solute.transport", "c", 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, vx="stresult_conf_vel_1_x", vy="stresult_conf_vel_1_y")
@@ -82,7 +82,7 @@
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",\
+grass.run_command("r.solute.transport", "c", 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_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, 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-23 13:47:02 UTC (rev 40627)
+++ grass/trunk/raster/r.solute.transport/seguin_verify_well.py 2010-01-23 13:49:09 UTC (rev 40628)
@@ -47,7 +47,7 @@
grass.message("First compute a steady state groundwater flow")
# Compute the steady state groundwater flow
-grass.run_command("r.gwflow", "s", solver="cg", top="top_conf_1", bottom="bottom_1", phead="phead_1",\
+grass.run_command("r.gwflow", solver="cg", top="top_conf_1", bottom="bottom_1", phead="phead_1",\
status="status_1", hc_x="hydcond_1", hc_y="hydcond_1", \
q="well_1", s="syield_1", r="recharge_1", output="gwresult_conf_1",\
dt=8640000000000, type="confined")
@@ -72,7 +72,7 @@
AT=5
# Compute the solute transport using the above defined dispersivity coefficients for a timestep of 250d
-grass.run_command("r.solute.transport", "c", "s", error=0.000000000000001, maxit=1000, solver="bicgstab",\
+grass.run_command("r.solute.transport", "c", 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)
@@ -83,7 +83,7 @@
AT=1
# Compute the solute transport using the above defined dispersivity coefficients for a timestep of 250d
-grass.run_command("r.solute.transport", "c", "s", error=0.000000000000001, maxit=1000, solver="bicgstab",\
+grass.run_command("r.solute.transport", "c", 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)
Modified: grass/trunk/raster3d/r3.gwflow/main.c
===================================================================
--- grass/trunk/raster3d/r3.gwflow/main.c 2010-01-23 13:47:02 UTC (rev 40627)
+++ grass/trunk/raster3d/r3.gwflow/main.c 2010-01-23 13:49:09 UTC (rev 40628)
@@ -139,7 +139,7 @@
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");
+ _("Calculate and store the groundwater filter velocity vector part in z direction [m/s]\n");
param.budget = G_define_option();
param.budget->key = "budget";
@@ -147,7 +147,7 @@
param.budget->required = NO;
param.budget->gisprompt = "new,grid3,3d-raster";
param.budget->description =
- _("Store the groundwater budget for each cell\n");
+ _("Store the groundwater budget for each cell [m^3/s]\n");
param.dt = N_define_standard_option(N_OPT_CALC_TIME);
param.maxit = N_define_standard_option(N_OPT_MAX_ITERATIONS);
Modified: grass/trunk/raster3d/r3.gwflow/r3.gwflow.html
===================================================================
--- grass/trunk/raster3d/r3.gwflow/r3.gwflow.html 2010-01-23 13:47:02 UTC (rev 40627)
+++ grass/trunk/raster3d/r3.gwflow/r3.gwflow.html 2010-01-23 13:49:09 UTC (rev 40628)
@@ -1,29 +1,37 @@
<H2>DESCRIPTION</H2>
-This numerical program calculates transient, confined groundwater flow in three dimensions
-based on volume maps and the current 3d region resolution.
-All initial- and boundary-conditions must be provided as
-volume maps.
+This numerical program calculates implicit transient and steady state,
+confined groundwater flow in three dimensions
+based on volume maps and the current 3d region settings.
+All initial- and boundary-conditions must be provided as volume maps.
+The unit in the location must be meters.
<br>
+<p>
+This module is sensitive to mask settings. All cells which are outside the mask
+are ignored and handled as no flow boundaries.
<br>
-This module calculates the piezometric head and optionally the
-groundwater velocity field.
+This module calculates the piezometric head and optionally the water balance for each cell
+and the groundwater velocity field in 3 dimensions.
The vector components can be visualized with paraview if they are exported
with <em>r3.out.vtk</em>.
<br>
<br>
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 raster maps to zero.
+For stady state computation set the timestep
+to a large number (billions of seconds) or set the
+specific yield raster map to zero.
-
<H2>NOTES</H2>
-The groundwater flow calculation is based on Darcy's law and a
-finite volume discretization. The groundwater flow partial
+The groundwater flow calculation is based on Darcy's law and a numerical implicit
+finite volume discretization. The discretization results in a symmetric and positive definit
+linear equation system in form of <i>Ax = b</i>, which must be solved. The groundwater flow partial
differential equation is of the following form:
<p>
+(dh/dt)*S = div (K grad h) + q
+<p>
+In detail for 3 dimensions:
+<p>
(dh/dt)*S = Kxx * (d^2h/dx^2) + Kyy * (d^2h/dy^2) + Kzz * (d^2h/dz^2) + q
<ul>
@@ -34,13 +42,14 @@
<li>Kxx -- the hydraulic conductivity tensor part in x direction in meter per second [m/s]</li>
<li>Kyy -- the hydraulic conductivity tensor part in y direction in meter per seconds [m/s]</li>
<li>Kzz -- the hydraulic conductivity tensor part in z direction in meter per seconds [m/s]</li>
-<li>q - inner source in [1/s]</li>
+<li>q - inner source/sinc in [1/s]</li>
</ul>
<br>
<br>
Two different boundary conditions are implemented,
-the Dirichlet and Neumann conditions. By default the calculation area is surrounded by homogeneous Neumann boundary conditions.
+the Dirichlet and Neumann conditions. By default the calculation area
+is surrounded by homogeneous Neumann boundary conditions.
The calculation and boundary status of single cells can be set with the status map,
the following cell states are supportet:
@@ -49,10 +58,14 @@
<li>1 == active - this cell is used for groundwater calculation, inner sources can be defined for those cells</li>
<li>2 == Dirichlet - cells of this type will have a fixed piezometric head value which do not change over time </li>
</ul>
-
<br>
<br>
-The groundwater flow equation can be solved with several solvers.
+Note that all required raster maps are read into main memory. Additionally the
+linear equation system will be allocated, so the memory consumption of this
+module rapidely grow with the size of the input maps.
+<br>
+<br>
+The resulting linear equation system <i>Ax = b</i> 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
@@ -60,7 +73,7 @@
(maps of size 10.000 cells will need more than one gigabyte of RAM).
Always prefer a sparse matrix solver.
-<H2>EXAMPLE</H2>
+<H2>EXAMPLE 1</H2>
Use this small script to create a working
groundwater flow area and data. Make sure you are not in a lat/lon projection.
@@ -83,14 +96,50 @@
r3.out.vtk -p in=gwresult,status,budget vector=vx,vy,vz out=/tmp/gwdata3d.vtk
#now load the data into paraview
+paraview --data=/tmp/gwdata3d.vtk
</pre></div>
+<H2>EXAMPLE 2</H2>
+This will create a nice 3d model with geolgical layer with different
+hydraulic conductivities. Make sure you are not in a lat/lon projection.
+
+<div class="code"><pre>
+# set the region accordingly
+g.region res=15 res3=15 t=500 b=0 n=1000 s=0 w=0 e=1000
+
+#now create the input raster maps for a confined aquifer
+r3.mapcalc --o expression="phead = if(col() == 1 && depth() == 33, 50, 40)"
+r3.mapcalc --o expression="status = if(col() == 1 && depth() == 33, 2, 1)"
+r3.mapcalc --o expression="well = if(row() == 20 && col() == 20 && depth() == 3, -0.25, 0)"
+r3.mapcalc --o expression="well = if(row() == 50 && col() == 50 && depth() == 3, -0.25, well)"
+r3.mapcalc --o expression="hydcond = 0.0025"
+r3.mapcalc --o expression="hydcond = if(depth() < 30 && depth() > 23 && col() < 60, 0.000025, hydcond)"
+r3.mapcalc --o expression="hydcond = if(depth() < 20 && depth() > 13 && col() > 7, 0.000025, hydcond)"
+r3.mapcalc --o expression="hydcond = if(depth() < 10 && depth() > 7 && col() < 60, 0.000025, hydcond)"
+r3.mapcalc --o expression="syield = 0.0001"
+
+r3.gwflow --o solver=cg phead=phead status=status hc_x=hydcond hc_y=hydcond \
+hc_z=hydcond q=well s=syield 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,budget,hydcond,well vector=vx,vy,vz out=/tmp/gwdata3d.vtk
+
+#now load the data into paraview
+paraview --data=/tmp/gwdata3d.vtk
+</pre></div>
+
<H2>SEE ALSO</H2>
<EM><A HREF="r.gwflow.html">r.gwflow</A></EM><br>
+<EM><A HREF="r.solute.transport.html">r.solute.transport</A></EM><br>
<EM><A HREF="r3.out.vtk.html">r3.out.vtk</A></EM><br>
-<H2>AUTHOR</H2>
-Soeren Gebbert
+<h2>AUTHOR</h2>
+Sören Gebbert
+<p>
+This work is based on the Diploma Thesis of Sören Gebbert available
+<a href="http://www.hydrogeologie.tu-berlin.de/fileadmin/fg66/_hydro/Diplomarbeiten/2007_Diplomarbeit_Soeren_Gebbert.pdf">here</a>
+at Technical University Berlin in Germany.
+
<p><i>Last changed: $Date$</i>
More information about the grass-commit
mailing list