[GRASS-SVN] r40097 - grass/trunk/raster/r.solute.transport
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Dec 20 18:43:16 EST 2009
Author: huhabla
Date: 2009-12-20 18:43:15 -0500 (Sun, 20 Dec 2009)
New Revision: 40097
Added:
grass/trunk/raster/r.solute.transport/example.py
grass/trunk/raster/r.solute.transport/seguin_verify.py
grass/trunk/raster/r.solute.transport/seguin_verify_well.py
Modified:
grass/trunk/raster/r.solute.transport/main.c
grass/trunk/raster/r.solute.transport/r.solute.transport.html
Log:
Documentation update of r.solute.transport
Removed redundant examples
Converted shell example scripts into python scripts
Copied: grass/trunk/raster/r.solute.transport/example.py (from rev 40082, grass/trunk/raster/r.solute.transport/example.sh)
===================================================================
--- grass/trunk/raster/r.solute.transport/example.py (rev 0)
+++ grass/trunk/raster/r.solute.transport/example.py 2009-12-20 23:43:15 UTC (rev 40097)
@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+# This is an example script how groundwater flow and solute transport are
+# computed within grass
+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 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)
+grass.run_command("r.mapcalc", expression="phead=if(col() == 1 , 50, 40)")
+grass.run_command("r.mapcalc", expression="phead=if(col() ==200 , 45 + row()/40, phead)")
+grass.run_command("r.mapcalc", expression="status=if(col() == 1 || col() == 200 , 2, 1)")
+grass.run_command("r.mapcalc", expression="well=if((row() == 50 && col() == 175) || (row() == 10 && col() == 135) , -0.001, 0)")
+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")
+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",\
+ status="status", hc_x="hydcond", hc_y="hydcond", q="well", s="syield",\
+ r="recharge", output="gwresult_conf", dt=8640000000000, type="confined")
+
+grass.message(_("generate the transport data"))
+grass.run_command("r.mapcalc", expression="c=if(col() == 15 && row() == 75 , 500.0, 0.0)")
+grass.run_command("r.mapcalc", expression="cs=if(col() == 15 && row() == 75 , 0.0, 0.0)")
+grass.run_command("r.mapcalc", expression="tstatus=if(col() == 1 || col() == 200 , 3, 1)")
+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",\
+ 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",\
+ 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)
+
Modified: grass/trunk/raster/r.solute.transport/main.c
===================================================================
--- grass/trunk/raster/r.solute.transport/main.c 2009-12-20 20:48:48 UTC (rev 40096)
+++ grass/trunk/raster/r.solute.transport/main.c 2009-12-20 23:43:15 UTC (rev 40097)
@@ -82,12 +82,12 @@
param.diff_x = G_define_standard_option(G_OPT_R_INPUT);
param.diff_x->key = "diff_x";
param.diff_x->description =
- _("The x-part of the diffusion tensor in [m�/s]");
+ _("The x-part of the diffusion tensor in [m^2/s]");
param.diff_y = G_define_standard_option(G_OPT_R_INPUT);
param.diff_y->key = "diff_y";
param.diff_y->description =
- _("The y-part of the diffusion tensor in [m�/s]");
+ _("The y-part of the diffusion tensor in [m^2/s]");
param.q = G_define_standard_option(G_OPT_R_INPUT);
param.q->key = "q";
@@ -97,14 +97,21 @@
param.cin = G_define_standard_option(G_OPT_R_INPUT);
param.cin->key = "cin";
param.cin->required = NO;
- param.cin->description = _("Concentration sources and sinks in [kg/m^3]");
+ param.cin->gisprompt = "old,raster,raster";
+ param.cin->description = _("concentration sources and sinks bounded to a "
+ "water source or sink in [kg/s]");
+
param.cs = G_define_standard_option(G_OPT_R_INPUT);
param.cs->key = "cs";
- param.cs->description = _("Concentration sources and sinks in [kg/m^3]");
+ param.cs->type = TYPE_STRING;
+ param.cs->required = YES;
+ param.cs->gisprompt = "old,raster,raster";
+ param.cs->description = _("concentration of inner sources and inner sinks in [kg/s] "
+ "(i.e. a chemical reaction)");
param.r = G_define_standard_option(G_OPT_R_INPUT);
- param.r->key = "R";
+ param.r->key = "r";
param.r->description = _("Retardation factor [-]");
param.nf = G_define_standard_option(G_OPT_R_INPUT);
@@ -120,16 +127,15 @@
param.bottom->description = _("Bottom surface of the aquifer in [m]");
param.output = G_define_standard_option(G_OPT_R_OUTPUT);
- param.output->description =
- _("The result of the numericalsolute transport calculation "
- "will be written to this map. [kg/m�]");
+ 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]");
+ _("Calculate the groundwater distance velocity vector field and "
+ "write the x, and y components to maps named name_(xy), [m/s]");
param.dt = N_define_standard_option(N_OPT_CALC_TIME);
param.maxit = N_define_standard_option(N_OPT_MAX_ITERATIONS);
@@ -168,7 +174,7 @@
param.stab->answer = "full";
param.stab->options = "full,exp";
param.stab->description =
- _("Set the flow stabilizing scheme.");
+ _("Set the flow stabilizing scheme (full or exponential upwinding).");
param.sparse = G_define_flag();
@@ -213,7 +219,8 @@
G_add_keyword(_("raster"));
G_add_keyword(_("solute transport"));
module->description =
- _("Numerical calculation program for transient, confined and unconfined solute transport in two dimensions.");
+ _("Numerical calculation program for transient, confined and unconfined "
+ "solute transport in two dimensions");
/* Get parameters from user */
set_params();
Modified: grass/trunk/raster/r.solute.transport/r.solute.transport.html
===================================================================
--- grass/trunk/raster/r.solute.transport/r.solute.transport.html 2009-12-20 20:48:48 UTC (rev 40096)
+++ grass/trunk/raster/r.solute.transport/r.solute.transport.html 2009-12-20 23:43:15 UTC (rev 40097)
@@ -1,12 +1,12 @@
<H2>DESCRIPTION</H2>
This numerical program calculates transient
-solute transport in two dimensions based on
+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 starting- and boundary-conditions must be provided as
+All initial- and boundary-conditions must be provided as
raster maps.
<br>
<br>
-This module calculates the concentration and optional the
+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.
@@ -14,11 +14,22 @@
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
+are used to compute the flow direction and the mean velocity of the groundwater.
+They are the base of the concentration computation.
+
+<br>
+<br>
The solute transport will always be calculated transient.
If you want to calculate stady state, 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.
-
<H2>NOTES</H2>
The solute transport partial
differential equation is of the following form:
@@ -27,14 +38,14 @@
(dc/dt)*R = div ( D grad c - uc) + cs -q/nf(c - c_in)
<ul>
-<li>c -- the concentration </li>
+<li>c -- the concentration [kg/m^3]</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</li>
-<li>cs -- inner concentration sources/sinks</li>
-<li>c_in -- the solute concentration of influent water</li>
-<li>q -- inner well sources/sinks</li>
-<li>nf -- the effective porosity </li>
+<li>R -- the linear retardation coefficient [-]</li>
+<li>D -- the diffusion and dispersion tensor [m^2/s]</li>
+<li>cs -- inner concentration sources/sinks [kg/m^3]</li>
+<li>c_in -- the solute concentration of influent water [kg/m^3]</li>
+<li>q -- inner well sources/sinks [m^3/s]</li>
+<li>nf -- the effective porosity [-] </li>
</ul>
<br>
@@ -48,11 +59,17 @@
<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>
<li>1 == active - this cell is used for sloute transport calculation, inner sources can be defined for those cells</li>
<li>2 == Dirichlet - cells of this type will have a fixed concentration value which do not change over time </li>
-<li>3 == Transmission - cells of this type should be paced on out-flow boundaries to assure the flow of the solute stream out </li>
+<li>3 == Transmission - cells of this type should be placed on out-flow boundaries to assure the flow of the solute stream out </li>
</ul>
<br>
<br>
+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 solute transport equation 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.
@@ -66,54 +83,65 @@
Make sure you are not in a lat/lon projection.
<div class="code"><pre>
-#!/bin/sh
-# set the region accordingly
-g.region res=10 res3=10 t=100 b=0 n=1000 s=0 w=0 e=2000
-r.mapcalc "phead=if(col() == 1 , 50, 40)"
-r.mapcalc "phead=if(col() ==200 , 20 + row()/5, phead)"
-r.mapcalc "status=if(col() == 1 || col() == 200 , 2, 1)"
-r.mapcalc "well=if((row() == 50 && col() == 175) || (row() == 10 && col() == 135) , -0.01, 0)"
-r.mapcalc "hydcond=0.00005"
-r.mapcalc "reacharge=0"
-r.mapcalc "top_conf=20"
-r.mapcalc "top_unconf=60"
-r.mapcalc "bottom=0"
-r.mapcalc "poros=0.17"
-r.mapcalc "syield=0.0001"
-r.mapcalc "null=0.0"
+#!/usr/bin/env python
+# This is an example script how groundwater flow and solute transport are
+# computed within grass
+import sys
+import os
+import grass.script as grass
-#First compute a steady state groundwater flow
-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=reacharge output=gwresult_conf dt=8640000000000 type=confined
+# Overwrite existing maps
+grass.run_command("g.gisenv", set="OVERWRITE=1")
-#create the flow direction vectors
-r.flow elevin=gwresult_conf flout=gwresult_conf_flow skip=15 --o
-# create contour lines
-r.contour input=gwresult_conf output=gwresult_conf_contour step=2 --o
+grass.message(_("Set the region"))
-#generate the transport data
-r.mapcalc "c=if(col() == 15 && row() == 75 , 500.0, 0.0)"
-r.mapcalc "cs=if(col() == 15 && row() == 75 , 0.0, 0.0)"
-r.mapcalc "tstatus=if(col() == 1 || col() == 200 , 3, 1)"
-r.mapcalc "diff=0.00001"
-r.mapcalc "R=1.0"
+# 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)
+grass.run_command("r.mapcalc", expression="phead=if(col() == 1 , 50, 40)")
+grass.run_command("r.mapcalc", expression="phead=if(col() ==200 , 45 + row()/40, phead)")
+grass.run_command("r.mapcalc", expression="status=if(col() == 1 || col() == 200 , 2, 1)")
+grass.run_command("r.mapcalc", expression="well=if((row() == 50 && col() == 175) || (row() == 10 && col() == 135) , -0.001, 0)")
+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")
+grass.run_command("r.mapcalc", expression="null=0.0")
-r.solute.transport --o -s solver=bicgstab top=top_conf bottom=bottom phead=gwresult_conf status=tstatus \
- hc_x=hydcond hc_y=hydcond R=R cs=cs q=null nf=nf output=stresult_conf dt=2592000 diff_x=diff diff_y=diff \
- c=c
+grass.message(_("Compute a steady state groundwater flow"))
+grass.run_command("r.gwflow", "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=8640000000000, type="confined")
-#compute in an infinite loop the transport with a timestept of 30 days
-while true ; do
-r.solute.transport --o -s solver=bicgstab top=top_conf bottom=bottom phead=gwresult_conf status=tstatus \
- hc_x=hydcond hc_y=hydcond R=R cs=cs q=null nf=nf output=stresult_conf dt=2592000 diff_x=diff diff_y=diff \
- c=stresult_conf
-done
+grass.message(_("generate the transport data"))
+grass.run_command("r.mapcalc", expression="c=if(col() == 15 && row() == 75 , 500.0, 0.0)")
+grass.run_command("r.mapcalc", expression="cs=if(col() == 15 && row() == 75 , 0.0, 0.0)")
+grass.run_command("r.mapcalc", expression="tstatus=if(col() == 1 || col() == 200 , 3, 1)")
+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",\
+ 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",\
+ 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)
+
+
</pre></div>
<H2>SEE ALSO</H2>
+<EM><A HREF="r.gwflow.html">r.gwflow</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>
Copied: grass/trunk/raster/r.solute.transport/seguin_verify.py (from rev 40082, grass/trunk/raster/r.solute.transport/seguin_verify.sh)
===================================================================
--- grass/trunk/raster/r.solute.transport/seguin_verify.py (rev 0)
+++ grass/trunk/raster/r.solute.transport/seguin_verify.py 2009-12-20 23:43:15 UTC (rev 40097)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+# Shellscript to verify r.solute.transport calculation, this calculation is based on
+# the example 1.1 and 1.2 at page 175 of the following book:
+# title = "Str{\"o}mungs und Transportmodellierung",
+# author = "Lege, T. and Kolditz, O. and Zielke W.",
+# publisher = "Springer (Berlin; Heidelberg; New York; Barcelona; Hongkong; London; Mailand; Paris; Singapur; Tokio)",
+# edition = "2. Auflage",
+# year = "1996",
+# series = "Handbuch zur Erkundung des Untergrundes von Deponien und Altlasten"
+
+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=25, res3=25, t=25, b=0, n=1000, s=0, w=0, e=2000)
+
+grass.message(_("Create all the input maps needed for groundwater flow computation"))
+
+# Initial condition of the piezometric head, we have a huge gradient from 280m to 50m
+# over a distance of 2000m
+grass.run_command("r.mapcalc", expression="phead_1=if(col() == 1 , 200, 50)")
+# Set the active cells and the dirichlet boundary condition
+grass.run_command("r.mapcalc", expression="status_1=if(col() == 1 || col() == 80 , 2, 1)")
+# We have a no wells
+grass.run_command("r.mapcalc", expression="well_1=0")
+# The hydraulic conductivity
+grass.run_command("r.mapcalc", expression="hydcond_1=0.0001")
+# The recharge, well we have no recharge at all
+grass.run_command("r.mapcalc", expression="recharge_1=0")
+# We have a confined aquifer with a height of 25m
+grass.run_command("r.mapcalc", expression="top_conf_1=25")
+# Bottom of the aquifer starts at 0m
+grass.run_command("r.mapcalc", expression="bottom_1=0")
+# This porosity of sand aquifer
+grass.run_command("r.mapcalc", expression="poros_1=0.17")
+# The specific yield of a confined aquifer
+grass.run_command("r.mapcalc", expression="syield_1=0.0001")
+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",\
+ 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")
+
+grass.message(_("generate the transport data"))
+
+# The initial concentration is zero
+grass.run_command("r.mapcalc", expression="c_1=if(col() == 20 && row() == 20 , 0.0, 0.0)")
+# 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)")
+# 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
+grass.run_command("r.mapcalc", expression="diff_1=0.0")
+# Normal retardation
+grass.run_command("r.mapcalc", expression="R_1=1.0")
+
+# Longitudinal and transversal dispersivity length
+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",\
+ 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")
+
+# The second computation uses different porosity for higher groundwater velocity
+
+# Used to compute a lower velocity
+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",\
+ top="top_conf_1", bottom="bottom_1", phead="gwresult_conf_2", 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")
Copied: grass/trunk/raster/r.solute.transport/seguin_verify_well.py (from rev 40082, grass/trunk/raster/r.solute.transport/seguin_verify_well.sh)
===================================================================
--- grass/trunk/raster/r.solute.transport/seguin_verify_well.py (rev 0)
+++ grass/trunk/raster/r.solute.transport/seguin_verify_well.py 2009-12-20 23:43:15 UTC (rev 40097)
@@ -0,0 +1,90 @@
+#!/usr/bin/env python
+# Shellscript to verify r.solute.transport calculation, this calculation is based on
+# the example 2.1 and 2.2 at page 181 of the following book:
+# title = "Str{\"o}mungs und Transportmodellierung",
+# author = "Lege, T. and Kolditz, O. and Zielke W.",
+# publisher = "Springer (Berlin; Heidelberg; New York; Barcelona; Hongkong; London; Mailand; Paris; Singapur; Tokio)",
+# edition = "2. Auflage",
+# year = "1996",
+# series = "Handbuch zur Erkundung des Untergrundes von Deponien und Altlasten"
+#
+
+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=25, res3=25, t=25, b=0, n=1000, s=0, w=0, e=2000)
+
+grass.message(_("Create all the input maps needed for groundwater flow computation"))
+
+# Initial condition of the piezometric head, we have a huge gradient from 200m to 50m
+# over a distance of 2000m
+grass.run_command("r.mapcalc", expression="phead_1=if(col() == 1 , 200, 50)")
+# Set the active cells and the dirichlet boundary condition
+grass.run_command("r.mapcalc", expression="status_1=if(col() == 1 || col() == 80 , 2, 1)")
+# We have a single well with a small influent pumping rate
+grass.run_command("r.mapcalc", expression="well_1=if((row() == 20 && col() == 20), 0.000002818287, 0)")
+# The hydraulic conductivity
+grass.run_command("r.mapcalc", expression="hydcond_1=0.0001")
+# The recharge, well we have no recharge at all
+grass.run_command("r.mapcalc", expression="recharge_1=0")
+# We have a confined aquifer with a height of 25m
+grass.run_command("r.mapcalc", expression="top_conf_1=25")
+# Bottom of the aquifer starts at 0m
+grass.run_command("r.mapcalc", expression="bottom_1=0")
+# This is the standard porosity of sand aquifer
+grass.run_command("r.mapcalc", expression="poros_1=0.11209524")
+# The specific yield of a confined aquifer
+grass.run_command("r.mapcalc", expression="syield_1=0.0001")
+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",\
+ 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")
+
+grass.message(_("generate the transport data"))
+
+# The initial concentration is zero
+grass.run_command("r.mapcalc", expression="c_1=if(col() == 20 && row() == 20 , 0.0, 0.0)")
+# No inner sources or sinks (result of chemical reaction)
+grass.run_command("r.mapcalc", expression="cs_1=0.0")
+# The pollution is inserted by a well
+grass.run_command("r.mapcalc", expression="cin_1=if(col() == 20 && row() == 20 ," + str(0.2/3600-0/24.0) + ", 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
+grass.run_command("r.mapcalc", expression="diff_1=0.0")
+# Normal retardation
+grass.run_command("r.mapcalc", expression="R_1=1.0")
+
+# Longitudinal and transversal dispersivity length
+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",\
+ 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")
+
+# The second computation uses different dispersivities
+
+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",\
+ 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")
+
More information about the grass-commit
mailing list