[GRASS-SVN] r68320 - in grass/trunk/raster/r.sim: . r.sim.sediment r.sim.water simlib test
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Apr 27 11:45:35 PDT 2016
Author: wenzeslaus
Date: 2016-04-27 11:45:35 -0700 (Wed, 27 Apr 2016)
New Revision: 68320
Added:
grass/trunk/raster/r.sim/simlib/simlib.h
grass/trunk/raster/r.sim/test/
grass/trunk/raster/r.sim/test/test.sh
Modified:
grass/trunk/raster/r.sim/r.sim.sediment/main.c
grass/trunk/raster/r.sim/r.sim.water/main.c
grass/trunk/raster/r.sim/simlib/Makefile
grass/trunk/raster/r.sim/simlib/hydro.c
grass/trunk/raster/r.sim/simlib/input.c
grass/trunk/raster/r.sim/simlib/observation_points.c
grass/trunk/raster/r.sim/simlib/output.c
grass/trunk/raster/r.sim/simlib/waterglobs.h
Log:
r.sim.water and r.sim.sediment: do not use lib globals outside of the lib
Using library global variables outside of a dynamically linked
library may not be portable (MS Windows). See #2966.
Use structure and a setter function to set all the library
variables. Some other calls and structures needed to be
moved around to not use the library globals from main
functions of individual modules. New header added which
is the interface for the SIMWE library.
Modified: grass/trunk/raster/r.sim/r.sim.sediment/main.c
===================================================================
--- grass/trunk/raster/r.sim/r.sim.sediment/main.c 2016-04-26 15:38:36 UTC (rev 68319)
+++ grass/trunk/raster/r.sim/r.sim.sediment/main.c 2016-04-27 18:45:35 UTC (rev 68320)
@@ -81,7 +81,7 @@
/********************************/
-#include <grass/waterglobs.h>
+#include <grass/simlib.h>
char fncdsm[32];
char filnam[10];
@@ -102,6 +102,9 @@
int ii;
int ret_val;
static int rand1 = 12345;
+ struct Cell_head cellhd;
+ struct WaterParams wp;
+ struct options parm;
G_gisinit(argv[0]);
@@ -261,122 +264,115 @@
G_get_set_window(&cellhd);
- conv = G_database_units_to_meters_factor();
+ WaterParams_init(&wp);
- mixx = cellhd.west * conv;
- maxx = cellhd.east * conv;
- miyy = cellhd.south * conv;
- mayy = cellhd.north * conv;
+ wp.conv = G_database_units_to_meters_factor();
- stepx = cellhd.ew_res * conv;
- stepy = cellhd.ns_res * conv;
- /* step = amin1(stepx,stepy); */
- step = (stepx + stepy) / 2.;
- mx = cellhd.cols;
- my = cellhd.rows;
- xmin = 0.;
- ymin = 0.;
- xp0 = xmin + stepx / 2.;
- yp0 = ymin + stepy / 2.;
- xmax = xmin + stepx * (float)mx;
- ymax = ymin + stepy * (float)my;
- hhc = hhmax = 0.;
+ wp.mixx = cellhd.west * wp.conv;
+ wp.maxx = cellhd.east * wp.conv;
+ wp.miyy = cellhd.south * wp.conv;
+ wp.mayy = cellhd.north * wp.conv;
+ wp.stepx = cellhd.ew_res * wp.conv;
+ wp.stepy = cellhd.ns_res * wp.conv;
+ /* wp.step = amin1(wp.stepx,wp.stepy); */
+ wp.step = (wp.stepx + wp.stepy) / 2.;
+ wp.mx = cellhd.cols;
+ wp.my = cellhd.rows;
+ wp.xmin = 0.;
+ wp.ymin = 0.;
+ wp.xp0 = wp.xmin + wp.stepx / 2.;
+ wp.yp0 = wp.ymin + wp.stepy / 2.;
+ wp.xmax = wp.xmin + wp.stepx * (float)wp.mx;
+ wp.ymax = wp.ymin + wp.stepy * (float)wp.my;
+ wp.hhc = wp.hhmax = 0.;
+
#if 0
- bxmi = 2093113. * conv;
- bymi = 731331. * conv;
- bxma = 2093461. * conv;
- byma = 731529. * conv;
- bresx = 2. * conv;
- bresy = 2. * conv;
- maxwab = 100000;
+ wp.bxmi = 2093113. * wp.conv;
+ wp.bymi = 731331. * wp.conv;
+ wp.bxma = 2093461. * wp.conv;
+ wp.byma = 731529. * wp.conv;
+ wp.bresx = 2. * wp.conv;
+ wp.bresy = 2. * wp.conv;
+ wp.maxwab = 100000;
- mx2o = (int)((bxma - bxmi) / bresx);
- my2o = (int)((byma - bymi) / bresy);
+ wp.mx2o = (int)((wp.bxma - wp.bxmi) / wp.bresx);
+ wp.my2o = (int)((wp.byma - wp.bymi) / wp.bresy);
/* relative small box coordinates: leave 1 grid layer for overlap */
- bxmi = bxmi - mixx + stepx;
- bymi = bymi - miyy + stepy;
- bxma = bxma - mixx - stepx;
- byma = byma - miyy - stepy;
- mx2 = mx2o - 2 * ((int)(stepx / bresx));
- my2 = my2o - 2 * ((int)(stepy / bresy));
+ wp.bxmi = wp.bxmi - wp.mixx + wp.stepx;
+ wp.bymi = wp.bymi - wp.miyy + wp.stepy;
+ wp.bxma = wp.bxma - wp.mixx - wp.stepx;
+ wp.byma = wp.byma - wp.miyy - wp.stepy;
+ wp.mx2 = wp.mx2o - 2 * ((int)(wp.stepx / wp.bresx));
+ wp.my2 = wp.my2o - 2 * ((int)(wp.stepy / wp.bresy));
#endif
- elevin = parm.elevin->answer;
- wdepth = parm.wdepth->answer;
- dxin = parm.dxin->answer;
- dyin = parm.dyin->answer;
- detin = parm.detin->answer;
- tranin = parm.tranin->answer;
- tauin = parm.tauin->answer;
- manin = parm.manin->answer;
- tc = parm.tc->answer;
- et = parm.et->answer;
- conc = parm.conc->answer;
- flux = parm.flux->answer;
- erdep = parm.erdep->answer;
- outwalk = parm.outwalk->answer;
+ wp.elevin = parm.elevin->answer;
+ wp.wdepth = parm.wdepth->answer;
+ wp.dxin = parm.dxin->answer;
+ wp.dyin = parm.dyin->answer;
+ wp.detin = parm.detin->answer;
+ wp.tranin = parm.tranin->answer;
+ wp.tauin = parm.tauin->answer;
+ wp.manin = parm.manin->answer;
+ wp.tc = parm.tc->answer;
+ wp.et = parm.et->answer;
+ wp.conc = parm.conc->answer;
+ wp.flux = parm.flux->answer;
+ wp.erdep = parm.erdep->answer;
+ wp.outwalk = parm.outwalk->answer;
- /* sscanf(parm.nwalk->answer, "%d", &maxwa); */
- sscanf(parm.niter->answer, "%d", ×ec);
- sscanf(parm.outiter->answer, "%d", &iterout);
-/* sscanf(parm.density->answer, "%d", &ldemo); */
- sscanf(parm.diffc->answer, "%lf", &frac);
- sscanf(parm.maninval->answer, "%lf", &manin_val);
+ /* sscanf(parm.nwalk->answer, "%d", &wp.maxwa); */
+ sscanf(parm.niter->answer, "%d", &wp.timesec);
+ sscanf(parm.outiter->answer, "%d", &wp.iterout);
+/* sscanf(parm.density->answer, "%d", &wp.ldemo); */
+ sscanf(parm.diffc->answer, "%lf", &wp.frac);
+ sscanf(parm.maninval->answer, "%lf", &wp.manin_val);
/* Recompute timesec from user input in minutes
* to real timesec in seconds */
- timesec = timesec * 60;
- iterout = iterout * 60;
- if ((timesec / iterout) > 100)
+ wp.timesec = wp.timesec * 60;
+ wp.iterout = wp.iterout * 60;
+ if ((wp.timesec / wp.iterout) > 100)
G_message(_("More than 100 files are going to be created !!!!!"));
/* compute how big the raster is and set this to appr 2 walkers per cell */
if (parm.nwalk->answer == NULL) {
- maxwa = mx * my * 2;
- rwalk = (double)(mx * my * 2.);
- G_message(_("default nwalk=%d, rwalk=%f"), maxwa, rwalk);
+ wp.maxwa = wp.mx * wp.my * 2;
+ wp.rwalk = (double)(wp.mx * wp.my * 2.);
+ G_message(_("default nwalk=%d, rwalk=%f"), wp.maxwa, wp.rwalk);
}
else {
- sscanf(parm.nwalk->answer, "%d", &maxwa);
- rwalk = (double)maxwa;
+ sscanf(parm.nwalk->answer, "%d", &wp.maxwa);
+ wp.rwalk = (double)wp.maxwa;
}
/*rwalk = (double) maxwa; */
- if (conv != 1.0)
- G_message(_("Using metric conversion factor %f, step=%f"), conv,
- step);
+ if (wp.conv != 1.0)
+ G_message(_("Using metric conversion factor %f, step=%f"), wp.conv,
+ wp.step);
+ init_library_globals(&wp);
- if ((tc == NULL) && (et == NULL) && (conc == NULL) && (flux == NULL) &&
- (erdep == NULL))
+ if ((wp.tc == NULL) && (wp.et == NULL) && (wp.conc == NULL) && (wp.flux == NULL) &&
+ (wp.erdep == NULL))
G_warning(_("You are not outputting any raster or site files"));
ret_val = input_data();
if (ret_val != 1)
G_fatal_error(_("Input failed"));
- /* mandatory for si,sigma */
+ alloc_grids_sediment();
- si = G_alloc_matrix(my, mx);
- sigma = G_alloc_matrix(my, mx);
-
- /* memory allocation for output grids */
-
- dif = G_alloc_fmatrix(my, mx);
- if (erdep != NULL || et != NULL)
- er = G_alloc_fmatrix(my, mx);
-
G_srand48(rand1);
grad_check();
-
- if (et != NULL)
- erod(si);
+ init_grids_sediment();
/* treba dat output pre topoerdep */
main_loop();
- if (tserie == NULL) {
+ /* always true for sediment? */
+ if (wp.tserie == NULL) {
ii = output_data(0, 1.);
if (ii != 1)
G_fatal_error(_("Cannot write raster maps"));
Modified: grass/trunk/raster/r.sim/r.sim.water/main.c
===================================================================
--- grass/trunk/raster/r.sim/r.sim.water/main.c 2016-04-26 15:38:36 UTC (rev 68319)
+++ grass/trunk/raster/r.sim/r.sim.water/main.c 2016-04-27 18:45:35 UTC (rev 68320)
@@ -88,7 +88,7 @@
/********************************/
-#include <grass/waterglobs.h>
+#include <grass/simlib.h>
/****************************************/
@@ -102,6 +102,10 @@
double x_orig, y_orig;
static int rand1 = 12345;
struct GModule *module;
+ struct Cell_head cellhd;
+ struct WaterParams wp;
+ struct options parm;
+ struct flags flag;
G_gisinit(argv[0]);
@@ -298,54 +302,56 @@
G_get_set_window(&cellhd);
- conv = G_database_units_to_meters_factor();
+ WaterParams_init(&wp);
- G_debug(3, "Conversion factor is set to: %f", conv);
+ wp.conv = G_database_units_to_meters_factor();
- mixx = conv * cellhd.west;
- maxx = conv * cellhd.east;
- miyy = conv * cellhd.south;
- mayy = conv * cellhd.north;
+ G_debug(3, "Conversion factor is set to: %f", wp.conv);
- stepx = cellhd.ew_res * conv;
- stepy = cellhd.ns_res * conv;
+ wp.mixx = wp.conv * cellhd.west;
+ wp.maxx = wp.conv * cellhd.east;
+ wp.miyy = wp.conv * cellhd.south;
+ wp.mayy = wp.conv * cellhd.north;
+
+ wp.stepx = cellhd.ew_res * wp.conv;
+ wp.stepy = cellhd.ns_res * wp.conv;
/* step = amin1(stepx,stepy); */
- step = (stepx + stepy) / 2.;
- mx = cellhd.cols;
- my = cellhd.rows;
- x_orig = cellhd.west * conv;
- y_orig = cellhd.south * conv; /* do we need this? */
- xmin = 0.;
- ymin = 0.;
- xp0 = xmin + stepx / 2.;
- yp0 = ymin + stepy / 2.;
- xmax = xmin + stepx * (float)mx;
- ymax = ymin + stepy * (float)my;
+ wp.step = (wp.stepx + wp.stepy) / 2.;
+ wp.mx = cellhd.cols;
+ wp.my = cellhd.rows;
+ x_orig = cellhd.west * wp.conv;
+ y_orig = cellhd.south * wp.conv; /* do we need this? */
+ wp.xmin = 0.;
+ wp.ymin = 0.;
+ wp.xp0 = wp.xmin + wp.stepx / 2.;
+ wp.yp0 = wp.ymin + wp.stepy / 2.;
+ wp.xmax = wp.xmin + wp.stepx * (float)wp.mx;
+ wp.ymax = wp.ymin + wp.stepy * (float)wp.my;
- G_debug(3, "xmax: %f, ymax: %f", xmax, ymax);
+ G_debug(3, "xmax: %f, ymax: %f", wp.xmax, wp.ymax);
- ts = flag.tserie->answer;
+ wp.ts = flag.tserie->answer;
- elevin = parm.elevin->answer;
- dxin = parm.dxin->answer;
- dyin = parm.dyin->answer;
- rain = parm.rain->answer;
- infil = parm.infil->answer;
- traps = parm.traps->answer;
- manin = parm.manin->answer;
- depth = parm.depth->answer;
- disch = parm.disch->answer;
- err = parm.err->answer;
- outwalk = parm.outwalk->answer;
+ wp.elevin = parm.elevin->answer;
+ wp.dxin = parm.dxin->answer;
+ wp.dyin = parm.dyin->answer;
+ wp.rain = parm.rain->answer;
+ wp.infil = parm.infil->answer;
+ wp.traps = parm.traps->answer;
+ wp.manin = parm.manin->answer;
+ wp.depth = parm.depth->answer;
+ wp.disch = parm.disch->answer;
+ wp.err = parm.err->answer;
+ wp.outwalk = parm.outwalk->answer;
G_debug(3, "Parsing numeric parameters");
- sscanf(parm.niter->answer, "%d", ×ec);
- sscanf(parm.outiter->answer, "%d", &iterout);
- sscanf(parm.diffc->answer, "%lf", &frac);
- sscanf(parm.hmax->answer, "%lf", &hhmax);
- sscanf(parm.halpha->answer, "%lf", &halpha);
- sscanf(parm.hbeta->answer, "%lf", &hbeta);
+ sscanf(parm.niter->answer, "%d", &wp.timesec);
+ sscanf(parm.outiter->answer, "%d", &wp.iterout);
+ sscanf(parm.diffc->answer, "%lf", &wp.frac);
+ sscanf(parm.hmax->answer, "%lf", &wp.hhmax);
+ sscanf(parm.halpha->answer, "%lf", &wp.halpha);
+ sscanf(parm.hbeta->answer, "%lf", &wp.hbeta);
G_debug(3, "Parsing rain parameters");
@@ -355,28 +361,28 @@
/* if no rain unique value input */
if (parm.rainval->answer == NULL) {
/*No rain input so use default */
- sscanf(RAINVAL, "%lf", &rain_val);
+ sscanf(RAINVAL, "%lf", &wp.rain_val);
/* if rain unique input exist, load it */
}
else {
/*Unique value input only */
- sscanf(parm.rainval->answer, "%lf", &rain_val);
+ sscanf(parm.rainval->answer, "%lf", &wp.rain_val);
}
/* if Rain map exists */
}
else {
/*Map input, so set rain_val to -999.99 */
if (parm.rainval->answer == NULL) {
- rain_val = -999.99;
+ wp.rain_val = -999.99;
}
else {
/*both map and unique value exist */
/*Choose the map, discard the unique value */
- rain_val = -999.99;
+ wp.rain_val = -999.99;
}
}
/* Report the final value of rain_val */
- G_debug(3, "rain_val is set to: %f\n", rain_val);
+ G_debug(3, "rain_val is set to: %f\n", wp.rain_val);
/* if no Mannings map, then: */
if (parm.manin->answer == NULL) {
@@ -384,28 +390,28 @@
/* if no Mannings unique value input */
if (parm.maninval->answer == NULL) {
/*No Mannings input so use default */
- sscanf(MANINVAL, "%lf", &manin_val);
+ sscanf(MANINVAL, "%lf", &wp.manin_val);
/* if Mannings unique input value exists, load it */
}
else {
/*Unique value input only */
- sscanf(parm.maninval->answer, "%lf", &manin_val);
+ sscanf(parm.maninval->answer, "%lf", &wp.manin_val);
}
/* if Mannings map exists */
}
else {
/* Map input, set manin_val to -999.99 */
if (parm.maninval->answer == NULL) {
- manin_val = -999.99;
+ wp.manin_val = -999.99;
}
else {
/*both map and unique value exist */
/*Choose map, discard the unique value */
- manin_val = -999.99;
+ wp.manin_val = -999.99;
}
}
/* Report the final value of manin_val */
- G_debug(1, "manin_val is set to: %f\n", manin_val);
+ G_debug(1, "manin_val is set to: %f\n", wp.manin_val);
/* if no infiltration map, then: */
if (parm.infil->answer == NULL) {
@@ -413,68 +419,63 @@
/*if no infiltration unique value input */
if (parm.infilval->answer == NULL) {
/*No infiltration unique value so use default */
- sscanf(INFILVAL, "%lf", &infil_val);
+ sscanf(INFILVAL, "%lf", &wp.infil_val);
/* if infiltration unique value exists, load it */
}
else {
/*unique value input only */
- sscanf(parm.infilval->answer, "%lf", &infil_val);
+ sscanf(parm.infilval->answer, "%lf", &wp.infil_val);
}
/* if infiltration map exists */
}
else {
/* Map input, set infil_val to -999.99 */
if (parm.infilval->answer == NULL) {
- infil_val = -999.99;
+ wp.infil_val = -999.99;
}
else {
/*both map and unique value exist */
/*Choose map, discard the unique value */
- infil_val = -999.99;
+ wp.infil_val = -999.99;
}
}
/* Report the final value of infil_val */
- G_debug(1, "infil_val is set to: %f\n", infil_val);
+ G_debug(1, "infil_val is set to: %f\n", wp.infil_val);
/* Recompute timesec from user input in minutes
* to real timesec in seconds */
- timesec = timesec * 60.0;
- iterout = iterout * 60.0;
- if ((timesec / iterout) > 100.0 && ts == 1)
+ wp.timesec = wp.timesec * 60.0;
+ wp.iterout = wp.iterout * 60.0;
+ if ((wp.timesec / wp.iterout) > 100.0 && wp.ts == 1)
G_message(_("More than 100 files are going to be created !!!!!"));
/* compute how big the raster is and set this to appr 2 walkers per cell */
if (parm.nwalk->answer == NULL) {
- maxwa = mx * my * 2;
- rwalk = (double)(mx * my * 2.);
- G_message(_("default nwalk=%d, rwalk=%f"), maxwa, rwalk);
+ wp.maxwa = wp.mx * wp.my * 2;
+ wp.rwalk = (double)(wp.mx * wp.my * 2.);
+ G_message(_("default nwalk=%d, rwalk=%f"), wp.maxwa, wp.rwalk);
}
else {
- sscanf(parm.nwalk->answer, "%d", &maxwa);
- rwalk = (double)maxwa;
+ sscanf(parm.nwalk->answer, "%d", &wp.maxwa);
+ wp.rwalk = (double)wp.maxwa;
}
/* rwalk = (double) maxwa; */
- if (conv != 1.0)
- G_message(_("Using metric conversion factor %f, step=%f"), conv,
- step);
+ if (wp.conv != 1.0)
+ G_message(_("Using metric conversion factor %f, step=%f"), wp.conv,
+ wp.step);
- if ((depth == NULL) && (disch == NULL) && (err == NULL))
+ init_library_globals(&wp);
+
+ if ((wp.depth == NULL) && (wp.disch == NULL) && (wp.err == NULL))
G_warning(_("You are not outputting any raster maps"));
ret_val = input_data();
if (ret_val != 1)
G_fatal_error(_("Input failed"));
+ alloc_grids_water();
- /* memory allocation for output grids */
- G_debug(1, "beginning memory allocation for output grids");
-
- gama = G_alloc_matrix(my, mx);
- if (err != NULL)
- gammas = G_alloc_matrix(my, mx);
- dif = G_alloc_fmatrix(my, mx);
-
G_debug(1, "seeding randoms");
G_srand48(rand1);
grad_check();
Modified: grass/trunk/raster/r.sim/simlib/Makefile
===================================================================
--- grass/trunk/raster/r.sim/simlib/Makefile 2016-04-26 15:38:36 UTC (rev 68319)
+++ grass/trunk/raster/r.sim/simlib/Makefile 2016-04-27 18:45:35 UTC (rev 68320)
@@ -9,4 +9,4 @@
include $(MODULE_TOPDIR)/include/Make/Lib.make
-default: $(ARCH_INCDIR)/waterglobs.h lib
+default: $(ARCH_INCDIR)/waterglobs.h $(ARCH_INCDIR)/simlib.h lib
Modified: grass/trunk/raster/r.sim/simlib/hydro.c
===================================================================
--- grass/trunk/raster/r.sim/simlib/hydro.c 2016-04-26 15:38:36 UTC (rev 68319)
+++ grass/trunk/raster/r.sim/simlib/hydro.c 2016-04-27 18:45:35 UTC (rev 68320)
@@ -24,6 +24,7 @@
#include <grass/glocale.h>
#include <grass/waterglobs.h>
+#include <grass/simlib.h>
/*
* Soeren 8. Mar 2011 TODO:
@@ -32,8 +33,6 @@
*
*/
-struct options parm;
-struct flags flag;
struct _points points;
char *elevin;
@@ -43,11 +42,12 @@
char *infil;
char *traps;
char *manin;
-/* char *observation; */
char *depth;
char *disch;
char *err;
-char *outwalk;
+char *outwalk;
+char *observation;
+char *logfile;
char *mapset;
char *mscale;
char *tserie;
@@ -68,8 +68,6 @@
struct seed seed;
-struct Cell_head cellhd;
-
double xmin, ymin, xmax, ymax;
double mayy, miyy, maxx, mixx;
int mx, my;
Modified: grass/trunk/raster/r.sim/simlib/input.c
===================================================================
--- grass/trunk/raster/r.sim/simlib/input.c 2016-04-26 15:38:36 UTC (rev 68319)
+++ grass/trunk/raster/r.sim/simlib/input.c 2016-04-27 18:45:35 UTC (rev 68320)
@@ -7,6 +7,7 @@
#include <grass/glocale.h>
#include <grass/linkm.h>
#include <grass/gmath.h>
+#include <grass/simlib.h>
#include <grass/waterglobs.h>
@@ -18,7 +19,255 @@
static void copy_matrix_undef_double_to_float_values(int rows, int cols, double **source, float **target);
static void copy_matrix_undef_float_values(int rows, int cols, float **source, float **target);
+/*!
+ * \brief Initialize WaterParams structure.
+ */
+void WaterParams_init(struct WaterParams *wp)
+{
+ /* this is little bit lengthy and perhaps error-prone
+ * but it simplifies initialization since then there is no
+ * difference in between initialization in water and sediment
+ * for the variables which are not used and would have been
+ * initialized if they were just global variables */
+ wp->ymin = 0;
+ wp->xmax = 0;
+ wp->ymax = 0;
+ wp->mayy = 0;
+ wp->miyy = 0;
+ wp->maxx = 0;
+ wp->mixx = 0;
+ wp->mx = 0;
+ wp->my = 0;
+ wp->mx2 = 0;
+ wp->my2 = 0;
+ wp->bxmi = 0;
+ wp->bymi = 0;
+ wp->bxma = 0;
+ wp->byma = 0;
+ wp->bresx = 0;
+ wp->bresy = 0;
+ wp->maxwab = 0;
+ wp->step = 0;
+ wp->conv = 0;
+
+ wp->frac = 0;
+ wp->bxmi = 0;
+ wp->bymi = 0;
+
+ wp->hbeta = 0;
+ wp->hhmax = 0;
+ wp->sisum = 0;
+ wp->vmean = 0;
+ wp->infsum = 0;
+ wp->infmean = 0;
+ wp->maxw = 0;
+ wp->maxwa = 0;
+ wp->nwalk = 0;
+ wp->rwalk = 0;
+ wp->bresx = 0;
+ wp->bresy = 0;
+ wp->xrand = 0;
+ wp->yrand = 0;
+ wp->stepx = 0;
+ wp->stepy = 0;
+ wp->xp0 = 0;
+ wp->yp0 = 0;
+ wp->chmean = 0;
+ wp->si0 = 0;
+ wp->deltap = 0;
+ wp->deldif = 0;
+ wp->cch = 0;
+ wp->hhc = 0;
+ wp->halpha = 0;
+ wp->eps = 0;
+ wp->maxwab = 0;
+ wp->nstack = 0;
+ wp->iterout = 0;
+ wp->mx2o = 0;
+ wp->my2o = 0;
+ wp->miter = 0;
+ wp->nwalka = 0;
+ wp->timec = 0;
+ wp->ts = 0;
+ wp->timesec = 0;
+
+ wp->rain_val = 0;
+ wp->manin_val = 0;
+ wp->infil_val = 0;
+
+ wp->elevin = NULL;
+ wp->dxin = NULL;
+ wp->dyin = NULL;
+ wp->rain = NULL;
+ wp->infil = NULL;
+ wp->traps = NULL;
+ wp->manin = NULL;
+ wp->depth = NULL;
+ wp->disch = NULL;
+ wp->err = NULL;
+ wp->outwalk = NULL;
+ wp->observation = NULL;
+ wp->logfile = NULL;
+ wp->mapset = NULL;
+ wp->mscale = NULL;
+ wp->tserie = NULL;
+
+ wp->wdepth = NULL;
+ wp->detin = NULL;
+ wp->tranin = NULL;
+ wp->tauin = NULL;
+ wp->tc = NULL;
+ wp->et = NULL;
+ wp->conc = NULL;
+ wp->flux = NULL;
+ wp->erdep = NULL;
+
+ wp->rainval = NULL;
+ wp->maninval = NULL;
+ wp->infilval = NULL;
+}
+
+/*!
+ * \brief Initialize global variables in the library.
+ */
+void init_library_globals(struct WaterParams *wp)
+{
+ /* this is little bit lengthy and perhaps error-prone
+ * but it separates library from its interface */
+ ymin = wp->ymin;
+ xmax = wp->xmax;
+ ymax = wp->ymax;
+ mayy = wp->mayy;
+ miyy = wp->miyy;
+ maxx = wp->maxx;
+ mixx = wp->mixx;
+ mx = wp->mx;
+ my = wp->my;
+ mx2 = wp->mx2;
+ my2 = wp->my2;
+
+ bxmi = wp->bxmi;
+ bymi = wp->bymi;
+ bxma = wp->bxma;
+ byma = wp->byma;
+ bresx = wp->bresx;
+ bresy = wp->bresy;
+ maxwab = wp->maxwab;
+ step = wp->step;
+ conv = wp->conv;
+
+ frac = wp->frac;
+ bxmi = wp->bxmi;
+ bymi = wp->bymi;
+
+ hbeta = wp->hbeta;
+ hhmax = wp->hhmax;
+ sisum = wp->sisum;
+ vmean = wp->vmean;
+ infsum = wp->infsum;
+ infmean = wp->infmean;
+ maxw = wp->maxw;
+ maxwa = wp->maxwa;
+ nwalk = wp->nwalk;
+ rwalk = wp->rwalk;
+ bresx = wp->bresx;
+ bresy = wp->bresy;
+ xrand = wp->xrand;
+ yrand = wp->yrand;
+ stepx = wp->stepx;
+ stepy = wp->stepy;
+ xp0 = wp->xp0;
+ yp0 = wp->yp0;
+ chmean = wp->chmean;
+ si0 = wp->si0;
+ deltap = wp->deltap;
+ deldif = wp->deldif;
+ cch = wp->cch;
+ hhc = wp->hhc;
+ halpha = wp->halpha;
+ eps = wp->eps;
+ maxwab = wp->maxwab;
+ nstack = wp->nstack;
+ iterout = wp->iterout;
+ mx2o = wp->mx2o;
+ my2o = wp->my2o;
+ miter = wp->miter;
+ nwalka = wp->nwalka;
+ timec = wp->timec;
+ ts = wp->ts;
+ timesec = wp->timesec;
+
+ rain_val = wp->rain_val;
+ manin_val = wp->manin_val;
+ infil_val = wp->infil_val;
+
+ elevin = wp->elevin;
+ dxin = wp->dxin;
+ dyin = wp->dyin;
+ rain = wp->rain;
+ infil = wp->infil;
+ traps = wp->traps;
+ manin = wp->manin;
+ depth = wp->depth;
+ disch = wp->disch;
+ err = wp->err;
+ outwalk = wp->outwalk;
+ observation = wp->observation;
+ logfile = wp->logfile;
+ mapset = wp->mapset;
+ mscale = wp->mscale;
+ tserie = wp->tserie;
+
+ wdepth = wp->wdepth;
+ detin = wp->detin;
+ tranin = wp->tranin;
+ tauin = wp->tauin;
+ tc = wp->tc;
+ et = wp->et;
+ conc = wp->conc;
+ flux = wp->flux;
+ erdep = wp->erdep;
+
+ rainval = wp->rainval;
+ maninval = wp->maninval;
+ infilval = wp->infilval;
+}
+
+/* we do the allocation inside because we anyway need to set the variables */
+
+void alloc_grids_water()
+{
+ /* memory allocation for output grids */
+ G_debug(1, "beginning memory allocation for output grids");
+
+ gama = G_alloc_matrix(my, mx);
+ if (err != NULL)
+ gammas = G_alloc_matrix(my, mx);
+ dif = G_alloc_fmatrix(my, mx);
+}
+
+void alloc_grids_sediment()
+{
+ /* mandatory for si,sigma */
+
+ si = G_alloc_matrix(my, mx);
+ sigma = G_alloc_matrix(my, mx);
+
+ /* memory allocation for output grids */
+
+ dif = G_alloc_fmatrix(my, mx);
+ if (erdep != NULL || et != NULL)
+ er = G_alloc_fmatrix(my, mx);
+}
+
+void init_grids_sediment()
+{
+ /* this should be fulfilled for sediment but not water */
+ if (et != NULL)
+ erod(si);
+}
+
/* ************************************************************** */
/* GRASS input procedures, allocations */
/* *************************************************************** */
Modified: grass/trunk/raster/r.sim/simlib/observation_points.c
===================================================================
--- grass/trunk/raster/r.sim/simlib/observation_points.c 2016-04-26 15:38:36 UTC (rev 68319)
+++ grass/trunk/raster/r.sim/simlib/observation_points.c 2016-04-27 18:45:35 UTC (rev 68320)
@@ -19,24 +19,26 @@
struct line_cats *cts;
double x, y;
int type, cat, i;
+ struct Cell_head cellhd;
- if(parm.observation->answer != NULL)
+ if(observation != NULL)
if_log += 1;
- if(parm.logfile->answer != NULL)
+ if(logfile != NULL)
if_log += 1;
/* Nothing to do */
if(if_log == 0)
return;
+ /* why both are required ? */
if(if_log == 1)
G_fatal_error("Observation vector map and logfile must be provided");
Vect_set_open_level(1);
- if (Vect_open_old(&Map, parm.observation->answer, "") < 0)
- G_fatal_error(_("Unable to open vector map <%s>"), parm.observation->answer);
+ if (Vect_open_old(&Map, observation, "") < 0)
+ G_fatal_error(_("Unable to open vector map <%s>"), observation);
Vect_rewind(&Map);
@@ -46,6 +48,9 @@
/* Initialize point structure */
init_points(&points, 128);
+ /* get the current computational region set in the runtime */
+ G_get_set_window(&cellhd);
+
/* Read all vector points */
while(1) {
type = Vect_read_next_line(&Map, pts, cts);
@@ -56,7 +61,7 @@
if(type == -1) {
Vect_close(&Map);
- G_fatal_error(_("Unable to read points from map %s"), parm.observation->answer);
+ G_fatal_error(_("Unable to read points from map %s"), observation);
}
if(type == GV_POINT) {
@@ -74,10 +79,10 @@
Vect_close(&Map);
/* Open the logfile */
- points.output = fopen(parm.logfile->answer, "w");
+ points.output = fopen(logfile, "w");
if(points.output == NULL)
- G_fatal_error(_("Unable to open observation logfile %s for writing"), parm.logfile->answer);
+ G_fatal_error(_("Unable to open observation logfile %s for writing"), logfile);
points.is_open = 1;
Modified: grass/trunk/raster/r.sim/simlib/output.c
===================================================================
--- grass/trunk/raster/r.sim/simlib/output.c 2016-04-26 15:38:36 UTC (rev 68319)
+++ grass/trunk/raster/r.sim/simlib/output.c 2016-04-27 18:45:35 UTC (rev 68320)
@@ -122,7 +122,7 @@
/* Write the output walkers */
output_walker_as_vector(tt_minutes, ndigit, ×tamp);
- Rast_set_window(&cellhd);
+ /* we write in the same region as we used for reading */
if (my != Rast_window_rows())
G_fatal_error("OOPS: rows changed from %d to %d", mx,
@@ -636,7 +636,7 @@
/* float a1,a2; */
- Rast_set_window(&cellhd);
+ /* we write in the same region as we used for reading */
if (et) {
et_cell = Rast_allocate_f_buf();
Added: grass/trunk/raster/r.sim/simlib/simlib.h
===================================================================
--- grass/trunk/raster/r.sim/simlib/simlib.h (rev 0)
+++ grass/trunk/raster/r.sim/simlib/simlib.h 2016-04-27 18:45:35 UTC (rev 68320)
@@ -0,0 +1,96 @@
+#ifndef __SIMLIB_H__
+#define __SIMLIB_H__
+
+/*! \file simlib.h
+ * \brief This is the interface for the simlib (SIMWE) library.
+ */
+
+struct WaterParams
+{
+ double xmin, ymin, xmax, ymax;
+ double mayy, miyy, maxx, mixx;
+ int mx, my;
+ int mx2, my2;
+
+ double bxmi, bymi, bxma, byma, bresx, bresy;
+ int maxwab;
+ double step, conv;
+
+ double frac;
+
+ double hbeta;
+ double hhmax, sisum, vmean;
+ double infsum, infmean;
+ int maxw, maxwa, nwalk;
+ double rwalk, xrand, yrand;
+ double stepx, stepy, xp0, yp0;
+ double chmean, si0, deltap, deldif, cch, hhc, halpha;
+ double eps;
+ int nstack;
+ int iterout, mx2o, my2o;
+ int miter, nwalka;
+ double timec;
+ int ts, timesec;
+
+ double rain_val;
+ double manin_val;
+ double infil_val;
+
+ char *elevin;
+ char *dxin;
+ char *dyin;
+ char *rain;
+ char *infil;
+ char *traps;
+ char *manin;
+ char *depth;
+ char *disch;
+ char *err;
+ char *outwalk;
+ char *observation;
+ char *logfile;
+ char *mapset;
+ char *mscale;
+ char *tserie;
+
+ char *wdepth;
+ char *detin;
+ char *tranin;
+ char *tauin;
+ char *tc;
+ char *et;
+ char *conc;
+ char *flux;
+ char *erdep;
+
+ char *rainval;
+ char *maninval;
+ char *infilval;
+};
+
+void WaterParams_init(struct WaterParams *wp);
+void init_library_globals(struct WaterParams *wp);
+void alloc_grids_water();
+void alloc_grids_sediment();
+void init_grids_sediment();
+
+int input_data(void);
+int grad_check(void);
+void main_loop(void);
+int output_data(int, double);
+
+struct options
+{
+ struct Option *elevin, *dxin, *dyin, *rain, *infil, *traps, *manin,
+ *observation, *depth, *disch, *err, *outwalk, *nwalk, *niter, *outiter,
+ *density, *diffc, *hmax, *halpha, *hbeta, *wdepth, *detin, *tranin,
+ *tauin, *tc, *et, *conc, *flux, *erdep, *rainval, *maninval,
+ *infilval, *logfile;
+};
+
+struct flags
+{
+ struct Flag *mscale, *tserie;
+};
+
+#endif /* __SIMLIB_H__ */
Property changes on: grass/trunk/raster/r.sim/simlib/simlib.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Modified: grass/trunk/raster/r.sim/simlib/waterglobs.h
===================================================================
--- grass/trunk/raster/r.sim/simlib/waterglobs.h 2016-04-26 15:38:36 UTC (rev 68319)
+++ grass/trunk/raster/r.sim/simlib/waterglobs.h 2016-04-27 18:45:35 UTC (rev 68320)
@@ -14,11 +14,12 @@
extern char *infil;
extern char *traps;
extern char *manin;
-/* extern char *observation; */
extern char *depth;
extern char *disch;
extern char *err;
extern char *outwalk;
+extern char *observation;
+extern char *logfile;
extern char *mapset;
extern char *mscale;
extern char *tserie;
@@ -37,24 +38,6 @@
extern char *maninval;
extern char *infilval;
-struct options
-{
- struct Option *elevin, *dxin, *dyin, *rain, *infil, *traps, *manin,
- *observation, *depth, *disch, *err, *outwalk, *nwalk, *niter, *outiter,
- *density, *diffc, *hmax, *halpha, *hbeta, *wdepth, *detin, *tranin,
- *tauin, *tc, *et, *conc, *flux, *erdep, *rainval, *maninval,
- *infilval, *logfile;
-};
-
-extern struct options parm;
-
-struct flags
-{
- struct Flag *mscale, *tserie;
-};
-
-extern struct flags flag;
-
struct seed
{
long int is1, is2;
@@ -62,9 +45,6 @@
extern struct seed seed;
-
-extern struct Cell_head cellhd;
-
struct _points
{
double *x; /* x coor for each point */
@@ -77,11 +57,7 @@
};
extern struct _points points;
-extern int input_data(void);
-extern int grad_check(void);
extern void erod(double **);
-extern void main_loop(void);
-extern int output_data(int, double);
extern int output_et(void);
extern double simwe_rand(void);
extern double gasdev(void);
Added: grass/trunk/raster/r.sim/test/test.sh
===================================================================
--- grass/trunk/raster/r.sim/test/test.sh (rev 0)
+++ grass/trunk/raster/r.sim/test/test.sh 2016-04-27 18:45:35 UTC (rev 68320)
@@ -0,0 +1,47 @@
+#!/bin/sh
+
+# assumes full NC SPM
+
+# execute the following in old version in a separate mapset
+# the switch to another mapset and execute it with a new version
+
+g.region rural_1m res=2 -p
+v.surf.rst -d input=elev_lid792_bepts elevation=elev_lid792_2m slope=dx_2m aspect=dy_2m tension=15 smooth=1.5 npmin=150
+r.sim.water -t elevation=elev_lid792_2m dx=dx_2m dy=dy_2m rain_value=50 infil_value=0 man_value=0.05 depth=wdp_2m discharge=disch_2m nwalkers=100000 niterations=120 output_step=20
+r.mapcalc "tranin = 0.001"
+r.mapcalc "detin = 0.001"
+r.mapcalc "tauin = 0.01"
+g.copy rast=wdp_2m.120,wdp_2m
+r.sim.sediment elevation=elev_lid792_2m dx=dx_2m dy=dy_2m water_depth=wdp_2m detachment_coeff=detin transport_coeff=tranin shear_stress=tauin man_value=0.05 nwalkers=1000000 niterations=120 transport_capacity=tcapacity tlimit_erosion_deposition=erdepmax sediment_flux=sedflow erosion_deposition=erdepsimwe
+
+g.copy rast=disch_2m.120,disch_2m
+
+# execute the following only in one of the mapsets
+
+REF_MAPSET=rsim_before
+
+r.mapcalc "elev_lid792_2m_diff = elev_lid792_2m@$REF_MAPSET - elev_lid792_2m"
+r.mapcalc "dx_2m_diff = dx_2m@$REF_MAPSET - "
+r.mapcalc "dy_2m_diff = dy_2m@$REF_MAPSET - dy_2m"
+
+r.mapcalc "wdp_2m_diff = wdp_2m@$REF_MAPSET - wdp_2m"
+r.mapcalc "disch_2m_diff = disch_2m@$REF_MAPSET - disch_2m"
+
+r.mapcalc "tcapacity_diff = tcapacity@$REF_MAPSET - tcapacity"
+r.mapcalc "erdepmax_diff = erdepmax@$REF_MAPSET - erdepmax"
+r.mapcalc "sedflow_diff = sedflow@$REF_MAPSET - sedflow"
+r.mapcalc "erdepsimwe_diff = erdepsimwe@$REF_MAPSET - erdepsimwe"
+
+# expecting zeros
+
+r.univar elev_lid792_2m_diff
+r.univar dx_2m_diff
+r.univar dy_2m_diff
+
+r.univar wdp_2m_diff
+r.univar disch_2m_diff
+
+r.univar tcapacity_diff
+r.univar erdepmax_diff
+r.univar sedflow_diff
+r.univar erdepsimwe_diff
Property changes on: grass/trunk/raster/r.sim/test/test.sh
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/x-sh
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list