[GRASS-SVN] r40078 - in grass/trunk/raster: . r.solute.transport
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Dec 19 16:25:49 EST 2009
Author: huhabla
Date: 2009-12-19 16:25:49 -0500 (Sat, 19 Dec 2009)
New Revision: 40078
Added:
grass/trunk/raster/r.solute.transport/
grass/trunk/raster/r.solute.transport/Disper_calc.sh
grass/trunk/raster/r.solute.transport/Makefile
grass/trunk/raster/r.solute.transport/example.sh
grass/trunk/raster/r.solute.transport/example2.sh
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.sh
grass/trunk/raster/r.solute.transport/seguin_verify_well.sh
Log:
New module to compute solute transport in porous media, including several examples and verification shell scripts.
Added: grass/trunk/raster/r.solute.transport/Disper_calc.sh
===================================================================
--- grass/trunk/raster/r.solute.transport/Disper_calc.sh (rev 0)
+++ grass/trunk/raster/r.solute.transport/Disper_calc.sh 2009-12-19 21:25:49 UTC (rev 40078)
@@ -0,0 +1,117 @@
+#!/bin/sh
+# Shellscript to verify r.solute.transport calculation
+# The difference between solute transport calculated with and without the CFL citeria is shown
+# set the region
+g.region res=1 res3=1 t=10 b=0 n=100 s=0 w=0 e=400
+r.mapcalc "phead_1=if(col() == 1 , 55, 50)"
+r.mapcalc "status_1=if(col() == 1 || col() == 400 , 2, 1)"
+r.mapcalc "well_1=if((row() == 50 && col() == 350), 0.0, 0)"
+#r.mapcalc "well_1=0.0"
+r.mapcalc "hydcond_1=0.0001"
+r.mapcalc "reacharge_1=0"
+r.mapcalc "top_conf_1=25"
+r.mapcalc "bottom_1=0"
+r.mapcalc "poros_1=0.2"
+r.mapcalc "poros_2=0.2"
+r.mapcalc "syield_1=0.0001"
+r.mapcalc "null_1=0.0"
+
+#generate the transport data
+r.mapcalc "c_1=if(col() == 50 && row() == 50 , 0, 0.0)"
+r.mapcalc "cs_1=if(col() == 50 && row() == 50 , 0.01, 0.0)"
+r.mapcalc "cin_1=if(col() == 50 && row() == 50 , 0.0, 0.0)"
+r.mapcalc "tstatus_1=if(col() == 1, 2, 1)"
+r.mapcalc "tstatus_1=if(col() == 400 , 3, tstatus_1)"
+r.mapcalc "diff_1=0.0"
+r.mapcalc "R_1=1.0"
+
+AL=1
+AT=0.1
+export GRASS_TRUECOLOR=TRUE
+export GRASS_WIDTH=1280
+export GRASS_HEIGHT=330
+#First compute a steady state groundwater flow
+r.gwflow --o -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=reacharge_1 output=gwresult_conf_1\
+ dt=8640000000000 type=confined
+
+# create contour lines
+r.contour input=gwresult_conf_1 output=gwresult_conf_contour_1 step=1 --o
+
+###
+### compute with low velocity
+###
+
+r.solute.transport --o -cs maxit=10000 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 error=0.000000000001 \
+ diff_x=diff_1 diff_y=diff_1 cin=cin_1 c=c_1 al=$AL at=$AT velocity=stresult_conf_vel_1
+
+eval `r.univar stresult_conf_1 -g`
+
+r.mapcalc "cont_map = stresult_conf_1/$max"
+r.mapcalc "stresult_conf_vel_1_x = stresult_conf_vel_1_x *24.0 * 3600.0"
+r.info stresult_conf_vel_1_x
+r.contour input=cont_map output=stresult_conf_contour_1 step=0.05 --o
+
+#export as postscript and convert into pdf
+export GRASS_PNGFILE=Disper_calc_${AL}_${AT}_a.png
+d.mon start=PNG
+d.mon select=PNG
+d.rast cont_map
+d.vect stresult_conf_contour_1 color=grey display=shape,attr attrcol=level
+d.vect gwresult_conf_contour_1 color=blue display=shape,attr attrcol=level lcolor=blue
+d.barscale at=1,85
+echo "Solute transport 1000d al=$AL at=$AT" | d.text size=6 color=black
+d.mon stop=PNG
+convert Disper_calc_${AL}_${AT}_a.png Disper_calc_${AL}_${AT}_a.eps
+convert Disper_calc_${AL}_${AT}_a.png Disper_calc_${AL}_${AT}_a.pdf
+
+###
+### compute without CFL criteria
+###
+
+r.solute.transport --o -s loops=1 maxit=10000 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 error=0.000000000001 \
+ diff_x=diff_1 diff_y=diff_1 cin=cin_1 c=c_1 al=$AL at=$AT velocity=stresult_conf_vel_2
+
+eval `r.univar stresult_conf_2 -g`
+
+r.mapcalc "cont_map_a = stresult_conf_2/$max"
+r.mapcalc "stresult_conf_vel_2_x = stresult_conf_vel_2_x *24.0 * 3600.0"
+r.info stresult_conf_vel_2_x
+r.contour input=cont_map_a output=stresult_conf_contour_2 step=0.05 --o
+
+#export as postscript and convert into pdf
+export GRASS_PNGFILE=Disper_calc_${AL}_${AT}_b.png
+d.mon start=PNG
+d.mon select=PNG
+d.rast cont_map_a
+d.vect stresult_conf_contour_2 color=grey display=shape,attr attrcol=level
+d.vect gwresult_conf_contour_1 color=blue display=shape,attr attrcol=level lcolor=blue
+d.barscale at=1,85
+echo "Solute transport 1000d al=$AL at=$AT" | d.text size=6 color=black
+d.mon stop=PNG
+convert Disper_calc_${AL}_${AT}_b.png Disper_calc_${AL}_${AT}_b.eps
+convert Disper_calc_${AL}_${AT}_b.png Disper_calc_${AL}_${AT}_b.pdf
+
+# compute the differences between the maps
+r.mapcalc "diff=stresult_conf_1 - stresult_conf_2 "
+
+#export as postscript and convert into pdf
+export GRASS_PNGFILE=Disper_calc_${AL}_${AT}_ab_diff.png
+d.mon start=PNG
+d.mon select=PNG
+d.rast diff
+d.vect stresult_conf_contour_1 color=red display=shape,attr attrcol=level lcolor=red
+d.vect stresult_conf_contour_2 color=black display=shape,attr attrcol=level lcolor=black
+d.vect gwresult_conf_contour_1 color=blue display=shape,attr attrcol=level lcolor=blue
+d.legend at=12,17,10,55 map=diff labelnum=5 color=black
+d.barscale at=1,10
+echo "Difference between maps" | d.text size=6 color=black
+d.mon stop=PNG
+convert Disper_calc_${AL}_${AT}_ab_diff.png Disper_calc_${AL}_${AT}_ab_diff.eps
+convert Disper_calc_${AL}_${AT}_ab_diff.png Disper_calc_${AL}_${AT}_ab_diff.pdf
+
Property changes on: grass/trunk/raster/r.solute.transport/Disper_calc.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: grass/trunk/raster/r.solute.transport/Makefile
===================================================================
--- grass/trunk/raster/r.solute.transport/Makefile (rev 0)
+++ grass/trunk/raster/r.solute.transport/Makefile 2009-12-19 21:25:49 UTC (rev 40078)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM=r.solute.transport
+
+LIBES = $(GISLIB) $(RASTERLIB) $(GPDELIB) $(GMATHLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP) $(GMATHDEP) $(GPDEDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass/trunk/raster/r.solute.transport/example.sh
===================================================================
--- grass/trunk/raster/r.solute.transport/example.sh (rev 0)
+++ grass/trunk/raster/r.solute.transport/example.sh 2009-12-19 21:25:49 UTC (rev 40078)
@@ -0,0 +1,80 @@
+#!/bin/sh
+# set the region accordingly
+g.region res=1 res3=1 t=10 b=0 n=100 s=0 w=0 e=200
+r.mapcalc "phead=if(col() == 1 , 50, 40)"
+r.mapcalc "phead=if(col() ==200 , 45 + row()/40, phead)"
+r.mapcalc "status=if(col() == 1 || col() == 200 , 2, 1)"
+r.mapcalc "well=if((row() == 50 && col() == 175) || (row() == 10 && col() == 135) , -0.001, 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"
+
+#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
+
+#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=0.5 --o
+
+
+#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.0000001"
+r.mapcalc "R=1.0"
+
+r.solute.transport --o -cs 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 dt=3600 diff_x=diff diff_y=diff \
+ c=c al=0.1 at=0.01
+
+days=0
+count=0
+
+d.erase
+d.rast stresult_conf
+d.vect gwresult_conf_flow render=l
+d.vect gwresult_conf_contour color=grey render=l
+d.legend at=8,12,15,85 map=stresult_conf
+d.barscale at=1,10
+echo "Solute transport timestep: 0d" | d.text size=6 color=black
+d.out.png res=1 output=/tmp/test_0000$count.png
+
+#compute in an infinite loop the transport
+
+while true ; do
+
+r.solute.transport --o -cs 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 dt=86400 diff_x=diff diff_y=diff \
+ c=stresult_conf al=0.1 at=0.01
+
+days=`expr $days + 1`
+count=`expr $count + 1`
+
+#show the map
+d.erase
+d.rast stresult_conf
+d.vect gwresult_conf_flow render=l
+d.vect gwresult_conf_contour color=grey render=l
+d.legend at=8,12,15,85 map=stresult_conf
+d.barscale at=1,10
+echo "Solute transport timestep: ${days}d" | d.text size=6 color=black
+if [ `expr $count \< 10` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/test_0000$count.png
+elif [ `expr $count \< 100` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/test_000$count.png
+elif [ `expr $count \< 1000` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/test_00$count.png
+elif [ `expr $count \< 10000` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/test_0$count.png
+fi
+
+done
+
Property changes on: grass/trunk/raster/r.solute.transport/example.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: grass/trunk/raster/r.solute.transport/example2.sh
===================================================================
--- grass/trunk/raster/r.solute.transport/example2.sh (rev 0)
+++ grass/trunk/raster/r.solute.transport/example2.sh 2009-12-19 21:25:49 UTC (rev 40078)
@@ -0,0 +1,121 @@
+#!/bin/sh
+# set the region accordingly
+g.region res=1 res3=1 t=10 b=0 n=100 s=0 w=0 e=200
+# create groundwater flow data
+r.mapcalc "phead=if(col() == 1 , 50, 45)"
+#r.mapcalc "phead=if(col() ==200 , 45 + row()/40, phead)"
+r.mapcalc "status=if(col() == 1 || col() == 200 , 2, 1)"
+r.mapcalc "well=if((row() == 50 && col() == 175) || (row() == 10 && col() == 135) , -0.005, 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.1"
+r.mapcalc "null=0.0"
+
+#generate the transport data
+r.mapcalc "c=if(col() == 15 && row() == 75 , 0.1, 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.0000001"
+r.mapcalc "R=1.0"
+
+#dispersivity lengths
+AL="0.1"
+AT="0.01"
+
+#First compute a steady state groundwater flow
+r.gwflow --o -s solver=cg top=top_unconf bottom=bottom phead=phead status=status hc_x=hydcond hc_y=hydcond \
+ q=null s=syield r=reacharge output=gwresult_conf dt=8640000000000 type=unconfined
+
+#Second initiate the solute transport
+r.solute.transport --o -cs solver=bicgstab top=gwresult_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 dt=30 diff_x=diff diff_y=diff \
+ c=c al=$AL at=$AT
+
+#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=0.5 --o
+
+d.erase
+d.rast stresult_conf
+d.vect gwresult_conf_flow
+d.vect gwresult_conf_contour color=grey
+d.legend at=8,12,15,85 map=stresult_conf
+d.barscale at=1,10
+echo "Solute transport timestep: 0d" | d.text size=6 color=black
+d.out.png res=1 output=/tmp/test_0000$count.png
+
+days=0
+count=0
+SINK=well
+
+#compute in an infinite loop the transport
+while true ; do
+
+
+#First compute groundwater flow
+r.gwflow --o -s solver=cg top=top_conf bottom=bottom phead=gwresult_conf status=status hc_x=hydcond \
+ hc_y=hydcond q=$SINK s=syield r=reacharge output=gwresult_conf dt=86400 type=unconfined
+
+#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=0.5 --o
+
+#Second compute soulte transport
+r.solute.transport --o -s stab=exp 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 dt=86400 diff_x=diff diff_y=diff \
+ c=stresult_conf al=$AL at=$AT
+
+days=`expr $days + 1`
+count=`expr $count + 1`
+
+if [ `expr $count \< 76` -eq 1 ] || [ `expr $count \> 151` -eq 1 ]; then
+SINK=well
+else
+SINK=null
+fi
+
+d.erase
+d.rast gwresult_conf
+d.vect gwresult_conf_flow
+d.vect gwresult_conf_contour color=grey
+d.legend at=7,11,15,85 map=gwresult_conf
+d.barscale at=1,10 &
+echo "Groundwater flow timestep: ${days}d" | d.text size=6 color=black
+if [ `expr $count \< 10` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/gwflow_0000$count.png
+elif [ `expr $count \< 100` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/gwflow_000$count.png
+elif [ `expr $count \< 1000` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/gwflow_00$count.png
+elif [ `expr $count \< 10000` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/gwflow_0$count.png
+fi
+
+eval `r.univar stresult_conf -g`
+
+d.erase
+d.rast stresult_conf
+d.vect gwresult_conf_flow
+d.vect gwresult_conf_contour color=grey
+d.legend at=7,11,15,85 map=stresult_conf
+d.barscale at=1,10
+echo "Solute transport timestep: ${days}d" | d.text size=6 color=black
+echo "Quantity: $sum min-max: $min - $max" | d.text size=4 color=black at=1,12
+if [ `expr $count \< 10` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/transport_0000$count.png
+elif [ `expr $count \< 100` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/transport_000$count.png
+elif [ `expr $count \< 1000` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/transport_00$count.png
+elif [ `expr $count \< 10000` -eq 1 ] ; then
+ d.out.png res=1 output=/tmp/transport_0$count.png
+fi
+
+done
+
Property changes on: grass/trunk/raster/r.solute.transport/example2.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: grass/trunk/raster/r.solute.transport/main.c
===================================================================
--- grass/trunk/raster/r.solute.transport/main.c (rev 0)
+++ grass/trunk/raster/r.solute.transport/main.c 2009-12-19 21:25:49 UTC (rev 40078)
@@ -0,0 +1,582 @@
+
+/****************************************************************************
+*
+* MODULE: r.solute.transport
+*
+* AUTHOR(S): Original author
+* Soeren Gebbert soerengebbert <at> gmx <dot> de
+* 27 11 2006 Berlin
+* PURPOSE: Calculates transient two dimensional solute transport
+* in porous media
+*
+* COPYRIGHT: (C) 2006 by the GRASS Development Team
+*
+* This program is free software under the GNU General Public
+* License (>=v2). Read the file COPYING that comes with GRASS
+* for details.
+*
+*****************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include <grass/gmath.h>
+#include <grass/N_pde.h>
+#include <grass/N_solute_transport.h>
+
+
+/*- Parameters and global variables -----------------------------------------*/
+typedef struct
+{
+ 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,
+ *al, *at, *loops, *stab;
+ struct Flag *sparse;
+ struct Flag *cfl;
+} paramType;
+
+paramType param; /*Parameters */
+
+/*- prototypes --------------------------------------------------------------*/
+void set_params(); /*Fill the paramType structure */
+void copy_result(N_array_2d * status, N_array_2d * c_start, double *result,
+ struct Cell_head *region, N_array_2d * target, int tflag);
+N_les *create_solve_les(N_geom_data * geom, N_solute_transport_data2d * data,
+ N_les_callback_2d * call, const char *solver, int maxit,
+ double error, double sor);
+
+/* ************************************************************************* */
+/* Set up the arguments we are expecting ********************************** */
+/* ************************************************************************* */
+void set_params()
+{
+ param.c = G_define_option();
+ param.c->key = "c";
+ param.c->type = TYPE_STRING;
+ param.c->required = YES;
+ param.c->gisprompt = "old,raster,raster";
+ param.c->description = _("The initial concentration in [kg/m^3]");
+
+ param.phead = G_define_option();
+ param.phead->key = "phead";
+ param.phead->type = TYPE_STRING;
+ param.phead->required = YES;
+ param.phead->gisprompt = "old,raster,raster";
+ param.phead->description = _("The piezometric head in [m]");
+
+ param.hc_x = G_define_option();
+ param.hc_x->key = "hc_x";
+ param.hc_x->type = TYPE_STRING;
+ param.hc_x->required = YES;
+ param.hc_x->gisprompt = "old,raster,raster";
+ param.hc_x->description =
+ _("The x-part of the hydraulic conductivity tensor in [m/s]");
+
+ param.hc_y = G_define_option();
+ param.hc_y->key = "hc_y";
+ param.hc_y->type = TYPE_STRING;
+ param.hc_y->required = YES;
+ param.hc_y->gisprompt = "old,raster,raster";
+ param.hc_y->description =
+ _("The y-part of the hydraulic conductivity tensor in [m/s]");
+
+
+ param.status = G_define_option();
+ param.status->key = "status";
+ param.status->type = TYPE_STRING;
+ param.status->required = YES;
+ param.status->gisprompt = "old,raster,raster";
+ param.status->description =
+ _
+ ("The status for each cell, = 0 - inactive cell, 1 - active cell, 2 - dirichlet- and 3 - transfer boundary condition");
+
+ param.diff_x = G_define_option();
+ param.diff_x->key = "diff_x";
+ param.diff_x->type = TYPE_STRING;
+ param.diff_x->required = YES;
+ param.diff_x->gisprompt = "old,raster,raster";
+ param.diff_x->description =
+ _("The x-part of the diffusion tensor in [m�/s]");
+
+ param.diff_y = G_define_option();
+ param.diff_y->key = "diff_y";
+ param.diff_y->type = TYPE_STRING;
+ param.diff_y->required = YES;
+ param.diff_y->gisprompt = "old,raster,raster";
+ param.diff_y->description =
+ _("The y-part of the diffusion tensor in [m�/s]");
+
+ param.q = G_define_option();
+ param.q->key = "q";
+ param.q->type = TYPE_STRING;
+ param.q->required = NO;
+ param.q->gisprompt = "old,raster,raster";
+ param.q->description = _("groundwater sources and sinks in [m^3/s]");
+
+ param.cin = G_define_option();
+ param.cin->key = "cin";
+ param.cin->type = TYPE_STRING;
+ param.cin->required = NO;
+ param.cin->gisprompt = "old,raster,raster";
+ param.cin->description = _("concentration sources and sinks in [kg/m^3]");
+
+ param.cs = G_define_option();
+ param.cs->key = "cs";
+ param.cs->type = TYPE_STRING;
+ param.cs->required = YES;
+ param.cs->gisprompt = "old,raster,raster";
+ param.cs->description = _("concentration sources and sinks in [kg/m^3]");
+
+ param.r = G_define_option();
+ param.r->key = "R";
+ param.r->type = TYPE_STRING;
+ param.r->required = YES;
+ param.r->gisprompt = "old,raster,raster";
+ param.r->description = _("Retardation factor [-]");
+
+ param.nf = G_define_option();
+ param.nf->key = "nf";
+ param.nf->type = TYPE_STRING;
+ param.nf->required = YES;
+ param.nf->gisprompt = "old,raster,raster";
+ param.nf->description = _("Effective porosity [-]");
+
+ param.top = G_define_option();
+ param.top->key = "top";
+ param.top->type = TYPE_STRING;
+ param.top->required = YES;
+ param.top->gisprompt = "old,raster,raster";
+ param.top->description = _("Top surface of the aquifer in [m]");
+
+ param.bottom = G_define_option();
+ param.bottom->key = "bottom";
+ param.bottom->type = TYPE_STRING;
+ param.bottom->required = YES;
+ param.bottom->gisprompt = "old,raster,raster";
+ param.bottom->description = _("Bottom surface of the aquifer in [m]");
+
+ param.output = G_define_option();
+ param.output->key = "output";
+ param.output->type = TYPE_STRING;
+ param.output->required = YES;
+ param.output->gisprompt = "new,raster,raster";
+ param.output->description =
+ _
+ ("The result of the numericalsolute transport calculation will be written to this map. [kg/m�]");
+
+ param.vector = G_define_option();
+ param.vector->key = "velocity";
+ param.vector->type = TYPE_STRING;
+ param.vector->required = NO;
+ param.vector->gisprompt = "new,raster,raster";
+ param.vector->description =
+ _
+ ("Calculate the groundwater 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);
+ param.error = N_define_standard_option(N_OPT_ITERATION_ERROR);
+ param.solver = N_define_standard_option(N_OPT_SOLVER_UNSYMM);
+ param.sor = N_define_standard_option(N_OPT_SOR_VALUE);
+
+ param.al = G_define_option();
+ param.al->key = "al";
+ param.al->type = TYPE_DOUBLE;
+ param.al->required = NO;
+ param.al->answer = "0.0";
+ param.al->description =
+ _("The longditudinal dispersivity length. [m]");
+
+ param.at = G_define_option();
+ param.at->key = "at";
+ param.at->type = TYPE_DOUBLE;
+ param.at->required = NO;
+ param.at->answer = "0.0";
+ param.at->description =
+ _("The transversal dispersivity length. [m]");
+
+ param.loops = G_define_option();
+ param.loops->key = "loops";
+ param.loops->type = TYPE_DOUBLE;
+ param.loops->required = NO;
+ param.loops->answer = "1";
+ param.loops->description =
+ _("Use this number of time loops if the CFL flag is off. The timestep will become dt/loops.");
+
+ param.stab = G_define_option();
+ param.stab->key = "stab";
+ param.stab->type = TYPE_STRING;
+ param.stab->required = NO;
+ param.stab->answer = "full";
+ param.stab->options = "full,exp";
+ param.stab->description =
+ _("Set the flow stabilizing scheme.");
+
+
+ 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 =
+ _
+ ("Use the Courant-Friedrichs-Lewy criteria for time step calculation");
+
+
+}
+
+/* ************************************************************************* */
+/* Main function *********************************************************** */
+/* ************************************************************************* */
+int main(int argc, char *argv[])
+{
+ struct GModule *module = NULL;
+ N_solute_transport_data2d *data = NULL;
+ N_geom_data *geom = NULL;
+ N_les *les = NULL;
+ N_les_callback_2d *call = NULL;
+ struct Cell_head region;
+ double error, sor;
+ char *solver;
+ int x, y, stat, i, maxit = 1;
+ double loops = 1;
+ N_array_2d *xcomp = NULL;
+ N_array_2d *ycomp = NULL;
+ 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;
+
+ /* Initialize GRASS */
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("solute transport"));
+ module->description =
+ _
+ ("Numerical calculation program for transient, confined and unconfined solute transport in two dimensions");
+
+ /* Get parameters from user */
+ set_params();
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /*Set the maximum iterations */
+ sscanf(param.maxit->answer, "%i", &(maxit));
+ /*Set the calculation error break criteria */
+ sscanf(param.error->answer, "%lf", &(error));
+ sscanf(param.sor->answer, "%lf", &(sor));
+ /*number of loops*/
+ sscanf(param.loops->answer, "%lf", &(loops));
+ /*Set the solver */
+ solver = param.solver->answer;
+
+ if (strcmp(solver, G_MATH_SOLVER_DIRECT_LU) == 0 && param.sparse->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)
+ G_fatal_error(_
+ ("The direct Gauss solver do not work with sparse matrices"));
+
+
+ /*get the current region */
+ G_get_set_window(®ion);
+
+ /*allocate the geometry structure for geometry and area calculation */
+ geom = N_init_geom_data_2d(®ion, geom);
+
+ /*Set the function callback to the groundwater flow function */
+ call = N_alloc_les_callback_2d();
+ N_set_les_callback_2d_func(call, (*N_callback_solute_transport_2d)); /*solute_transport 2d */
+
+ /*Allocate the groundwater flow data structure */
+ data = N_alloc_solute_transport_data2d(geom->cols, geom->rows);
+
+ /*Set the stabilizing scheme*/
+ if (strncmp("full", param.stab->answer, 4) == 0) {
+ data->stab = N_UPWIND_FULL;
+ }
+ if (strncmp("exp", param.stab->answer, 3) == 0) {
+ data->stab = N_UPWIND_EXP;
+ }
+
+ /*the dispersivity lengths*/
+ sscanf(param.al->answer, "%lf", &(data->al));
+ sscanf(param.at->answer, "%lf", &(data->at));
+
+ /*Set the calculation time */
+ sscanf(param.dt->answer, "%lf", &(data->dt));
+
+ /*read all input maps into the memory and take care of the
+ * null values.*/
+ N_read_rast_to_array_2d(param.c->answer, data->c);
+ N_convert_array_2d_null_to_zero(data->c);
+ N_read_rast_to_array_2d(param.c->answer, data->c_start);
+ N_convert_array_2d_null_to_zero(data->c_start);
+ N_read_rast_to_array_2d(param.status->answer, data->status);
+ N_convert_array_2d_null_to_zero(data->status);
+ N_read_rast_to_array_2d(param.diff_x->answer, data->diff_x);
+ N_convert_array_2d_null_to_zero(data->diff_x);
+ N_read_rast_to_array_2d(param.diff_y->answer, data->diff_y);
+ N_convert_array_2d_null_to_zero(data->diff_y);
+ N_read_rast_to_array_2d(param.q->answer, data->q);
+ N_convert_array_2d_null_to_zero(data->q);
+ N_read_rast_to_array_2d(param.nf->answer, data->nf);
+ N_convert_array_2d_null_to_zero(data->nf);
+ N_read_rast_to_array_2d(param.cs->answer, data->cs);
+ N_convert_array_2d_null_to_zero(data->cs);
+ N_read_rast_to_array_2d(param.top->answer, data->top);
+ N_convert_array_2d_null_to_zero(data->top);
+ N_read_rast_to_array_2d(param.bottom->answer, data->bottom);
+ N_convert_array_2d_null_to_zero(data->bottom);
+ N_read_rast_to_array_2d(param.r->answer, data->R);
+ N_convert_array_2d_null_to_zero(data->R);
+
+ if(param.cin->answer) {
+ N_read_rast_to_array_2d(param.cin->answer, data->cin);
+ N_convert_array_2d_null_to_zero(data->cin);
+ }
+
+ /*initiate the values for velocity calculation*/
+ hc_x = N_alloc_array_2d(geom->cols, geom->rows, 1, DCELL_TYPE);
+ hc_x = N_read_rast_to_array_2d(param.hc_x->answer, hc_x);
+ N_convert_array_2d_null_to_zero(hc_x);
+ hc_y = N_alloc_array_2d(geom->cols, geom->rows, 1, DCELL_TYPE);
+ hc_y = N_read_rast_to_array_2d(param.hc_y->answer, hc_y);
+ N_convert_array_2d_null_to_zero(hc_y);
+ phead = N_alloc_array_2d(geom->cols, geom->rows, 1, DCELL_TYPE);
+ phead = N_read_rast_to_array_2d(param.phead->answer, phead);
+ N_convert_array_2d_null_to_zero(phead);
+
+ /* Set the inactive values to zero, to assure a no flow boundary */
+ for (y = 0; y < geom->rows; y++) {
+ for (x = 0; x < geom->cols; x++) {
+ stat = (int)N_get_array_2d_d_value(data->status, x, y);
+ if (stat == N_CELL_INACTIVE) { /*only inactive cells */
+ N_put_array_2d_d_value(data->diff_x, x, y, 0);
+ N_put_array_2d_d_value(data->diff_y, x, y, 0);
+ N_put_array_2d_d_value(data->cs, x, y, 0);
+ N_put_array_2d_d_value(data->q, x, y, 0);
+ }
+ }
+ }
+
+ /*compute the velocities */
+ N_math_array_2d(hc_x, data->nf, hc_x, N_ARRAY_DIV);
+ 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);
+
+ /***************************************/
+ /*the Courant-Friedrichs-Lewy criteria */
+ /*Compute the correct time step */
+ if (geom->dx > geom->dy)
+ length = geom->dx;
+ else
+ length = geom->dy;
+
+ if (fabs(data->grad->max) > fabs(data->grad->min)) {
+ cfl = (double)data->dt * fabs(data->grad->max) / length;
+ time_step = 1*length / fabs(data->grad->max);
+ }
+ else {
+ cfl = (double)data->dt * fabs(data->grad->min) / length;
+ time_step = 1*length / fabs(data->grad->min);
+ }
+
+ G_message("The Courant-Friedrichs-Lewy criteria is %g it should be within [0:1]", cfl);
+ G_message("The largest stable time step is %g", time_step);
+
+ /*Set the number of inner loops and the time step*/
+ if (data->dt > time_step && param.cfl->answer) {
+ /*safe the user time step */
+ time_sum = data->dt;
+ time_loops = data->dt / time_step;
+ time_loops = floor(time_loops) + 1;
+ data->dt = data->dt / time_loops;
+ G_message("Number of inner loops is %g", time_loops);
+ G_message("Time step for each loop %g", data->dt);
+ }
+ else {
+ if(data->dt > time_step)
+ G_warning("The time step is to large: %gs. The largest time step should be of size %gs.", data->dt, time_step);
+
+ time_loops = loops;
+ data->dt = data->dt / loops;
+ }
+
+ N_free_array_2d(phead);
+ N_free_array_2d(hc_x);
+ N_free_array_2d(hc_y);
+
+ /*Compute for each time step*/
+ for (i = 0; i < time_loops; i++) {
+ G_message("Time step %i with time sum %g", i + 1, (i + 1)*data->dt);
+
+ /*assemble the linear equation system and solve it */
+ les = create_solve_les(geom, data, call, solver, maxit, error, sor);
+
+ /* copy the result into the c array for output */
+ copy_result(data->status, data->c_start, les->x, ®ion, data->c, 1);
+ N_convert_array_2d_null_to_zero(data->c_start);
+
+ if (les)
+ N_free_les(les);
+
+ /*Set the start array*/
+ N_copy_array_2d(data->c, data->c_start);
+ /*Set the transmission boundary*/
+ N_calc_solute_transport_transmission_2d(data);
+
+ }
+
+ /*write the result to the output file */
+ 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) {
+ 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 (xcomp)
+ N_free_array_2d(xcomp);
+ if (ycomp)
+ N_free_array_2d(ycomp);
+ }
+
+
+ if (data)
+ N_free_solute_transport_data2d(data);
+ if (geom)
+ N_free_geom_data(geom);
+ if (call)
+ G_free(call);
+
+ return (EXIT_SUCCESS);
+}
+
+
+/* ************************************************************************* */
+/* this function copies the result from the x vector to a N_array_2d array * */
+/* ************************************************************************* */
+void
+copy_result(N_array_2d * status, N_array_2d * c_start, double *result,
+ struct Cell_head *region, N_array_2d * target, int tflag)
+{
+ 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(status, x, y);
+ if (stat == N_CELL_ACTIVE) { /*only active cells */
+ d1 = result[count];
+ val = (DCELL) d1;
+ count++;
+ }
+ else if (stat == N_CELL_DIRICHLET) { /*dirichlet cells */
+ d1 = N_get_array_2d_d_value(c_start, x, y);
+ val = (DCELL) d1;
+ }
+ else if (tflag == 1 && stat == N_CELL_TRANSMISSION) {/*transmission cells*/
+ d1 = N_get_array_2d_d_value(c_start, x, y);
+ val = (DCELL) d1;
+ }
+ else {
+ Rast_set_null_value(&val, 1, DCELL_TYPE);
+ }
+ N_put_array_2d_d_value(target, x, y, val);
+ }
+ }
+
+ return;
+}
+
+/* *************************************************************** */
+/* ***** create and solve the linear equation system ************* */
+/* *************************************************************** */
+N_les *create_solve_les(N_geom_data * geom, N_solute_transport_data2d * data,
+ N_les_callback_2d * call, const char *solver, int maxit,
+ double error, double sor)
+{
+
+ N_les *les;
+
+ /*assemble the linear equation system */
+ if (param.sparse->answer)
+ les =
+ N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c,
+ (void *)data, call);
+ else
+ les =
+ N_assemble_les_2d(N_NORMAL_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)
+ 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);
+ }
+
+ if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_SOR) == 0)
+ {
+ if (param.sparse->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);
+ }
+
+ if (strcmp(solver, G_MATH_SOLVER_ITERATIVE_BICGSTAB) == 0)
+ {
+ if (param.sparse->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);
+ }
+
+ if (strcmp(solver, G_MATH_SOLVER_DIRECT_LU) == 0)
+ G_math_solver_lu(les->A, les->x, les->b, les->rows);
+
+ if (strcmp(solver, G_MATH_SOLVER_DIRECT_GAUSS) == 0)
+ G_math_solver_gauss(les->A, les->x, les->b, les->rows);
+
+ if (les == NULL)
+ G_fatal_error(_
+ ("Could not create and solve the linear equation system"));
+
+ return les;
+}
+
Added: grass/trunk/raster/r.solute.transport/r.solute.transport.html
===================================================================
--- grass/trunk/raster/r.solute.transport/r.solute.transport.html (rev 0)
+++ grass/trunk/raster/r.solute.transport/r.solute.transport.html 2009-12-19 21:25:49 UTC (rev 40078)
@@ -0,0 +1,123 @@
+<H2>DESCRIPTION</H2>
+This numerical program calculates transient
+solute transport in two dimensions based on
+raster maps and the current region resolution.
+All starting- and boundary-conditions must be provided as
+raster maps.
+<br>
+<br>
+This module calculates the concentration and optional the
+velocity field, based on the hydraulic conductivity,
+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>
+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).
+
+
+<H2>NOTES</H2>
+The solute transport partial
+differential equation is of the following form:
+
+<p>
+(dc/dt)*R = div ( D grad c - uc) + cs -q/nf(c - c_in)
+
+<ul>
+<li>c -- the concentration </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>
+</ul>
+
+<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:
+
+<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>
+<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>
+</ul>
+
+<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.
+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).
+
+<H2>EXAMPLE</H2>
+Use this small script to create a working
+groundwater flow / solute transport area and data.
+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"
+
+#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
+
+#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
+
+#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"
+
+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
+
+
+#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
+
+</pre></div>
+
+<H2>SEE ALSO</H2>
+
+<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
+
+<p><i>Last changed: $Date: 2007/02/15 23:27:58 $</i>
Added: grass/trunk/raster/r.solute.transport/seguin_verify.sh
===================================================================
--- grass/trunk/raster/r.solute.transport/seguin_verify.sh (rev 0)
+++ grass/trunk/raster/r.solute.transport/seguin_verify.sh 2009-12-19 21:25:49 UTC (rev 40078)
@@ -0,0 +1,149 @@
+#!/bin/sh
+# Shellscript to verify r.solute.transport calculation
+# The difference between solute transport calculated with and without the CFL citeria is shown
+# set the region
+g.region res=1 res3=1 t=10 b=0 n=100 s=0 w=0 e=400
+r.mapcalc "phead_1=if(col() == 1 , 53, 50)"
+r.mapcalc "status_1=if(col() == 1 || col() == 400 , 2, 1)"
+r.mapcalc "well_1=if((row() == 50 && col() == 350), 0.0, 0)"
+#r.mapcalc "well_1=0.0"
+r.mapcalc "hydcond_1=0.0001"
+r.mapcalc "reacharge_1=0"
+r.mapcalc "top_conf_1=25"
+r.mapcalc "bottom_1=0"
+r.mapcalc "poros_1=0.2"
+r.mapcalc "poros_2=0.2"
+r.mapcalc "syield_1=0.0001"
+r.mapcalc "null_1=0.0"
+
+#generate the transport data
+r.mapcalc "c_1=if(col() == 50 && row() == 50 , 0, 0.0)"
+r.mapcalc "cs_1=if(col() == 50 && row() == 50 , 0.01, 0.0)"
+r.mapcalc "cin_1=if(col() == 50 && row() == 50 , 0.0, 0.0)"
+r.mapcalc "tstatus_1=if(col() == 1, 2, 1)"
+r.mapcalc "tstatus_1=if(col() == 400 , 3, tstatus_1)"
+r.mapcalc "diff_1=0.0"
+r.mapcalc "R_1=1.0"
+
+AL=10
+AT=1
+LOOPS=2
+
+#First compute a steady state groundwater flow
+r.gwflow --o -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=reacharge_1 output=gwresult_conf_1\
+ dt=8640000000000 type=confined
+
+# create contour lines
+r.contour input=gwresult_conf_1 output=gwresult_conf_contour_1 step=1 --o
+
+###
+### Compute with full upwind
+###
+
+r.solute.transport --o -s loops=$LOOPS stab=full maxit=10000 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 error=0.000000000001 \
+ diff_x=diff_1 diff_y=diff_1 cin=cin_1 c=c_1 al=$AL at=$AT velocity=stresult_conf_vel_1
+
+eval `r.univar stresult_conf_1 -g`
+
+r.mapcalc "cont_map = stresult_conf_1/$max"
+r.mapcalc "stresult_conf_vel_1_x = stresult_conf_vel_1_x *24.0 * 3600.0"
+r.info stresult_conf_vel_1_x
+r.contour input=cont_map output=stresult_conf_contour_1 step=0.02 --o
+
+#export as postscript and convert into pdf
+export GRASS_TRUECOLOR=TRUE
+export GRASS_PSFILE=Disper_calc_${AL}_${AT}_a.ps
+d.mon start=PS
+d.mon select=PS
+d.rast cont_map
+d.vect stresult_conf_contour_1 color=grey display=shape,attr attrcol=level
+d.vect gwresult_conf_contour_1 color=blue display=shape,attr attrcol=level lcolor=blue
+d.barscale at=1,85
+echo "Solute transport 1000d al=$AL at=$AT" | d.text size=4 color=black
+d.mon stop=PS
+ps2epsi Disper_calc_${AL}_${AT}_a.ps
+epstopdf Disper_calc_${AL}_${AT}_a.epsi
+
+###
+### compute with exp upwind
+###
+
+r.solute.transport --o -s loops=$LOOPS stab=exp maxit=10000 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 error=0.000000000001 \
+ diff_x=diff_1 diff_y=diff_1 cin=cin_1 c=c_1 al=$AL at=$AT velocity=stresult_conf_vel_2
+
+eval `r.univar stresult_conf_2 -g`
+
+r.mapcalc "cont_map_a = stresult_conf_2/$max"
+r.mapcalc "stresult_conf_vel_2_x = stresult_conf_vel_2_x *24.0 * 3600.0"
+r.info stresult_conf_vel_2_x
+r.contour input=cont_map_a output=stresult_conf_contour_2 step=0.02 --o
+
+#export as postscript and convert into pdf
+export GRASS_TRUECOLOR=TRUE
+export GRASS_PSFILE=Disper_calc_${AL}_${AT}_b.ps
+d.mon start=PS
+d.mon select=PS
+d.rast cont_map_a
+d.vect stresult_conf_contour_2 color=grey display=shape,attr attrcol=level
+d.vect gwresult_conf_contour_1 color=blue display=shape,attr attrcol=level lcolor=blue
+d.barscale at=1,85
+echo "Solute transport 1000d al=$AL at=$AT" | d.text size=4 color=black
+d.mon stop=PS
+ps2epsi Disper_calc_${AL}_${AT}_b.ps
+epstopdf Disper_calc_${AL}_${AT}_b.epsi
+
+###
+### compute with weight upwind
+###
+
+r.solute.transport --o -s loops=$LOOPS stab=weight maxit=10000 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_3 dt=86400000 error=0.000000000001 \
+ diff_x=diff_1 diff_y=diff_1 cin=cin_1 c=c_1 al=$AL at=$AT velocity=stresult_conf_vel_3
+
+eval `r.univar stresult_conf_2 -g`
+
+r.mapcalc "cont_map_a = stresult_conf_3/$max"
+r.mapcalc "stresult_conf_vel_3_x = stresult_conf_vel_3_x *24.0 * 3600.0"
+r.info stresult_conf_vel_3_x
+r.contour input=cont_map_a output=stresult_conf_contour_3 step=0.02 --o
+
+#export as postscript and convert into pdf
+export GRASS_TRUECOLOR=TRUE
+export GRASS_PSFILE=Disper_calc_${AL}_${AT}_b.ps
+d.mon start=PS
+d.mon select=PS
+d.rast cont_map_a
+d.vect stresult_conf_contour_3 color=grey display=shape,attr attrcol=level
+d.vect gwresult_conf_contour_1 color=blue display=shape,attr attrcol=level lcolor=blue
+d.barscale at=1,85
+echo "Solute transport with weight upwind 1000d al=$AL at=$AT" | d.text size=4 color=black
+d.mon stop=PS
+ps2epsi Disper_calc_${AL}_${AT}_c.ps
+epstopdf Disper_calc_${AL}_${AT}_c.epsi
+
+# compute the differences between the maps
+r.mapcalc "diff=stresult_conf_2 - stresult_conf_3 "
+
+#export as postscript and convert into pdf
+export GRASS_TRUECOLOR=TRUE
+export GRASS_PSFILE=Disper_calc_${AL}_${AT}_c.ps
+d.mon start=PS
+d.mon select=PS
+d.rast diff
+d.vect stresult_conf_contour_2 color=red display=shape,attr attrcol=level lcolor=red
+d.vect stresult_conf_contour_3 color=black display=shape,attr attrcol=level lcolor=black
+d.vect gwresult_conf_contour_1 color=blue display=shape,attr attrcol=level lcolor=blue
+d.legend at=8,12,15,85 map=diff labelnum=5 color=black
+d.barscale at=1,10
+echo "Difference between maps" | d.text size=6 color=black
+d.mon stop=PS
+ps2epsi Disper_calc_${AL}_${AT}_bc_diff.ps
+epstopdf Disper_calc_${AL}_${AT}_bc_diff.epsi
+
Property changes on: grass/trunk/raster/r.solute.transport/seguin_verify.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: grass/trunk/raster/r.solute.transport/seguin_verify_well.sh
===================================================================
--- grass/trunk/raster/r.solute.transport/seguin_verify_well.sh (rev 0)
+++ grass/trunk/raster/r.solute.transport/seguin_verify_well.sh 2009-12-19 21:25:49 UTC (rev 40078)
@@ -0,0 +1,105 @@
+#!/bin/sh
+# 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"
+#
+# set the region
+g.region res=50 res3=1 t=10 b=0 n=1000 s=0 w=0 e=2000
+r.mapcalc "phead_1=if(col() == 1 , 200, 50)"
+r.mapcalc "status_1=if(col() == 1 || col() == 40 , 2, 1)"
+r.mapcalc "well_1=if((row() == 10 && col() == 10), 0.000002818287, 0)"
+r.mapcalc "well_1=0.0"
+r.mapcalc "hydcond_1=0.0001"
+r.mapcalc "reacharge_1=0"
+r.mapcalc "top_conf_1=25"
+r.mapcalc "bottom_1=0"
+r.mapcalc "poros_1=0.11209524"
+r.mapcalc "syield_1=0.0001"
+r.mapcalc "null_1=0.0"
+
+#generate the transport data
+r.mapcalc "c_1=if(col() == 10 && row() == 10 , 0.0, 0.0)"
+r.mapcalc "cs_1=if(col() == 10 && row() == 10 , 0.000000000000000000001, 0.0)"
+r.mapcalc "cin_1=if(col() == 10 && row() == 10 , 0.00000231481, 0.0)"
+r.mapcalc "tstatus_1=if(col() == 1 || col() == 40 , 3, 1)"
+r.mapcalc "diff_1=0.0"
+r.mapcalc "R_1=1.0"
+
+AL=50
+AT=5
+
+#First compute a steady state groundwater flow
+r.gwflow --o -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=reacharge_1 output=gwresult_conf_1\
+ dt=8640000000000 type=confined
+
+# create contour lines
+r.contour input=gwresult_conf_1 output=gwresult_conf_contour_1 step=10 --o
+
+###
+### compute with low velocity
+###
+
+r.solute.transport --o -cs 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
+
+eval `r.univar stresult_conf_1 -g`
+
+#export as postscript and convert into pdf
+export GRASS_TRUECOLOR=TRUE
+export GRASS_PSFILE=seguin_verify_c.ps
+d.mon start=PS
+d.mon select=PS
+d.rast cont_map
+#d.rast.num cont_map dp=2
+d.vect stresult_conf_contour_1 color=grey display=shape,attr attrcol=level
+d.legend at=8,12,15,85 map=cont_map labelnum=5 color=black thin=500
+d.barscale at=1,10
+echo "Solute transport 250d al=$AL at=$AT" | d.text size=6 color=black
+d.mon stop=PS
+ps2epsi seguin_verify_c.ps
+epstopdf seguin_verify_c.epsi
+
+###
+### compute different dispersion
+###
+
+AL=10
+AT=1
+
+r.solute.transport --o -cs 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
+
+eval `r.univar stresult_conf_2 -g`
+
+r.mapcalc "cont_map_a = stresult_conf_2/$max"
+r.mapcalc "stresult_conf_vel_2_x = stresult_conf_vel_2_x *24.0 * 3600.0"
+r.info stresult_conf_vel_2_x
+r.contour input=cont_map_a output=stresult_conf_contour_2 step=0.1 --o
+
+
+#export as postscript and convert into pdf
+export GRASS_TRUECOLOR=TRUE
+export GRASS_PSFILE=seguin_verify_d.ps
+d.mon start=PS
+d.mon select=PS
+d.rast cont_map_a
+#d.rast.num cont_map dp=2
+d.vect stresult_conf_contour_2 color=grey display=shape,attr attrcol=level
+d.legend at=8,12,15,85 map=cont_map_a labelnum=5 color=black thin=500
+d.barscale at=1,10
+echo "Solute transport 250d al=$AL at=$AT" | d.text size=6 color=black
+d.mon stop=PS
+ps2epsi seguin_verify_d.ps
+epstopdf seguin_verify_d.epsi
+
Property changes on: grass/trunk/raster/r.solute.transport/seguin_verify_well.sh
___________________________________________________________________
Added: svn:executable
+ *
More information about the grass-commit
mailing list