Major documentation update.

Added link to diploma thesis.
Small code cleanup.
Examples updated and extended.

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

raster/r.gwflow/main.c

@@ -148,7 +148,7 @@ void set_params(void)
     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 @@ int main(int argc, char *argv[])
 
     /*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.*/

raster/r.gwflow/r.gwflow.html

@@ -1,10 +1,12 @@
 <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>
@@ -24,9 +26,9 @@ with <em>r.out.vtk</em>.
 <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 @@ 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 
-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 @@ differential equation is of the following form:
 <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 @@ the following states are supportet:
 <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 @@ paraview --data=/tmp/gwdata_unconf2d.vtk &
 
 <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&ouml;ren Gebbert
+<p>
+This work is based on the Diploma Thesis of S&ouml;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>

raster/r.solute.transport/example.py

@@ -19,7 +19,6 @@ grass.run_command("r.mapcalc", expression="well=if((row() == 50 && col() == 175)
 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.run_command("r.mapcalc", expression="null=0.0")
 #
 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="diff=0.0000001")
 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")

raster/r.solute.transport/main.c

@@ -35,7 +35,7 @@ typedef struct
 	*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,11 +186,10 @@ void set_params()
     param.stab->description =
 	_("Set the flow stabilizing scheme (full or exponential upwinding).");
 
-
-    param.sparse = G_define_flag();
-    param.sparse->key = 's';
-    param.sparse->description =
-	_("Use a sparse linear equation system, only available with iterative solvers");
+    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.cfl = G_define_flag();
     param.cfl->key = 'c';
@@ -247,9 +246,9 @@ int main(int argc, char *argv[])
     /*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 *create_solve_les(N_geom_data * geom, N_solute_transport_data2d * data,
     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 @@ N_les *create_solve_les(N_geom_data * geom, N_solute_transport_data2d * data,
 
     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 @@ N_les *create_solve_les(N_geom_data * geom, N_solute_transport_data2d * data,
 
     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);

raster/r.solute.transport/r.solute.transport.html

@@ -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 @@ differential equation is of the following form:
 
 <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 @@ differential equation is of the following form:
 
 <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 @@ module rapidely grow with the size of the input maps.
 
 <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="well=if((row() == 50 && col() == 175)
 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.run_command("r.mapcalc", expression="null=0.0")
 
 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="diff=0.0000001")
 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 @@ for dt in range(30):
 <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&ouml;ren Gebbert
+<p>
+This work is based on the Diploma Thesis of S&ouml;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>

raster/r.solute.transport/seguin_verify.py

@@ -46,7 +46,7 @@ grass.run_command("r.mapcalc", expression="null_1=0.0")
 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 @@ AL=100
 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.solute.transport", "c", "s", error=0.000000000000001, maxit
 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")

raster/r.solute.transport/seguin_verify_well.py

@@ -47,7 +47,7 @@ grass.run_command("r.mapcalc", expression="null_1=0.0")
 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 @@ AL=50
 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 @@ AL=10
 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)

raster3d/r3.gwflow/main.c

@@ -139,7 +139,7 @@ void set_params(void)
     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 @@ void set_params(void)
     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);

raster3d/r3.gwflow/r3.gwflow.html

@@ -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 @@ differential equation is of the following form:
 <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 @@ the following cell states are supportet:
 <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 @@ only work with normal quadratic matrices, so be careful using them with large ma
 (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 @@ hc_z=hydcond q=well s=syield r=recharge output=gwresult dt=8640000 vx=vx vy=vy v
 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&ouml;ren Gebbert
+<p>
+This work is based on the Diploma Thesis of S&ouml;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>